Handle dependencies between reactions, only recalculating propensity when affected...

diff --git a/reaction.h b/reaction.h
index 2e5e52c..d5eb1d8 100644 (file)
--- a/reaction.h
+++ b/reaction.h
@@ -7,13 +7,17 @@
 
 typedef std::set<unsigned int> MutatedSiteSequence;
 
+enum class RxnType {Elementary, Virus, VirusDeath, Killing, CTLActivation, NTVirus, NTVirusDeath, AAKilling, CTLaaActivation};
+
 class Reaction {
 public:
     Reaction() {}
-    Reaction(double k) : rate_constant(k) {}
+    Reaction(double k) : rate_constant(k), affects() {}
 
     double rate_constant;  // must manually ensure it includes division by statistical degeneracy factor (in case of reaction order >1 for any species)
     double propensity;
+    std::string name;
+    std::vector<Reaction*> affects; // reactions whose propensity must be updated when this one fires
     virtual double recalc() = 0;
     virtual void fire() = 0;
 };
@@ -150,5 +154,31 @@ public:
     void fire();
 };
 
+
+inline RxnType get_rxn_type(const Reaction* const rxn)
+{
+    if (dynamic_cast<const ElementaryReaction*>(rxn))
+        return RxnType::Elementary;
+    else if (dynamic_cast<const VirusReaction*>(rxn))
+        return RxnType::Virus;
+    else if (dynamic_cast<const VirusDeathReaction*>(rxn))
+        return RxnType::VirusDeath;
+    else if (dynamic_cast<const KillingReaction*>(rxn))
+        return RxnType::Killing;
+    else if (dynamic_cast<const CTLActivationReaction*>(rxn))
+        return RxnType::CTLActivation;
+    else if (dynamic_cast<const NTVirusReaction*>(rxn))
+        return RxnType::NTVirus;
+    else if (dynamic_cast<const NTVirusDeathReaction*>(rxn))
+        return RxnType::NTVirusDeath;
+    else if (dynamic_cast<const AAKillingReaction*>(rxn))
+        return RxnType::AAKilling;
+    else if (dynamic_cast<const CTLaaActivationReaction*>(rxn))
+        return RxnType::CTLaaActivation;
+    else
+        throw "unknown reaction type?!";
+}
+
+
 #endif  // REACTION_H
 

diff --git a/ss.cpp b/ss.cpp
index e10176d..d92bd58 100644 (file)
--- a/ss.cpp
+++ b/ss.cpp
@@ -39,6 +39,182 @@ static unsigned sim_random_seed() {
 }
 
 
