Full system, including parameters.

diff --git a/reaction.cpp b/reaction.cpp
index c4ac28a..3ee908f 100644 (file)
--- a/reaction.cpp
+++ b/reaction.cpp
@@ -14,7 +14,7 @@ double ElementaryReaction::recalc()
     //double denom = 1.0;
     this->propensity = this->rate_constant;
 
-    std::vector<SimpleSpecies*>::iterator reac = this->reactants.begin(),
+    std::vector<Species*>::iterator reac = this->reactants.begin(),
                                           reac_end = this->reactants.end();
     std::vector<int>::iterator reacN = this->reactant_stoich.begin(),
                                reacN_end = this->reactant_stoich.end();
@@ -33,7 +33,7 @@ double ElementaryReaction::recalc()
 
 void ElementaryReaction::fire()
 {
-    std::vector<SimpleSpecies*>::iterator reac = this->reactants.begin(),
+    std::vector<Species*>::iterator reac = this->reactants.begin(),
                                           reac_end = this->reactants.end();
     std::vector<int>::iterator reacN = this->reactant_stoich.begin(),
                                reacN_end = this->reactant_stoich.end();
@@ -46,7 +46,7 @@ void ElementaryReaction::fire()
     }
     */
 
-    std::vector<SimpleSpecies*>::iterator prod = this->products.begin(),
+    std::vector<Species*>::iterator prod = this->products.begin(),
                                           prod_end = this->products.end();
     std::vector<int>::iterator prodN = this->product_stoich.begin(),
                                prodN_end = this->product_stoich.end();
@@ -98,7 +98,8 @@ void VirusReaction::fire()
     // pick how many sites to mutate in the new virus
     // adapted from WF Virus::mutate()
     //double mu = 6.0e-5;  // mut prob upon reproduction
-    double mu = 1.0e-1;
+    //double mu = 1.0e-1;
+    double mu = iter->first.mu;
     Virus v(H,mu,iter->first.mutated_sites);
     unsigned int n = gsl_ran_binomial(rnd,mu,H.size);
     MutatedSiteSequence sites_to_mutate;

diff --git a/reaction.h b/reaction.h
index 3ea501a..2f1ab57 100644 (file)
--- a/reaction.h
+++ b/reaction.h
@@ -20,13 +20,14 @@ public:
 
 // a A + b B + ... -> c C + d D + ...
 // mass action kinetics
