Counter (dummy species) for viruses that have escaped T cell recognition.

diff --git a/reaction.cpp b/reaction.cpp
index 84bd466..857a985 100644 (file)
--- 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<Reaction*>(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<Reaction*>(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 (file)
--- a/reaction.h
+++ b/reaction.h
@@ -7,7 +7,7 @@
 
 typedef std::set<unsigned int> 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<const CTLaaActivationReaction*>(rxn))
         return RxnType::CTLaaActivation;
+    else if (dynamic_cast<const EscapeCountingReaction*>(rxn))
+        return RxnType::EscapeCounting;
+    else if (dynamic_cast<const EscapeAACountingReaction*>(rxn))
+        return RxnType::EscapeAACounting;
     else
         throw "unknown reaction type?!";
 }

diff --git a/ss.cpp b/ss.cpp
index a9ee543..ea4696a 100644 (file)
--- 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<CTLActivationReaction*>(rxn)->V || s == dynamic_cast<CTLActivationReaction*>(rxn)->Tfrom) { return true; }  // rxn->Tto is irrelevant to rate
         break;
+    case RxnType::EscapeCounting:
+        if (s == dynamic_cast<EscapeCountingReaction*>(rxn)->V) { return true; }
+        break;
     case RxnType::NTVirus:
         if (s == dynamic_cast<NTVirusReaction*>(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<CTLaaActivationReaction*>(rxn)->V || s == dynamic_cast<CTLaaActivationReaction*>(rxn)->Tfrom) { return true; }  // rxn->Tto is irrelevant to rate
         break;
+    case RxnType::EscapeAACounting:
+        if (s == dynamic_cast<EscapeAACountingReaction*>(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 (file)
--- 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<std::set<unsigned int> > initPop;   // Initial population sequences
     std::vector<double>                  initFrac;  // Initial population fraction
@@ -226,6 +230,8 @@ public:
 
         usePotts = false;
 
+        trackEscaped = false;
+
     }