+bool does_species_appear_in_reactants(Species *s, Reaction *rxn)
+{
+    bool appears = false;
+
+    RxnType T = get_rxn_type(rxn);
+    switch (T) {
+    case RxnType::Elementary:
+        for (Species_parray::iterator t = dynamic_cast<ElementaryReaction*>(rxn)->reactants.begin(),
+                                      f = dynamic_cast<ElementaryReaction*>(rxn)->reactants.end();
+                                      t != f; ++t) {
+            if (s == *t){std::cout << "elem: " << rxn->name << '\n'; return true;}
+        }
+        break;
+    case RxnType::Virus:
+        if (s == dynamic_cast<VirusReaction*>(rxn)->V){std::cout << "virus: " << rxn->name << '\n'; return true;}
+        break;
+    case RxnType::VirusDeath:
+        if (s == dynamic_cast<VirusDeathReaction*>(rxn)->V){std::cout << "vdeath: " << rxn->name << '\n'; return true;}
+        break;
+    case RxnType::Killing:
+        if (s == dynamic_cast<KillingReaction*>(rxn)->V || s == dynamic_cast<KillingReaction*>(rxn)->T){std::cout << "killing: " << rxn->name << '\n'; return true;}
+        break;
+    case RxnType::CTLActivation:
+        if (s == dynamic_cast<CTLActivationReaction*>(rxn)->V || s == dynamic_cast<CTLActivationReaction*>(rxn)->Tfrom){std::cout << "ctlact: " << rxn->name << '\n'; return true;}
+        break;
+    case RxnType::NTVirus:
+        if (s == dynamic_cast<NTVirusReaction*>(rxn)->V){std::cout << "ntvirus: " << rxn->name << '\n'; return true;}
+        break;
+    case RxnType::NTVirusDeath:
+        if (s == dynamic_cast<NTVirusDeathReaction*>(rxn)->V){std::cout << "ntvdeath: " << rxn->name << '\n'; return true;}
+        break;
+    case RxnType::AAKilling:
+        if (s == dynamic_cast<AAKillingReaction*>(rxn)->V || s == dynamic_cast<AAKillingReaction*>(rxn)->T){std::cout << "aakilling: " << rxn->name << '\n'; return true;}
+        break;
+    case RxnType::CTLaaActivation:
+        if (s == dynamic_cast<CTLaaActivationReaction*>(rxn)->V || s == dynamic_cast<CTLaaActivationReaction*>(rxn)->Tfrom){std::cout << "ctlaaact: " << rxn->name << '\n'; return true;}
+        break;
+    }
+
+    return appears;
+}
+
+
+// does newrxn affect rxn?
+// i.e. if newrxn fires, will rxn's propensity change?
+bool does_affect(Reaction *newrxn, Reaction *rxn)
+{
+    bool affects = false;
+
+    RxnType newT = get_rxn_type(newrxn);
+
+    switch (newT) {
+
+    case RxnType::Elementary: {
+        // If any of newrxn's "products" (species whose count changes upon firing)
+        // are in rxn's "reactants" (species needed to calculate rate),
+        // then newrxn does indeed affect rxn.
+        for (Species_parray::iterator s = dynamic_cast<ElementaryReaction*>(newrxn)->products.begin(),
+                                      e = dynamic_cast<ElementaryReaction*>(newrxn)->products.end();
+                                      s != e; ++s) {
+            std::cout << "newrxn = Elementary: " << newrxn->name << '\n';
+            affects |= does_species_appear_in_reactants(*s,rxn);
+        }
+        break;
+    }
+
+    case RxnType::Virus: {
+        std::cout << "newrxn = Virus: " << newrxn->name << '\n';
+        Species *s = dynamic_cast<VirusReaction*>(newrxn)->V;
+        affects = does_species_appear_in_reactants(s,rxn);
+        return affects;
+        break;
+    }
+
+    case RxnType::VirusDeath: {
+        std::cout << "newrxn = VirusDeath: " << newrxn->name << '\n';
+        Species *s = dynamic_cast<VirusDeathReaction*>(newrxn)->V;
+        affects = does_species_appear_in_reactants(s,rxn);
+        return affects;
+        break;
+    }
+
+    case RxnType::Killing: {
+        std::cout << "newrxn = Killing: " << newrxn->name << '\n';
+        Species *s = dynamic_cast<KillingReaction*>(newrxn)->V;
+        affects = does_species_appear_in_reactants(s,rxn);
+        //return affects;
+        Species *s2 = dynamic_cast<KillingReaction*>(newrxn)->T;
+        affects |= does_species_appear_in_reactants(s2,rxn);
+        return affects;
+        break;
+    }
+
+    case RxnType::CTLActivation: {
+        std::cout << "newrxn = CTLActivation: " << newrxn->name << '\n';
+        Species *s;
+        //s = dynamic_cast<CTLActivationReaction*>(newrxn)->V;
+        //affects = does_species_appear_in_reactants(s,rxn);
+        //return affects;
+        s = dynamic_cast<CTLActivationReaction*>(newrxn)->Tfrom;
+        affects = does_species_appear_in_reactants(s,rxn);
+        //return affects;
+        s = dynamic_cast<CTLActivationReaction*>(newrxn)->Tto;
+        affects |= does_species_appear_in_reactants(s,rxn);
+        return affects;
+        break;
+    }
+
+    case RxnType::NTVirus: {
+        std::cout << "newrxn = NTVirus: " << newrxn->name << '\n';
+        Species *s = dynamic_cast<NTVirusReaction*>(newrxn)->V;
+        affects = does_species_appear_in_reactants(s,rxn);
+        return affects;
+        break;
+    }
+
+    case RxnType::NTVirusDeath: {
+        std::cout << "newrxn = NTVirusDeath: " << newrxn->name << '\n';
+        Species *s = dynamic_cast<NTVirusDeathReaction*>(newrxn)->V;
+        affects = does_species_appear_in_reactants(s,rxn);
+        break;
+    }
+
+    case RxnType::AAKilling: {
+        std::cout << "newrxn = AAKilling: " << newrxn->name << '\n';
+        Species *s = dynamic_cast<AAKillingReaction*>(newrxn)->V;
+        affects = does_species_appear_in_reactants(s,rxn);
+        return affects;
+        //Species *s2 = dynamic_cast<AAKillingReaction*>(newrxn)->T;
+        //affects = does_species_appear_in_reactants(s2,rxn);
+        //return affects;
+        break;
+    }
+
+    case RxnType::CTLaaActivation: {
+        std::cout << "newrxn = CTLaaActivation: " << newrxn->name << '\n';
+        Species *s = dynamic_cast<CTLaaActivationReaction*>(newrxn)->V;
+        affects = does_species_appear_in_reactants(s,rxn);
+        //return affects;
+        Species *s2 = dynamic_cast<CTLaaActivationReaction*>(newrxn)->Tfrom;
+        affects |= does_species_appear_in_reactants(s2,rxn);
+        return affects;
+        break;
+    }
+
+
+    } // switch (newT)
+
+    return affects;
+}
+
+
+// call after any reaction is appended to the reaction array
+// in order to update the existing Reaction objects' "affects" member
+void update_rxn_dependents(Rxn_parray &reactions)
+{
+    printf("----------\n");
+    Reaction *newrxn = reactions.back();
+
+    Rxn_parray::iterator rxn = reactions.begin(),
+                         end = std::prev(reactions.end());    // don't include self
+
+    for (; rxn != end; ++rxn) {
+        if (does_affect(newrxn,*rxn)) { // newrxn affects old?
+            printf("new affects old\n\n");
+            newrxn->affects.push_back(*rxn);
+        }
+        if (does_affect(*rxn,newrxn)) { // old reaction affects new?
+            printf("old affects new\n\n");
+            (*rxn)->affects.push_back(newrxn);
+        }
+    }
+}
+
+
+
 void initialize_Ising(RunParameters_SS &r, Species_parray &species, Rxn_parray &reactions, Species_parray &print_spec)
 {
 
@@ -71,15 +247,18 @@ void initialize_Ising(RunParameters_SS &r, Species_parray &species, Rxn_parray &
     r1->mu = r.mu;
     r1->H = H;
     r1->V = s1;
-    reactions.push_back(r1);
+    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;
-    reactions.push_back(r2);
+    r2->name = "V --(u)--> 0";
+    reactions.push_back(r2);    update_rxn_dependents(reactions);
 
 }
 
+
 void initialize_Potts(RunParameters_SS &r, Species_parray &species, Rxn_parray &reactions, Species_parray &print_spec)
 {
 
@@ -114,12 +293,14 @@ void initialize_Potts(RunParameters_SS &r, Species_parray &species, Rxn_parray &
     r1->mu = r.mu;
     r1->H = H;
     r1->V = s1;
-    reactions.push_back(r1);
+    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;
-    reactions.push_back(r2);
+    r2->name = "V --(u)--> 0";
+    reactions.push_back(r2);    update_rxn_dependents(reactions);
 
 }
 
@@ -260,12 +441,15 @@ void simulate(Rxn_parray &reactions, Species_parray &species, double t_end, long
     print_species_traj(t,step,print_spec,verbose);   // XXX
 
     Rxn_parray::iterator rxn, end;
+    Rxn_parray::iterator rxn2, end2;
 
     for (rxn = reactions.begin(), end = reactions.end(); rxn != end; ++rxn) {
         (*rxn)->recalc();
         total_propensity += (*rxn)->propensity;
     }
 
+    std::cout << "total propensity = " << total_propensity << '\n';
+
 
     // stop if time runs out or too many steps or no more reactions will occur
     for (; (max_steps == 0 || step < max_steps) && (t < t_end) && (total_propensity > 0.0); ++step) {
@@ -291,10 +475,27 @@ void simulate(Rxn_parray &reactions, Species_parray &species, double t_end, long
         if (rxn == end)
             break;
 
+        (*rxn)->recalc();
+        std::cout << "fired: " << (*rxn)->name << '\n';
+        for (rxn2 = (*rxn)->affects.begin(), end2 = (*rxn)->affects.end(); rxn2 != end2; ++rxn2) {
+            std::cout << "affects: " << (*rxn2)->name << '\n';
+            double oldA = (*rxn2)->propensity;
+            //std::cout << step << "\thit it\n";
+            //std::cout << step << "\told = " << (*rxn2)->propensity;
+            printf("%ld\told = %.4f",step,oldA);
+            (*rxn2)->recalc();
+            //std::cout << "\tnew = " << (*rxn2)->propensity << '\n';
+            double newA = (*rxn2)->propensity;
+            printf("\tnew = %.4f",newA);
+            if (oldA == newA)
+                printf("\tSAME");
+            printf("\n");
+        }
+
         // recalculate rates
         total_propensity = 0.0;
         for (rxn = reactions.begin(), end = reactions.end(); rxn != end; ++rxn) {
-            (*rxn)->recalc();
+            //(*rxn)->recalc();     // dependencies handled above
             total_propensity += (*rxn)->propensity;
         }
 
@@ -307,6 +508,7 @@ void simulate(Rxn_parray &reactions, Species_parray &species, double t_end, long
         if ((step_sample_interval != 0) && (step_next_sample <= step)) {
             step_next_sample += step_sample_interval;
             print_species_traj(t,step,print_spec,verbose);
+            std::cout << "total propensity = " << total_propensity << '\n';
         }
 
     }
@@ -507,25 +709,29 @@ void importEpitope(RunParameters_SS &r, Species_parray &species, Rxn_parray &rea
         for (it = Tgen.begin(); it != end; ++it) {   // p
             KillingReaction* rx = new KillingReaction(r.rate_p);
             rx->T = *it;  rx->V = V;
-            reactions.push_back(rx);
+            rx->name = (*it)->name + " + V --(p)--> " + (*it)->name;
+            reactions.push_back(rx);    update_rxn_dependents(reactions);
         }
         {   // a
             CTLActivationReaction* rx = new CTLActivationReaction(r.rate_a);
             rx->Tfrom = N;  rx->Tto = T;  rx->V = V;
-            reactions.push_back(rx);
+            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
             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);
             rx->products.push_back(Tgen[i+1]); rx->product_stoich.push_back(2);
-            reactions.push_back(rx);
+            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
             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);
-            reactions.push_back(rx);
+            rx->name = Tgen.back()->name + " --(r)--> 0";
+            reactions.push_back(rx);    update_rxn_dependents(reactions);
         }
         // un-commenting the following reaction, and removing the above reaction,
         // causes there to be no difference between the effector cell generations,
@@ -535,48 +741,55 @@ void importEpitope(RunParameters_SS &r, Species_parray &species, Rxn_parray &rea
         //    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);
-        //    reactions.push_back(rx);
+        //    reactions.push_back(rx);  update_rxn_dependents(reactions);
         //}
         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);
-            reactions.push_back(rx);
+            rx->name = (*it)->name + " --(d)--> 0";
+            reactions.push_back(rx);    update_rxn_dependents(reactions);
         }
         {   // 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);
-            reactions.push_back(rx);
+            rx->name = "N --(b)--> 2 N";
+            reactions.push_back(rx);    update_rxn_dependents(reactions);
         }
         {   // 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);
-            reactions.push_back(rx);
+            rx->name = "N --(e)--> 0";
+            reactions.push_back(rx);    update_rxn_dependents(reactions);
         }
         {   // w
             ElementaryReaction* rx = new ElementaryReaction(r.rate_w);
             rx->products.push_back(N); rx->product_stoich.push_back(1);
-            reactions.push_back(rx);
+            rx->name = "0 --(w)--> N";
+            reactions.push_back(rx);    update_rxn_dependents(reactions);
         }
         {   // a'
             CTLActivationReaction* rx = new CTLActivationReaction(r.rate_aprime);
             rx->Tfrom = M;  rx->Tto = T;  rx->V = V;
-            reactions.push_back(rx);
+            rx->name = "M + V --(a')--> T1 + V";
+            reactions.push_back(rx);    update_rxn_dependents(reactions);
         }
         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);
             rx->products.push_back(M); rx->product_stoich.push_back(1);
-            reactions.push_back(rx);
+            rx->name = (*it)->name + " --(g)--> M";
+            reactions.push_back(rx);    update_rxn_dependents(reactions);
         }
         {   // 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);
-            reactions.push_back(rx);
+            rx->name = "M --(h)--> 0";
+            reactions.push_back(rx);    update_rxn_dependents(reactions);
         }
 
         for (size_t i = 0; i < T->num_ep; ++i) {  // XXX
@@ -682,25 +895,28 @@ void importEpitope_Potts(RunParameters_SS &r, Species_parray &species, Rxn_parra
         for (it = Tgen.begin(); it != end; ++it) {   // p
             AAKillingReaction* rx = new AAKillingReaction(r.rate_p);
             rx->T = *it;  rx->V = V;
-            reactions.push_back(rx);
+            reactions.push_back(rx);    update_rxn_dependents(reactions);
         }
         {   // a
             CTLaaActivationReaction* rx = new CTLaaActivationReaction(r.rate_a);
             rx->Tfrom = N;  rx->Tto = T;  rx->V = V;
-            reactions.push_back(rx);
+            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
             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);
             rx->products.push_back(Tgen[i+1]); rx->product_stoich.push_back(2);
-            reactions.push_back(rx);
+            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
             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);
-            reactions.push_back(rx);
+            rx->name = Tgen.back()->name + " --(r)--> 0";
+            reactions.push_back(rx);    update_rxn_dependents(reactions);
         }
         // un-commenting the following reaction, and removing the above reaction,
         // causes there to be no difference between the effector cell generations,
@@ -710,48 +926,55 @@ void importEpitope_Potts(RunParameters_SS &r, Species_parray &species, Rxn_parra
         //    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);
-        //    reactions.push_back(rx);
+        //    reactions.push_back(rx);  update_rxn_dependents(reactions);
         //}
         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);
-            reactions.push_back(rx);
+            rx->name = (*it)->name + " --(d)--> 0";
+            reactions.push_back(rx);    update_rxn_dependents(reactions);
         }
         {   // 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);
-            reactions.push_back(rx);
+            rx->name = "N --(b)--> 2 N";
+            reactions.push_back(rx);    update_rxn_dependents(reactions);
         }
         {   // 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);
-            reactions.push_back(rx);
+            rx->name = "N --(e)--> 0";
+            reactions.push_back(rx);    update_rxn_dependents(reactions);
         }
         {   // w
             ElementaryReaction* rx = new ElementaryReaction(r.rate_w);
             rx->products.push_back(N); rx->product_stoich.push_back(1);
-            reactions.push_back(rx);
+            rx->name = "0 --(w)--> N";
+            reactions.push_back(rx);    update_rxn_dependents(reactions);
         }
         {   // a'
             CTLaaActivationReaction* rx = new CTLaaActivationReaction(r.rate_aprime);
             rx->Tfrom = M;  rx->Tto = T;  rx->V = V;
-            reactions.push_back(rx);
+            rx->name = "M + V --(a')--> T1 + V";
+            reactions.push_back(rx);    update_rxn_dependents(reactions);
         }
         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);
             rx->products.push_back(M); rx->product_stoich.push_back(1);
-            reactions.push_back(rx);
+            rx->name = (*it)->name + " --(g)--> M";
+            reactions.push_back(rx);    update_rxn_dependents(reactions);
         }
         {   // 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);
-            reactions.push_back(rx);
+            rx->name = "M --(h)--> 0";
+            reactions.push_back(rx);    update_rxn_dependents(reactions);
         }
 
         for (size_t e=0; e<M->num_ep; e++) {    // XXX