From ba10f693bdd41de5095f9ac7b281cf5a74918e90 Mon Sep 17 00:00:00 2001 From: Dariusz Murakowski Date: Thu, 25 Jun 2015 22:55:21 -0400 Subject: [PATCH] Counter (dummy species) for viruses that have escaped T cell recognition. --- reaction.cpp | 153 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++- reaction.h | 34 ++++++++++++- ss.cpp | 45 ++++++++++++++++++ ss.h | 6 +++ 4 files changed, 235 insertions(+), 3 deletions(-) diff --git a/reaction.cpp b/reaction.cpp index 84bd466..857a985 100644 --- a/reaction.cpp +++ b/reaction.cpp @@ -240,7 +240,10 @@ double KillingReaction::recalc(Reaction *rxn_in) case RxnType::NTVirus: case RxnType::NTVirusDeath: case RxnType::AAKilling: - case RxnType::CTLaaActivation: { // did not change V->count + case RxnType::CTLaaActivation: + case RxnType::EscapeCounting: + case RxnType::EscapeAACounting: + { // did not change V->count propensity = this->recalc(); // TODO //std::cerr << "shouldn't have hit here in KillingReaction::recalc(Reaction*)....\n"; @@ -509,7 +512,10 @@ double AAKillingReaction::recalc(Reaction *rxn_in) case RxnType::Virus: case RxnType::VirusDeath: case RxnType::Killing: - case RxnType::CTLaaActivation: { // did not change V->count + case RxnType::CTLaaActivation: + case RxnType::EscapeCounting: + case RxnType::EscapeAACounting: + { // did not change V->count propensity = this->recalc(); // TODO //std::cerr << "shouldn't have hit here in KillingReaction::recalc(Reaction*)....\n"; @@ -579,3 +585,146 @@ void CTLaaActivationReaction::fire() } +double EscapeCountingReaction::recalc() +{ + this->propensity = 0.0; + + esc->count = 0; + virus_map::iterator iter = V->pop.begin(), + end = V->pop.end(); + for(; iter != end; ++iter) { + if (esc->recognized(iter->first) == 0.0) + esc->count += iter->second; + } + + return this->propensity; +} + +// adapted from KillingReaction::recalc(Reaction *rxn_in) +double EscapeCountingReaction::recalc(Reaction *rxn_in) +{ + RxnType type = get_rxn_type(rxn_in); + + switch (type) { + + case RxnType::Virus: { // increased V->count + //VirusReaction *rxn = dynamic_cast(rxn_in); + virus_map::iterator iter = V->lastModified; + if (esc->recognized(iter->first) == 0.0) + esc->count += 1; + break; + } + + case RxnType::VirusDeath: + case RxnType::Killing: { // decreased V->count + virus_map::iterator iter = V->lastModified; + if (iter == V->pop.end()) { + if (esc->recognized(V->lastDeleted) == 0.0) + esc->count -= 1; + } + else { + if (esc->recognized(iter->first) == 0.0) + esc->count -= 1; + } + + break; + } + + case RxnType::Elementary: + case RxnType::CTLActivation: + case RxnType::NTVirus: + case RxnType::NTVirusDeath: + case RxnType::AAKilling: + case RxnType::CTLaaActivation: + case RxnType::EscapeCounting: + case RxnType::EscapeAACounting: + { // did not change V->count + break; + } + + } // switch (type) + + if (V->lastModified == V->pop.end()) { // deleted such that none is left + // TODO + } + + //return this->recalc(); + return 0.0; + +} + +void EscapeCountingReaction::fire() +{ + // do nothing +} + +double EscapeAACountingReaction::recalc() +{ + this->propensity = 0.0; + + esc->count = 0; + NTVirus_map::iterator iter = V->pop.begin(), + end = V->pop.end(); + for(; iter != end; ++iter) { + if (esc->recognized(iter->first) == 0.0) + esc->count += iter->second; + } + + return this->propensity; +} + +double EscapeAACountingReaction::recalc(Reaction *rxn_in) +{ + RxnType type = get_rxn_type(rxn_in); + + switch (type) { + + case RxnType::NTVirus: { // increased V->count + //VirusReaction *rxn = dynamic_cast(rxn_in); + NTVirus_map::iterator iter = V->lastModified; + if (esc->recognized(iter->first) == 0.0) + esc->count += 1; + break; + } + + case RxnType::NTVirusDeath: + case RxnType::AAKilling: { // decreased V->count + NTVirus_map::iterator iter = V->lastModified; + if (iter == V->pop.end()) { + if (esc->recognized(V->lastDeleted) == 0.0) + esc->count -= 1; + } + else { + if (esc->recognized(iter->first) == 0.0) + esc->count -= 1; + } + + break; + } + + case RxnType::Elementary: + case RxnType::CTLActivation: + case RxnType::Virus: + case RxnType::VirusDeath: + case RxnType::Killing: + case RxnType::CTLaaActivation: + case RxnType::EscapeCounting: + case RxnType::EscapeAACounting: + { // did not change V->count + break; + } + + } // switch (type) + + if (V->lastModified == V->pop.end()) { // deleted such that none is left + // TODO + } + + //return this->recalc(); + return 0.0; +} + +void EscapeAACountingReaction::fire() +{ + // do nothing +} diff --git a/reaction.h b/reaction.h index f8a29b7..d95d274 100644 --- a/reaction.h +++ b/reaction.h @@ -7,7 +7,7 @@ typedef std::set MutatedSiteSequence; -enum class RxnType {Elementary, Virus, VirusDeath, Killing, CTLActivation, NTVirus, NTVirusDeath, AAKilling, CTLaaActivation}; +enum class RxnType {Elementary, Virus, VirusDeath, Killing, CTLActivation, NTVirus, NTVirusDeath, AAKilling, CTLaaActivation, EscapeCounting, EscapeAACounting}; class Reaction { public: @@ -164,6 +164,34 @@ public: void fire(); }; +// fake reaction to track number of escape variants +// propensity is always 0, but esc->count gets updated +class EscapeCountingReaction : public Reaction { +public: + EscapeCountingReaction() : Reaction() {} + EscapeCountingReaction(double k) : Reaction(k) {} + + VirusSpecies* V; + CTLSpecies* esc; + double recalc(); + double recalc(Reaction *rxn); + void fire(); +}; + +// fake reaction to track number of escape variants +// propensity is always 0, but esc->count gets updated +class EscapeAACountingReaction : public Reaction { +public: + EscapeAACountingReaction() : Reaction() {} + EscapeAACountingReaction(double k) : Reaction(k) {} + + NTVirusSpecies* V; + CTLaaSpecies* esc; + double recalc(); + double recalc(Reaction *rxn); + void fire(); +}; + inline RxnType get_rxn_type(const Reaction* const rxn) { @@ -185,6 +213,10 @@ inline RxnType get_rxn_type(const Reaction* const rxn) return RxnType::AAKilling; else if (dynamic_cast(rxn)) return RxnType::CTLaaActivation; + else if (dynamic_cast(rxn)) + return RxnType::EscapeCounting; + else if (dynamic_cast(rxn)) + return RxnType::EscapeAACounting; else throw "unknown reaction type?!"; } diff --git a/ss.cpp b/ss.cpp index a9ee543..ea4696a 100644 --- a/ss.cpp +++ b/ss.cpp @@ -64,6 +64,9 @@ bool does_species_appear_in_reactants(Species *s, Reaction *rxn) case RxnType::CTLActivation: if (s == dynamic_cast(rxn)->V || s == dynamic_cast(rxn)->Tfrom) { return true; } // rxn->Tto is irrelevant to rate break; + case RxnType::EscapeCounting: + if (s == dynamic_cast(rxn)->V) { return true; } + break; case RxnType::NTVirus: if (s == dynamic_cast(rxn)->V) { return true; } break; @@ -76,6 +79,9 @@ bool does_species_appear_in_reactants(Species *s, Reaction *rxn) case RxnType::CTLaaActivation: if (s == dynamic_cast(rxn)->V || s == dynamic_cast(rxn)->Tfrom) { return true; } // rxn->Tto is irrelevant to rate break; + case RxnType::EscapeAACounting: + if (s == dynamic_cast(rxn)->V) { return true; } + break; } return appears; @@ -163,6 +169,14 @@ bool does_affect(Reaction *newrxn, Reaction *rxn) break; } + // these reactions never fire + case RxnType::EscapeCounting: + case RxnType::EscapeAACounting: + { + affects = false; + break; + } + } // switch (newT) @@ -657,6 +671,20 @@ void importEpitope(RunParameters_SS &r, Species_parray &species, Rxn_parray &rea } + if (r.trackEscaped) { // escape counter + CTLSpecies* esc = new CTLSpecies(*M); + esc->name = "Iesc"; esc->count = 0; + species.push_back(esc); + print_spec.push_back(esc); + + EscapeCountingReaction* rx = new EscapeCountingReaction(0.0); + rx->name = "Iesc -(count)-> Iesc"; + rx->V = V; rx->esc = esc; + reactions.push_back(rx); update_rxn_dependents(reactions); + + rx->recalc(); // initialize with the correct count + } + species.push_back(M); // memory CTLSpecies* N = new CTLSpecies(*M); // naive @@ -843,6 +871,20 @@ void importEpitope_Potts(RunParameters_SS &r, Species_parray &species, Rxn_parra } + if (r.trackEscaped) { // escape counter + CTLaaSpecies* esc = new CTLaaSpecies(*M); + esc->name = "Iesc"; esc->count = 0; + species.push_back(esc); + print_spec.push_back(esc); + + EscapeAACountingReaction* rx = new EscapeAACountingReaction(0.0); + rx->name = "Iesc -(count)-> Iesc"; + rx->V = V; rx->esc = esc; + reactions.push_back(rx); update_rxn_dependents(reactions); + + rx->recalc(); // initialize with the correct count + } + species.push_back(M); // memory CTLaaSpecies* N = new CTLaaSpecies(*M); // naive @@ -1000,6 +1042,7 @@ void usage() " -rate_X (double) X in [s,u,p,a,r,d,dprime,b,e,w,aprime,g,h]\n" " -const_ND (int) number of effector CTL divisions (generations) [default=9]\n" " -vol (double) volume scaling factor; affects rate constants, not initial values [default=1.0]\n" +" -track_esc flag to count escaped viruses (in Iesc species)\n" ; std::cout.flush(); } @@ -1077,6 +1120,8 @@ int main(int argc, char *argv[]) { else if (strcmp(argv[i],"-vol")==0) { if (++i==argc) break; else r.vol=strtodouble(argv[i]); } + else if (strcmp(argv[i],"-track_esc")==0) { r.trackEscaped = true; } + else if (strcmp(argv[i],"-h")==0 || strcmp(argv[i],"--help")==0) { usage(); return 0; } else printf("Unrecognized argument! '%s'\n",argv[i]); diff --git a/ss.h b/ss.h index d84e542..df9e0a3 100644 --- a/ss.h +++ b/ss.h @@ -138,6 +138,10 @@ public: bool useEpitope; // Include selective pressure on an epitope bool useVerbose; // If true, print extra information while program is running + bool trackEscaped; // Maintain counter to track number of escaped viruses? + // If true, use a dummy T cell species "Iesc" with special reactions, + // updating amount when viruses aren't recognized. + std::vector > initPop; // Initial population sequences std::vector initFrac; // Initial population fraction @@ -226,6 +230,8 @@ public: usePotts = false; + trackEscaped = false; + } -- 2.7.4