Initialize multiple CTL clones+epitopes. Output at intervals. Pretty-print the sequence.

diff --git a/pop_ss.h b/pop_ss.h
index 764cf9b..12e21ab 100644 (file)
--- a/pop_ss.h
+++ b/pop_ss.h
@@ -56,22 +56,32 @@ public:
     bool escaped(const MutatedSiteSequence &mutated_sites) const;
 };
 
+
 class Species {
 public:
     long count;
     std::string name;
+    Species() : count(0) {}
+    Species(long c) : count(c) {}
 };
 
 class SimpleSpecies : public Species {
+    SimpleSpecies() : Species() {}
+    SimpleSpecies(long c) : Species(c) {}
 };
 
 class VirusSpecies : public Species {
 public:
+    VirusSpecies() : Species() {}
+    VirusSpecies(long c) : Species(c) {}
     virus_map pop;
 };
 
 class CTLSpecies : public Species {
 public:
+    CTLSpecies() : Species() {}
+    CTLSpecies(long c) : Species(c) {}
+
     std::vector<EpitopeRecognizer> epitopes;
     std::vector<double> affinity;
     size_t num_ep;

diff --git a/reaction.cpp b/reaction.cpp
index 20f55e3..c4ac28a 100644 (file)
--- a/reaction.cpp
+++ b/reaction.cpp
@@ -38,11 +38,13 @@ void ElementaryReaction::fire()
     std::vector<int>::iterator reacN = this->reactant_stoich.begin(),
                                reacN_end = this->reactant_stoich.end();
 
+    /*
     for (; (reac != reac_end) && (reacN != reacN_end); ++reac, ++reacN) {
         //if (reac->count - reacN >= 0) {
             (*reac)->count = (*reac)->count - *reacN;
         //}
     }
+    */
 
     std::vector<SimpleSpecies*>::iterator prod = this->products.begin(),
                                           prod_end = this->products.end();

diff --git a/reaction.h b/reaction.h
index 9865256..3ea501a 100644 (file)
--- a/reaction.h
+++ b/reaction.h
@@ -9,6 +9,9 @@ typedef std::set<unsigned int> MutatedSiteSequence;
 
 class Reaction {
 public:
+    Reaction() {}
+    Reaction(double k) : rate_constant(k) {}
+
     double rate_constant;  // must manually ensure it includes division by statistical degeneracy factor (in case of reaction order >1 for any species)
     double propensity;
     virtual double recalc() = 0;
@@ -19,10 +22,13 @@ public:
 // mass action kinetics
 class ElementaryReaction : public Reaction {
 public:
-    std::vector<SimpleSpecies*> reactants;
-    std::vector<SimpleSpecies*> products;
+    ElementaryReaction() : Reaction() {}
+    ElementaryReaction(double k) : Reaction(k) {}
+
+    std::vector<SimpleSpecies*> reactants;  // anything that participates in reaction rate
+    std::vector<SimpleSpecies*> products;  // effects of reaction
     std::vector<int> reactant_stoich;   // reaction order
-    std::vector<int> product_stoich;
+    std::vector<int> product_stoich;  // must be negative for reactants that are consumed!
     double recalc();
     void fire();
 };
@@ -32,6 +38,9 @@ public:
 // depending on energy, with mutation
 class VirusReaction : public Reaction {
 public:
+    VirusReaction() : Reaction() {}
+    VirusReaction(double k) : Reaction(k) {}
+
     VirusSpecies* V;
     Hamiltonian H;
     double recalc();
@@ -43,6 +52,9 @@ public:
 // in principle energy-dependent, but not here
 class VirusDeathReaction : public Reaction {
 public:
+    VirusDeathReaction() : Reaction() {}
+    VirusDeathReaction(double k) : Reaction(k) {}
+
     VirusSpecies* V;
     Hamiltonian H;
     double recalc();
@@ -53,6 +65,9 @@ public:
 // depending on epitope affinity
 class KillingReaction : public Reaction {
 public:
+    KillingReaction() : Reaction() {}
+    KillingReaction(double k) : Reaction(k) {}
+
     VirusSpecies* V;
     CTLSpecies* T;
     double recalc();
@@ -62,6 +77,9 @@ public:
 // V_s + T_k -> V_s + T_k'
 class CTLActivationReaction : public Reaction {
 public:
+    CTLActivationReaction() : Reaction() {}
+    CTLActivationReaction(double k) : Reaction(k) {}
+
     VirusSpecies* V;
     CTLSpecies* Tfrom;
     CTLSpecies* Tto;

diff --git a/ss.cpp b/ss.cpp
index 7f08189..96841bb 100644 (file)
--- a/ss.cpp
+++ b/ss.cpp
@@ -52,20 +52,24 @@ void run(RunParameters_SS &r, unsigned seed) {
     gsl_rng_set(rnd,seed);
 
     r.setFiles();
-    if (r.importState) importState(r);      // <infile>.st
-    if (r.useEpitope)  importEpitope(r);    // <infile>.ep
+
+    // master arrays of (pointers to) species and reactions
+    Species_parray species;
+    Rxn_parray reactions;
 
     std::string couplingsFile = "neutral_2site.j";
     Hamiltonian H(couplingsFile);
     //double mu = 6.0e-5;
-    double mu = 1.0e-1;
+    //double mu = 1.0e-1;
 
-    std::vector<Species*> species;
-    VirusSpecies s1; s1.count = 100;
-    s1.pop[Virus(H,mu)] = s1.count;
+    // initialize virus
+    VirusSpecies s1; s1.count = r.n;
+    s1.pop[Virus(H,r.mu)] = s1.count;
     species.push_back(&s1);
 
-    Rxn_parray reactions;
+    if (r.importState) importState(r);      // <infile>.st
+    if (r.useEpitope)  importEpitope(r,species,reactions,&s1);    // <infile>.ep
+
     // V -> V + V'
     VirusReaction* r1 = new VirusReaction();
     r1->rate_constant = 3.0;
@@ -82,9 +86,10 @@ void run(RunParameters_SS &r, unsigned seed) {
     Species_parray print_spec;
     print_spec.push_back(&s1);
 
-    double t_end = 1.0;
+    double t_end = 10.0;
     long max_steps = LONG_MAX;
-    simulate(reactions, species, t_end, max_steps, print_spec);
+    double sample_interval = 1.0;
+    simulate(reactions, species, t_end, max_steps, sample_interval, print_spec);
 
 
 #if 0
@@ -497,7 +502,7 @@ void run(RunParameters_SS &r, unsigned seed) {
 }
 
 
-void simulate(Rxn_parray &reactions, Species_parray &species, double t_end, long max_steps, const Species_parray &print_spec)
+void simulate(Rxn_parray &reactions, Species_parray &species, double t_end, long max_steps, double sample_interval, const Species_parray &print_spec)
 {
     long step = 0;
 
@@ -505,6 +510,8 @@ void simulate(Rxn_parray &reactions, Species_parray &species, double t_end, long
     double t_remain = t_end;
     double dt = 0.0;
 
+    double t_next_sample = t + sample_interval;
+
     double total_propensity = 0.0;
 
     Rxn_parray::iterator rxn, end;
@@ -555,15 +562,21 @@ void simulate(Rxn_parray &reactions, Species_parray &species, double t_end, long
 
         ++step;
 
-        std::cout << t << '\t';
-        for (Species_parray::const_iterator spec = print_spec.begin(),
-                                       end = print_spec.end();
-          spec != end; ++spec)
-        {
-            // TODO: prettier representation
-            std::cout << (*spec)->count << '\t';
+
+        // print copy numbers
+        if (t_next_sample <= t) {
+            t_next_sample += sample_interval;
+
+            std::cout << t << '\t';
+            for (Species_parray::const_iterator spec = print_spec.begin(),
+                                           end = print_spec.end();
+              spec != end; ++spec)
+            {
+                // TODO: prettier representation
+                std::cout << (*spec)->count << '\t';
+            }
+            std::cout << '\n';
         }
-        std::cout << '\n';
 
 #if 0
         // print copy numbers
@@ -652,64 +665,110 @@ void importState(RunParameters_SS &r) {
 
 // load epitope definitions from file
 
-void importEpitope(RunParameters_SS &r) {
+void importEpitope(RunParameters_SS &r, Species_parray &species, Rxn_parray &reactions, VirusSpecies *V)
+{
+
+    if ((r.num_CTL_clones != r.epitopefiles.size()) || 
+        (r.num_CTL_clones != r.init_CTL_num.size())) {
+        perror((std::string("ERROR in importEpitope: invalid number of CTL clones").c_str()));
+        exit(1);
+    }
 
-    std::ifstream input(r.epitopeInfile.c_str());   // <infile>.ep
-    if (!input) { perror((std::string("ERROR in importEpitope: ") + r.epitopeInfile).c_str()); exit(1); }
+    for (size_t i = 0; i != r.num_CTL_clones; ++i) {
+        CTLSpecies* T = new CTLSpecies(r.init_CTL_num[i]);
 
-    std::string readStr;
-    while (std::getline(input,readStr)) {
+        std::ifstream input(r.epitopefiles[i].c_str());   // <infile>.ep
+        if (!input) { perror((std::string("ERROR in importEpitope: ") + r.epitopeInfile).c_str()); exit(1); }
 
-        std::string word;
-        unsigned int site;
-        double aff;
-        std::vector<unsigned int> tmpEp;
-        size_t pos0 = 0;
-        size_t posBar = readStr.find('|',pos0);
-        size_t posEnd = std::string::npos; //readStr.size();
+        std::string readStr;
+        while (std::getline(input,readStr)) {
 
-        // could use std::noskipws to keep tab ('\t') characters?
+            T->num_ep += 1;
+            T->epitopes.push_back(EpitopeRecognizer());
 
-        std::stringstream readStrStrm(std::string(readStr,pos0,posBar-pos0));
+            std::string word;
+            unsigned int site;
+            double aff;
+            std::vector<unsigned int> tmpEp;
+            size_t pos0 = 0;
+            size_t posBar = readStr.find('|',pos0);
+            size_t posEnd = std::string::npos; //readStr.size();
 
-        // first entry is affinity
-        readStrStrm >> word;
-        aff = strtodouble(word);
-        printf("%.6e ",aff);
+            // could use std::noskipws to keep tab ('\t') characters?
 
-        // next entries are WT sites
-        while (readStrStrm >> word) {
-            std::cout << word << " "; // XXX
-            std::istringstream i(word);
-            if (i >> site)
-                tmpEp.push_back(site);
-            else    // must have encountered '|'
-                break;
-        }
-        r.eWT.push_back(tmpEp);
+            std::stringstream readStrStrm(std::string(readStr,pos0,posBar-pos0));
 
-        tmpEp.clear();  // reset temp epitope
-        readStrStrm.str("");  readStrStrm.clear();  // reset stream
+            // first entry is affinity
+            readStrStrm >> word;
+            aff = strtodouble(word);
+            printf("%.6e ",aff);
 
-        std::cout << "| "; // XXX
+            T->affinity.push_back(aff);
+
+            // next entries are WT sites
+            while (readStrStrm >> word) {
+                std::cout << word << " "; // XXX
+                std::istringstream i(word);
+                if (i >> site)
+                    tmpEp.push_back(site);
+                else    // must have encountered '|'
+                    break;
+            }
+            //r.eWT.push_back(tmpEp);
+            T->epitopes.back().epitopeWT = tmpEp;
+
+            tmpEp.clear();  // reset temp epitope
+            readStrStrm.str("");  readStrStrm.clear();  // reset stream
+
+            std::cout << "| "; // XXX
+
+            readStrStrm.str(std::string(readStr,posBar+1,posEnd/*-posBar-1*/));
+            while (readStrStrm >> word) {
+                std::cout << word << " "; // XXX
+                std::istringstream i(word);
+                if (i >> site)
+                    tmpEp.push_back(site);
+                else
+                    break;
+            }
+            //r.eMut.push_back(tmpEp);
+            T->epitopes.back().epitopeMut = tmpEp;
+
+            //r.numEpitopes++;
+
+            std::cout << "\n"; // XXX
 
-        readStrStrm.str(std::string(readStr,posBar+1,posEnd/*-posBar-1*/));
-        while (readStrStrm >> word) {
-            std::cout << word << " "; // XXX
-            std::istringstream i(word);
-            if (i >> site)
-                tmpEp.push_back(site);
-            else
-                break;
         }
-        r.eMut.push_back(tmpEp);
 
-        r.numEpitopes++;
+        species.push_back(T);   // effector
 
-        std::cout << "\n"; // XXX
+        CTLSpecies* M = new CTLSpecies(*T);     // memory
+        species.push_back(M);
 
-    }
+        CTLSpecies* N = new CTLSpecies(*T);
+        species.push_back(N);   // naive
+
+        KillingReaction* r1 = new KillingReaction(1e-5);
+        r1->T = T;  r1->V = V;
+        reactions.push_back(r1);
+
+        CTLActivationReaction* r2 = new CTLActivationReaction(1e-6);
+        r2->Tfrom = M;  r2->Tto = T;  r2->V = V;
+        reactions.push_back(r2);
+
+        CTLActivationReaction* r3 = new CTLActivationReaction(1e-7);
+        r3->Tfrom = M;  r3->Tto = T;  r3->V = V;
+        reactions.push_back(r3);
 
+
+        for (size_t i = 0; i < T->num_ep; ++i) {  // XXX
+            std::cout << T->count << '\t'
+                      << T->affinity[i] << '\t'
+                      << T->epitopes[i].epitopeWT << '|'
+                      << T->epitopes[i].epitopeMut << '\n';
+        }
+
+    }
 }
 
 
@@ -737,6 +796,7 @@ void usage()
 " -o  (string)        output file stem\n"
 " -s  (string)        input state file, containing initial population fraction\n"
 " -e  (string)        input file containing targeted epitope\n"
+" -en (int)           corresponding number of that CTL (ordered!)\n"
 " -n  (int/double)    population size\n"
 " -g  (int/double)    number of generations\n"
 " -mu (double)        mutation rate\n"
@@ -794,7 +854,8 @@ int main(int argc, char *argv[]) {
         else if (strcmp(argv[i],"-i")==0) { if (++i==argc) break; else r.infile=argv[i];                            }
         else if (strcmp(argv[i],"-o")==0) { if (++i==argc) break; else r.outfile=argv[i];                           }
         else if (strcmp(argv[i],"-s")==0) { if (++i==argc) break; else { r.statefile=argv[i]; r.importState=true; } }
-        else if (strcmp(argv[i],"-e")==0) { if (++i==argc) break; else { r.epitopefiles.push_back(argv[i]); r.useEpitope=true; } }
+        else if (strcmp(argv[i],"-e")==0) { if (++i==argc) break; else { r.epitopefiles.push_back(argv[i]); r.useEpitope=true; ++r.num_CTL_clones; } }
+        else if (strcmp(argv[i],"-en")==0) { if (++i==argc) break; else { r.init_CTL_num.push_back((long)strtodouble(argv[i])); r.useEpitope=true; } }
         
         else if (strcmp(argv[i],"-n")==0)  { if (++i==argc) break; else r.n=(unsigned int)strtodouble(argv[i]); }
         else if (strcmp(argv[i],"-g")==0)  { if (++i==argc) break; else r.g=(unsigned int)strtodouble(argv[i]); }

diff --git a/ss.h b/ss.h
index 747094a..be805c4 100644 (file)
--- a/ss.h
+++ b/ss.h
@@ -70,6 +70,28 @@ inline double strtodouble(const std::string &s) {
 
 }
 
+// print a set of mutations to stream as a CSV-like string
+std::ostream& operator<<(std::ostream& os, std::vector<unsigned int> s)
+{
+    std::vector<unsigned int>::iterator i = s.begin(), end = s.end();
+    if (i == end)
+        return os;  // stop immediately
+    os << *i;  // print first element
+    for (++i; i != end; ++i)  // then with comma as prefix
+        os << ',' << *i;
+    return os;
+}
+std::ostream& operator<<(std::ostream& os, MutatedSiteSequence s)
+{
+    MutatedSiteSequence::iterator i = s.begin(), end = s.end();
+    if (i == end)
+        return os;  // stop immediately
+    os << *i;  // print first element
+    for (++i; i != end; ++i)  // then with comma as prefix
+        os << ',' << *i;
+    return os;
+}
+
 
 // PROGRAM SETTINGS
 
@@ -128,6 +150,8 @@ public:
     
 
     std::vector<std::string> epitopefiles;
+    std::vector<long> init_CTL_num;
+    unsigned int num_CTL_clones;
     
     
     RunParameters_SS() {
@@ -162,8 +186,11 @@ public:
         write_mod = 1;
 
         numEpitopes = 0;
+
+        num_CTL_clones = 0;
         
     }
+
     void setFiles() {
         
         if (strcmp(directory.c_str(),"")!=0) {
@@ -204,8 +231,8 @@ public:
 
 void run(RunParameters_SS &r);
 void importState(RunParameters_SS &r);
-void importEpitope(RunParameters_SS &r);
-void simulate(Rxn_parray &reactions, Species_parray &species, double t_end, long max_steps, const Species_parray &print_spec);
+void importEpitope(RunParameters_SS &r, Species_parray &species, Rxn_parray &reactions, VirusSpecies *V);
+void simulate(Rxn_parray &reactions, Species_parray &species, double t_end, long max_steps, double sample_interval, const Species_parray &print_spec);
 
 
 #endif