#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<epitopeWT.size();i++) {
+ if (mutated_sites.count(epitopeWT[i])==1) { escape=true; break; }
+ }
+ if (!escape) {
+ for (unsigned i=0;i<epitopeMut.size();i++) {
+ if (mutated_sites.count(epitopeMut[i])==0) { escape=true; break; }
+ }
+ }
+ return escape;
+}
+
+// returns the affinity to which a virus is recognized
+// or 0 if it is not recognized (i.e. has escaped from all epitopes)
+// If multiple epitopes are present and recognized, returns the total affinity
+double CTLSpecies::recognized(const Virus &v) const {
+ return CTLSpecies::recognized(v.mutated_sites);
+}
+double CTLSpecies::recognized(const MutatedSiteSequence &mutated_sites) const {
+ double bind = 0.0;
+ for (size_t i = 0; i < num_ep; ++i) {
+ if (!epitopes[i].escaped(mutated_sites)) {
+ bind += affinity[i];
+ }
+ }
+ return bind;
+}
//Constructor that assembles a population of N viruses of wild type
};
#endif
+typedef std::set<unsigned int> MutatedSiteSequence;
+
class EpitopeRecognizer {
public:
- std::vector<unsigned int> epitopeWT;
- std::vector<unsigned int> epitopeMut;
+ std::vector<unsigned int> epitopeWT; // sites with state=0 are recognized
+ std::vector<unsigned int> 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 {
class CTLSpecies : public Species {
public:
std::vector<EpitopeRecognizer> epitopes;
+ std::vector<double> affinity;
+ size_t num_ep;
+
+ double recognized(const Virus &v) const;
+ double recognized(const MutatedSiteSequence &mutated_sites) const;
};
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;
}
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<unsigned int> 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<unsigned int>::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);
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;
}
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);
}
#include "pop_ss.h"
+typedef std::set<unsigned int> 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)
virtual void fire() = 0;
};
-
+// a A + b B + ... -> c C + d D + ...
+// mass action kinetics
class ElementaryReaction : public Reaction {
public:
std::vector<SimpleSpecies*> reactants;
};
+// V_s -> V_s + V_s'
+// depending on energy, with mutation
class VirusReaction : public Reaction {
public:
VirusSpecies* V;
};
+// V_s -> 0
+// in principle energy-dependent, but not here
class VirusDeathReaction : public Reaction {
public:
VirusSpecies* V;
void fire();
};
-
+// V_s + T_k -> T_k
+// depending on epitope affinity
class KillingReaction : public Reaction {
public:
VirusSpecies* V;
void fire();
};
+class CTLTransitionReaction : public Reaction {
+public:
+ VirusSpecies* V;
+ CTLSpecies* Tfrom;
+ CTLSpecies* Tto;
+ double recalc();
+ void fire();
+};
+
#endif // REACTION_H
#include "ss.h"
#include "reaction.h"
+typedef std::vector<Reaction*> Rxn_parray;
+
gsl_rng *rnd; //pointer to random number generator state
// generate a random number generator seed from /dev/urandom
s1.pop[Virus(H,mu)] = s1.count;
species.push_back(&s1);
- std::vector<Reaction*> reactions;
+ Rxn_parray reactions;
// V -> V + V'
VirusReaction* r1 = new VirusReaction();
r1->rate_constant = 3.0;
*/
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<Reaction*>::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<Reaction*>::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();
}
// TODO: handle case with total propensity == 0.0
+ if (rxn == end)
+ break;
// recalculate rates
total_propensity = 0.0;
- for (std::vector<Reaction*>::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;
}
std::cout << (*spec)->count << '\t';
}
- std::set<unsigned int> ms;
+ MutatedSiteSequence ms;
FILE* output = stdout;
for (virus_map::iterator iter = s1.pop.begin(),
end = s1.pop.end();
std::cout << iter->second << '(';
ms = iter->first.mutated_sites;
// print sequence as CSV (comma-separated)
- std::set<unsigned int>::iterator ms_iter=ms.begin();
+ MutatedSiteSequence::iterator ms_iter=ms.begin();
if (ms_iter != ms.end()) {
fprintf(output,"%d",*ms_iter);
++ms_iter;
projects / VirEvoDyn.git / commitdiff