Define all the reactions for Potts virus, handling nucleotide-level mutations.

diff --git a/ham_ss.h b/ham_ss.h
index ea0a946..76d0c62 100644 (file)
--- 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 (file)
--- 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 (file)
--- 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