import java.util.HashMap; /** * Created with IntelliJ IDEA. * User: rebeccamancy * Date: 26/10/2012 * Time: 14:50 * To change this template use File | Settings | File Templates. */ public class OH2000SIRSModel extends Model { private double eParam; private double cParam; private double alpha; private double rParam; // new recovery rate parameter /** * Constructor sets up model with three states (S-I-R) * @param alpha * @param cParam * @param eParam */ public OH2000SIRSModel(double eParam, double cParam, double alpha, double rParam) { super.setModelStates(3,1); // states are S-I-R or susceptible-infectious-recovered (immune) (S=0, I=1, R=2) this.eParam = eParam; this.cParam = cParam; this.alpha = alpha; this.rParam = rParam; } /** * Updates state at focal patch and rates at focal and all other patches * @param focalPatch patch that is changing state * @param patchStatus array of all patch status values (with focal patch having NEW state) * @param rates array of transition rates at all sites * @param hab habitat giving K and distances between sites * */ public double[] updateRates(int focalPatch, int[] patchStatus, double[] rates, Habitat hab) { double contrib = 0; // Set rate for focal patch and all other patches given new status if (patchStatus[focalPatch] == 0) { // new state is S = 0 // Update rate at focal patch rates[focalPatch] = getRate01(focalPatch, patchStatus, hab); // new colonisation rate at focal patch // Currently S ... no change to other patch rates as was previously R so couldn't infect then either /*for (int otherPatch = 0; otherPatch < hab.getNumPatches(); otherPatch++) { // loop over all patches if ((otherPatch != focalPatch) && (patchStatus[otherPatch] == 0)) { // all unoccupied other patches contrib = colonRate(alpha, hab.getK(focalPatch), hab.lookupDist(focalPatch, otherPatch) ); rates[otherPatch] -= contrib; //System.out.println("Decreasing rate at site " + otherPatch + " by " + contrib); //System.out.println(" with parameters dist=" + hab.lookupDist(focalPatch, otherPatch)); } }*/ } else if (patchStatus[focalPatch] == 1) { // new state is state I = 1 // Update rate at focal patch rates[focalPatch] = getRate12(hab.getK(focalPatch)); // new I-R transition rate at focal patch // Currently I, so INCREASE colonisation rate at other S patches by contribution of focal for (int otherPatch = 0; otherPatch < hab.getNumPatches(); otherPatch++) { // loop over all patches if ((otherPatch != focalPatch) && (canBeInfected(patchStatus[otherPatch])==1) ) { // all other S patches contrib = colonRate(alpha, hab.getK(focalPatch), hab.lookupDist(focalPatch, otherPatch) ); rates[otherPatch] += contrib; //System.out.println("Increasing rate at site " + otherPatch + " by " + contrib); //System.out.println(" with parameters dist=" + hab.lookupDist(focalPatch, otherPatch)); } } } else if (patchStatus[focalPatch] == 2) { // new state is R = 2 // Update rate at focal patch rates[focalPatch] = getRate20(hab.getK(focalPatch)); // No longer infectious (R) so DECREASE colonisation rate at other S patches by contribution of focal for (int otherPatch = 0; otherPatch < hab.getNumPatches(); otherPatch++) { // loop over all patches if ((otherPatch != focalPatch) && (canBeInfected(patchStatus[otherPatch])==1) ) { // all other S patches contrib = colonRate(alpha, hab.getK(focalPatch), hab.lookupDist(focalPatch, otherPatch) ); rates[otherPatch] -= contrib; //System.out.println("Decreasing rate at site " + otherPatch + " by " + contrib); //System.out.println(" with parameters dist=" + hab.lookupDist(focalPatch, otherPatch)); } } } return rates; } /** * Returns rate S-I (infection) for a specified patch * @param focPatch identifier of patch * @param hab habitat * @return (extinction) rate */ private double getRate01(int focPatch, int[] patchStatus, Habitat hab) { // Rate of S-I is given by sum over all I patches * ( size of patch * e^(alpha * dist) ) double rate = 0; for (int otherPatch = 0; otherPatch < hab.getNumPatches(); otherPatch++) { rate = rate + (canInfect(patchStatus[otherPatch]) * colonRate(this.alpha, hab.getK(otherPatch), hab.lookupDist(focPatch, otherPatch))); } rate = this.cParam * rate; return rate; } /** * Returns rate of I-R (recovery) for a specified patch given by eParam / K; * @param K * @return */ private double getRate12(double K) { // extinction of the disease (patch recovery rate) still size dependent todo?? return this.eParam / K; } /** * Returns rate of R-S (loss of patch immunity) for a specified patch given by rParam * K * @param K Size of the patch * @return rate of transition R-S */ private double getRate20(double K) { // transition back to susceptible state from recovered return this.rParam * K; // todo - what should this be? Here larger patches recover more quickly as higher number of births?? } /** * Returns 1 if can infect other sites (if Infectious); 0 otherwise * @param state * @return */ private int canInfect(int state) { if (state == 1) { return 1; // in OH2000SIRModel, can only colonise / infect if in state 1 } else { return 0; } } /** * Returns 1 if can be infected (i.e. if susceptible); 0 otherwise (so state 2=R returns 0 here as not infectious) * @param state * @return */ private int canBeInfected(int state) { if (state == 0) { return 1; // state 0 = S so can be infected } else { return 0; // otherwise can't be infected } } /** * Returns the rate of a patch given its ID, status and the habitat (called from Simulator) * @param focalPatch patch for which rate is to be calculated * @param patchStatus current patch status * @param hab habitat * @return */ public double getRate(int focalPatch, int[] patchStatus, Habitat hab) { int focalStatus = patchStatus[focalPatch]; if (focalStatus == 0) { return getRate01(focalPatch, patchStatus, hab); } else if (focalStatus == 1) { return getRate12(hab.getK(focalPatch)); } else { return getRate20(hab.getK(focalPatch)); } } /** * Calculates the contribution of the colonisation rate of a specific site to a site at a given distance (when occupied) * @param alpha alpha parameter * @param K carrying capacity of patch * @param dist distance between patches * @return contribution to colonisation rate */ private double colonRate(double alpha, double K, double dist) { return K * Math.exp(-alpha * dist); } /** * Calculate contribution rate from one site to another (colonisation) - for matC * @param hab * @param from * @param to * @return */ public double contribRate (Habitat hab, int from, int to) { double K = hab.getK(from); // it's the size of the from patch that counts in OH2000 double dist = hab.lookupDist(from,to); return K * Math.exp(-alpha * dist); } // todo this is probably wrong as doesn't include parameter c ... ! /** * Calculate extinction rate from a site (todo - to is not used here?) * @param hab * @param from * @param to * @return */ public double extRate (Habitat hab, int from, int to) { double K = hab.getK(from); return this.eParam/K; } private double getConnex(double dist) { return Math.exp(-alpha * dist); } /** * Calculates the appropriate m_{ij} value * @param focalPatch first patch * @param otherPatch other patch * @return the mij value for the landscape matrix M */ public double mijValue(int focalPatch, int otherPatch, Habitat hab) { if (focalPatch == otherPatch) { return 0.0; } else { return Math.exp(-alpha * hab.lookupDist(focalPatch,otherPatch)) * hab.getK(otherPatch) * hab.getK(focalPatch); } } public void seteParam(double eParam) { this.eParam = eParam; } public void setcParam(double cParam) { this.cParam = cParam; } public void setAlpha(double alpha) { this.alpha = alpha; } public double geteParam() { return this.eParam; } public double getcParam() { return this.cParam; } public double getAlpha() { return this.alpha; } }