+// the Species element should really be a SimpleSpecies, but really anything other than VirusSpecies is okay
 class ElementaryReaction : public Reaction {
 public:
     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<Species*> reactants;  // anything that participates in reaction rate
+    std::vector<Species*> products;  // effects of reaction
     std::vector<int> reactant_stoich;   // reaction order
     std::vector<int> product_stoich;  // must be negative for reactants that are consumed!
     double recalc();

diff --git a/ss.cpp b/ss.cpp
index 96841bb..858e4af 100644 (file)
--- a/ss.cpp
+++ b/ss.cpp
@@ -59,34 +59,36 @@ void run(RunParameters_SS &r, unsigned seed) {
 
     std::string couplingsFile = "neutral_2site.j";
     Hamiltonian H(couplingsFile);
+    H.set_temp(0.01);
     //double mu = 6.0e-5;
     //double mu = 1.0e-1;
 
     // initialize virus
     VirusSpecies s1; s1.count = r.n;
-    s1.pop[Virus(H,r.mu)] = s1.count;
+    s1.pop[Virus(H,r.mu)] = s1.count;   // default mu = 6.0e-5
     species.push_back(&s1);
 
+    Species_parray print_spec;
+    print_spec.push_back(&s1);
+
     if (r.importState) importState(r);      // <infile>.st
-    if (r.useEpitope)  importEpitope(r,species,reactions,&s1);    // <infile>.ep
+    if (r.useEpitope)  importEpitope(r,species,reactions,&s1,print_spec);    // <infile>.ep
 
     // V -> V + V'
-    VirusReaction* r1 = new VirusReaction();
-    r1->rate_constant = 3.0;
+    VirusReaction* r1 = new VirusReaction(1.5);  // s_k
+    //r1->rate_constant = 1.5;
     r1->H = H;
     r1->V = &s1;
     reactions.push_back(r1);
     // V -> 0
-    VirusDeathReaction* r2 = new VirusDeathReaction();
-    r2->rate_constant = 1.5;
+    VirusDeathReaction* r2 = new VirusDeathReaction(0.5);  // u
+    //r2->rate_constant = 0.5;
     r2->H = H;
     r2->V = &s1;
     reactions.push_back(r2);
 
-    Species_parray print_spec;
-    print_spec.push_back(&s1);
 
-    double t_end = 10.0;
+    double t_end = 1000.0;
     long max_steps = LONG_MAX;
     double sample_interval = 1.0;
     simulate(reactions, species, t_end, max_steps, sample_interval, print_spec);
@@ -665,7 +667,7 @@ void importState(RunParameters_SS &r) {
 
 // load epitope definitions from file
 
-void importEpitope(RunParameters_SS &r, Species_parray &species, Rxn_parray &reactions, VirusSpecies *V)
+void importEpitope(RunParameters_SS &r, Species_parray &species, Rxn_parray &reactions, VirusSpecies *V, Species_parray &print_spec)
 {
 
     if ((r.num_CTL_clones != r.epitopefiles.size()) || 
@@ -742,24 +744,73 @@ void importEpitope(RunParameters_SS &r, Species_parray &species, Rxn_parray &rea
 
         species.push_back(T);   // effector
 
+        CTLSpecies* N = new CTLSpecies(*T);     // naive
+        species.push_back(N);
+
         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);
+        print_spec.push_back(T);
+        print_spec.push_back(N);
+        print_spec.push_back(M);
 
+        {   // p
+            KillingReaction* r = new KillingReaction(1e-5);
+            r->T = T;  r->V = V;
+            reactions.push_back(r);
+        }
+        {   // a
+            CTLActivationReaction* r = new CTLActivationReaction(1e-7);
+            r->Tfrom = N;  r->Tto = T;  r->V = V;
+            reactions.push_back(r);
+        }
+        {   // r
+            ElementaryReaction* r = new ElementaryReaction(4.0);
+            r->reactants.push_back(T); r->reactant_stoich.push_back(1);
+            r->products.push_back(T); r->product_stoich.push_back(1);
+            reactions.push_back(r);
+        }
+        {   // d = 0.5  // note d' = 3.0
+            ElementaryReaction* r = new ElementaryReaction(3.0);
+            r->reactants.push_back(T); r->reactant_stoich.push_back(1);
+            r->products.push_back(T); r->product_stoich.push_back(-1);
+            reactions.push_back(r);
+        }
+        {   // b
+            ElementaryReaction* r = new ElementaryReaction(1e-4);
+            r->reactants.push_back(N); r->reactant_stoich.push_back(1);
+            r->products.push_back(N); r->product_stoich.push_back(1);
+            reactions.push_back(r);
+        }
+        {   // e
+            ElementaryReaction* r = new ElementaryReaction(3e-4);
+            r->reactants.push_back(N); r->reactant_stoich.push_back(1);
+            r->products.push_back(N); r->product_stoich.push_back(-1);
+            reactions.push_back(r);
+        }
+        {   // w
+            ElementaryReaction* r = new ElementaryReaction(7.5);
+            r->products.push_back(N); r->product_stoich.push_back(1);
+            reactions.push_back(r);
+        }
+        {   // a'
+            CTLActivationReaction* r = new CTLActivationReaction(1e-6);
+            r->Tfrom = M;  r->Tto = T;  r->V = V;
+            reactions.push_back(r);
+        }
+        {   // g
+            ElementaryReaction* r = new ElementaryReaction(0.03);
+            r->reactants.push_back(T); r->reactant_stoich.push_back(1);
+            r->products.push_back(T); r->product_stoich.push_back(-1);
+            r->products.push_back(M); r->product_stoich.push_back(1);
+            reactions.push_back(r);
+        }
+        {   // h
+            ElementaryReaction* r = new ElementaryReaction(0.03);
+            r->reactants.push_back(M); r->reactant_stoich.push_back(1);
+            r->products.push_back(M); r->product_stoich.push_back(-1);
+            reactions.push_back(r);
+        }
 
         for (size_t i = 0; i < T->num_ep; ++i) {  // XXX
             std::cout << T->count << '\t'

diff --git a/ss.h b/ss.h
index be805c4..2c26221 100644 (file)
--- a/ss.h
+++ b/ss.h
@@ -231,7 +231,7 @@ public:
 
 void run(RunParameters_SS &r);
 void importState(RunParameters_SS &r);
-void importEpitope(RunParameters_SS &r, Species_parray &species, Rxn_parray &reactions, VirusSpecies *V);
+void importEpitope(RunParameters_SS &r, Species_parray &species, Rxn_parray &reactions, VirusSpecies *V, 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);
 
 

diff --git a/virus.cpp b/virus.cpp
index fd69023..f523ca4 100644 (file)
--- a/virus.cpp
+++ b/virus.cpp
@@ -101,7 +101,3 @@ bool operator<(const Virus& lhs, const Virus& rhs) {
 
 }
 
-
-
-       
-