pExit, param.pExit, timeDep.pExit, diff.pExit and hostCount.pExit

• pExit is the first required fundamental parameter and provides a daily probability for a host to leave the simulation (either cured, died, etc.).
• param.pExit is the list of functions needed to individually parameterize pExit (see Get started). The name of each function in the list has to match the name of the parameter it is sampling for pExit.
• timeDep.pExit allows for pExit to be dependent on the “absolute” time of the simulation, to account - for example - for seasonality or other external time-related covariates. By default, timeDep.pExit is set to FALSE.
• diff.pExit allows pExit to differ according to the current discrete state of the host. This can be useful, for example, if one state has a higher mortality rate (or better cures!) for the infection, in that case the probability to exit the simulation is higher. By default, diff.pExit is set to FALSE. Be careful, every state should give back a result for pExit.
• hostCount.pExit allows pExit to differ according to the number of hosts currently in a state. By default, hostCount.pExit is set to FALSE. To use hostCount.pExit, diff.pExit has to be set to TRUE too.

pMove, param.pMove, timeDep.pMove, diff.pMove and hostCount.pMove

• pMove is the probability (per unit of time) for a host to move on the map. It should not be confused with the distance the host will travel linked to sdMove (see below).
• param.pMove is the list of functions needed to individually parameterize pMove (see Get started). The name of each function in the list has to match the name of the parameter it is sampling for pMove.
• timeDep.pMove allows for pMove to be dependent of the “absolute” time of the simulation, to account, for example, for seasonality or other external time related covariates. By default, timeDep.pMove is set to FALSE.
• diff.pMove allows pMove to be different according to the current environmental value of the location where the host is located. By default, diff.pMove is set to FALSE.
• hostCount.pMove allows pMove to differ according to the number of hosts currently in a raster cell. By default, hostCount.pMove is set to FALSE. To use hostCount.pMove, diff.pMove has to be set to TRUE too.

sdMove, param.sdMove, timeDep.sdMove, diff.sdMove and hostCount.sdMove

• sdMove represents the standard deviation of the Brownian motion movement when a host undergo a move (set by pMove).
• param.sdMove is the list of function need to individually parameterize sdMove (see Get started). The name of each function in the list has to have the same name as the parameter it is sampling for sdMove.
• timeDep.sdMove allows for sdMove to be dependent of the “absolute” time of the simulation, to account, for example, for seasonality or other external time related covariates. By default, timeDep.sdMove is set to FALSE.
• diff.sdMove allows sdMove to be different according to the current environmental value of the location where the host is located. By default, diff.sdMove is set to FALSE.
• hostCount.sdMove allows sdMove to differ according to the number of hosts currently in a raster cell. By default, hostCount.sdMove is set to FALSE. To use hostCount.sdMove, diff.sdMove has to be set to TRUE too.

nContact, param.nContact, timeDep.nContact, diff.nContact and hostCount.nContact

• nContact represents the number (expressed as a positive integer) of potentially infectious contacts an infected hosts can encounter per unit of time. At each time point, a number of contacts will be determined for each active host in the simulation. The number of contacts (i.e. the output of your function) has to be an integer and can be set to zero.
• param.nContact is the list of functions needed to individually parameterize nContact (see Get started). The name of each function in the list has to match the name of the parameter it is sampling for nContact.
• timeDep.nContact allows for nContact to be dependent on the “absolute” time of the simulation, to account - for example - for seasonality or other external time-related covariates. By default, timeDep.nContact is set to FALSE.
• diff.nContact allows for nContact to differ according to the current discrete state of the host. By default, diff.nContact is set to FALSE.
• hostCount.nContact allows nContact to differ according to the number of hosts currently in a raster cell. This can be useful to adjust the number of contact to the number of potentially susceptible hosts if the infected population is close to the maximum size of the population in a raster cell. By default, hostCount.nContact is set to FALSE. To use hostCount.nContact, diff.nContact has to be set to TRUE too.

pTrans, param.pTrans, timeDep.pTrans,diff.pTrans and hostCount.pTrans

