From f4f552ccd4cdac5bd61515962f2576b39c9021e5 Mon Sep 17 00:00:00 2001 From: Dariusz Murakowski Date: Thu, 9 Apr 2015 01:10:37 -0400 Subject: [PATCH] Flesh out reactions for virus reproduction (with mutation) and death. --- reaction.cpp | 115 ++++++++++++++++++++++++++++++++++++++++++++++++++- reaction.h | 10 +++++ ss.cpp | 132 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++--- virus.h | 2 +- 4 files changed, 252 insertions(+), 7 deletions(-) diff --git a/reaction.cpp b/reaction.cpp index 077b33c..37c34e5 100644 --- a/reaction.cpp +++ b/reaction.cpp @@ -1,6 +1,12 @@ +#include +#include +#include +#include #include "reaction.h" +extern gsl_rng *rnd; //pointer to random number generator state + //static long factorial[] = {1, 1, 2, 6, 24, 120, 720}; double ElementaryReaction::recalc() @@ -49,16 +55,123 @@ void ElementaryReaction::fire() } } +double RC_fun(double E) +{ + // sigmoid function + //double z = exp(-E); + //return z/(1.0+z); + return (E>0.0) ? 1.0/(1.0+exp(-E)) : exp(E)/(1.0+exp(E)); +} double VirusReaction::recalc() { - // sum stuff + V->count = 0; + this->propensity = this->rate_constant; + + double partial_sum = 0.0; + virus_map::iterator iter = this->V->pop.begin(), + end = this->V->pop.end(); + for(; iter != end; ++iter) { + V->count += iter->second; + partial_sum += RC_fun(iter->first.energy) * iter->second; + } + + this->propensity *= partial_sum; return this->propensity; } void VirusReaction::fire() { + // pick which virus reproduces + double rand = gsl_rng_uniform(rnd) * this->propensity; + double a_sum = 0.0; + virus_map::iterator iter = V->pop.begin(), + end = V->pop.begin(); + for(; iter != end; ++iter) { + a_sum += RC_fun(iter->first.energy) * iter->second; + if (rand <= a_sum) + break; + } + + // 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; + Virus v(H,mu,iter->first.mutated_sites); + unsigned int n = gsl_ran_binomial(rnd,mu,H.size); + std::set sites_to_mutate; + while (sites_to_mutate.size() < n) { + sites_to_mutate.insert(gsl_rng_uniform_int(rnd,H.size)); + } + for (std::set::iterator it = sites_to_mutate.begin(), + 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); + } + /* + // this version is no good when mu is high or length is too small: + // it can propose duplicates, under-estimating the true extent + for (unsigned i = 0; i < n; ++i) { + unsigned int x = gsl_rng_uniform_int(rnd,H.size); + if (v.mutated_sites.count(x)==0) v.mutated_sites.insert(x); + else v.mutated_sites.erase(x); + } + */ + v.update_energy(H); + // update copy numbers + if (V->pop.count(v)==0) V->pop[v] = 1; // new variant + else V->pop[v] += 1; // add to existing pool + +} + +double VirusDeathReaction::recalc() +{ + V->count = 0; + propensity = rate_constant; + + double partial_sum = 0.0; + virus_map::iterator iter = this->V->pop.begin(), + end = this->V->pop.end(); + for(; iter != end; ++iter) { + V->count += iter->second; + partial_sum += iter->second; + } + + propensity *= partial_sum; + return propensity; +} + +void VirusDeathReaction::fire() +{ + // pick which virus dies + double rand = gsl_rng_uniform(rnd) * this->propensity; + double a_sum = 0.0; + virus_map::iterator iter = V->pop.begin(), + end = V->pop.begin(); + for (; iter != end; ++iter) { + a_sum += iter->second; + if (rand <= a_sum) + break; + } + + // update copy number + V->pop[iter->first] -= 1; + + // delete zero-population strains + if (V->pop[iter->first] == 0) + V->pop.erase(iter); +} + + +double KillingReaction::recalc() +{ + return this->propensity; +} +void KillingReaction::fire() +{ + // update copy numbers } diff --git a/reaction.h b/reaction.h index 90039ae..b7aae00 100644 --- a/reaction.h +++ b/reaction.h @@ -28,6 +28,16 @@ public: class VirusReaction : public Reaction { public: VirusSpecies* V; + Hamiltonian H; + double recalc(); + void fire(); +}; + + +class VirusDeathReaction : public Reaction { +public: + VirusSpecies* V; + Hamiltonian H; double recalc(); void fire(); }; diff --git a/ss.cpp b/ss.cpp index c64ee93..b40ac8f 100644 --- a/ss.cpp +++ b/ss.cpp @@ -14,6 +14,7 @@ #include "ss.h" #include "reaction.h" +gsl_rng *rnd; //pointer to random number generator state // generate a random number generator seed from /dev/urandom // borrowed from SSC src/runtime/ssc_rts.c @@ -42,13 +43,135 @@ void run(RunParameters_SS &r, unsigned seed) { // Initialize RNG and set initial state, if importing from file - static gsl_rng *rnd = gsl_rng_alloc(gsl_rng_taus2); //pointer to random number generator state + //static gsl_rng *rnd = gsl_rng_alloc(gsl_rng_taus2); //pointer to random number generator state + rnd = gsl_rng_alloc(gsl_rng_taus2); //pointer to random number generator state //srand((unsigned)time(0)); //gsl_rng_set(rnd,rand()); gsl_rng_set(rnd,seed); #if 1 + std::string couplingsFile = "neutral_2site.j"; + Hamiltonian H(couplingsFile); + //double mu = 6.0e-5; + double mu = 1.0e-1; + + std::vector species; + VirusSpecies s1; s1.count = 100; + s1.pop[Virus(H,mu)] = s1.count; + species.push_back(&s1); + + std::vector reactions; + // V -> V + V' + VirusReaction* r1 = new VirusReaction(); + r1->rate_constant = 3.0; + r1->H = H; + r1->V = &s1; + reactions.push_back(r1); + // V -> 0 + VirusDeathReaction* r2 = new VirusDeathReaction(); + r2->rate_constant = 1.5; + r2->H = H; + r2->V = &s1; + reactions.push_back(r2); + +/* + double total_propensity = 0.0; + for (std::vector::iterator rxn = reactions.begin(), + end = reactions.end(); + rxn != end; ++rxn) { + (*rxn)->recalc(); + total_propensity += (*rxn)->propensity; + } + std::cout << total_propensity << std::endl; + + std::cout << ((VirusSpecies*)species.back())->pop.begin()->second << std::endl; + std::cout << ((VirusReaction*)reactions.back())->V->pop.begin()->second << std::endl; + + reactions.back()->fire(); + + std::cout << ((VirusSpecies*)species.back())->pop.begin()->second << std::endl; + std::cout << ((VirusReaction*)reactions.back())->V->pop.begin()->second << std::endl; +*/ + + double t = 0.0; + double t_end = 100.0; + double dt = 0.0; + + double total_propensity = 0.0; + for (std::vector::iterator rxn = reactions.begin(), + end = reactions.end(); + rxn != end; ++rxn) { + (*rxn)->recalc(); + total_propensity += (*rxn)->propensity; + } + + + while (t < t_end) { + + // determine next reaction + + double rand = gsl_rng_uniform(rnd) * total_propensity; // uniform on [0,R) + double a_sum = 0.0; + + dt = -log(gsl_rng_uniform(rnd)) / total_propensity; + t += dt; + + for (std::vector::iterator rxn = reactions.begin(), + end = reactions.end(); + rxn != end; ++rxn) { + a_sum += (*rxn)->propensity; + if (rand <= a_sum) { + (*rxn)->fire(); + break; + } + } + + // TODO: handle case with total propensity == 0.0 + + // recalculate rates + total_propensity = 0.0; + for (std::vector::iterator rxn = reactions.begin(), + end = reactions.end(); + rxn != end; ++rxn) { + (*rxn)->recalc(); + total_propensity += (*rxn)->propensity; + } + + // print copy numbers + std::cout << t << '\t'; + for (std::vector::iterator spec = species.begin(), + end = species.end(); + spec != end; ++spec) { + std::cout << (*spec)->count << '\t'; + } + + std::set ms; + FILE* output = stdout; + for (virus_map::iterator iter = s1.pop.begin(), + end = s1.pop.end(); + iter != end; ++iter) { + std::cout << iter->second << '('; + ms = iter->first.mutated_sites; + // print sequence as CSV (comma-separated) + std::set::iterator ms_iter=ms.begin(); + if (ms_iter != ms.end()) { + fprintf(output,"%d",*ms_iter); + ++ms_iter; + } + for (; ms_iter!=ms.end(); ++ms_iter) + fprintf(output,",%d",*ms_iter); + std::cout << ')'; + std::cout << '\t'; + } + std::cout << '\n'; + + } + + +#endif + +#if 0 std::vector species; SimpleSpecies s1; s1.count = 100; s1.name = "A"; SimpleSpecies s2; s2.count = 50; s2.name = "B"; @@ -77,7 +200,6 @@ void run(RunParameters_SS &r, unsigned seed) { double t_end = 100.0; double dt = 0.0; - // initialize double total_propensity = 0.0; for (std::vector::iterator rxn = reactions.begin(), end = reactions.end(); @@ -94,6 +216,9 @@ void run(RunParameters_SS &r, unsigned seed) { double rand = gsl_rng_uniform(rnd) * total_propensity; // uniform on [0,R) double a_sum = 0.0; + dt = -log(gsl_rng_uniform(rnd)) / total_propensity; + t += dt; + for (std::vector::iterator rxn = reactions.begin(), end = reactions.end(); rxn != end; ++rxn) { @@ -115,9 +240,6 @@ void run(RunParameters_SS &r, unsigned seed) { total_propensity += (*rxn)->propensity; } - dt = -log(gsl_rng_uniform(rnd)) / total_propensity; - t += dt; - // print copy numbers std::cout << t << '\t'; diff --git a/virus.h b/virus.h index b254d80..f985684 100644 --- a/virus.h +++ b/virus.h @@ -28,7 +28,7 @@ public: friend class Population; -private: +//private: double mu; unsigned int L; -- 2.7.4