Separate function for running the dynamics.

diff --git a/reaction.cpp b/reaction.cpp
index 3eca7ee..20f55e3 100644 (file)
--- a/reaction.cpp
+++ b/reaction.cpp
@@ -202,3 +202,27 @@ void KillingReaction::fire()
         V->pop.erase(iter);
 }
 
+
+double CTLActivationReaction::recalc()
+{
+    propensity = rate_constant * Tfrom->count;
+
+    double partial_sum = 0.0;
+    virus_map::iterator iter = V->pop.begin(),
+                        end  = V->pop.end();
+    for(; iter != end; ++iter) {
+        partial_sum += iter->second * Tfrom->recognized(iter->first);
+    }
+
+    propensity *= partial_sum;
+    return propensity;
+
+}
+
+void CTLActivationReaction::fire()
+{
+    Tfrom->count -= 1;
+    Tto->count += 1;
+}
+
+

diff --git a/reaction.h b/reaction.h
index 7501941..9865256 100644 (file)
--- a/reaction.h
+++ b/reaction.h
@@ -59,7 +59,8 @@ public:
     void fire();
 };
 
-class CTLTransitionReaction : public Reaction {
+// V_s + T_k -> V_s + T_k'
+class CTLActivationReaction : public Reaction {
 public:
     VirusSpecies* V;
     CTLSpecies* Tfrom;

diff --git a/ss.cpp b/ss.cpp
index d7b09b8..7f08189 100644 (file)
--- a/ss.cpp
+++ b/ss.cpp
@@ -5,6 +5,7 @@
 #include <cstdlib>
 #include <cmath>
 #include <numeric>  // for partial_sum
+#include <climits>
 
 #include <gsl/gsl_rng.h>
 
@@ -14,8 +15,6 @@
 #include "ss.h"
 #include "reaction.h"
 
-typedef std::vector<Reaction*> Rxn_parray;
-
 gsl_rng *rnd;  //pointer to random number generator state
 
 // generate a random number generator seed from /dev/urandom
@@ -52,7 +51,10 @@ void run(RunParameters_SS &r, unsigned seed) {
     //gsl_rng_set(rnd,rand());
     gsl_rng_set(rnd,seed);
 
-#if 1
+    r.setFiles();
+    if (r.importState) importState(r);      // <infile>.st
+    if (r.useEpitope)  importEpitope(r);    // <infile>.ep
+
     std::string couplingsFile = "neutral_2site.j";
     Hamiltonian H(couplingsFile);
     //double mu = 6.0e-5;
@@ -77,6 +79,16 @@ void run(RunParameters_SS &r, unsigned seed) {
     r2->V = &s1;
     reactions.push_back(r2);
 
+    Species_parray print_spec;
+    print_spec.push_back(&s1);
+
+    double t_end = 1.0;
+    long max_steps = LONG_MAX;
+    simulate(reactions, species, t_end, max_steps, print_spec);
+
+
+#if 0
+
 /*
     double total_propensity = 0.0;
     for (std::vector<Reaction*>::iterator rxn = reactions.begin(),
@@ -485,6 +497,111 @@ 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)
+{
+    long step = 0;
+
+    double t = 0.0;
+    double t_remain = t_end;
+    double dt = 0.0;
+
+    double total_propensity = 0.0;
+
+    Rxn_parray::iterator rxn, end;
+
+    for (rxn = reactions.begin(), end = reactions.end(); rxn != end; ++rxn) {
+        (*rxn)->recalc();
+        total_propensity += (*rxn)->propensity;
+    }
+
+
+    while (1) {
+
+        double rand = (1.0-gsl_rng_uniform_pos(rnd)) * total_propensity;  // uniform on (0,R]
+        double a_sum = 0.0;
+
+        // time to reaction
+        dt = -log(gsl_rng_uniform(rnd)) / total_propensity;
+        t += dt;
+        t_remain -= dt;
+
+        // stop if time ran out
+        if (t_remain < 0.0)
+            break;
+
+        // determine next reaction
+        for (rxn = reactions.begin(), end = reactions.end(); rxn != end; ++rxn) {
+            a_sum += (*rxn)->propensity;
+            if (rand <= a_sum) {
+                (*rxn)->fire();
+                break;
+            }
+        }
+        
+        // TODO: better handle case with total propensity == 0.0
+
+        if (rxn == end)
+            break;
+
+        // recalculate rates
+        total_propensity = 0.0;
+        for (rxn = reactions.begin(), end = reactions.end(); rxn != end; ++rxn) {
+            (*rxn)->recalc();
+            total_propensity += (*rxn)->propensity;
+        }
+
+        if (step >= max_steps)
+            break;
+
+        ++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';
+        }
+        std::cout << '\n';
+
+#if 0
+        // print copy numbers
+        std::cout << t << '\t';
+        for (std::vector<Species*>::iterator spec = species.begin(),
+                                        end = species.end();
+         spec != end; ++spec) {
+            std::cout << (*spec)->count << '\t';
+        }
+
+        MutatedSiteSequence ms;
+        FILE* output = stdout;
+        for (virus_map::iterator iter = s1.pop.begin(),
+                                 end = s1.pop.end();
+          iter != end; ++iter) {
+            std::cout << iter->second << '(';
+            ms = iter->first.mutated_sites;
+            // print sequence as CSV (comma-separated)
+            MutatedSiteSequence::iterator ms_iter=ms.begin();
+            if (ms_iter != ms.end()) {
+                fprintf(output,"%d",*ms_iter);
+                ++ms_iter;
+            }
+            for (; ms_iter!=ms.end(); ++ms_iter)
+                fprintf(output,",%d",*ms_iter);
+            std::cout << ')';
+            std::cout << '\t';
+        }
+        std::cout << '\n';
+#endif
+
+    }
+
+
+
+}
+
+
 // Import initial state from a state file
 
 void importState(RunParameters_SS &r) {
@@ -545,6 +662,7 @@ void importEpitope(RunParameters_SS &r) {
 
         std::string word;
         unsigned int site;
+        double aff;
         std::vector<unsigned int> tmpEp;
         size_t pos0 = 0;
         size_t posBar = readStr.find('|',pos0);
@@ -553,6 +671,13 @@ void importEpitope(RunParameters_SS &r) {
         // could use std::noskipws to keep tab ('\t') characters?
 
         std::stringstream readStrStrm(std::string(readStr,pos0,posBar-pos0));
+
+        // first entry is affinity
+        readStrStrm >> word;
+        aff = strtodouble(word);
+        printf("%.6e ",aff);
+
+        // next entries are WT sites
         while (readStrStrm >> word) {
             std::cout << word << " "; // XXX
             std::istringstream i(word);
@@ -669,7 +794,7 @@ 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.epitopefile=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; } }
         
         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 ad52952..747094a 100644 (file)
--- a/ss.h
+++ b/ss.h
@@ -10,6 +10,12 @@
 #include <iostream>
 #include <stdio.h>
 
+#include "reaction.h"
+#include "pop_ss.h"
+
+typedef std::vector<Reaction*> Rxn_parray;
+typedef std::vector<Species*> Species_parray;
+
 
 // STRING MANIPULATION
 
@@ -118,6 +124,10 @@ public:
     enum PenaltyTYPE { PenaltyEACH, PenaltyTOTAL };
 
     PenaltyTYPE penaltyType;
+
+    
+
+    std::vector<std::string> epitopefiles;
     
     
     RunParameters_SS() {
@@ -163,6 +173,10 @@ public:
             supplementaryOutfile=directory+"/"+outfile+".sum";
             stateInfile=directory+"/"+statefile+".st";
             epitopeInfile=directory+"/"+epitopefile+".ep";
+
+            for (size_t i = 0; i != epitopefiles.size(); ++i) {
+                epitopefiles[i] = directory+"/"+epitopefiles[i]+".ep";
+            }
             
         }
         
@@ -173,6 +187,10 @@ public:
             supplementaryOutfile=outfile+".sum";
             stateInfile=statefile+".st";
             epitopeInfile=epitopefile+".ep";
+
+            for (size_t i = 0; i != epitopefiles.size(); ++i) {
+                epitopefiles[i] = epitopefiles[i]+".ep";
+            }
         
         }
         
@@ -187,6 +205,7 @@ 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);
 
 
 #endif