Detect when escape from *all* epitopes (vs *any* immune pressure).

diff --git a/hamiltonian.cpp b/hamiltonian.cpp
index bcedac6..40a3697 100644 (file)
--- a/hamiltonian.cpp
+++ b/hamiltonian.cpp
@@ -244,7 +244,7 @@ bool EpitopeHamiltonian::escaped(const std::set<unsigned int> &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<unsigned int> &eWT, const std::vector<unsigned int> &eMut) {
 
@@ -253,7 +253,7 @@ bool EpitopeHamiltonian::escaped(const Virus &v, const std::vector<unsigned int>
 }
 
 
-// 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<unsigned int> &mutated_sites, const std::vector<unsigned int> &eWT, const std::vector<unsigned int> &eMut) const {
 
@@ -280,6 +280,23 @@ bool EpitopeHamiltonian::escaped(const std::set<unsigned int> &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<unsigned int> &mutated_sites) const {
+    for (unsigned ep=0; ep<penalty.size(); ++ep) {
+        if (!escaped(mutated_sites, epitopeWT[ep], epitopeMut[ep])) {
+            return false;
+        }
+    }
+    return true;
+}
+
+
+
 // Return the energy for a given set of mutated sites
 // default: do include contribution from epitopes
 

diff --git a/hamiltonian.h b/hamiltonian.h
index 3d2faf4..a0ba314 100644 (file)
--- a/hamiltonian.h
+++ b/hamiltonian.h
@@ -31,6 +31,8 @@ public:
 
     virtual bool escaped(const Virus &v) { return false; }
     virtual bool escaped(const std::set<unsigned int> &mutated_sites) const { return false; }
+    virtual bool escaped_all(const Virus &v) { return false; }
+    virtual bool escaped_all(const std::set<unsigned int> &mutated_sites) const { return false; }
     virtual double get_energy(const std::set<unsigned int> &mutated_sites) const;
     
     void set_temp(double x)           { bh=x; bJ=x; }
@@ -55,6 +57,10 @@ public:
     bool escaped(const std::set<unsigned int> &mutated_sites) const;
     bool escaped(const Virus &v, const std::vector<unsigned int> &eWT, const std::vector<unsigned int> &eMut);
     bool escaped(const std::set<unsigned int> &mutated_sites, const std::vector<unsigned int> &eWT, const std::vector<unsigned int> &eMut) const;
+
+    bool escaped_all(const Virus &v);
+    bool escaped_all(const std::set<unsigned int> &mutated_sites) const;
+
     double get_energy(const std::set<unsigned int> &mutated_sites) const;
     template<bool include_epitope> double get_energy(const std::set<unsigned int> &mutated_sites) const;
 

diff --git a/population.cpp b/population.cpp
index 6844d11..0bc1835 100644 (file)
--- a/population.cpp
+++ b/population.cpp
@@ -302,6 +302,52 @@ unsigned int Population::escape_variant(Hamiltonian &H, std::set<unsigned int> &
 
 
 
+// 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<unsigned int> &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 (file)
--- 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<unsigned int> &mutant);
+
+    unsigned int compute_num_escaped_all(Hamiltonian &H);
+    bool escaped_all(Hamiltonian &H);
+    unsigned int escape_variant_all(Hamiltonian &H, std::set<unsigned int> &mutant);
         
 };
 

diff --git a/wf.cpp b/wf.cpp
index 0064dea..6d2a594 100644 (file)
--- 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<unsigned int> 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<unsigned int>::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 (file)
--- 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;