Only add reactions with non-zero rate constants.

diff --git a/ss.cpp b/ss.cpp
index b602a7e..a9ee543 100644 (file)
--- a/ss.cpp
+++ b/ss.cpp
@@ -219,18 +219,22 @@ void initialize_Ising(RunParameters_SS &r, Species_parray &species, Rxn_parray &
     if (r.useEpitope)  importEpitope(r,species,reactions,s1,print_spec);    // <infile>.ep
 
     // V -> V + V'
-    VirusReaction* r1 = new VirusReaction(r.rate_s);  // s_k
-    r1->mu = r.mu;
-    r1->H = H;
-    r1->V = s1;
-    r1->name = "V --(s)--> V + V'";
-    reactions.push_back(r1);    update_rxn_dependents(reactions);
+    if (r.rate_s != 0.0) {
+        VirusReaction* r1 = new VirusReaction(r.rate_s);  // s_k
+        r1->mu = r.mu;
+        r1->H = H;
+        r1->V = s1;
+        r1->name = "V --(s)--> V + V'";
+        reactions.push_back(r1);    update_rxn_dependents(reactions);
+    }
     // V -> 0
-    VirusDeathReaction* r2 = new VirusDeathReaction(r.rate_u);  // u
-    r2->H = H;
-    r2->V = s1;
-    r2->name = "V --(u)--> 0";
-    reactions.push_back(r2);    update_rxn_dependents(reactions);
+    if (r.rate_u != 0.0) {
+        VirusDeathReaction* r2 = new VirusDeathReaction(r.rate_u);  // u
+        r2->H = H;
+        r2->V = s1;
+        r2->name = "V --(u)--> 0";
+        reactions.push_back(r2);    update_rxn_dependents(reactions);
+    }
 
 }
 
@@ -265,18 +269,22 @@ void initialize_Potts(RunParameters_SS &r, Species_parray &species, Rxn_parray &
     if (r.useEpitope)  importEpitope_Potts(r,species,reactions,s1,print_spec);    // <infile>.ep
 
     // V -> V + V'
-    NTVirusReaction* r1 = new NTVirusReaction(r.rate_s);  // s_k
-    r1->mu = r.mu;
-    r1->H = H;
-    r1->V = s1;
-    r1->name = "V --(s)--> V + V'";
-    reactions.push_back(r1);    update_rxn_dependents(reactions);
+    if (r.rate_s != 0.0) {
+        NTVirusReaction* r1 = new NTVirusReaction(r.rate_s);  // s_k
+        r1->mu = r.mu;
+        r1->H = H;
+        r1->V = s1;
+        r1->name = "V --(s)--> V + V'";
+        reactions.push_back(r1);    update_rxn_dependents(reactions);
+    }
     // V -> 0
-    NTVirusDeathReaction* r2 = new NTVirusDeathReaction(r.rate_u);  // u
-    r2->H = H;
-    r2->V = s1;
-    r2->name = "V --(u)--> 0";
-    reactions.push_back(r2);    update_rxn_dependents(reactions);
+    if (r.rate_u != 0.0) {
+        NTVirusDeathReaction* r2 = new NTVirusDeathReaction(r.rate_u);  // u
+        r2->H = H;
+        r2->V = s1;
+        r2->name = "V --(u)--> 0";
+        reactions.push_back(r2);    update_rxn_dependents(reactions);
+    }
 
 }
 
