From c18beda2cadf901d86c41c61e03dce8124d5b8bb Mon Sep 17 00:00:00 2001
From: Dariusz Murakowski <murakdar@mit.edu>
Date: Tue, 21 Apr 2015 14:24:40 -0400
Subject: [PATCH] Virus reproduction reaction stores its own mutation rate.

---
 reaction.cpp | 11 +++++------
 reaction.h   |  4 +++-
 ss.cpp       |  1 +
 3 files changed, 9 insertions(+), 7 deletions(-)

diff --git a/reaction.cpp b/reaction.cpp
index 3ee908f..9d3ad2a 100644
--- a/reaction.cpp
+++ b/reaction.cpp
@@ -57,6 +57,7 @@ void ElementaryReaction::fire()
     }
 }
 
+
 double RC_fun(double E)
 {
     // sigmoid function
@@ -97,17 +98,15 @@ 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 = iter->first.mu;
-    Virus v(H,mu,iter->first.mutated_sites);
+    //double mu = iter->first.mu;
+    Virus v(H,this->mu,iter->first.mutated_sites);
     unsigned int n = gsl_ran_binomial(rnd,mu,H.size);
     MutatedSiteSequence sites_to_mutate;
     while (sites_to_mutate.size() < n) {
         sites_to_mutate.insert(gsl_rng_uniform_int(rnd,H.size));
     }
     for (MutatedSiteSequence::iterator it  = sites_to_mutate.begin(),
-                                          end = sites_to_mutate.end();
+                                       end = sites_to_mutate.end();
       it != end; ++it) {
         if (v.mutated_sites.count(*it)==0) v.mutated_sites.insert(*it);
         else                               v.mutated_sites.erase(*it);
@@ -129,6 +128,7 @@ void VirusReaction::fire()
 
 }
 
+
 double VirusDeathReaction::recalc()
 {
     V->count = 0;
@@ -228,4 +228,3 @@ void CTLActivationReaction::fire()
     Tto->count += 1;
 }
 
-
diff --git a/reaction.h b/reaction.h
index 2f1ab57..b74051a 100644
--- a/reaction.h
+++ b/reaction.h
@@ -40,10 +40,12 @@ public:
 class VirusReaction : public Reaction {
 public:
     VirusReaction() : Reaction() {}
-    VirusReaction(double k) : Reaction(k) {}
+    VirusReaction(double k) : Reaction(k), mu(0) {}
 
+    double mu;
     VirusSpecies* V;
     Hamiltonian H;
+
     double recalc();
     void fire();
 };
diff --git a/ss.cpp b/ss.cpp
index f1487a6..ba25e86 100644
--- a/ss.cpp
+++ b/ss.cpp
@@ -96,6 +96,7 @@ void run(RunParameters_SS &r, unsigned seed) {
     // V -> V + V'
     VirusReaction* r1 = new VirusReaction(r.rate_s);  // s_k
     //r1->rate_constant = 1.5;
+    r1->mu = r.mu;
     r1->H = H;
     r1->V = s1;
     reactions.push_back(r1);
-- 
2.7.4