modelMers.r

In rocheben/ibmR: Individual-Based Model in R for phylodynamics studies

rm(list=ls())
library(ibmR)

times=seq(0,20,by=0.001);

model=new("scheduler")
#Adults camel population
model=addPopulation(model,nbIndividuals=1000,evtCtctRate=0.1,intraSpecCtctRate=0.1,lifespan=1/70,name="Camel adults");
#Juvenile camel population
model=addPopulation(model,nbIndividuals=1000,evtCtctRate=0.1,intraSpecCtctRate=0.1,lifespan=1/70,name="Camel juvenile");
#Environment
model=addEnvrionment(model,seasonality=sin(times*2*pi))
#Add pathogen
model=addPathogen(model,mutationRate=0.1,crossProtection=0,probaInfection=0.01,latencyPeriod=365/7,infectiousPeriod=365/7,recoveryPeriod=365/7)

#Modulation of probability infection according to infection history
modulationProbaInfectionHistory<-function(proba,pPathogen){
	return(proba);
}

#Individual step function
step<-function(){
	infArray
	popArray
	patArray
	indexPat=sample(1:length(patArray))
	for(k in indexPat){
		#Infection from direct contact with other populations
		proba=patArray[indexPat[k]]$probaInfection*(intraSpecCtctRate*length(infArray[infArray$pop==currentPop & infArray$pathogen==patArray[indexPart[k]]$id]))+sum((interSpecCtctRate*length(infArray[infArray$pop~=currentPop & infArray$pathogen==patArray[indexPart[k]]$id])));
		proba=modulationProbaInfectionHistory(proba,patArray[k]);
		if(rand()<proba){
			#Infection
			addPathogen(patArray[indexPat[k]]);
		}
		#Infection from Environment
		loadEvt=getLoad(patArray[k]);
		probaInfection=evtCtctRate*(loadEvt/(loadEvt+patArray[k]$loadNeeded))
		proba=modulationProbaInfectionHistory(proba,patArray[k]);
		if(rand()<proba){
			#Infection
			addPathogen(patArray[indexPat[k]]);
		}
	}
	#Moving from age class
	indexPop=sample(1:length(popArray))
	for(k in indexPop){
		proba=transitionRate[currentPop,indexPop[k]];
		if(rand()<proba){
			#Infection
			changePop(indexPop[k]);
		}
	}

	if(rand()<mutationRate){
		mutation();
	}
}

#Environmental step function
step<-function(){
	pathArray
	for(k in 1:length(pathArray)){
		pathArray[k]$Load=pathArray[k]$Load*(pathArray[k]$growth*seasonality(currentTime)-pathArray[k]$decay*seasonalityseasonality(currentTime));
	}
	if(rand()<mutationRate){
		mutation();
	}
}