From 2095bbaa73c5b79c3737c8e2bab105d30ffe69b4 Mon Sep 17 00:00:00 2001 From: Dariusz Murakowski Date: Fri, 24 Apr 2015 04:05:17 -0400 Subject: [PATCH] Define all the reactions for Potts virus, handling nucleotide-level mutations. --- ham_ss.h | 2 +- reaction.cpp | 197 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++- reaction.h | 60 ++++++++++++++++++ 3 files changed, 257 insertions(+), 2 deletions(-) diff --git a/ham_ss.h b/ham_ss.h index ea0a946..76d0c62 100644 --- a/ham_ss.h +++ b/ham_ss.h @@ -77,7 +77,7 @@ inline int index(int i, int j, size_t length) { class PottsHamiltonian { public: - unsigned int size; + unsigned int size; // amino acid length double bh; // fields "inverse temperature" multiplier double bJ; // couplings "inverse temperature" multiplier Coupling J; diff --git a/reaction.cpp b/reaction.cpp index e8a2e25..99f0daf 100644 --- a/reaction.cpp +++ b/reaction.cpp @@ -98,7 +98,7 @@ void VirusReaction::fire() // pick how many sites to mutate in the new virus // adapted from WF Virus::mutate() - Virus v(H,iter->first.mutated_sites); + Virus v(H,iter->first.mutated_sites); // TODO: replace with copy constructor unsigned int n = gsl_ran_binomial(rnd,mu,H.size); MutatedSiteSequence sites_to_mutate; while (sites_to_mutate.size() < n) { @@ -227,3 +227,198 @@ void CTLActivationReaction::fire() Tto->count += 1; } + +//////////////////////////////////////// +// Potts code below +//////////////////////////////////////// + + + +double NTVirusReaction::recalc() +{ + V->count = 0; + this->propensity = this->rate_constant; + + double partial_sum = 0.0; + NTVirus_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 NTVirusReaction::fire() +{ + // pick which virus reproduces + double rand = gsl_rng_uniform_pos(rnd) * this->propensity; + double a_sum = 0.0; + NTVirus_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 * rate_constant) + break; + } + + // pick how many sites (nucleotides) to mutate in the new virus + // adapted from WF Virus::mutate() + NTVirus v(H,iter->first.nt_seq); // TODO: replace with copy constructor + unsigned int n = gsl_ran_binomial(rnd,mu,3*H.size); + MutatedSiteSequence sites_to_mutate; + while (sites_to_mutate.size() < n) { + sites_to_mutate.insert(gsl_rng_uniform_int(rnd,3*H.size)); + } + + bool potts_state_changed = false; + for (MutatedSiteSequence::iterator it = sites_to_mutate.begin(), + end = sites_to_mutate.end(); + it != end; ++it) { + unsigned siteNT = *it; + + // nucleotide to codon + int oldNT = str2nt.at(v.nt_seq[siteNT]); // silent conversion from enum nt + int newNT; + do { // pick rand on {0,1,2,3} until we get a new one + newNT = gsl_rng_uniform_int(rnd,4-1); + } while (newNT == oldNT); + v.nt_seq[siteNT] = newNT; + + // codon to amino acid + unsigned siteAA = siteNT / 3; + aa oldAA = v.aa_seq[siteAA]; + aa newAA = codon2aa.at(str2codon.at(v.nt_seq.substr(siteAA,3))); + + // amino acid to Potts state + if (newAA != oldAA) { + v.aa_seq[siteAA] = newAA; + unsigned oldState = v.config[siteAA]; + unsigned newState = H.aa2state[siteAA][newAA]; + if (newState != oldState) { + v.config[siteAA] = newState; + potts_state_changed = true; + } + } + } + /* + // 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); + if (v.mutated_sites.count(x)==0) v.mutated_sites.insert(x); + else v.mutated_sites.erase(x); + } + */ + if (potts_state_changed) + 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 NTVirusDeathReaction::recalc() +{ + V->count = 0; + propensity = rate_constant; + + double partial_sum = 0.0; + NTVirus_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 NTVirusDeathReaction::fire() +{ + // pick which virus dies + double rand = gsl_rng_uniform_pos(rnd) * this->propensity; + double a_sum = 0.0; + NTVirus_map::iterator iter = V->pop.begin(), + end = V->pop.begin(); + for (; iter != end; ++iter) { + a_sum += iter->second; + if (rand <= a_sum * rate_constant) + break; + } + + // update copy number + V->pop[iter->first] -= 1; + + // delete zero-population strains + if (V->pop[iter->first] == 0) + V->pop.erase(iter); +} + + +double AAKillingReaction::recalc() +{ + V->count = 0; + propensity = rate_constant * T->count; + + double partial_sum = 0.0; + NTVirus_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 AAKillingReaction::fire() +{ + // pick which virus to kill + double rand = gsl_rng_uniform_pos(rnd) * this->propensity; + double a_sum = 0.0; + NTVirus_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); +} + + +double CTLaaActivationReaction::recalc() +{ + propensity = rate_constant * Tfrom->count; + + double partial_sum = 0.0; + NTVirus_map::iterator iter = V->pop.begin(), + end = V->pop.end(); + for(; iter != end; ++iter) { + partial_sum += iter->second * Tfrom->recognized(iter->first); + } + + propensity *= partial_sum; + return propensity; + +} + +void CTLaaActivationReaction::fire() +{ + Tfrom->count -= 1; + Tto->count += 1; +} + + diff --git a/reaction.h b/reaction.h index b74051a..2e5e52c 100644 --- a/reaction.h +++ b/reaction.h @@ -90,5 +90,65 @@ public: void fire(); }; + +//////////////////////////////////////// +// Potts code below +//////////////////////////////////////// + +// V_s -> V_s + V_s' +// depending on energy, with mutation +class NTVirusReaction : public Reaction { +public: + NTVirusReaction() : Reaction() {} + NTVirusReaction(double k) : Reaction(k), mu(0) {} + + double mu; + NTVirusSpecies* V; + PottsHamiltonian H; + + double recalc(); + void fire(); +}; + + +// V_s -> 0 +// in principle energy-dependent, but not here +class NTVirusDeathReaction : public Reaction { +public: + NTVirusDeathReaction() : Reaction() {} + NTVirusDeathReaction(double k) : Reaction(k) {} + + NTVirusSpecies* V; + PottsHamiltonian H; + double recalc(); + void fire(); +}; + +// V_s + T_k -> T_k +// depending on epitope affinity +class AAKillingReaction : public Reaction { +public: + AAKillingReaction() : Reaction() {} + AAKillingReaction(double k) : Reaction(k) {} + + NTVirusSpecies* V; + CTLaaSpecies* T; + double recalc(); + void fire(); +}; + +// V_s + T_k -> V_s + T_k' +class CTLaaActivationReaction : public Reaction { +public: + CTLaaActivationReaction() : Reaction() {} + CTLaaActivationReaction(double k) : Reaction(k) {} + + NTVirusSpecies* V; + CTLaaSpecies* Tfrom; + CTLaaSpecies* Tto; + double recalc(); + void fire(); +}; + #endif // REACTION_H -- 2.7.4