Framework for multiple divisions (i.e. generations) of effector T cells.

diff --git a/ss.cpp b/ss.cpp
index 0fc519c..95ed444 100644 (file)
--- a/ss.cpp
+++ b/ss.cpp
@@ -386,6 +386,14 @@ void importState_Potts(RunParameters_SS &r)
 }
 
 
+// helper function adapted from /usr/include/c++/4.4.7/bits/basic_string.h
+  inline std::string
+  to_string_uint(unsigned __val)
+  { return __gnu_cxx::__to_xstring<std::string>(&std::vsnprintf,
+                       4 * sizeof(unsigned),
+                       "%u", __val); }
+
+
 // load epitope definitions from file
 
 void importEpitope(RunParameters_SS &r, Species_parray &species, Rxn_parray &reactions, VirusSpecies *V, Species_parray &print_spec)
@@ -397,7 +405,7 @@ void importEpitope(RunParameters_SS &r, Species_parray &species, Rxn_parray &rea
     }
 
     for (size_t i = 0; i != r.num_CTL_clones; ++i) {
-        CTLSpecies* M = new CTLSpecies("E",r.init_CTL_numM[i]);
+        CTLSpecies* M = new CTLSpecies("M",r.init_CTL_numM[i]);
 
         std::ifstream input(r.epitopefiles[i].c_str());   // <infile>.ep
         if (!input) { perror((std::string("ERROR in importEpitope: ") + r.epitopeInfile).c_str()); exit(1); }
@@ -469,12 +477,24 @@ void importEpitope(RunParameters_SS &r, Species_parray &species, Rxn_parray &rea
         species.push_back(N);
 
         CTLSpecies* T = new CTLSpecies(*M);     // effector
-        T->name = "T";  T->count = r.init_CTL_numT[i];
+        T->name = "T1";  T->count = r.init_CTL_numT[i];
         species.push_back(T);
 
-        print_spec.push_back(T);
-        print_spec.push_back(N);
         print_spec.push_back(M);
+        print_spec.push_back(N);
+        print_spec.push_back(T);
+
+        std::vector<CTLSpecies*> Tgen(r.num_T_gen);
+        Tgen[0] = T;
+        for (unsigned i=1; i<r.num_T_gen; i++) {
+            CTLSpecies* Tn = new CTLSpecies(*T);
+            //Tn->name = "T" + std::to_string(i+1);
+            Tn->name = "T" + to_string_uint(i+1);
+            Tn->count = 0;
+            Tgen[i] = Tn;
+            species.push_back(Tn);
+            print_spec.push_back(Tn);
+        }
 
         {   // p
             KillingReaction* rx = new KillingReaction(r.rate_p);
@@ -613,12 +633,24 @@ void importEpitope_Potts(RunParameters_SS &r, Species_parray &species, Rxn_parra
         species.push_back(N);
 
         CTLaaSpecies* T = new CTLaaSpecies(*M);     // effector
-        T->name = "T";  T->count = r.init_CTL_numT[i];
+        T->name = "T1";  T->count = r.init_CTL_numT[i];
         species.push_back(T);
 
-        print_spec.push_back(T);
-        print_spec.push_back(N);
         print_spec.push_back(M);
+        print_spec.push_back(N);
+        print_spec.push_back(T);
+
+        std::vector<CTLaaSpecies*> Tgen(r.num_T_gen);
+        Tgen[0] = T;
+        for (unsigned i=1; i<r.num_T_gen; i++) {
+            CTLaaSpecies* Tn = new CTLaaSpecies(*T);
+            //Tn->name = "T" + std::to_string(i+1);
+            Tn->name = "T" + to_string_uint(i+1);
+            Tn->count = 0;
+            Tgen[i] = Tn;
+            species.push_back(Tn);
+            print_spec.push_back(Tn);
+        }
 
         {   // p
             AAKillingReaction* rx = new AAKillingReaction(r.rate_p);
@@ -719,6 +751,7 @@ void usage()
 " -potts              flag for nucleotide/amino-acid/Potts states\n"
 " -numruns (int)      number of trajectories to simulate\n"
 " -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"
 ;   std::cout.flush();
 }
 
@@ -791,6 +824,8 @@ int main(int argc, char *argv[]) {
         else if (strcmp(argv[i],"-rate_g")==0) { if (++i==argc) break; else r.rate_g=strtodouble(argv[i]); }
         else if (strcmp(argv[i],"-rate_h")==0) { if (++i==argc) break; else r.rate_h=strtodouble(argv[i]); }
 
+        else if (strcmp(argv[i],"-const_ND")==0) { if (++i==argc) break; else r.num_T_gen=(unsigned int)strtodouble(argv[i]); }
+
         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 eb0c85b..d971f19 100644 (file)
--- a/ss.h
+++ b/ss.h
@@ -173,6 +173,8 @@ public:
     double rate_aprime;
     double rate_g;
     double rate_h;
+
+    unsigned int num_T_gen;    // ND, number of effector cell generations (i.e. divisions)
     
     bool usePotts;
     std::string seq2stateInfile;
@@ -214,6 +216,8 @@ public:
         rate_g = 0.03;
         rate_h = 0.03;
 
+        num_T_gen = 9;
+
         usePotts = false;
 
     }