Merge pull request #19 from trvrb/waning

Waning

Host.java

@@ -98,6 +98,23 @@ public void mutate() {
 		infection = mutV;
 	}
 
+	// remove random phenotype from host's immune profile, do nothing if empty
+	public void waneImmunity() {
+		int length = immuneHistory.length;
+		if (length > 0) {
+			int remove = Random.nextInt(0, length-1);
+			Phenotype[] newHistory = new Phenotype[length - 1];
+			int currentIndex = 0;
+			for (int i = 0; i < length; i++) {
+				if (i != remove) {
+					newHistory[currentIndex] = immuneHistory[i];
+					currentIndex++;
+				}
+			}
+			immuneHistory = newHistory;
+		}
+	}
+
 	// history methods
 	public Phenotype[] getHistory() {
 		return immuneHistory;

HostPopulation.java

@@ -243,6 +243,9 @@ public void stepForward() {
 		if (Parameters.transcendental) { 
 			loseImmunity(); 
 		}
+		if (Parameters.waning) {
+			waneImmunity();
+		}
 		mutate();
 		sample();
 
@@ -436,9 +439,20 @@ public void loseImmunity() {
 			}
 		}			
 	}	
+
+	// draw a Poisson distributed number of waning events
+	public void waneImmunity() {
+		// each host wanes at a per-day rate of waningRate
+		double totalWaningRate = getN() * Parameters.waningRate * Parameters.deltaT;		
+		int wanings = Random.nextPoisson(totalWaningRate);
+		for (int i = 0; i < wanings; i++) {
+			Host h = getRandomHost();
+			h.waneImmunity();
+		}	
+	}	
 
 	// draw a Poisson distributed number of mutations and mutate based upon this
-	// mutate should not impact other Virus's Phenotypes through reference
+	// mutation should not impact other Virus's Phenotypes through reference
 	public void mutate() {
 		// each infected mutates at a per-day rate of mu
 		double totalMutationRate = getI() * Parameters.muPhenotype * Parameters.deltaT;

Parameters.java

@@ -39,14 +39,14 @@ public class Parameters {
 	public static int[] initialNs = {1000000,1000000,1000000};	
 
 	// host parameters
-	public static double birthRate = 0.000091;				// in births per individual per day, 1/30 years = 0.000091
-	public static double deathRate = 0.000091;				// in deaths per individual per day, 1/30 years = 0.000091
+	public static double birthRate = 0.000091;					// in births per individual per day, 1/30 years = 0.000091
+	public static double deathRate = 0.000091;					// in deaths per individual per day, 1/30 years = 0.000091
 	public static boolean swapDemography = true;				// whether to keep overall population size constant
 
 	// epidemiological parameters
 	public static int initialI = 10;							// in individuals
-	public static int initialDeme = 2;						// index of deme where infection starts, 1..n
-	public static double initialPrR = 0.5; 					// as proportion of population
+	public static int initialDeme = 2;							// index of deme where infection starts, 1..n
+	public static double initialPrR = 0.5; 						// as proportion of population
 	public static double beta = 0.36; // 0.3					// in contacts per individual per day
 	public static double nu = 0.2; //0.2						// in recoveries per individual per day
 	public static double betweenDemePro = 0.0005;				// relative to within-deme beta	
@@ -64,6 +64,8 @@ public class Parameters {
 	// phenotype parameters
 	public static String phenotypeSpace = "geometric";			// options include: "geometric", "geometric3d", "geometric10d"
 	public static double muPhenotype = 0.005; 					// in mutations per individual per day
+	public static boolean waning = false;						// whether to allow waning of host immunity
+	public static double waningRate = 0.01;						// rate per day of a host removing a random phenotype from their immune history
 
 	// parameters specific to GeometricPhenotype
 	public static double smithConversion = 0.1;					// multiplier to distance to give cross-immunity
@@ -224,6 +226,12 @@ public static void load() {
 			if (map.get("muPhenotype") != null) {				
 				muPhenotype = (double) map.get("muPhenotype");
 			}
+			if (map.get("waning") != null) {			
+				waning = (boolean) map.get("waning");
+			}
+			if (map.get("waningRate") != null) {			
+				waningRate = (double) map.get("waningRate");
+			}		
 			if (map.get("smithConversion") != null) {				
 				smithConversion = (double) map.get("smithConversion");
 			}

parameters.yml

@@ -50,9 +50,12 @@ demeBaselines: [1., 1., 1.]                 # baseline of seasonality
 demeAmplitudes: [0.1, 0., 0.1]              # amplitude of seasonality
 demeOffsets: [0., 0., 0.5]                  # seasonal offset relative to the year
 
-                                            # phenotype parameters  
+                                            # immune parameters  
 phenotypeSpace: "geometric"                 # options include: geometric, geometric3d, geometric10d
 muPhenotype: 0.005                          # in mutations per individual per day
+waning: false                               # whether to allow waning of host immunity
+waningRate: 0.01                            # rate per day of a host removing a random phenotype from their immune history
+
 smithConversion: 0.1                        # multiplier to distance to give cross-immunity
 homologousImmunity: 0.95                    # immunity raised to antigenically identical virus
 initialTraitA: -6.                          # value in dimension 1 for initial host immunity