• pTrans is the heart of the transmission process and represents the probability of transmission over time (when a contact occurs).
• param.pTrans is the list of functions needed to individually parameterize pTrans (see Get started). The name of each function in the list has to match the name of the parameter it is sampling for pTrans.
• timeDep.pTrans allows for pTrans to be dependent on the “absolute” time of the simulation, to account - for example - for seasonality or other external time-related covariates. By default, timeDep.pTrans is set to FALSE.
• diff.pTrans allows for pTrans to be different according to the current discrete state of the host. This can be used to account of different dynamics linked to external factors common within a state, such as temperature for example. By default, diff.pTrans is set to FALSE.
• hostCount.pTrans allows pTrans to differ according to the number of hosts currently in a raster cell. By default, hostCount.pTrans is set to FALSE. To use hostCount.pTrans, diff.pTrans has to be set to TRUE too.

pExit

For pExit, we choose a probability that only depends on the environmental value of the location where the host currently is. The maximum value of our environmental raster is close to 70, and the higher the environmental value is, the less likely the host is to die:

p_Exit_fct  <- function(t, current.env.value){
  if(current.env.value > 60){p=0.02}
  if(current.env.value < 60 && current.env.value > 30){p=0.04}
  if(current.env.value < 30){p=0.08}
  return(p)
}

Taken together, the value for pExit according to the environmental value will be the following:

Remember that pExit, like the other core functions has to be function of t, even if t is not used. Since pExit is dependent on the location, diff.pMove=TRUE. However, there is no use of the “absolute” time of the simulation nor individual-based parameters, so param.pExit=NA, and timeDep.pExit=FALSE.

pMove

We choose a constant value for pMove, namely 0.1, i.e. an infected host has 10% chance to change its location for each unit of time.

p_Move_fct  <- function(t){return(0.1)}

Remember that pMove, like the other core functions has to be function of t, even if t is not used. Since pMove is not dependent on the location, diff.pMove=FALSE. Similarly, there is no use of the “absolute” time of the simulation nor individual-based parameters, so param.pMove=NA, and timeDep.pMove=FALSE.

sdMove

We choose a constant value for sdMove, namely 0.25. sdMove is the standard deviation of the Brownian motion underlying the movement in nosoi.

sd_Move_fct  <- function(t){return(0.25)}

Remember that sdMove, like the other core functions has to be function of t, even if t is not used. Since sdMove is not dependent on the location, diff.pMove=FALSE. Similarly, there is no use of the “absolute” time of the simulation nor individual-based parameters, so param.sdMove=NA, and timeDep.sdMove=FALSE.

nContact

For nContact, we choose a constant function that will draw a value in a normal distribution of mean = 0.5 and sd = 1, round it, and take its absolute value.

n_contact_fct <- function(t){abs(round(rnorm(1, 0.5, 1), 0))}

The distribution of nContact looks as follows:

At each time and for each infected host, nContact will be drawn anew. Remember that nContact, like the other core functions has to be function of t, even if t is not used. Since nContact is constant here, there is no use of the “absolute” time of the simulation, the location of the host, nor individual-based parameters. So param.nContact=NA, timeDep.nContact=FALSE and diff.nContact=FALSE.

pTrans

We choose pTrans in the form of a threshold function: before a certain amount of time since initial infection, the host does not transmit (incubation time, which we call t_incub), and after that time, it will transmit with a certain (constant) probability (which we call p_max). This function is dependent of the time since the host’s infection t.

p_Trans_fct <- function(t, p_max, t_incub){
    if(t < t_incub){p=0}
    if(t >= t_incub){p=p_max}
    return(p)
}

Because each host is different (slightly different biotic and abiotic factors), you can expect each host to exhibit differences in the dynamics of infection, and hence the probability of transmission over time. Thus, t_incub and p_max will be sampled for each host individually according to a certain distribution. t_incub will be sampled from a normal distribution of mean = 7 and sd = 1, while p_max will be sampled from a beta distribution with shape parameters α = 5 and β = 2:

