From 922d997cc023f1c4631c26173171880aa24abccf Mon Sep 17 00:00:00 2001 From: Dariusz Murakowski Date: Thu, 9 Apr 2015 08:03:00 -0400 Subject: [PATCH] Framework for killing by CTL. --- pop_ss.cpp | 35 +++++++++++++++++++++++++++++++++++ pop_ss.h | 14 ++++++++++++-- reaction.cpp | 45 ++++++++++++++++++++++++++++++++++++--------- reaction.h | 21 +++++++++++++++++++-- ss.cpp | 41 ++++++++++++++++++++++++----------------- 5 files changed, 126 insertions(+), 30 deletions(-) diff --git a/pop_ss.cpp b/pop_ss.cpp index dc480d8..e0a197b 100644 --- a/pop_ss.cpp +++ b/pop_ss.cpp @@ -8,6 +8,41 @@ #include "pop_ss.h" +// Did virus escape from immune pressure? +bool EpitopeRecognizer::escaped(const Virus &v) const { + return EpitopeRecognizer::escaped(v.mutated_sites); +} + +// Return true if the sequence has escaped from immune pressure at particular epitope +// i.e. if any sites listed in epitopeWT are mutated, or sites in epitopeMut are wildtype +bool EpitopeRecognizer::escaped(const MutatedSiteSequence &mutated_sites) const { + bool escape=false; + for (unsigned i=0;i MutatedSiteSequence; + class EpitopeRecognizer { public: - std::vector epitopeWT; - std::vector epitopeMut; + std::vector epitopeWT; // sites with state=0 are recognized + std::vector epitopeMut; // sites with state=1 are recognized + //unsigned int ep_len; + bool escaped(const Virus &v) const; + bool escaped(const MutatedSiteSequence &mutated_sites) const; }; class Species { @@ -68,6 +73,11 @@ public: class CTLSpecies : public Species { public: std::vector epitopes; + std::vector affinity; + size_t num_ep; + + double recognized(const Virus &v) const; + double recognized(const MutatedSiteSequence &mutated_sites) const; }; diff --git a/reaction.cpp b/reaction.cpp index 37c34e5..3eca7ee 100644 --- a/reaction.cpp +++ b/reaction.cpp @@ -83,13 +83,13 @@ double VirusReaction::recalc() void VirusReaction::fire() { // pick which virus reproduces - double rand = gsl_rng_uniform(rnd) * this->propensity; + double rand = gsl_rng_uniform_pos(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) + if (rand <= a_sum * rate_constant) break; } @@ -99,18 +99,18 @@ void VirusReaction::fire() 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; + MutatedSiteSequence 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(), + for (MutatedSiteSequence::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: + // this version might be 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); @@ -146,13 +146,13 @@ double VirusDeathReaction::recalc() void VirusDeathReaction::fire() { // pick which virus dies - double rand = gsl_rng_uniform(rnd) * this->propensity; + double rand = gsl_rng_uniform_pos(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) + if (rand <= a_sum * rate_constant) break; } @@ -167,11 +167,38 @@ void VirusDeathReaction::fire() double KillingReaction::recalc() { - return this->propensity; + V->count = 0; + propensity = rate_constant * T->count; + + double partial_sum = 0.0; + virus_map::iterator iter = V->pop.begin(), + end = V->pop.end(); + for(; iter != end; ++iter) { + V->count += iter->second; + partial_sum += iter->second * T->recognized(iter->first); + } + + propensity *= partial_sum; + return propensity; } void KillingReaction::fire() { - // update copy numbers + // pick which virus to kill + double rand = gsl_rng_uniform_pos(rnd) * this->propensity; + double a_sum = 0.0; + virus_map::iterator iter = V->pop.begin(), + end = V->pop.end(); + for(; iter != end; ++iter) { + a_sum += iter->second * T->recognized(iter->first); + if (rand <= a_sum * rate_constant * T->count) + break; + } + + // update copy number + V->pop[iter->first] -= 1; + // delete zero-population strains + if (V->pop[iter->first] == 0) + V->pop.erase(iter); } diff --git a/reaction.h b/reaction.h index b7aae00..7501941 100644 --- a/reaction.h +++ b/reaction.h @@ -5,6 +5,8 @@ #include "pop_ss.h" +typedef std::set MutatedSiteSequence; + class Reaction { public: double rate_constant; // must manually ensure it includes division by statistical degeneracy factor (in case of reaction order >1 for any species) @@ -13,7 +15,8 @@ public: virtual void fire() = 0; }; - +// a A + b B + ... -> c C + d D + ... +// mass action kinetics class ElementaryReaction : public Reaction { public: std::vector reactants; @@ -25,6 +28,8 @@ public: }; +// V_s -> V_s + V_s' +// depending on energy, with mutation class VirusReaction : public Reaction { public: VirusSpecies* V; @@ -34,6 +39,8 @@ public: }; +// V_s -> 0 +// in principle energy-dependent, but not here class VirusDeathReaction : public Reaction { public: VirusSpecies* V; @@ -42,7 +49,8 @@ public: void fire(); }; - +// V_s + T_k -> T_k +// depending on epitope affinity class KillingReaction : public Reaction { public: VirusSpecies* V; @@ -51,5 +59,14 @@ public: void fire(); }; +class CTLTransitionReaction : public Reaction { +public: + VirusSpecies* V; + CTLSpecies* Tfrom; + CTLSpecies* Tto; + double recalc(); + void fire(); +}; + #endif // REACTION_H diff --git a/ss.cpp b/ss.cpp index b40ac8f..d7b09b8 100644 --- a/ss.cpp +++ b/ss.cpp @@ -14,6 +14,8 @@ #include "ss.h" #include "reaction.h" +typedef std::vector Rxn_parray; + gsl_rng *rnd; //pointer to random number generator state // generate a random number generator seed from /dev/urandom @@ -61,7 +63,7 @@ void run(RunParameters_SS &r, unsigned seed) { s1.pop[Virus(H,mu)] = s1.count; species.push_back(&s1); - std::vector reactions; + Rxn_parray reactions; // V -> V + V' VirusReaction* r1 = new VirusReaction(); r1->rate_constant = 3.0; @@ -95,31 +97,36 @@ void run(RunParameters_SS &r, unsigned seed) { */ double t = 0.0; - double t_end = 100.0; + double t_end = 1.0; + double t_remain = t_end; double dt = 0.0; double total_propensity = 0.0; - for (std::vector::iterator rxn = reactions.begin(), - end = reactions.end(); - rxn != end; ++rxn) { + + Rxn_parray::iterator rxn, end; + + for (rxn = reactions.begin(), end = reactions.end(); rxn != end; ++rxn) { (*rxn)->recalc(); total_propensity += (*rxn)->propensity; } - while (t < t_end) { - - // determine next reaction + while (1) { - double rand = gsl_rng_uniform(rnd) * total_propensity; // uniform on [0,R) + double rand = (1.0-gsl_rng_uniform_pos(rnd)) * total_propensity; // uniform on (0,R] double a_sum = 0.0; + // time to reaction dt = -log(gsl_rng_uniform(rnd)) / total_propensity; t += dt; + t_remain -= dt; - for (std::vector::iterator rxn = reactions.begin(), - end = reactions.end(); - rxn != end; ++rxn) { + // stop if time ran out + if (t_remain < 0.0) + break; + + // determine next reaction + for (rxn = reactions.begin(), end = reactions.end(); rxn != end; ++rxn) { a_sum += (*rxn)->propensity; if (rand <= a_sum) { (*rxn)->fire(); @@ -128,12 +135,12 @@ void run(RunParameters_SS &r, unsigned seed) { } // TODO: handle case with total propensity == 0.0 + if (rxn == end) + break; // recalculate rates total_propensity = 0.0; - for (std::vector::iterator rxn = reactions.begin(), - end = reactions.end(); - rxn != end; ++rxn) { + for (rxn = reactions.begin(), end = reactions.end(); rxn != end; ++rxn) { (*rxn)->recalc(); total_propensity += (*rxn)->propensity; } @@ -146,7 +153,7 @@ void run(RunParameters_SS &r, unsigned seed) { std::cout << (*spec)->count << '\t'; } - std::set ms; + MutatedSiteSequence ms; FILE* output = stdout; for (virus_map::iterator iter = s1.pop.begin(), end = s1.pop.end(); @@ -154,7 +161,7 @@ void run(RunParameters_SS &r, unsigned seed) { std::cout << iter->second << '('; ms = iter->first.mutated_sites; // print sequence as CSV (comma-separated) - std::set::iterator ms_iter=ms.begin(); + MutatedSiteSequence::iterator ms_iter=ms.begin(); if (ms_iter != ms.end()) { fprintf(output,"%d",*ms_iter); ++ms_iter; -- 2.7.4