From: Dariusz Murakowski Date: Thu, 9 Apr 2015 14:31:39 +0000 (-0400) Subject: Initialize multiple CTL clones+epitopes. Output at intervals. Pretty-print the sequence. X-Git-Url: http://src.murakowski.org/?a=commitdiff_plain;h=19be763c5b68eae7194d6ab7e03b4e57aa0734ce;p=VirEvoDyn.git Initialize multiple CTL clones+epitopes. Output at intervals. Pretty-print the sequence. --- diff --git a/pop_ss.h b/pop_ss.h index 764cf9b..12e21ab 100644 --- a/pop_ss.h +++ b/pop_ss.h @@ -56,22 +56,32 @@ public: bool escaped(const MutatedSiteSequence &mutated_sites) const; }; + class Species { public: long count; std::string name; + Species() : count(0) {} + Species(long c) : count(c) {} }; class SimpleSpecies : public Species { + SimpleSpecies() : Species() {} + SimpleSpecies(long c) : Species(c) {} }; class VirusSpecies : public Species { public: + VirusSpecies() : Species() {} + VirusSpecies(long c) : Species(c) {} virus_map pop; }; class CTLSpecies : public Species { public: + CTLSpecies() : Species() {} + CTLSpecies(long c) : Species(c) {} + std::vector epitopes; std::vector affinity; size_t num_ep; diff --git a/reaction.cpp b/reaction.cpp index 20f55e3..c4ac28a 100644 --- a/reaction.cpp +++ b/reaction.cpp @@ -38,11 +38,13 @@ void ElementaryReaction::fire() std::vector::iterator reacN = this->reactant_stoich.begin(), reacN_end = this->reactant_stoich.end(); + /* for (; (reac != reac_end) && (reacN != reacN_end); ++reac, ++reacN) { //if (reac->count - reacN >= 0) { (*reac)->count = (*reac)->count - *reacN; //} } + */ std::vector::iterator prod = this->products.begin(), prod_end = this->products.end(); diff --git a/reaction.h b/reaction.h index 9865256..3ea501a 100644 --- a/reaction.h +++ b/reaction.h @@ -9,6 +9,9 @@ typedef std::set MutatedSiteSequence; class Reaction { public: + Reaction() {} + Reaction(double k) : rate_constant(k) {} + double rate_constant; // must manually ensure it includes division by statistical degeneracy factor (in case of reaction order >1 for any species) double propensity; virtual double recalc() = 0; @@ -19,10 +22,13 @@ public: // mass action kinetics class ElementaryReaction : public Reaction { public: - std::vector reactants; - std::vector products; + ElementaryReaction() : Reaction() {} + ElementaryReaction(double k) : Reaction(k) {} + + 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; + std::vector product_stoich; // must be negative for reactants that are consumed! double recalc(); void fire(); }; @@ -32,6 +38,9 @@ public: // depending on energy, with mutation class VirusReaction : public Reaction { public: + VirusReaction() : Reaction() {} + VirusReaction(double k) : Reaction(k) {} + VirusSpecies* V; Hamiltonian H; double recalc(); @@ -43,6 +52,9 @@ public: // in principle energy-dependent, but not here class VirusDeathReaction : public Reaction { public: + VirusDeathReaction() : Reaction() {} + VirusDeathReaction(double k) : Reaction(k) {} + VirusSpecies* V; Hamiltonian H; double recalc(); @@ -53,6 +65,9 @@ public: // depending on epitope affinity class KillingReaction : public Reaction { public: + KillingReaction() : Reaction() {} + KillingReaction(double k) : Reaction(k) {} + VirusSpecies* V; CTLSpecies* T; double recalc(); @@ -62,6 +77,9 @@ public: // V_s + T_k -> V_s + T_k' class CTLActivationReaction : public Reaction { public: + CTLActivationReaction() : Reaction() {} + CTLActivationReaction(double k) : Reaction(k) {} + VirusSpecies* V; CTLSpecies* Tfrom; CTLSpecies* Tto; diff --git a/ss.cpp b/ss.cpp index 7f08189..96841bb 100644 --- a/ss.cpp +++ b/ss.cpp @@ -52,20 +52,24 @@ void run(RunParameters_SS &r, unsigned seed) { gsl_rng_set(rnd,seed); r.setFiles(); - if (r.importState) importState(r); // .st - if (r.useEpitope) importEpitope(r); // .ep + + // master arrays of (pointers to) species and reactions + Species_parray species; + Rxn_parray reactions; std::string couplingsFile = "neutral_2site.j"; Hamiltonian H(couplingsFile); //double mu = 6.0e-5; - double mu = 1.0e-1; + //double mu = 1.0e-1; - std::vector species; - VirusSpecies s1; s1.count = 100; - s1.pop[Virus(H,mu)] = s1.count; + // initialize virus + VirusSpecies s1; s1.count = r.n; + s1.pop[Virus(H,r.mu)] = s1.count; species.push_back(&s1); - Rxn_parray reactions; + if (r.importState) importState(r); // .st + if (r.useEpitope) importEpitope(r,species,reactions,&s1); // .ep + // V -> V + V' VirusReaction* r1 = new VirusReaction(); r1->rate_constant = 3.0; @@ -82,9 +86,10 @@ void run(RunParameters_SS &r, unsigned seed) { Species_parray print_spec; print_spec.push_back(&s1); - double t_end = 1.0; + double t_end = 10.0; long max_steps = LONG_MAX; - simulate(reactions, species, t_end, max_steps, print_spec); + double sample_interval = 1.0; + simulate(reactions, species, t_end, max_steps, sample_interval, print_spec); #if 0 @@ -497,7 +502,7 @@ void run(RunParameters_SS &r, unsigned seed) { } -void simulate(Rxn_parray &reactions, Species_parray &species, double t_end, long max_steps, const 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) { long step = 0; @@ -505,6 +510,8 @@ void simulate(Rxn_parray &reactions, Species_parray &species, double t_end, long double t_remain = t_end; double dt = 0.0; + double t_next_sample = t + sample_interval; + double total_propensity = 0.0; Rxn_parray::iterator rxn, end; @@ -555,15 +562,21 @@ void simulate(Rxn_parray &reactions, Species_parray &species, double t_end, long ++step; - std::cout << t << '\t'; - for (Species_parray::const_iterator spec = print_spec.begin(), - end = print_spec.end(); - spec != end; ++spec) - { - // TODO: prettier representation - std::cout << (*spec)->count << '\t'; + + // print copy numbers + if (t_next_sample <= t) { + t_next_sample += sample_interval; + + std::cout << t << '\t'; + for (Species_parray::const_iterator spec = print_spec.begin(), + end = print_spec.end(); + spec != end; ++spec) + { + // TODO: prettier representation + std::cout << (*spec)->count << '\t'; + } + std::cout << '\n'; } - std::cout << '\n'; #if 0 // print copy numbers @@ -652,64 +665,110 @@ void importState(RunParameters_SS &r) { // load epitope definitions from file -void importEpitope(RunParameters_SS &r) { +void importEpitope(RunParameters_SS &r, Species_parray &species, Rxn_parray &reactions, VirusSpecies *V) +{ + + if ((r.num_CTL_clones != r.epitopefiles.size()) || + (r.num_CTL_clones != r.init_CTL_num.size())) { + perror((std::string("ERROR in importEpitope: invalid number of CTL clones").c_str())); + exit(1); + } - std::ifstream input(r.epitopeInfile.c_str()); // .ep - if (!input) { perror((std::string("ERROR in importEpitope: ") + r.epitopeInfile).c_str()); exit(1); } + for (size_t i = 0; i != r.num_CTL_clones; ++i) { + CTLSpecies* T = new CTLSpecies(r.init_CTL_num[i]); - std::string readStr; - while (std::getline(input,readStr)) { + std::ifstream input(r.epitopefiles[i].c_str()); // .ep + if (!input) { perror((std::string("ERROR in importEpitope: ") + r.epitopeInfile).c_str()); exit(1); } - std::string word; - unsigned int site; - double aff; - std::vector tmpEp; - size_t pos0 = 0; - size_t posBar = readStr.find('|',pos0); - size_t posEnd = std::string::npos; //readStr.size(); + std::string readStr; + while (std::getline(input,readStr)) { - // could use std::noskipws to keep tab ('\t') characters? + T->num_ep += 1; + T->epitopes.push_back(EpitopeRecognizer()); - std::stringstream readStrStrm(std::string(readStr,pos0,posBar-pos0)); + std::string word; + unsigned int site; + double aff; + std::vector tmpEp; + size_t pos0 = 0; + size_t posBar = readStr.find('|',pos0); + size_t posEnd = std::string::npos; //readStr.size(); - // first entry is affinity - readStrStrm >> word; - aff = strtodouble(word); - printf("%.6e ",aff); + // could use std::noskipws to keep tab ('\t') characters? - // next entries are WT sites - while (readStrStrm >> word) { - std::cout << word << " "; // XXX - std::istringstream i(word); - if (i >> site) - tmpEp.push_back(site); - else // must have encountered '|' - break; - } - r.eWT.push_back(tmpEp); + std::stringstream readStrStrm(std::string(readStr,pos0,posBar-pos0)); - tmpEp.clear(); // reset temp epitope - readStrStrm.str(""); readStrStrm.clear(); // reset stream + // first entry is affinity + readStrStrm >> word; + aff = strtodouble(word); + printf("%.6e ",aff); - std::cout << "| "; // XXX + T->affinity.push_back(aff); + + // next entries are WT sites + while (readStrStrm >> word) { + std::cout << word << " "; // XXX + std::istringstream i(word); + if (i >> site) + tmpEp.push_back(site); + else // must have encountered '|' + break; + } + //r.eWT.push_back(tmpEp); + T->epitopes.back().epitopeWT = tmpEp; + + tmpEp.clear(); // reset temp epitope + readStrStrm.str(""); readStrStrm.clear(); // reset stream + + std::cout << "| "; // XXX + + readStrStrm.str(std::string(readStr,posBar+1,posEnd/*-posBar-1*/)); + while (readStrStrm >> word) { + std::cout << word << " "; // XXX + std::istringstream i(word); + if (i >> site) + tmpEp.push_back(site); + else + break; + } + //r.eMut.push_back(tmpEp); + T->epitopes.back().epitopeMut = tmpEp; + + //r.numEpitopes++; + + std::cout << "\n"; // XXX - readStrStrm.str(std::string(readStr,posBar+1,posEnd/*-posBar-1*/)); - while (readStrStrm >> word) { - std::cout << word << " "; // XXX - std::istringstream i(word); - if (i >> site) - tmpEp.push_back(site); - else - break; } - r.eMut.push_back(tmpEp); - r.numEpitopes++; + species.push_back(T); // effector - std::cout << "\n"; // XXX + 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); + + for (size_t i = 0; i < T->num_ep; ++i) { // XXX + std::cout << T->count << '\t' + << T->affinity[i] << '\t' + << T->epitopes[i].epitopeWT << '|' + << T->epitopes[i].epitopeMut << '\n'; + } + + } } @@ -737,6 +796,7 @@ void usage() " -o (string) output file stem\n" " -s (string) input state file, containing initial population fraction\n" " -e (string) input file containing targeted epitope\n" +" -en (int) corresponding number of that CTL (ordered!)\n" " -n (int/double) population size\n" " -g (int/double) number of generations\n" " -mu (double) mutation rate\n" @@ -794,7 +854,8 @@ int main(int argc, char *argv[]) { else if (strcmp(argv[i],"-i")==0) { if (++i==argc) break; else r.infile=argv[i]; } else if (strcmp(argv[i],"-o")==0) { if (++i==argc) break; else r.outfile=argv[i]; } else if (strcmp(argv[i],"-s")==0) { if (++i==argc) break; else { r.statefile=argv[i]; r.importState=true; } } - else if (strcmp(argv[i],"-e")==0) { if (++i==argc) break; else { r.epitopefiles.push_back(argv[i]); r.useEpitope=true; } } + else if (strcmp(argv[i],"-e")==0) { if (++i==argc) break; else { r.epitopefiles.push_back(argv[i]); r.useEpitope=true; ++r.num_CTL_clones; } } + else if (strcmp(argv[i],"-en")==0) { if (++i==argc) break; else { r.init_CTL_num.push_back((long)strtodouble(argv[i])); r.useEpitope=true; } } else if (strcmp(argv[i],"-n")==0) { if (++i==argc) break; else r.n=(unsigned int)strtodouble(argv[i]); } else if (strcmp(argv[i],"-g")==0) { if (++i==argc) break; else r.g=(unsigned int)strtodouble(argv[i]); } diff --git a/ss.h b/ss.h index 747094a..be805c4 100644 --- a/ss.h +++ b/ss.h @@ -70,6 +70,28 @@ inline double strtodouble(const std::string &s) { } +// print a set of mutations to stream as a CSV-like string +std::ostream& operator<<(std::ostream& os, std::vector s) +{ + std::vector::iterator i = s.begin(), end = s.end(); + if (i == end) + return os; // stop immediately + os << *i; // print first element + for (++i; i != end; ++i) // then with comma as prefix + os << ',' << *i; + return os; +} +std::ostream& operator<<(std::ostream& os, MutatedSiteSequence s) +{ + MutatedSiteSequence::iterator i = s.begin(), end = s.end(); + if (i == end) + return os; // stop immediately + os << *i; // print first element + for (++i; i != end; ++i) // then with comma as prefix + os << ',' << *i; + return os; +} + // PROGRAM SETTINGS @@ -128,6 +150,8 @@ public: std::vector epitopefiles; + std::vector init_CTL_num; + unsigned int num_CTL_clones; RunParameters_SS() { @@ -162,8 +186,11 @@ public: write_mod = 1; numEpitopes = 0; + + num_CTL_clones = 0; } + void setFiles() { if (strcmp(directory.c_str(),"")!=0) { @@ -204,8 +231,8 @@ public: void run(RunParameters_SS &r); void importState(RunParameters_SS &r); -void importEpitope(RunParameters_SS &r); -void simulate(Rxn_parray &reactions, Species_parray &species, double t_end, long max_steps, const Species_parray &print_spec); +void importEpitope(RunParameters_SS &r, Species_parray &species, Rxn_parray &reactions, VirusSpecies *V); +void simulate(Rxn_parray &reactions, Species_parray &species, double t_end, long max_steps, double sample_interval, const Species_parray &print_spec); #endif