@@ -677,19 +685,19 @@ void importEpitope(RunParameters_SS &r, Species_parray &species, Rxn_parray &rea
 
         std::vector<CTLSpecies*>::iterator it = Tgen.begin(), end = Tgen.end();
 
-        for (it = Tgen.begin(); it != end; ++it) {   // p
+        if (r.rate_p != 0.0) for (it = Tgen.begin(); it != end; ++it) {   // p
             KillingReaction* rx = new KillingReaction(r.rate_p / r.vol);
             rx->T = *it;  rx->V = V;
             rx->name = (*it)->name + " + V --(p)--> " + (*it)->name;
             reactions.push_back(rx);    update_rxn_dependents(reactions);
         }
-        {   // a
+        if (r.rate_a != 0.0) {   // a
             CTLActivationReaction* rx = new CTLActivationReaction(r.rate_a / r.vol);
             rx->Tfrom = N;  rx->Tto = T;  rx->V = V;
             rx->name = "N + V --(a)--> T1 + V";
             reactions.push_back(rx);    update_rxn_dependents(reactions);
         }
-        for (size_t i=0; i<r.num_T_gen-1; i++) {   // r
+        if (r.rate_r != 0.0) for (size_t i=0; i<r.num_T_gen-1; i++) {   // r
             ElementaryReaction* rx = new ElementaryReaction(r.rate_r);
             rx->reactants.push_back(Tgen[i]); rx->reactant_stoich.push_back(1);
             rx->products.push_back(Tgen[i]); rx->product_stoich.push_back(-1);
@@ -697,7 +705,7 @@ void importEpitope(RunParameters_SS &r, Species_parray &species, Rxn_parray &rea
             rx->name = Tgen[i]->name + " --(r)--> 2 " + Tgen[i+1]->name;
             reactions.push_back(rx);    update_rxn_dependents(reactions);
         }
-        {   // r    // terminally differentiated effector cell just dies
+        if (r.rate_r != 0.0) {   // r    // terminally differentiated effector cell just dies
             ElementaryReaction* rx = new ElementaryReaction(r.rate_r);
             rx->reactants.push_back(Tgen.back()); rx->reactant_stoich.push_back(1);
             rx->products.push_back(Tgen.back()); rx->product_stoich.push_back(-1);
@@ -714,40 +722,40 @@ void importEpitope(RunParameters_SS &r, Species_parray &species, Rxn_parray &rea
         //    rx->products.push_back(Tgen.back()); rx->product_stoich.push_back(1);
         //    reactions.push_back(rx);  update_rxn_dependents(reactions);
         //}
-        for (it = Tgen.begin(); it != end; ++it) {   // d = 0.5  // note d' = 3.0
+        if (r.rate_d != 0.0) for (it = Tgen.begin(); it != end; ++it) {   // d = 0.5  // note d' = 3.0
             ElementaryReaction* rx = new ElementaryReaction(r.rate_d);
             rx->reactants.push_back(*it); rx->reactant_stoich.push_back(1);
             rx->products.push_back(*it); rx->product_stoich.push_back(-1);
             rx->name = (*it)->name + " --(d)--> 0";
             reactions.push_back(rx);    update_rxn_dependents(reactions);
         }
-        {   // b
+        if (r.rate_b != 0.0) {   // b
             ElementaryReaction* rx = new ElementaryReaction(r.rate_b);
             rx->reactants.push_back(N); rx->reactant_stoich.push_back(1);
             rx->products.push_back(N); rx->product_stoich.push_back(1);
             rx->name = "N --(b)--> 2 N";
             reactions.push_back(rx);    update_rxn_dependents(reactions);
         }
-        {   // e
+        if (r.rate_e != 0.0) {   // e
             ElementaryReaction* rx = new ElementaryReaction(r.rate_e);
             rx->reactants.push_back(N); rx->reactant_stoich.push_back(1);
             rx->products.push_back(N); rx->product_stoich.push_back(-1);
             rx->name = "N --(e)--> 0";
             reactions.push_back(rx);    update_rxn_dependents(reactions);
         }
-        {   // w
+        if (r.rate_w != 0.0) {   // w
             ElementaryReaction* rx = new ElementaryReaction(r.rate_w * r.vol);
             rx->products.push_back(N); rx->product_stoich.push_back(1);
             rx->name = "0 --(w)--> N";
             reactions.push_back(rx);    update_rxn_dependents(reactions);
         }
-        {   // a'
+        if (r.rate_aprime != 0.0) {   // a'
             CTLActivationReaction* rx = new CTLActivationReaction(r.rate_aprime / r.vol);
             rx->Tfrom = M;  rx->Tto = T;  rx->V = V;
             rx->name = "M + V --(a')--> T1 + V";
             reactions.push_back(rx);    update_rxn_dependents(reactions);
         }
-        for (it = Tgen.begin(); it != end; ++it) {   // g
+        if (r.rate_g != 0.0) for (it = Tgen.begin(); it != end; ++it) {   // g
             ElementaryReaction* rx = new ElementaryReaction(r.rate_g);
             rx->reactants.push_back(*it); rx->reactant_stoich.push_back(1);
             rx->products.push_back(*it); rx->product_stoich.push_back(-1);
@@ -755,7 +763,7 @@ void importEpitope(RunParameters_SS &r, Species_parray &species, Rxn_parray &rea
             rx->name = (*it)->name + " --(g)--> M";
             reactions.push_back(rx);    update_rxn_dependents(reactions);
         }
-        {   // h
+        if (r.rate_h != 0.0) {   // h
             ElementaryReaction* rx = new ElementaryReaction(r.rate_h);
             rx->reactants.push_back(M); rx->reactant_stoich.push_back(1);
             rx->products.push_back(M); rx->product_stoich.push_back(-1);
@@ -863,18 +871,18 @@ void importEpitope_Potts(RunParameters_SS &r, Species_parray &species, Rxn_parra
 
         std::vector<CTLaaSpecies*>::iterator it = Tgen.begin(), end = Tgen.end();
 
-        for (it = Tgen.begin(); it != end; ++it) {   // p
+        if (r.rate_p != 0.0) for (it = Tgen.begin(); it != end; ++it) {   // p
             AAKillingReaction* rx = new AAKillingReaction(r.rate_p / r.vol);
             rx->T = *it;  rx->V = V;
             reactions.push_back(rx);    update_rxn_dependents(reactions);
         }
-        {   // a
+        if (r.rate_a != 0.0) {   // a
             CTLaaActivationReaction* rx = new CTLaaActivationReaction(r.rate_a / r.vol);
             rx->Tfrom = N;  rx->Tto = T;  rx->V = V;
             rx->name = "N + V --(a)--> T1 + V";
             reactions.push_back(rx);    update_rxn_dependents(reactions);
         }
-        for (size_t i=0; i<r.num_T_gen-1; i++) {   // r
+        if (r.rate_r != 0.0) for (size_t i=0; i<r.num_T_gen-1; i++) {   // r
             ElementaryReaction* rx = new ElementaryReaction(r.rate_r);
             rx->reactants.push_back(Tgen[i]); rx->reactant_stoich.push_back(1);
             rx->products.push_back(Tgen[i]); rx->product_stoich.push_back(-1);
@@ -882,7 +890,7 @@ void importEpitope_Potts(RunParameters_SS &r, Species_parray &species, Rxn_parra
             rx->name = Tgen[i]->name + " --(r)--> 2 " + Tgen[i+1]->name;
             reactions.push_back(rx);    update_rxn_dependents(reactions);
         }
-        {   // r   // terminally differentiated effector cell just dies
+        if (r.rate_r != 0.0) {   // r   // terminally differentiated effector cell just dies
             ElementaryReaction* rx = new ElementaryReaction(r.rate_r);
             rx->reactants.push_back(Tgen.back()); rx->reactant_stoich.push_back(1);
             rx->products.push_back(Tgen.back()); rx->product_stoich.push_back(-1);
@@ -899,40 +907,40 @@ void importEpitope_Potts(RunParameters_SS &r, Species_parray &species, Rxn_parra
         //    rx->products.push_back(Tgen.back()); rx->product_stoich.push_back(1);
         //    reactions.push_back(rx);  update_rxn_dependents(reactions);
         //}
-        for (it = Tgen.begin(); it != end; ++it) {   // d = 0.5  // note d' = 3.0
+        if (r.rate_d != 0.0) for (it = Tgen.begin(); it != end; ++it) {   // d = 0.5  // note d' = 3.0
             ElementaryReaction* rx = new ElementaryReaction(r.rate_d);
             rx->reactants.push_back(*it); rx->reactant_stoich.push_back(1);
             rx->products.push_back(*it); rx->product_stoich.push_back(-1);
             rx->name = (*it)->name + " --(d)--> 0";
             reactions.push_back(rx);    update_rxn_dependents(reactions);
         }
-        {   // b
+        if (r.rate_b != 0.0) {   // b
             ElementaryReaction* rx = new ElementaryReaction(r.rate_b);
             rx->reactants.push_back(N); rx->reactant_stoich.push_back(1);
             rx->products.push_back(N); rx->product_stoich.push_back(1);
             rx->name = "N --(b)--> 2 N";
             reactions.push_back(rx);    update_rxn_dependents(reactions);
         }
-        {   // e
+        if (r.rate_e != 0.0) {   // e
             ElementaryReaction* rx = new ElementaryReaction(r.rate_e);
             rx->reactants.push_back(N); rx->reactant_stoich.push_back(1);
             rx->products.push_back(N); rx->product_stoich.push_back(-1);
             rx->name = "N --(e)--> 0";
             reactions.push_back(rx);    update_rxn_dependents(reactions);
         }
-        {   // w
+        if (r.rate_w != 0.0) {   // w
             ElementaryReaction* rx = new ElementaryReaction(r.rate_w * r.vol);
             rx->products.push_back(N); rx->product_stoich.push_back(1);
             rx->name = "0 --(w)--> N";
             reactions.push_back(rx);    update_rxn_dependents(reactions);
         }
-        {   // a'
+        if (r.rate_aprime != 0.0) {   // a'
             CTLaaActivationReaction* rx = new CTLaaActivationReaction(r.rate_aprime / r.vol);
             rx->Tfrom = M;  rx->Tto = T;  rx->V = V;
             rx->name = "M + V --(a')--> T1 + V";
             reactions.push_back(rx);    update_rxn_dependents(reactions);
         }
-        for (it = Tgen.begin(); it != end; ++it) {   // g
+        if (r.rate_g != 0.0) for (it = Tgen.begin(); it != end; ++it) {   // g
             ElementaryReaction* rx = new ElementaryReaction(r.rate_g);
             rx->reactants.push_back(*it); rx->reactant_stoich.push_back(1);
             rx->products.push_back(*it); rx->product_stoich.push_back(-1);
@@ -940,7 +948,7 @@ void importEpitope_Potts(RunParameters_SS &r, Species_parray &species, Rxn_parra
             rx->name = (*it)->name + " --(g)--> M";
             reactions.push_back(rx);    update_rxn_dependents(reactions);
         }
-        {   // h
+        if (r.rate_h != 0.0) {   // h
             ElementaryReaction* rx = new ElementaryReaction(r.rate_h);
             rx->reactants.push_back(M); rx->reactant_stoich.push_back(1);
             rx->products.push_back(M); rx->product_stoich.push_back(-1);