Nothing
1 - Running a simple simulation
In landsepi: Landscape Epidemiology and Evolution
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
library(landsepi)
General presentation of the package
See ?lansdepi for a complete description of the model, assumptions and available functions.
See also the vignette list of parameters for a detailed description of all model parameters.
Take a quick overview of what landsepi can do.
Initialisation of simulation parameters
The function createSimulParams() will create a directory to store simulation outputs,
and return an object of class LandsepiParams that will further contain all simulation parameters.
simul_params <- createSimulParams(outputDir = getwd())
In the following, the object simul_params is to be updated with all simulation parameters using set*() functions.
To help parameterisation, built-in parameters are available and may be loaded via load*() functions.
Setting the seed and time parameters
A seed (for random number generator) is randomly generated when simul_params is initialised by
createSimulParams(). If a specific seed is required, it can be set using the function setSeed():
simul_params@Seed
simul_params <- setSeed(simul_params, seed = 1)
simul_params@Seed
The number of cropping seasons to simulate (e.g. 6 years) and the number of time steps per
cropping season (e.g. 120 days/year) can be set using setTime():
simul_params <- setTime(simul_params, Nyears = 6, nTSpY = 120)
simul_params@TimeParam
Setting pathogen parameters
Pathogen parameters must be stored in a list of aggressiveness components defined on a susceptible host for a
pathogen genotype not adapted to resistance.
Buit-in parameterisation for different pathogens (e.g. rusts of cereal crops, mildew of grapevine,
black sigatoka of banana) are available using function loadPathogen():
basic_patho_param <- loadPathogen(disease = "rust")
basic_patho_param
basic_patho_param <- loadPathogen(disease = "mildew")
basic_patho_param
basic_patho_param <- loadPathogen(disease = "sigatoka")
basic_patho_param
This list may be updated to change a specific parameter. For instance, to change the infection rate to 50%:
basic_patho_param <- loadPathogen("rust")
basic_patho_param$infection_rate <- 0.5
basic_patho_param
Alternatively, the list may be generated manually to control all parameters:
basic_patho_param <- list(infection_rate = 0.4
, latent_period_mean = 10
, latent_period_var = 9
, propagule_prod_rate = 3.125
, infectious_period_mean = 24
, infectious_period_var = 105
, survival_prob = 1e-4
, repro_sex_prob = 0
, sigmoid_kappa = 5.333
, sigmoid_sigma = 3
, sigmoid_plateau = 1
, sex_propagule_viability_limit = 1
, sex_propagule_release_mean = 1
, clonal_propagule_gradual_release = 0)
Then, the list is used to fuel the object simul_params via the function setPathogen():
simul_params <- setPathogen(simul_params, patho_params = basic_patho_param)
simul_params@Pathogen
Setting inoculum
The function setInoculum() updates simul_params with the probability for an individual to be infectious
(i.e. state I) at the beginning of the simulation. Note that the inoculum is present on the first cultivar
only (usually defined as the susceptible cultivar) and is not adapted to any plant resistance gene.
simul_params <- setInoculum(simul_params, val = 5e-4)
simul_params@PI0
Setting the landscape and pathogen dispersal
The landscape (i.e. boundaries of fields) must be in shapefile format. Five built-in landscapes of about 150
fields are available using the function loadLandscape():
landscape <- loadLandscape(id = 1)
length(landscape)
plot(landscape, main = "Landscape structure")
See also tutorial on how to parameterise landscape and dispersal
to use your own landscape and compute your own dispersal matrices.
Dispersal is given by a vectorised matrix giving the probability of dispersal from any field of the
landscape to any other field. The size of the matrix must be the square of the number of fields in the landscape.
It is thus specific to both the pathogen and the landscape. For rusts pathogens,
a built-in dispersal matrix is available for each landscape using the function loadDispersalPathogen():
disp_patho_clonal <- loadDispersalPathogen(id = 1)[[1]]
head(disp_patho_clonal)
length(landscape)^2 == length(disp_patho_clonal)
Then, the object simul_params is updated with the landscape and dispersal matrix via the
functions setLandscape() and setDispersalPathogen(), respectively:
simul_params <- setLandscape(simul_params, land = landscape)
simul_params <- setDispersalPathogen(simul_params, disp_patho_clonal)
Setting croptypes, cultivars and resistance genes
Fields of the landscape are cultivated with different croptypes that can rotate through time;
each croptype is composed of a pure cultivar or a mixture; and each cultivar may carry one or
several resistance genes.
Cultivars
Characteristics of each host genotype (i.e. cultivar) are summarised in a dataframe, which
contains parameters representing the cultivar as if it was cultivated in a pure crop.
Note that the name of the cultivar cannot accept spaces.
A buit-in parameterisation for classical types of cultivars is available using loadCultivar()
to generate each line of the dataframe. The only difference between the types "growingHost" and
"nongrowingHost" is the absence of growth in the later. Type "nonCrop" allows the simulation of
forest, fallows, etc. i.e. everything that is not planted, does not cost anything and does not
yield anything (with regard to the subject of the study).
It is advised to implement a susceptible cultivar at first line of the dataframe to allow
pathogen initial inoculation.
cultivar1 <- loadCultivar(name = "Susceptible", type = "growingHost")
cultivar2 <- loadCultivar(name = "Resistant1", type = "growingHost")
cultivar3 <- loadCultivar(name = "Resistant2", type = "growingHost")
cultivar4 <- loadCultivar(name = "Resistant3", type = "nongrowingHost")
cultivar5 <- loadCultivar(name = "Forest", type = "nonCrop")
cultivars <- data.frame(rbind(cultivar1, cultivar2, cultivar3, cultivar4, cultivar5)
, stringsAsFactors = FALSE)
cultivars
Similarly as pathogen parameters, characteristics of the cultivars may be updated as required.
For example, to change the growth rate of the susceptible cultivar:
cultivars[cultivars$cultivarName == "Susceptible", "growth_rate"] <- 0.2
cultivars
Finally, the dataframe cultivars can also be generated entirely from scratch:
cultivars_new <- data.frame(cultivarName = c("Susceptible", "Resistant"),
initial_density = c(0.1, 0.2),
max_density = c(2.0, 3.0),
growth_rate = c(0.1, 0.2),
reproduction_rate = c(0.0, 0.0),
yield_H = c(2.5, 2.0),
yield_L = c(0.0, 0.0),
yield_I = c(0.0, 0.0),
yield_R = c(0.0, 0.0),
planting_cost = c(225, 300),
market_value = c(200, 150),
stringsAsFactors = FALSE)
cultivars_new
Note, assuming that only healthy hosts (state H) contribute to host growth, the production of healthy biomass between $t$ and $t+1$ is computed using the following logistic equation:
$$H_{t+1} = H_{t} \times \left[1 + growth_rate \times (1-\frac{N_{t}}{K})\right]$$
with $N_t$ the total number of individual hosts at time-step $t$ and $K$ the carrying capacity.
Resistance genes
This part can be skipped if no resistance gene is to be simulated.
Characteristics of each plant resistance gene and each corresponding pathogenicity gene are summarised in a dataframe.
A built-in parameterisation for classical resistance sources is available using loadGene() to generate each line of the
dataframe. Type "majorGene" is for completely efficient resistance to which pathogen may adapt via a single mutation. Type "APR"
stands for Adult Plant Resistance, i.e. a major gene which activates after a delay after planting, and type "QTL" is a
partially efficient resistance source to which the pathogen may adapt gradually. The type "immunity" code for a completely
efficient resistance gene that cannot be overcome, in such a way that infection is totally impossible; it is helpful
to parameterise non-host cultivars.
gene1 <- loadGene(name = "MG 1", type = "majorGene")
gene2 <- loadGene(name = "Lr34", type = "APR")
gene3 <- loadGene(name = "gene 3", type = "QTL")
gene4 <- loadGene(name = "nonhost resistance", type = "immunity")
genes <- data.frame(rbind(gene1, gene2, gene3, gene4), stringsAsFactors = FALSE)
genes
Similarly as pathogen parameters, characteristics of the genes may be updated as required. For example, to change the
mutation probability of the pathogen with regard to its adaptation to "MG 1":
genes[genes$geneName == "MG 1", "mutation_prob"] <- 1e-3
genes
Alternatively, the dataframe genes can also be generated entirely from scratch:
genes_new <- data.frame(geneName = c("MG1", "MG2"),
efficiency = c(1.0 , 0.8 ),
age_of_activ_mean = c(0.0 , 0.0 ),
age_of_activ_var = c(0.0 , 0.0 ),
mutation_prob = c(1E-7 , 1E-4),
Nlevels_aggressiveness = c(2 , 2 ),
adaptation_cost = c(0.50 , 0.75 ),
tradeoff_strength = c(1.0 , 1.0 ),
target_trait = c("IR" , "LAT"),
recombination_sd = c(1.0,1.0),
stringsAsFactors = FALSE)
genes_new
Note that the column "efficiency" refers to the percentage of reduction of the targeted aggressiveness
component on hosts carrying the gene.
For example, a 80% efficient resistance against infection rate (target_trait = "IR")
means that the infection rate of a non-adapted pathogen is reduced by 80% on a resistant host compared
to what it is on a susceptible host. For resistances targetting the latent period,
the percentage of reduction is applied to the inverse of the latent period in such a way that latent
period is higher in resistant than in susceptible hosts.
Allocating genes to cultivars
The object simul_params can be updated with setGenes() and setCultivars():
simul_params <- setGenes(simul_params, dfGenes = genes)
simul_params <- setCultivars(simul_params, dfCultivars = cultivars)
simul_params@Genes
simul_params@Cultivars
Then the function allocateCultivarGenes() allows the attribution of resistance genes to cultivars:
simul_params <- allocateCultivarGenes(simul_params
, cultivarName = "Resistant1"
, listGenesNames = c("MG 1"))
simul_params <- allocateCultivarGenes(simul_params
, cultivarName = "Resistant2"
, listGenesNames = c("Lr34", "gene 3"))
simul_params <- allocateCultivarGenes(simul_params
, cultivarName = "Resistant3"
, listGenesNames = c("nonhost resistance"))
simul_params@Cultivars
With this example of parameterisation:
- "Susceptible" is a susceptible cultivar (initially infected by a wild-type pathogen)
- "Resistant1" is a mono-resistant cultivar
- "Resistant2" is a pyramided cultivar
- "Resistant3" is a nonhost cultivar
- "Forest" is not a crop.
Infection is impossible on both "Resistance3" and "Forest", but the former will
be considered for host growth, yield and planting costs, whereas the latter won't.
Allocating cultivars to croptypes
Characteristics of each croptype (a croptype is a set of hosts cultivated in a field with specific proportions)
are summarised in a dataframe.
A buit-in parameterisation for classical croptypes is available
using loadCroptypes() to generate the whole table filled with zeros:
croptypes <- loadCroptypes(simul_params, names = c("Susceptible crop"
, "Pure resistant crop"
, "Mixture"
, "Other"))
croptypes
Then croptypes is updated by allocateCroptypeCultivars()
to specify the composition of every croptype with regard to cultivars (and proportion for mixtures):
croptypes <- allocateCroptypeCultivars(croptypes
, croptypeName = "Susceptible crop"
, cultivarsInCroptype = "Susceptible")
croptypes <- allocateCroptypeCultivars(croptypes
, croptypeName = "Pure resistant crop"
, cultivarsInCroptype = "Resistant1")
croptypes <- allocateCroptypeCultivars(croptypes
, croptypeName = "Mixture"
, cultivarsInCroptype = c("Resistant2","Resistant3")
, prop = c(0.4, 0.6))
croptypes <- allocateCroptypeCultivars(croptypes
, croptypeName = "Other"
, cultivarsInCroptype = "Forest")
croptypes
Finally the object simul_params is updated using setCroptypes():
simul_params <- setCroptypes(simul_params, dfCroptypes = croptypes)
simul_params@Croptypes
Alternatively, the dataframe croptypes can be generated from scratch:
croptypes <- data.frame(croptypeID = c(0, 1, 2, 3) ## must start at 0 and match with values from landscape "croptypeID" layer
, croptypeName = c("Susceptible crop"
, "Pure resistant crop"
, "Mixture"
, "Other")
, Susceptible = c(1,0,0 ,0)
, Resistant1 = c(0,1,0 ,0)
, Resistant2 = c(0,0,0.5,0)
, Resistant3 = c(0,0,0.5,0)
, Forest = c(0,0,0 ,1)
, stringsAsFactors = FALSE)
simul_params <- setCroptypes(simul_params, croptypes)
Allocating croptypes to fields of the landscape
The function allocateLandscapeCroptypes() manages the allocation of croptypes in time (for rotations
of different croptypes on the same fields) and space (for mosaic of fields cultivated with different croptypes).
See ?allocateLandscapeCroptypes for help in parameters.
Briefly, a rotation sequence is defined by a list. Each
element of this list is a vector containing the indices of the croptypes that are cultivated simultaneously in the
landscape. The "rotation_period" parameter defines the duration before switching from one element of "rotation_sequence"
to the next. For each element of "rotation_sequence" is associated a vector of proportions of each croptype in the
landscape (parameter "prop"). The parameter "aggreg" controls the level of spatial aggregation between croptypes.
For example, to generate a landscape whose surface is composed of
1/3 of forests, 1/3 of susceptible crop, and 1/3 of fields where a pure resistant crop is alternated with a mixture
every two years:
# croptypeIDs cultivated in each element of the rotation sequence:
rotation_sequence <- list(c(0,1,3), c(0,2,3))
rotation_period <- 2 # number of years before rotation of the landscape
prop <- list(rep(1/3, 3), rep(1/3, 3)) # proportion (in surface) of each croptype
aggreg <-1 # level of spatial aggregation
simul_params <- allocateLandscapeCroptypes(simul_params
, rotation_period = rotation_period
, rotation_sequence = rotation_sequence
, prop = prop
, aggreg = aggreg
, graphic = TRUE)
# plot(simul_params@Landscape)
Setting chemical treatments
This part can be skipped if no chemical treatment is to be simulated.
Cultivars may be treated with chemicals which reduce the pathogen infection rate.
Treatment efficiency is maximum (i.e. equal to the parameter treatment_efficiency)
at the time of treatment application (noted $t^*$); then it decreases with time (i.e. natural pesticide degradation)
and host growth (i.e. new biomass is not protected by treatments):
$$IR(t) = IR_{max} \times \left(1 - \frac{treatment_efficiency}{1+e^{4.0- 8.5 \times C(t)}}\right)$$
$C(t) = C_1 \times C_2$ is the treatment concentration at $t$, which depends on:
- the timelag $\Delta t = t - t^$ passed since the time of treatment application $t^$
(i.e. natural pesticide degradation via the parameter
treatment_degradation_rate)
$$C_1 = e^{-treatment_degradation_rate \times \Delta t}$$
- the newly produced host biomass, which is not protected by treatments ($N(t^*)$ and $N(t)$ are host biomass at the
time of treatment application and at time $t$, respectively)
$$C_2 = min(1.0, N(t^)/N(t))$$
Cultivars to be treated with chemicals and dates of (possible) applications are defined with parameters treatment_cultivars and treatment_timesteps, and are the same for all polygons cultivated with the cultivars to be treated.
However, the chemicals are applied in a polygon only if the disease severity (i.e. $I(t^) / N(t^*)$) at the application date is higher than a given threshold, defined by treatment_application_threshold.
The cost of a single treatment application (monetary units/ha) is defined by treatment_cost and will impact the economic outputs.
Treatment parameters can be set via the function setTreatment:
treatment <- list(treatment_degradation_rate = 0.1,
treatment_efficiency = 0.8,
treatment_timesteps = seq(1,120,14) ,
treatment_cultivars = c(0),
treatment_cost = 0,
treatment_application_threshold = c(0.0))
simul_params <- setTreatment(simul_params, treatment)
Choosing output variables
Several epidemiological, evolutionary and economic outputs can be generated by landsepi and represented in
text files (writeTXT=TRUE), graphics (graphic=TRUE) and a video (videoMP4=TRUE).
See equation below, as well as ?epid_output and ?evol_output for details on the different output variables.
Possible epidemiological outputs include:
- "audpc" : Area Under Disease Progress Curve (average number of diseased hosts per square meter, $AUDPC_{v,y}$)
- "audpc_rel" : relative Area Under Disease Progress Curve (average proportion of diseased hosts
relative to the total number of existing hosts, $AUDPC_{v,y}^r$)
- "gla" : Green Leaf Area (average number of healthy hosts per square meter, $GLA_{v,y}$)
- "gla_rel" : relative Green Leaf Area (average proportion of healthy hosts relative to the
total number of existing hosts, $GLA_{v,y}^r$)
- "eco_yield" : total crop yield (in weight or volume units per ha, $Yield_{v,y}$)
- "eco_cost" : operational crop costs (in monetary units per ha, $Operational_cost_{v,y}$)
- "eco_product" : total crop products (in monetary units per ha, $Product_{v,y}$)
- "eco_margin" : margin (products - operational costs, in monetary units per ha, $Margin_{v,y}$)
- "contrib": contribution of pathogen genotypes to LIR dynamics ($Contrib_{p,v,y}$)
- "HLIR_dynamics", "H_dynamics", "L_dynamics", "IR_dynamics", "HLI_dynamics", etc.: Epidemic dynamics
related to the specified sanitary status (H, L, I or R and all their combinations). Graphics only,
works only if graphic=TRUE.
- "all" : compute all these outputs (default)
- "" : none of these outputs will be generated.
Possible evolutionary outputs, based on the computation of genotype frequencies ($freq(I)_{p,t}$) include:
- "evol_patho": evolution of pathogen genotypes
- "evol_aggr": evolution of pathogen aggressiveness (i.e. phenotype)
- "durability": durability of resistance genes
- "all": compute all these outputs (default)
- "": none of these outputs will be generated.
Equations
In the following equations, $H_{i,v,t}$, $L_{i,v,p,t}$, $I_{i,v,p,t}$ and $R_{i,v,p,t}$ respectively denote the number of healthy, latent, infectious and removed host individuals, respectively, in field $i$ ($i=1,…,J$), for cultivar $v$ ($v=1,…,V$) at time step $t$ ($t=1,…,T \times Y$ with $Y$ the total number of cropping seasons and $T$ the number of time-steps per season).
$AUDPC_{v,y} = \frac{\sum_{t=t^0(y)}^{t^f(y)} \sum_{i=1}^{J} \sum_{p=1}^{P} {I_{i,v,p,t}+R_{i,v,p,t}}}{T \times \sum_{i=1}^{J} A_i}$
$AUDPC^r_{v,y} = \frac{\sum_{t=t^0(y)}^{t^f(y)} \sum_{i=1}^{J} \sum_{p=1}^{P} {I_{i,v,p,t}+R_{i,v,p,t}}}{\sum_{t=t^0(y)}^{t^f(y)} \sum_{i=1}^{J} \left(H_{i,v,t}+\sum_{p=1}^{P} {L_{i,v,p,t}+I_{i,v,p,t}+R_{i,v,p,t}}\right)}$
$GLA_{v,y} = \frac{\sum_{t=t^0(y)}^{t^f(y)} \sum_{i=1}^{J} H_{i,v,t}}{T \times \sum_{i=1}^{J} A_i}$
$GLA^r_{v,y} = \frac{\sum_{t=t^0(y)}^{t^f(y)} \sum_{i=1}^{J} H_{i,v,t}}{\sum_{t=t^0(y)}^{t^f(y)} \sum_{i=1}^{J} \left(H_{i,v,t}+\sum_{p=1}^{P} {L_{i,v,p,t}+I_{i,v,p,t}+R_{i,v,p,t}}\right)}$
$Yield_{v,y} = \frac{\sum_{t=t^0(y)}^{t^f(y)} \sum_{i=1}^{J} \left( y_{H,v} \times H_{i,v,t} + \sum_{p=1}^{P} {y_{L,v} \times L_{i,v,p,t} + y_{I,v} \times I_{i,v,p,t} + y_{R,v} \times R_{i,v,p,t}} \right)} {\sum_{t=t^0(y)}^{t^f(y)} \sum_{i=1}^{J} H^{i,v,t}}$
with $y{H,v}$, $y_{L,v}$, $y_{I,v}$ and $y_{R,v}$ the theoretical yield of cultivar v in pure crop (in weight or volume unit/ha/season) associated with the sanitary statuses ‘H’, ‘L’, ‘I’ and ‘R’, respectively. $H^_{i,v,t}$ is the number of healthy hosts in a pure crop and in absence of disease.
$Operational_cost_{v,y}=planting_cost_v \times \frac{\sum_{i \in \Omega_{v,y}}(area_i \times prop_{v,i})} {\sum_{i \in \Omega_{v,y}}area_i} + treatment_cost \times \frac{\sum_{i \in \Omega_{v,y}}(TFI_{i,v,y} \times area_i \times prop_{v,i})} {\sum_{i \in \Omega_{v,y}}area_i}$
with $\Omega_{v,y}$ the set of polygon indices where cultivar $v$ is cultivated at year $y$, and $prop_{v,i}$ the proportion of cultivar $v$ in polygon $i$ (for mixtures). $TFI$ stands for the Treatment Frequency Indicator (number of treatment applications per ha).
$Product_{v,y}=yield_{v,y} \times market_value$
$Margin_{v,y} = Product_{v,y} - Operational_Cost_{v,y}$
$Contrib_{p,v,y} = \frac{\sum_{t=t^0(y)}^{t^f(y)} \sum_{i=1}^{J} {L_{i,v,p,t}+I_{i,v,p,t}+R_{i,v,p,t}}} {\sum_{t=t^0(y)}^{t^f(y)} \sum_{i=1}^{J} \sum_{p=1}^{P} {L_{i,v,p,t}+I_{i,v,p,t}+R_{i,v,p,t}}}$
$freq(I){p,t} = \frac{\sum{i=1}^{J} \sum_{v=1}^{V} I_{i,v,p,t}} {\sum_{p=1}^{P} \sum_{i=1}^{J} \sum_{v=1}^{V} I_{i,v,p,t}}$
With respect to evolutionary outputs, for each pathogen genotype (evol_patho) or phenotype (evol_aggr, note that different pathogen genotypes may lead to the same phenotype on a resistant host, i.e. level of aggressiveness), several computations are performed:
- appearance: time to first appearance (as propagule);
- R_infection: time to first true infection of a resistant host;
- R_invasion: time to invasion, when the number of infections of resistant hosts reaches \code{thres_breakdown}, above which the genotype or phenotype is unlikely to go extinct.
The value Nyears + 1 time step is used if the genotype or phenotype never appeared/infected/invaded.
Durability is defined as the time to invasion of completely adapted pathogen individuals.
Parameterisation
A list of outputs can be generated using loadOutputs():
outputlist <- loadOutputs(epid_outputs = "all", evol_outputs = "all")
outputlist
Among the elements of outputList, "audpc100S" is the audpc in a fully susceptible landscape (used as reference value for graphics).
If necessary, the function compute_audpc100S helps compute this value in a single 1-km^2 field:
audpc100S <- compute_audpc100S("rust", "growingHost", area=1E6)
audpc100S <- compute_audpc100S("mildew", "grapevine", area=1E6)
audpc100S <- compute_audpc100S("sigatoka", "banana", area=1E6, nTSpY=182)
Then simul_params can be updated via setOutputs():
simul_params <- setOutputs(simul_params, outputlist)
See also tutorial on how to run a numerical experimental design
to compute your own output variables and to run several simulations within an experimental design
Running the simulation
The functions checkSimulParams() and saveDeploymentStrategy() check simulation parameters and
save the object simul_params (which contains all parameters associated with the deployment strategy)
into a GPKG file, respectively.
checkSimulParams(simul_params)
simul_params <- saveDeploymentStrategy(simul_params)
Then, the function runSimul() launches the simulation. Use ?runSimul to get all available options.
runSimul(simul_params, graphic = TRUE, videoMP4 = FALSE)
system(paste("rm -rf ", simul_params@OutputDir))
Try the landsepi package in your browser
Any scripts or data that you put into this service are public.
landsepi documentation built on July 26, 2023, 5:36 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set( collapse = TRUE, comment = "#>" )
library(landsepi)
General presentation of the package
See ?lansdepi for a complete description of the model, assumptions and available functions.
See also the vignette list of parameters for a detailed description of all model parameters.
Take a quick overview of what landsepi can do.
Initialisation of simulation parameters
The function createSimulParams() will create a directory to store simulation outputs,
and return an object of class LandsepiParams that will further contain all simulation parameters.
simul_params <- createSimulParams(outputDir = getwd())
In the following, the object simul_params is to be updated with all simulation parameters using set*() functions.
To help parameterisation, built-in parameters are available and may be loaded via load*() functions.
Setting the seed and time parameters
A seed (for random number generator) is randomly generated when simul_params is initialised by
createSimulParams(). If a specific seed is required, it can be set using the function setSeed():
simul_params@Seed simul_params <- setSeed(simul_params, seed = 1) simul_params@Seed
The number of cropping seasons to simulate (e.g. 6 years) and the number of time steps per
cropping season (e.g. 120 days/year) can be set using setTime():
simul_params <- setTime(simul_params, Nyears = 6, nTSpY = 120) simul_params@TimeParam
Setting pathogen parameters
Pathogen parameters must be stored in a list of aggressiveness components defined on a susceptible host for a pathogen genotype not adapted to resistance.
Buit-in parameterisation for different pathogens (e.g. rusts of cereal crops, mildew of grapevine,
black sigatoka of banana) are available using function loadPathogen():
basic_patho_param <- loadPathogen(disease = "rust") basic_patho_param basic_patho_param <- loadPathogen(disease = "mildew") basic_patho_param basic_patho_param <- loadPathogen(disease = "sigatoka") basic_patho_param
This list may be updated to change a specific parameter. For instance, to change the infection rate to 50%:
basic_patho_param <- loadPathogen("rust") basic_patho_param$infection_rate <- 0.5 basic_patho_param
Alternatively, the list may be generated manually to control all parameters:
basic_patho_param <- list(infection_rate = 0.4 , latent_period_mean = 10 , latent_period_var = 9 , propagule_prod_rate = 3.125 , infectious_period_mean = 24 , infectious_period_var = 105 , survival_prob = 1e-4 , repro_sex_prob = 0 , sigmoid_kappa = 5.333 , sigmoid_sigma = 3 , sigmoid_plateau = 1 , sex_propagule_viability_limit = 1 , sex_propagule_release_mean = 1 , clonal_propagule_gradual_release = 0)
Then, the list is used to fuel the object simul_params via the function setPathogen():
simul_params <- setPathogen(simul_params, patho_params = basic_patho_param) simul_params@Pathogen
Setting inoculum
The function setInoculum() updates simul_params with the probability for an individual to be infectious
(i.e. state I) at the beginning of the simulation. Note that the inoculum is present on the first cultivar
only (usually defined as the susceptible cultivar) and is not adapted to any plant resistance gene.
simul_params <- setInoculum(simul_params, val = 5e-4) simul_params@PI0
Setting the landscape and pathogen dispersal
The landscape (i.e. boundaries of fields) must be in shapefile format. Five built-in landscapes of about 150
fields are available using the function loadLandscape():
landscape <- loadLandscape(id = 1) length(landscape) plot(landscape, main = "Landscape structure")
See also tutorial on how to parameterise landscape and dispersal to use your own landscape and compute your own dispersal matrices.
Dispersal is given by a vectorised matrix giving the probability of dispersal from any field of the
landscape to any other field. The size of the matrix must be the square of the number of fields in the landscape.
It is thus specific to both the pathogen and the landscape. For rusts pathogens,
a built-in dispersal matrix is available for each landscape using the function loadDispersalPathogen():
disp_patho_clonal <- loadDispersalPathogen(id = 1)[[1]] head(disp_patho_clonal) length(landscape)^2 == length(disp_patho_clonal)
Then, the object simul_params is updated with the landscape and dispersal matrix via the
functions setLandscape() and setDispersalPathogen(), respectively:
simul_params <- setLandscape(simul_params, land = landscape) simul_params <- setDispersalPathogen(simul_params, disp_patho_clonal)
Setting croptypes, cultivars and resistance genes
Fields of the landscape are cultivated with different croptypes that can rotate through time; each croptype is composed of a pure cultivar or a mixture; and each cultivar may carry one or several resistance genes.
Cultivars
Characteristics of each host genotype (i.e. cultivar) are summarised in a dataframe, which
contains parameters representing the cultivar as if it was cultivated in a pure crop.
Note that the name of the cultivar cannot accept spaces.
A buit-in parameterisation for classical types of cultivars is available using loadCultivar()
to generate each line of the dataframe. The only difference between the types "growingHost" and
"nongrowingHost" is the absence of growth in the later. Type "nonCrop" allows the simulation of
forest, fallows, etc. i.e. everything that is not planted, does not cost anything and does not
yield anything (with regard to the subject of the study).
It is advised to implement a susceptible cultivar at first line of the dataframe to allow
pathogen initial inoculation.
cultivar1 <- loadCultivar(name = "Susceptible", type = "growingHost") cultivar2 <- loadCultivar(name = "Resistant1", type = "growingHost") cultivar3 <- loadCultivar(name = "Resistant2", type = "growingHost") cultivar4 <- loadCultivar(name = "Resistant3", type = "nongrowingHost") cultivar5 <- loadCultivar(name = "Forest", type = "nonCrop") cultivars <- data.frame(rbind(cultivar1, cultivar2, cultivar3, cultivar4, cultivar5) , stringsAsFactors = FALSE) cultivars
Similarly as pathogen parameters, characteristics of the cultivars may be updated as required. For example, to change the growth rate of the susceptible cultivar:
cultivars[cultivars$cultivarName == "Susceptible", "growth_rate"] <- 0.2 cultivars
Finally, the dataframe cultivars can also be generated entirely from scratch:
cultivars_new <- data.frame(cultivarName = c("Susceptible", "Resistant"), initial_density = c(0.1, 0.2), max_density = c(2.0, 3.0), growth_rate = c(0.1, 0.2), reproduction_rate = c(0.0, 0.0), yield_H = c(2.5, 2.0), yield_L = c(0.0, 0.0), yield_I = c(0.0, 0.0), yield_R = c(0.0, 0.0), planting_cost = c(225, 300), market_value = c(200, 150), stringsAsFactors = FALSE) cultivars_new
Note, assuming that only healthy hosts (state H) contribute to host growth, the production of healthy biomass between $t$ and $t+1$ is computed using the following logistic equation: $$H_{t+1} = H_{t} \times \left[1 + growth_rate \times (1-\frac{N_{t}}{K})\right]$$ with $N_t$ the total number of individual hosts at time-step $t$ and $K$ the carrying capacity.
Resistance genes
This part can be skipped if no resistance gene is to be simulated.
Characteristics of each plant resistance gene and each corresponding pathogenicity gene are summarised in a dataframe.
A built-in parameterisation for classical resistance sources is available using loadGene() to generate each line of the
dataframe. Type "majorGene" is for completely efficient resistance to which pathogen may adapt via a single mutation. Type "APR"
stands for Adult Plant Resistance, i.e. a major gene which activates after a delay after planting, and type "QTL" is a
partially efficient resistance source to which the pathogen may adapt gradually. The type "immunity" code for a completely
efficient resistance gene that cannot be overcome, in such a way that infection is totally impossible; it is helpful
to parameterise non-host cultivars.
gene1 <- loadGene(name = "MG 1", type = "majorGene") gene2 <- loadGene(name = "Lr34", type = "APR") gene3 <- loadGene(name = "gene 3", type = "QTL") gene4 <- loadGene(name = "nonhost resistance", type = "immunity") genes <- data.frame(rbind(gene1, gene2, gene3, gene4), stringsAsFactors = FALSE) genes
Similarly as pathogen parameters, characteristics of the genes may be updated as required. For example, to change the mutation probability of the pathogen with regard to its adaptation to "MG 1":
genes[genes$geneName == "MG 1", "mutation_prob"] <- 1e-3 genes
Alternatively, the dataframe genes can also be generated entirely from scratch:
genes_new <- data.frame(geneName = c("MG1", "MG2"), efficiency = c(1.0 , 0.8 ), age_of_activ_mean = c(0.0 , 0.0 ), age_of_activ_var = c(0.0 , 0.0 ), mutation_prob = c(1E-7 , 1E-4), Nlevels_aggressiveness = c(2 , 2 ), adaptation_cost = c(0.50 , 0.75 ), tradeoff_strength = c(1.0 , 1.0 ), target_trait = c("IR" , "LAT"), recombination_sd = c(1.0,1.0), stringsAsFactors = FALSE) genes_new
Note that the column "efficiency" refers to the percentage of reduction of the targeted aggressiveness component on hosts carrying the gene. For example, a 80% efficient resistance against infection rate (target_trait = "IR") means that the infection rate of a non-adapted pathogen is reduced by 80% on a resistant host compared to what it is on a susceptible host. For resistances targetting the latent period, the percentage of reduction is applied to the inverse of the latent period in such a way that latent period is higher in resistant than in susceptible hosts.
Allocating genes to cultivars
The object simul_params can be updated with setGenes() and setCultivars():
simul_params <- setGenes(simul_params, dfGenes = genes) simul_params <- setCultivars(simul_params, dfCultivars = cultivars) simul_params@Genes simul_params@Cultivars
Then the function allocateCultivarGenes() allows the attribution of resistance genes to cultivars:
simul_params <- allocateCultivarGenes(simul_params , cultivarName = "Resistant1" , listGenesNames = c("MG 1")) simul_params <- allocateCultivarGenes(simul_params , cultivarName = "Resistant2" , listGenesNames = c("Lr34", "gene 3")) simul_params <- allocateCultivarGenes(simul_params , cultivarName = "Resistant3" , listGenesNames = c("nonhost resistance")) simul_params@Cultivars
With this example of parameterisation:
- "Susceptible" is a susceptible cultivar (initially infected by a wild-type pathogen)
- "Resistant1" is a mono-resistant cultivar
- "Resistant2" is a pyramided cultivar
- "Resistant3" is a nonhost cultivar
- "Forest" is not a crop.
Infection is impossible on both "Resistance3" and "Forest", but the former will
be considered for host growth, yield and planting costs, whereas the latter won't.
Allocating cultivars to croptypes
Characteristics of each croptype (a croptype is a set of hosts cultivated in a field with specific proportions)
are summarised in a dataframe.
A buit-in parameterisation for classical croptypes is available
using loadCroptypes() to generate the whole table filled with zeros:
croptypes <- loadCroptypes(simul_params, names = c("Susceptible crop" , "Pure resistant crop" , "Mixture" , "Other")) croptypes
Then croptypes is updated by allocateCroptypeCultivars()
to specify the composition of every croptype with regard to cultivars (and proportion for mixtures):
croptypes <- allocateCroptypeCultivars(croptypes , croptypeName = "Susceptible crop" , cultivarsInCroptype = "Susceptible") croptypes <- allocateCroptypeCultivars(croptypes , croptypeName = "Pure resistant crop" , cultivarsInCroptype = "Resistant1") croptypes <- allocateCroptypeCultivars(croptypes , croptypeName = "Mixture" , cultivarsInCroptype = c("Resistant2","Resistant3") , prop = c(0.4, 0.6)) croptypes <- allocateCroptypeCultivars(croptypes , croptypeName = "Other" , cultivarsInCroptype = "Forest") croptypes
Finally the object simul_params is updated using setCroptypes():
simul_params <- setCroptypes(simul_params, dfCroptypes = croptypes) simul_params@Croptypes
Alternatively, the dataframe croptypes can be generated from scratch:
croptypes <- data.frame(croptypeID = c(0, 1, 2, 3) ## must start at 0 and match with values from landscape "croptypeID" layer , croptypeName = c("Susceptible crop" , "Pure resistant crop" , "Mixture" , "Other") , Susceptible = c(1,0,0 ,0) , Resistant1 = c(0,1,0 ,0) , Resistant2 = c(0,0,0.5,0) , Resistant3 = c(0,0,0.5,0) , Forest = c(0,0,0 ,1) , stringsAsFactors = FALSE) simul_params <- setCroptypes(simul_params, croptypes)
Allocating croptypes to fields of the landscape
The function allocateLandscapeCroptypes() manages the allocation of croptypes in time (for rotations
of different croptypes on the same fields) and space (for mosaic of fields cultivated with different croptypes).
See ?allocateLandscapeCroptypes for help in parameters.
Briefly, a rotation sequence is defined by a list. Each
element of this list is a vector containing the indices of the croptypes that are cultivated simultaneously in the
landscape. The "rotation_period" parameter defines the duration before switching from one element of "rotation_sequence"
to the next. For each element of "rotation_sequence" is associated a vector of proportions of each croptype in the
landscape (parameter "prop"). The parameter "aggreg" controls the level of spatial aggregation between croptypes.
For example, to generate a landscape whose surface is composed of 1/3 of forests, 1/3 of susceptible crop, and 1/3 of fields where a pure resistant crop is alternated with a mixture every two years:
# croptypeIDs cultivated in each element of the rotation sequence: rotation_sequence <- list(c(0,1,3), c(0,2,3)) rotation_period <- 2 # number of years before rotation of the landscape prop <- list(rep(1/3, 3), rep(1/3, 3)) # proportion (in surface) of each croptype aggreg <-1 # level of spatial aggregation simul_params <- allocateLandscapeCroptypes(simul_params , rotation_period = rotation_period , rotation_sequence = rotation_sequence , prop = prop , aggreg = aggreg , graphic = TRUE) # plot(simul_params@Landscape)
Setting chemical treatments
This part can be skipped if no chemical treatment is to be simulated.
Cultivars may be treated with chemicals which reduce the pathogen infection rate.
Treatment efficiency is maximum (i.e. equal to the parameter treatment_efficiency)
at the time of treatment application (noted $t^*$); then it decreases with time (i.e. natural pesticide degradation)
and host growth (i.e. new biomass is not protected by treatments):
$$IR(t) = IR_{max} \times \left(1 - \frac{treatment_efficiency}{1+e^{4.0- 8.5 \times C(t)}}\right)$$ $C(t) = C_1 \times C_2$ is the treatment concentration at $t$, which depends on:
- the timelag $\Delta t = t - t^$ passed since the time of treatment application $t^$
(i.e. natural pesticide degradation via the parameter
treatment_degradation_rate)
$$C_1 = e^{-treatment_degradation_rate \times \Delta t}$$
- the newly produced host biomass, which is not protected by treatments ($N(t^*)$ and $N(t)$ are host biomass at the time of treatment application and at time $t$, respectively)
$$C_2 = min(1.0, N(t^)/N(t))$$
Cultivars to be treated with chemicals and dates of (possible) applications are defined with parameters treatment_cultivars and treatment_timesteps, and are the same for all polygons cultivated with the cultivars to be treated.
However, the chemicals are applied in a polygon only if the disease severity (i.e. $I(t^) / N(t^*)$) at the application date is higher than a given threshold, defined by treatment_application_threshold.
The cost of a single treatment application (monetary units/ha) is defined by treatment_cost and will impact the economic outputs.
Treatment parameters can be set via the function setTreatment:
treatment <- list(treatment_degradation_rate = 0.1, treatment_efficiency = 0.8, treatment_timesteps = seq(1,120,14) , treatment_cultivars = c(0), treatment_cost = 0, treatment_application_threshold = c(0.0)) simul_params <- setTreatment(simul_params, treatment)
Choosing output variables
Several epidemiological, evolutionary and economic outputs can be generated by landsepi and represented in
text files (writeTXT=TRUE), graphics (graphic=TRUE) and a video (videoMP4=TRUE).
See equation below, as well as ?epid_output and ?evol_output for details on the different output variables.
Possible epidemiological outputs include:
- "audpc" : Area Under Disease Progress Curve (average number of diseased hosts per square meter, $AUDPC_{v,y}$)
- "audpc_rel" : relative Area Under Disease Progress Curve (average proportion of diseased hosts relative to the total number of existing hosts, $AUDPC_{v,y}^r$)
- "gla" : Green Leaf Area (average number of healthy hosts per square meter, $GLA_{v,y}$)
- "gla_rel" : relative Green Leaf Area (average proportion of healthy hosts relative to the total number of existing hosts, $GLA_{v,y}^r$)
- "eco_yield" : total crop yield (in weight or volume units per ha, $Yield_{v,y}$)
- "eco_cost" : operational crop costs (in monetary units per ha, $Operational_cost_{v,y}$)
- "eco_product" : total crop products (in monetary units per ha, $Product_{v,y}$)
- "eco_margin" : margin (products - operational costs, in monetary units per ha, $Margin_{v,y}$)
- "contrib": contribution of pathogen genotypes to LIR dynamics ($Contrib_{p,v,y}$) - "HLIR_dynamics", "H_dynamics", "L_dynamics", "IR_dynamics", "HLI_dynamics", etc.: Epidemic dynamics related to the specified sanitary status (H, L, I or R and all their combinations). Graphics only, works only if graphic=TRUE.
- "all" : compute all these outputs (default)
- "" : none of these outputs will be generated.Possible evolutionary outputs, based on the computation of genotype frequencies ($freq(I)_{p,t}$) include:
- "evol_patho": evolution of pathogen genotypes
- "evol_aggr": evolution of pathogen aggressiveness (i.e. phenotype)
- "durability": durability of resistance genes
- "all": compute all these outputs (default)
- "": none of these outputs will be generated.
Equations
In the following equations, $H_{i,v,t}$, $L_{i,v,p,t}$, $I_{i,v,p,t}$ and $R_{i,v,p,t}$ respectively denote the number of healthy, latent, infectious and removed host individuals, respectively, in field $i$ ($i=1,…,J$), for cultivar $v$ ($v=1,…,V$) at time step $t$ ($t=1,…,T \times Y$ with $Y$ the total number of cropping seasons and $T$ the number of time-steps per season).
$AUDPC_{v,y} = \frac{\sum_{t=t^0(y)}^{t^f(y)} \sum_{i=1}^{J} \sum_{p=1}^{P} {I_{i,v,p,t}+R_{i,v,p,t}}}{T \times \sum_{i=1}^{J} A_i}$
$AUDPC^r_{v,y} = \frac{\sum_{t=t^0(y)}^{t^f(y)} \sum_{i=1}^{J} \sum_{p=1}^{P} {I_{i,v,p,t}+R_{i,v,p,t}}}{\sum_{t=t^0(y)}^{t^f(y)} \sum_{i=1}^{J} \left(H_{i,v,t}+\sum_{p=1}^{P} {L_{i,v,p,t}+I_{i,v,p,t}+R_{i,v,p,t}}\right)}$
$GLA_{v,y} = \frac{\sum_{t=t^0(y)}^{t^f(y)} \sum_{i=1}^{J} H_{i,v,t}}{T \times \sum_{i=1}^{J} A_i}$
$GLA^r_{v,y} = \frac{\sum_{t=t^0(y)}^{t^f(y)} \sum_{i=1}^{J} H_{i,v,t}}{\sum_{t=t^0(y)}^{t^f(y)} \sum_{i=1}^{J} \left(H_{i,v,t}+\sum_{p=1}^{P} {L_{i,v,p,t}+I_{i,v,p,t}+R_{i,v,p,t}}\right)}$
$Yield_{v,y} = \frac{\sum_{t=t^0(y)}^{t^f(y)} \sum_{i=1}^{J} \left( y_{H,v} \times H_{i,v,t} + \sum_{p=1}^{P} {y_{L,v} \times L_{i,v,p,t} + y_{I,v} \times I_{i,v,p,t} + y_{R,v} \times R_{i,v,p,t}} \right)} {\sum_{t=t^0(y)}^{t^f(y)} \sum_{i=1}^{J} H^{i,v,t}}$ with $y{H,v}$, $y_{L,v}$, $y_{I,v}$ and $y_{R,v}$ the theoretical yield of cultivar v in pure crop (in weight or volume unit/ha/season) associated with the sanitary statuses ‘H’, ‘L’, ‘I’ and ‘R’, respectively. $H^_{i,v,t}$ is the number of healthy hosts in a pure crop and in absence of disease.
$Operational_cost_{v,y}=planting_cost_v \times \frac{\sum_{i \in \Omega_{v,y}}(area_i \times prop_{v,i})} {\sum_{i \in \Omega_{v,y}}area_i} + treatment_cost \times \frac{\sum_{i \in \Omega_{v,y}}(TFI_{i,v,y} \times area_i \times prop_{v,i})} {\sum_{i \in \Omega_{v,y}}area_i}$ with $\Omega_{v,y}$ the set of polygon indices where cultivar $v$ is cultivated at year $y$, and $prop_{v,i}$ the proportion of cultivar $v$ in polygon $i$ (for mixtures). $TFI$ stands for the Treatment Frequency Indicator (number of treatment applications per ha).
$Product_{v,y}=yield_{v,y} \times market_value$
$Margin_{v,y} = Product_{v,y} - Operational_Cost_{v,y}$
$Contrib_{p,v,y} = \frac{\sum_{t=t^0(y)}^{t^f(y)} \sum_{i=1}^{J} {L_{i,v,p,t}+I_{i,v,p,t}+R_{i,v,p,t}}} {\sum_{t=t^0(y)}^{t^f(y)} \sum_{i=1}^{J} \sum_{p=1}^{P} {L_{i,v,p,t}+I_{i,v,p,t}+R_{i,v,p,t}}}$
$freq(I){p,t} = \frac{\sum{i=1}^{J} \sum_{v=1}^{V} I_{i,v,p,t}} {\sum_{p=1}^{P} \sum_{i=1}^{J} \sum_{v=1}^{V} I_{i,v,p,t}}$
With respect to evolutionary outputs, for each pathogen genotype (evol_patho) or phenotype (evol_aggr, note that different pathogen genotypes may lead to the same phenotype on a resistant host, i.e. level of aggressiveness), several computations are performed:
- appearance: time to first appearance (as propagule);
- R_infection: time to first true infection of a resistant host;
- R_invasion: time to invasion, when the number of infections of resistant hosts reaches \code{thres_breakdown}, above which the genotype or phenotype is unlikely to go extinct.
The value Nyears + 1 time step is used if the genotype or phenotype never appeared/infected/invaded.
Durability is defined as the time to invasion of completely adapted pathogen individuals.
Parameterisation
A list of outputs can be generated using loadOutputs():
outputlist <- loadOutputs(epid_outputs = "all", evol_outputs = "all") outputlist
Among the elements of outputList, "audpc100S" is the audpc in a fully susceptible landscape (used as reference value for graphics).
If necessary, the function compute_audpc100S helps compute this value in a single 1-km^2 field:
audpc100S <- compute_audpc100S("rust", "growingHost", area=1E6) audpc100S <- compute_audpc100S("mildew", "grapevine", area=1E6) audpc100S <- compute_audpc100S("sigatoka", "banana", area=1E6, nTSpY=182)
Then simul_params can be updated via setOutputs():
simul_params <- setOutputs(simul_params, outputlist)
See also tutorial on how to run a numerical experimental design to compute your own output variables and to run several simulations within an experimental design
Running the simulation
The functions checkSimulParams() and saveDeploymentStrategy() check simulation parameters and
save the object simul_params (which contains all parameters associated with the deployment strategy)
into a GPKG file, respectively.
checkSimulParams(simul_params) simul_params <- saveDeploymentStrategy(simul_params)
Then, the function runSimul() launches the simulation. Use ?runSimul to get all available options.
runSimul(simul_params, graphic = TRUE, videoMP4 = FALSE)
system(paste("rm -rf ", simul_params@OutputDir))
Try the landsepi package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.