t_incub_fct <- function(x){rnorm(x,mean = 7,sd=1)}
p_max_fct <- function(x){rbeta(x,shape1 = 5,shape2=2)}

Note that here t_incub and p_max are functions of x and not t (they are not core functions but individual-based parameters), and x enters the function as the number of draws to make.

Taken together, the profile for pTrans for a subset of 200 individuals in the population will look as follows:

pTrans is not dependent on the “absolute” time of the simulation, nor on the hosts location, so timeDep.pTrans=FALSE and diff.pTrans=FALSE. However, since we make use of individual-based parameters, we have to provide a param.pTrans as a list of functions. The name of each element within this list should have the same name that the core function (here pTrans) uses as argument, e.g.:

t_incub_fct <- function(x){rnorm(x,mean = 7,sd=1)}
p_max_fct <- function(x){rbeta(x,shape1 = 5,shape2=2)}

param_pTrans = list(p_max=p_max_fct,t_incub=t_incub_fct)

Running

Once nosoiSim is set up, you can run the simulation (here the “seed” ensures that you will get the same results as in this tutorial).

library(nosoi)

#Raster is test.raster

#Starting position will be
start.pos <- c(0,0) # c(x,y)

#pExit
p_Exit_fct  <- function(t, current.env.value){
  if(current.env.value > 60){p=0.02}
  if(current.env.value < 60 && current.env.value > 30){p=0.04}
  if(current.env.value < 30){p=0.08}
  return(p)
}

#pMove
p_Move_fct  <- function(t){return(0.1)}

#sdMove
sd_Move_fct  <- function(t){return(0.25)}

#nContact
n_contact_fct = function(t){abs(round(rnorm(1, 0.5, 1), 0))}

#pTrans
proba <- function(t,p_max,t_incub){
  if(t <= t_incub){p=0}
  if(t >= t_incub){p=p_max}
  return(p)
}

t_incub_fct <- function(x){rnorm(x,mean = 5,sd=1)}
p_max_fct <- function(x){rbeta(x,shape1 = 5,shape2=2)}

param_pTrans = list(p_max=p_max_fct,t_incub=t_incub_fct)

# Starting the simulation ------------------------------------

set.seed(846)
SimulationSingle <- nosoiSim(type="single", popStructure="continuous",
                             length.sim=300, max.infected=300, init.individuals=1, 
                             
                             init.structure=start.pos, 
                             
                             structure.raster=test.raster,
                             
                             pExit = p_Exit_fct,
                             param.pExit = NA,
                             timeDep.pExit=FALSE,
                             diff.pExit=TRUE,
                             
                             pMove = p_Move_fct,
                             param.pMove = NA,
                             timeDep.pMove=FALSE,
                             diff.pMove=FALSE,
                             
                             sdMove = sd_Move_fct,
                             param.sdMove = NA,
                             timeDep.sdMove=FALSE,
                             diff.sdMove=FALSE,
                             
                             attracted.by.raster=FALSE,
                             
                             nContact=n_contact_fct,
                             param.nContact=NA,
                             timeDep.nContact=FALSE,
                             diff.nContact=FALSE,
                             
                             pTrans = proba,
                             param.pTrans = list(p_max=p_max_fct,t_incub=t_incub_fct),
                             timeDep.pTrans=FALSE,
                             diff.pTrans=FALSE,
                             
                             prefix.host="H",
                             print.progress=FALSE,
                             print.step=10)
#> Starting the simulation
#> Initializing ... running ...
#> done. 
#> The simulation has run for 40 units of time and a total of 350 hosts have been infected.

Once the simulation has finished, it reports the number of time units for which the simulation has run (40), and the maximum number of infected hosts (350). Note that the simulation has stopped here before reaching length.sim as it has crossed the max.infected threshold set at 300.

pExit.B

For pExit.B, we choose a value that depends on the “absolute” time of the simulation, for example cyclic climatic conditions (temperature). In that case, the function’s arguments should be t and prestime (the “absolute” time of the simulation), in that order:

p_Exit_fctB  <- function(t,prestime){(sin(prestime/(2*pi*10))+1)/16} #for a periodic function

