From: Dariusz Murakowski Date: Fri, 17 Apr 2015 17:18:25 +0000 (-0400) Subject: Full system, including parameters. X-Git-Url: http://src.murakowski.org/?a=commitdiff_plain;h=a7896b781628c26c508ec7e9e12d86ef696d1243;p=VirEvoDyn.git Full system, including parameters. --- diff --git a/reaction.cpp b/reaction.cpp index c4ac28a..3ee908f 100644 --- a/reaction.cpp +++ b/reaction.cpp @@ -14,7 +14,7 @@ double ElementaryReaction::recalc() //double denom = 1.0; this->propensity = this->rate_constant; - std::vector::iterator reac = this->reactants.begin(), + std::vector::iterator reac = this->reactants.begin(), reac_end = this->reactants.end(); std::vector::iterator reacN = this->reactant_stoich.begin(), reacN_end = this->reactant_stoich.end(); @@ -33,7 +33,7 @@ double ElementaryReaction::recalc() void ElementaryReaction::fire() { - std::vector::iterator reac = this->reactants.begin(), + std::vector::iterator reac = this->reactants.begin(), reac_end = this->reactants.end(); std::vector::iterator reacN = this->reactant_stoich.begin(), reacN_end = this->reactant_stoich.end(); @@ -46,7 +46,7 @@ void ElementaryReaction::fire() } */ - std::vector::iterator prod = this->products.begin(), + std::vector::iterator prod = this->products.begin(), prod_end = this->products.end(); std::vector::iterator prodN = this->product_stoich.begin(), prodN_end = this->product_stoich.end(); @@ -98,7 +98,8 @@ void VirusReaction::fire() // 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; + //double mu = 1.0e-1; + double mu = iter->first.mu; Virus v(H,mu,iter->first.mutated_sites); unsigned int n = gsl_ran_binomial(rnd,mu,H.size); MutatedSiteSequence sites_to_mutate; diff --git a/reaction.h b/reaction.h index 3ea501a..2f1ab57 100644 --- a/reaction.h +++ b/reaction.h @@ -20,13 +20,14 @@ public: // a A + b B + ... -> c C + d D + ... // mass action kinetics +// the Species element should really be a SimpleSpecies, but really anything other than VirusSpecies is okay class ElementaryReaction : public Reaction { public: ElementaryReaction() : Reaction() {} ElementaryReaction(double k) : Reaction(k) {} - std::vector reactants; // anything that participates in reaction rate - std::vector products; // effects of reaction + std::vector reactants; // anything that participates in reaction rate + std::vector products; // effects of reaction std::vector reactant_stoich; // reaction order std::vector product_stoich; // must be negative for reactants that are consumed! double recalc(); diff --git a/ss.cpp b/ss.cpp index 96841bb..858e4af 100644 --- a/ss.cpp +++ b/ss.cpp @@ -59,34 +59,36 @@ void run(RunParameters_SS &r, unsigned seed) { std::string couplingsFile = "neutral_2site.j"; Hamiltonian H(couplingsFile); + H.set_temp(0.01); //double mu = 6.0e-5; //double mu = 1.0e-1; // initialize virus VirusSpecies s1; s1.count = r.n; - s1.pop[Virus(H,r.mu)] = s1.count; + s1.pop[Virus(H,r.mu)] = s1.count; // default mu = 6.0e-5 species.push_back(&s1); + Species_parray print_spec; + print_spec.push_back(&s1); + if (r.importState) importState(r); // .st - if (r.useEpitope) importEpitope(r,species,reactions,&s1); // .ep + if (r.useEpitope) importEpitope(r,species,reactions,&s1,print_spec); // .ep // V -> V + V' - VirusReaction* r1 = new VirusReaction(); - r1->rate_constant = 3.0; + VirusReaction* r1 = new VirusReaction(1.5); // s_k + //r1->rate_constant = 1.5; r1->H = H; r1->V = &s1; reactions.push_back(r1); // V -> 0 - VirusDeathReaction* r2 = new VirusDeathReaction(); - r2->rate_constant = 1.5; + VirusDeathReaction* r2 = new VirusDeathReaction(0.5); // u + //r2->rate_constant = 0.5; r2->H = H; r2->V = &s1; reactions.push_back(r2); - Species_parray print_spec; - print_spec.push_back(&s1); - double t_end = 10.0; + double t_end = 1000.0; long max_steps = LONG_MAX; double sample_interval = 1.0; simulate(reactions, species, t_end, max_steps, sample_interval, print_spec); @@ -665,7 +667,7 @@ void importState(RunParameters_SS &r) { // load epitope definitions from file -void importEpitope(RunParameters_SS &r, Species_parray &species, Rxn_parray &reactions, VirusSpecies *V) +void importEpitope(RunParameters_SS &r, Species_parray &species, Rxn_parray &reactions, VirusSpecies *V, Species_parray &print_spec) { if ((r.num_CTL_clones != r.epitopefiles.size()) || @@ -742,24 +744,73 @@ void importEpitope(RunParameters_SS &r, Species_parray &species, Rxn_parray &rea species.push_back(T); // effector + CTLSpecies* N = new CTLSpecies(*T); // naive + species.push_back(N); + CTLSpecies* M = new CTLSpecies(*T); // memory species.push_back(M); - CTLSpecies* N = new CTLSpecies(*T); - species.push_back(N); // naive - - KillingReaction* r1 = new KillingReaction(1e-5); - r1->T = T; r1->V = V; - reactions.push_back(r1); - - CTLActivationReaction* r2 = new CTLActivationReaction(1e-6); - r2->Tfrom = M; r2->Tto = T; r2->V = V; - reactions.push_back(r2); - - CTLActivationReaction* r3 = new CTLActivationReaction(1e-7); - r3->Tfrom = M; r3->Tto = T; r3->V = V; - reactions.push_back(r3); + print_spec.push_back(T); + print_spec.push_back(N); + print_spec.push_back(M); + { // p + KillingReaction* r = new KillingReaction(1e-5); + r->T = T; r->V = V; + reactions.push_back(r); + } + { // a + CTLActivationReaction* r = new CTLActivationReaction(1e-7); + r->Tfrom = N; r->Tto = T; r->V = V; + reactions.push_back(r); + } + { // r + ElementaryReaction* r = new ElementaryReaction(4.0); + r->reactants.push_back(T); r->reactant_stoich.push_back(1); + r->products.push_back(T); r->product_stoich.push_back(1); + reactions.push_back(r); + } + { // d = 0.5 // note d' = 3.0 + ElementaryReaction* r = new ElementaryReaction(3.0); + r->reactants.push_back(T); r->reactant_stoich.push_back(1); + r->products.push_back(T); r->product_stoich.push_back(-1); + reactions.push_back(r); + } + { // b + ElementaryReaction* r = new ElementaryReaction(1e-4); + r->reactants.push_back(N); r->reactant_stoich.push_back(1); + r->products.push_back(N); r->product_stoich.push_back(1); + reactions.push_back(r); + } + { // e + ElementaryReaction* r = new ElementaryReaction(3e-4); + r->reactants.push_back(N); r->reactant_stoich.push_back(1); + r->products.push_back(N); r->product_stoich.push_back(-1); + reactions.push_back(r); + } + { // w + ElementaryReaction* r = new ElementaryReaction(7.5); + r->products.push_back(N); r->product_stoich.push_back(1); + reactions.push_back(r); + } + { // a' + CTLActivationReaction* r = new CTLActivationReaction(1e-6); + r->Tfrom = M; r->Tto = T; r->V = V; + reactions.push_back(r); + } + { // g + ElementaryReaction* r = new ElementaryReaction(0.03); + r->reactants.push_back(T); r->reactant_stoich.push_back(1); + r->products.push_back(T); r->product_stoich.push_back(-1); + r->products.push_back(M); r->product_stoich.push_back(1); + reactions.push_back(r); + } + { // h + ElementaryReaction* r = new ElementaryReaction(0.03); + r->reactants.push_back(M); r->reactant_stoich.push_back(1); + r->products.push_back(M); r->product_stoich.push_back(-1); + reactions.push_back(r); + } for (size_t i = 0; i < T->num_ep; ++i) { // XXX std::cout << T->count << '\t' diff --git a/ss.h b/ss.h index be805c4..2c26221 100644 --- a/ss.h +++ b/ss.h @@ -231,7 +231,7 @@ public: void run(RunParameters_SS &r); void importState(RunParameters_SS &r); -void importEpitope(RunParameters_SS &r, Species_parray &species, Rxn_parray &reactions, VirusSpecies *V); +void importEpitope(RunParameters_SS &r, Species_parray &species, Rxn_parray &reactions, VirusSpecies *V, Species_parray &print_spec); void simulate(Rxn_parray &reactions, Species_parray &species, double t_end, long max_steps, double sample_interval, const Species_parray &print_spec); diff --git a/virus.cpp b/virus.cpp index fd69023..f523ca4 100644 --- a/virus.cpp +++ b/virus.cpp @@ -101,7 +101,3 @@ bool operator<(const Virus& lhs, const Virus& rhs) { } - - - -