From 25b99d8c99852ef63af1d0f4a14f8f056fb97211 Mon Sep 17 00:00:00 2001
From: Dariusz Murakowski <murakdar@mit.edu>
Date: Fri, 12 Jun 2015 17:26:30 -0400
Subject: [PATCH] Rescale rates by system volume ratio (parameter).

---
 ss.cpp | 19 +++++++++++--------
 ss.h   |  4 ++++
 2 files changed, 15 insertions(+), 8 deletions(-)

diff --git a/ss.cpp b/ss.cpp
index 4a8e438..b602a7e 100644
--- a/ss.cpp
+++ b/ss.cpp
@@ -678,13 +678,13 @@ 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
-            KillingReaction* rx = new KillingReaction(r.rate_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
-            CTLActivationReaction* rx = new CTLActivationReaction(r.rate_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);
@@ -736,13 +736,13 @@ void importEpitope(RunParameters_SS &r, Species_parray &species, Rxn_parray &rea
             reactions.push_back(rx);    update_rxn_dependents(reactions);
         }
         {   // w
-            ElementaryReaction* rx = new ElementaryReaction(r.rate_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'
-            CTLActivationReaction* rx = new CTLActivationReaction(r.rate_aprime);
+            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);
@@ -864,12 +864,12 @@ 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
-            AAKillingReaction* rx = new AAKillingReaction(r.rate_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
-            CTLaaActivationReaction* rx = new CTLaaActivationReaction(r.rate_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);
@@ -921,13 +921,13 @@ void importEpitope_Potts(RunParameters_SS &r, Species_parray &species, Rxn_parra
             reactions.push_back(rx);    update_rxn_dependents(reactions);
         }
         {   // w
-            ElementaryReaction* rx = new ElementaryReaction(r.rate_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'
-            CTLaaActivationReaction* rx = new CTLaaActivationReaction(r.rate_aprime);
+            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);
@@ -991,6 +991,7 @@ void usage()
 " -numruns (int)      number of trajectories to simulate\n"
 " -rate_X (double)    X in [s,u,p,a,r,d,dprime,b,e,w,aprime,g,h]\n"
 " -const_ND (int)     number of effector CTL divisions (generations) [default=9]\n"
+" -vol (double)       volume scaling factor; affects rate constants, not initial values [default=1.0]\n"
 ;   std::cout.flush();
 }
 
@@ -1066,6 +1067,8 @@ int main(int argc, char *argv[]) {
 
         else if (strcmp(argv[i],"-const_ND")==0) { if (++i==argc) break; else r.num_T_gen=(unsigned int)strtodouble(argv[i]); }
 
+        else if (strcmp(argv[i],"-vol")==0) { if (++i==argc) break; else r.vol=strtodouble(argv[i]); }
+
         else if (strcmp(argv[i],"-h")==0 || strcmp(argv[i],"--help")==0)   { usage(); return 0;    }
 
         else printf("Unrecognized argument! '%s'\n",argv[i]);
diff --git a/ss.h b/ss.h
index 3a8f93c..d84e542 100644
--- a/ss.h
+++ b/ss.h
@@ -175,6 +175,8 @@ public:
     double rate_g;
     double rate_h;
 
+    double vol;
+
     unsigned int num_T_gen;    // ND, number of effector cell generations (i.e. divisions)
     
     bool usePotts;
@@ -218,6 +220,8 @@ public:
         rate_g = 0.03;
         rate_h = 0.03;
 
+        vol = 1.0;  // corresponds to real volume basis of XXX
+
         num_T_gen = 9;
 
         usePotts = false;
-- 
2.7.4