The values of pExit.B across the “absolute time” of the simulation will be the following:

Since pExit.B is dependent on the simulation’s time, do not forget to set timeDep.pExit.B to TRUE. Since there are no individual-based parameters nor is there any influence of the host’s location, param.pExit.B=NA and diff.pExit.B=NA.

pMove.B

We assume here that the hosts B do not move. pMove.B will hence be set to NA.

p_Move_fct.B  <- NA

Since pMove.B is not dependent on the location, diff.pMove.B=FALSE. Similarly, there is no use of the “absolute” time of the simulation nor individual-based parameters, so param.pMove.B=NA, and timeDep.pMove.B=FALSE.

sdMove.B

Since hosts B do not move, sdMove.B is irrelevant and will be set to NA.

sd_Move_fct.B  <- NA

Since sdMove.B is not dependent on the location, diff.sdMove.B=FALSE. Similarly, there is no use of the “absolute” time of the simulation nor individual-based parameters, so param.sdMove.B=NA, and timeDep.sdMove.B=FALSE.

nContact.B

For nContact.B, we choose a constant function that will sample a value from a provided list of values, with a certain probability:

n_contact_fct.B = function(t){sample(c(0,1,2),1,prob=c(0.6,0.3,0.1))}

The distribution of nContact.B looks as follows:

At each time and for each infected host, nContact.B will be drawn anew. Remember that nContact.B, like the other core functions has to be function of t, even if t is not used. Since nContact.B is constant here, there is no use of the “absolute” time of the simulation, the host’s location, nor individual-based parameters. So param.nContact.B=NA, timeDep.nContact.B=FALSE and diff.nContact.B=FALSE.

pTrans.B

We choose pTrans.B in the form of a Gaussian function. It will reach its maximum value at a certain time point (mean) after initial infection and will subsequently decrease until it reaches 0.

The function dnorm used here will achieve this objective: for each time t after infection, its return value will reach its maximum value (< 1) at its mean (here max.time, a parameter that will be set for each individual, see below) and then decline back to 0. Its increase and decline is parameterized by the standard deviation sd of the dnorm function.

p_Trans_fct.B <- function(t, max.time){
  dnorm(t, mean=max.time, sd=2)*5
}

Because each host is different (slightly different biotic and abiotic factors), you can expect each host to exhibit differences in the dynamics of infection, and hence the probability of transmission over time. Thus, max.time will be sampled for each host individually according to a certain distribution. max.time will be sampled from a normal distribution of parameters mean = 5 and sd = 1:

max.time_fct <- function(x){rnorm(x,mean = 5,sd=1)}

Note again that here max.time is a function of x and not t (i.e. not a core function but individual-based parameters), and x enters the function as the number of draws to make.

Taken together, the profile for pTrans for a subset of 200 individuals in the population will look as follow:

Since pTrans.B is not dependent of the “absolute” time of the simulation nor the host’s location, timeDep.pTrans.B=FALSE and diff.pTrans.B=FALSE. However, since we have made use of individual-based parameters, we have to provide a param.pTrans as a list of functions. The names of each element of the list should have the same name that the core function (pTrans.B here) uses as its argument, as shown here:

max.time_fct <- function(x){rnorm(x,mean = 5,sd=1)}

param_pTrans.B = list(max.time=max.time_fct)

Running

Once nosoiSim is set up, you can run the simulation (here the “seed” ensures that you will get the same results as in this tutorial).

library(nosoi)

#Raster is test.raster

#Starting position will be
start.pos <- c(0,0) # c(x,y)

#Host A -----------------------------------

#pExit
p_Exit_fct  <- function(t, current.env.value){
  if(current.env.value > 60){p=0.02}
  if(current.env.value < 60 && current.env.value > 30){p=0.04}
  if(current.env.value < 30){p=0.08}
  return(p)
}

#pMove
p_Move_fct  <- function(t){return(0.1)}

#sdMove
sd_Move_fct  <- function(t){return(0.25)}

#nContact
n_contact_fct = function(t){abs(round(rnorm(1, 0.5, 1), 0))}

