From: Dariusz Murakowski Date: Thu, 17 Apr 2014 01:47:20 +0000 (-0400) Subject: Detect when escape from *all* epitopes (vs *any* immune pressure). X-Git-Url: http://src.murakowski.org/?a=commitdiff_plain;h=5e20b7c5f9edad5561367731a3d3b6b6cf0024a9;p=VirEvoDyn.git Detect when escape from *all* epitopes (vs *any* immune pressure). --- diff --git a/hamiltonian.cpp b/hamiltonian.cpp index bcedac6..40a3697 100644 --- a/hamiltonian.cpp +++ b/hamiltonian.cpp @@ -244,7 +244,7 @@ bool EpitopeHamiltonian::escaped(const std::set &mutated_sites) co } -// Return true if the virus has escaped from immune pressure +// Return true if the virus has escaped from immune pressure at particular epitope bool EpitopeHamiltonian::escaped(const Virus &v, const std::vector &eWT, const std::vector &eMut) { @@ -253,7 +253,7 @@ bool EpitopeHamiltonian::escaped(const Virus &v, const std::vector } -// Return true if the sequence has escaped from immune pressure +// Return true if the sequence has escaped from immune pressure at particular epitope bool EpitopeHamiltonian::escaped(const std::set &mutated_sites, const std::vector &eWT, const std::vector &eMut) const { @@ -280,6 +280,23 @@ bool EpitopeHamiltonian::escaped(const std::set &mutated_sites, co } +// Did virus escape from *all* immune pressure everywhere? +bool EpitopeHamiltonian::escaped_all(const Virus &v) { + return EpitopeHamiltonian::escaped_all(v.mutated_sites); +} + +// Did sequence escape from *all* immune pressure everywhere? +bool EpitopeHamiltonian::escaped_all(const std::set &mutated_sites) const { + for (unsigned ep=0; ep &mutated_sites) const { return false; } + virtual bool escaped_all(const Virus &v) { return false; } + virtual bool escaped_all(const std::set &mutated_sites) const { return false; } virtual double get_energy(const std::set &mutated_sites) const; void set_temp(double x) { bh=x; bJ=x; } @@ -55,6 +57,10 @@ public: bool escaped(const std::set &mutated_sites) const; bool escaped(const Virus &v, const std::vector &eWT, const std::vector &eMut); bool escaped(const std::set &mutated_sites, const std::vector &eWT, const std::vector &eMut) const; + + bool escaped_all(const Virus &v); + bool escaped_all(const std::set &mutated_sites) const; + double get_energy(const std::set &mutated_sites) const; template double get_energy(const std::set &mutated_sites) const; diff --git a/population.cpp b/population.cpp index 6844d11..0bc1835 100644 --- a/population.cpp +++ b/population.cpp @@ -302,6 +302,52 @@ unsigned int Population::escape_variant(Hamiltonian &H, std::set & +// Compute the number of escaped viruses in the population + +unsigned int Population::compute_num_escaped_all(Hamiltonian &H) { + + unsigned int i=0; + + for (virus_map::iterator iter=pop.begin(); iter!=pop.end(); ++iter) { if (H.escaped_all(iter->first)) i+=iter->second; } + + return i; + +} + + +// Check whether most of the population has escaped from *all* immune pressure + +bool Population::escaped_all(Hamiltonian &H) { + + if (2*compute_num_escaped_all(H)>N) return true; + else return false; + +} + + +// Compute the number of escaped viruses in the population + +unsigned int Population::escape_variant_all(Hamiltonian &H, std::set &mutant) { + + unsigned int i=0; + + for (virus_map::iterator iter=pop.begin(); iter!=pop.end(); ++iter) { + + if (H.escaped_all(iter->first) && (iter->second)>i) { + + i=iter->second; + mutant=iter->first.mutated_sites; + + } + + } + + return i; + +} + + + // Output population to file // two-site two-allele diff --git a/population.h b/population.h index 9374b9a..ebac257 100644 --- a/population.h +++ b/population.h @@ -36,6 +36,10 @@ public: unsigned int compute_num_escaped(Hamiltonian &H); bool escaped(Hamiltonian &H); unsigned int escape_variant(Hamiltonian &H, std::set &mutant); + + unsigned int compute_num_escaped_all(Hamiltonian &H); + bool escaped_all(Hamiltonian &H); + unsigned int escape_variant_all(Hamiltonian &H, std::set &mutant); }; diff --git a/wf.cpp b/wf.cpp index 0064dea..6d2a594 100644 --- a/wf.cpp +++ b/wf.cpp @@ -130,11 +130,18 @@ void run(RunParameters &r, unsigned seed) { } if (r.runUntilEscape && P.escaped(H)) break; + + if (r.runUntilEscape_all && P.escaped_all(H)) break; } std::set esc_var; - unsigned int esc_num=P.escape_variant(H, esc_var); + unsigned int esc_num = 0; + if (!r.runUntilEscape && !r.runUntilEscape_all) esc_num = P.escape_variant(H, esc_var); // default behavior for backwards compatibility + if (r.runUntilEscape_all) esc_num = P.escape_variant_all(H, esc_var); + if (r.runUntilEscape) esc_num = P.escape_variant(H, esc_var); + // should behave better if both are true + // but more likely to escape one than all, so that one is output fprintf(supout,"%d\t%d",i,esc_num); for (std::set::iterator iter=esc_var.begin(); iter!=esc_var.end(); ++iter) fprintf(supout,"\t%d",*iter); @@ -314,6 +321,7 @@ void usage() " -b (double) \"inverse temperature\" multiplier\n" " -r flag to use relative energy condition for survival probability\n" " -esc flag to run until escape observed (or max number of generations reached)\n" +" -esc_all flag to run until escape from *all* immune pressure (or max num gen)\n" " -v flag for verbose output\n" "\n" " -2site flag for two-site two-allele system\n" @@ -370,6 +378,7 @@ int main(int argc, char *argv[]) { else if (strcmp(argv[i],"-r")==0) { r.useRelative=true; } else if (strcmp(argv[i],"-esc")==0) { r.runUntilEscape=true; } + else if (strcmp(argv[i],"-esc_all")==0) { r.runUntilEscape_all=true; } else if (strcmp(argv[i],"-v")==0) { r.useVerbose=true; } else if (strcmp(argv[i],"-numruns")==0) { if (++i==argc) break; else r.num_runs=(unsigned int)strtodouble(argv[i]); } diff --git a/wf.h b/wf.h index a787b4f..d849c63 100644 --- a/wf.h +++ b/wf.h @@ -89,6 +89,7 @@ public: int estart; // Epitope start position int eend; // Epitope end position bool runUntilEscape; // If true, run the simulation until the population escapes + bool runUntilEscape_all; // If true, run the simulation until the population escapes *all* immune pressure bool importState; // If true, import state from a state file bool useEpitope; // Include selective pressure on an epitope bool useRelative; // If true, use relative energy for survival probability @@ -132,6 +133,7 @@ public: penalty = 10.0; runUntilEscape=false; + runUntilEscape_all=false; importState=false; useEpitope=false; useRelative=false;