#pTrans
proba <- function(t,p_max,t_incub){
  if(t <= t_incub){p=0}
  if(t >= t_incub){p=p_max}
  return(p)
}

t_incub_fct <- function(x){rnorm(x,mean = 5,sd=1)}
p_max_fct <- function(x){rbeta(x,shape1 = 5,shape2=2)}

param_pTrans = list(p_max=p_max_fct,t_incub=t_incub_fct)

#Host B -----------------------------------

#pExit
p_Exit_fct.B  <- function(t,prestime){(sin(prestime/(2*pi*10))+1)/16}

#pMove
p_Move_fct.B  <- NA

#pMove
sd_Move_fct.B  <- NA

#nContact
n_contact_fct.B = function(t){sample(c(0,1,2),1,prob=c(0.6,0.3,0.1))}

#pTrans
p_Trans_fct.B <- function(t, max.time){
  dnorm(t, mean=max.time, sd=2)*5
}
 
max.time_fct <- function(x){rnorm(x,mean = 5,sd=1)}

param_pTrans.B = list(max.time=max.time_fct)

# Starting the simulation ------------------------------------

set.seed(60)
SimulationDual <- nosoiSim(type="dual", popStructure="continuous",
                        length.sim=300, 
                        max.infected.A=100,
                        max.infected.B=200,
                        init.individuals.A=1,
                        init.individuals.B=0,
                        init.structure.A=start.pos,
                        init.structure.B=NA,
                        structure.raster.A=test.raster,
                        structure.raster.B=test.raster,

                        pExit.A = p_Exit_fct,
                        param.pExit.A = NA,
                        timeDep.pExit.A=FALSE,
                        diff.pExit.A=TRUE,
                        
                        pMove.A = p_Move_fct,
                        param.pMove.A = NA,
                        timeDep.pMove.A=FALSE,
                        diff.pMove.A=FALSE,
                        
                        sdMove.A = sd_Move_fct,
                        param.sdMove.A = NA,
                        timeDep.sdMove.A=FALSE,
                        diff.sdMove.A=FALSE,
                        attracted.by.raster.A=FALSE,
                        
                        nContact.A=n_contact_fct,
                        param.nContact.A=NA,
                        timeDep.nContact.A=FALSE,
                        diff.nContact.A=FALSE,

                        pTrans.A = proba,
                        param.pTrans.A = list(p_max=p_max_fct,t_incub=t_incub_fct),
                        timeDep.pTrans.A=FALSE,
                        diff.pTrans.A=FALSE,
                        prefix.host.A="H",
                        
                        pExit.B = p_Exit_fct.B,
                        param.pExit.B = NA,
                        timeDep.pExit.B=TRUE,
                        diff.pExit.B=FALSE,
                        
                        pMove.B = p_Move_fct.B,
                        param.pMove.B = NA,
                        timeDep.pMove.B=FALSE,
                        diff.pMove.B=FALSE,
                        
                        sdMove.B = sd_Move_fct.B,
                        param.sdMove.B = NA,
                        timeDep.sdMove.B=FALSE,
                        diff.sdMove.B=FALSE,
                        attracted.by.raster.B=FALSE,
                        
                        nContact.B=n_contact_fct.B,
                        param.nContact.B=NA,
                        timeDep.nContact.B=FALSE,
                        diff.nContact.B=FALSE,

                        pTrans.B = p_Trans_fct.B,
                        param.pTrans.B = param_pTrans.B,
                        timeDep.pTrans.B=FALSE,
                        diff.pTrans.B=FALSE,
                        prefix.host.B="V",
                        
                        print.progress=FALSE)
#> Starting the simulation
#> Initializing ... running ...
#> done. 
#> The simulation has run for 33 units of time and a total of 102 (A) and 137 (B) hosts have been infected.

Once the simulation has finished, it reports the number of time units for which the simulation has run (33), and the maximum number of infected hosts A (102) and hosts B (137). Note that the simulation has stopped here before reaching length.sim as it has crossed the max.infected.A threshold set at 100.