library(landsepi)
See ?lansdepi
for a complete description of the model, assumptions and available functions. See also the detailed list of parameters for a detailed description of all model parameters.
Take a quick overview of what landsepi can do.
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.
createSimulParams(outputDir = getwd())
simul_params <-#> Created output directory : /tmp/RtmpB2hmhu/Rbuild2a0227a8993/landsepi/vignettes/simul_landsepi_2022-03-01_16-38-45
#> Seed for RGN set to : 1646149125
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.
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()
:
@Seed
simul_params#> [1] 1646149125
setSeed(simul_params, seed = 1)
simul_params <-#> Seed for RGN set to : 1
@Seed
simul_params#> [1] 1
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()
:
setTime(simul_params, Nyears = 6, nTSpY = 120)
simul_params <-@TimeParam
simul_params#> $Nyears
#> [1] 6
#>
#> $nTSpY
#> [1] 120
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 rust pathogens of cereal crops is available using function loadPathogen()
:
loadPathogen(disease = "rust")
basic_patho_param <-
basic_patho_param#> $name
#> [1] "rust"
#>
#> $survival_prob
#> [1] 1e-04
#>
#> $repro_sex_prob
#> [1] 0
#>
#> $infection_rate
#> [1] 0.4
#>
#> $propagule_prod_rate
#> [1] 3.125
#>
#> $latent_period_exp
#> [1] 10
#>
#> $latent_period_var
#> [1] 9
#>
#> $infectious_period_exp
#> [1] 24
#>
#> $infectious_period_var
#> [1] 105
#>
#> $sigmoid_kappa
#> [1] 5.333
#>
#> $sigmoid_sigma
#> [1] 3
#>
#> $sigmoid_plateau
#> [1] 1
This list may be updated to change a specific parameter. For instance, to change the infection rate to 50%:
loadPathogen("rust")
basic_patho_param <-$infection_rate <- 0.5
basic_patho_param
basic_patho_param#> $name
#> [1] "rust"
#>
#> $survival_prob
#> [1] 1e-04
#>
#> $repro_sex_prob
#> [1] 0
#>
#> $infection_rate
#> [1] 0.5
#>
#> $propagule_prod_rate
#> [1] 3.125
#>
#> $latent_period_exp
#> [1] 10
#>
#> $latent_period_var
#> [1] 9
#>
#> $infectious_period_exp
#> [1] 24
#>
#> $infectious_period_var
#> [1] 105
#>
#> $sigmoid_kappa
#> [1] 5.333
#>
#> $sigmoid_sigma
#> [1] 3
#>
#> $sigmoid_plateau
#> [1] 1
Finally, the list may be generated manually to control all parameters:
list(infection_rate = 0.4
basic_patho_param <-latent_period_exp = 10
, latent_period_var = 9
, propagule_prod_rate = 3.125
, infectious_period_exp = 24
, infectious_period_var = 105
, survival_prob = 1e-4
, repro_sex_prob = 0
, sigmoid_kappa = 5.333, sigmoid_sigma = 3, sigmoid_plateau = 1) ,
Then, the list is used to fuel the object simul_params
via the function setPathogen()
:
setPathogen(simul_params, patho_params = basic_patho_param)
simul_params <-@Pathogen
simul_params#> $infection_rate
#> [1] 0.4
#>
#> $latent_period_exp
#> [1] 10
#>
#> $latent_period_var
#> [1] 9
#>
#> $propagule_prod_rate
#> [1] 3.125
#>
#> $infectious_period_exp
#> [1] 24
#>
#> $infectious_period_var
#> [1] 105
#>
#> $survival_prob
#> [1] 1e-04
#>
#> $repro_sex_prob
#> [1] 0
#>
#> $sigmoid_kappa
#> [1] 5.333
#>
#> $sigmoid_sigma
#> [1] 3
#>
#> $sigmoid_plateau
#> [1] 1
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.
setInoculum(simul_params, val = 5e-4)
simul_params <-@PI0
simul_params#> [1] 5e-04
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()
:
loadLandscape(id = 1)
landscape <-length(landscape)
#> [1] 155
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()
:
loadDispersalPathogen(id = 1)
disp_patho <-head(disp_patho)
#> [1] 8.814254e-01 9.525884e-04 7.079895e-10 1.594379e-10 3.285800e-06
#> [6] 3.634297e-11
length(landscape)^2 == length(disp_patho)
#> [1] TRUE
Then, the object simul_params
is updated with the landscape and dispersal matrix via the functions setLandscape()
and setDispersalPathogen()
, respectively:
setLandscape(simul_params, land = landscape)
simul_params <- setDispersalPathogen(simul_params, mat = disp_patho) simul_params <-
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.
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.
loadCultivar(name = "Susceptible", type = "growingHost")
cultivar1 <- loadCultivar(name = "Resistant1", type = "growingHost")
cultivar2 <- loadCultivar(name = "Resistant2", type = "growingHost")
cultivar3 <- loadCultivar(name = "Resistant3", type = "nongrowingHost")
cultivar4 <- loadCultivar(name = "Forest", type = "nonCrop")
cultivar5 <- data.frame(rbind(cultivar1, cultivar2, cultivar3, cultivar4, cultivar5)
cultivars <-stringsAsFactors = FALSE)
,
cultivars#> cultivarName initial_density max_density growth_rate reproduction_rate
#> 1 Susceptible 0.1 2 0.1 0
#> 2 Resistant1 0.1 2 0.1 0
#> 3 Resistant2 0.1 2 0.1 0
#> 4 Resistant3 2.0 2 0.0 0
#> 5 Forest 0.0 2 0.0 0
#> death_rate yield_H yield_L yield_I yield_R planting_cost market_value
#> 1 0 2.5 0 0 0 225 200
#> 2 0 2.5 0 0 0 225 200
#> 3 0 2.5 0 0 0 225 200
#> 4 0 2.5 0 0 0 225 200
#> 5 0 0.0 0 0 0 0 0
Similarly as pathogen parameters, characteristics of the cultivars may be updated as required. For example, to change the growth rate of the susceptible cultivar:
$cultivarName == "Susceptible", "growth_rate"] <- 0.2
cultivars[cultivars
cultivars#> cultivarName initial_density max_density growth_rate reproduction_rate
#> 1 Susceptible 0.1 2 0.2 0
#> 2 Resistant1 0.1 2 0.1 0
#> 3 Resistant2 0.1 2 0.1 0
#> 4 Resistant3 2.0 2 0.0 0
#> 5 Forest 0.0 2 0.0 0
#> death_rate yield_H yield_L yield_I yield_R planting_cost market_value
#> 1 0 2.5 0 0 0 225 200
#> 2 0 2.5 0 0 0 225 200
#> 3 0 2.5 0 0 0 225 200
#> 4 0 2.5 0 0 0 225 200
#> 5 0 0.0 0 0 0 0 0
Finally, the dataframe cultivars
can also be generated entirely from scratch:
data.frame(cultivarName = c("Susceptible", "Resistant"),
cultivars_new <-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),
death_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#> cultivarName initial_density max_density growth_rate reproduction_rate
#> 1 Susceptible 0.1 2 0.1 0
#> 2 Resistant 0.2 3 0.2 0
#> death_rate yield_H yield_L yield_I yield_R planting_cost market_value
#> 1 0 2.5 0 0 0 225 200
#> 2 0 2.0 0 0 0 300 150
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.
loadGene(name = "MG 1", type = "majorGene")
gene1 <- loadGene(name = "Lr34", type = "APR")
gene2 <- loadGene(name = "gene 3", type = "QTL")
gene3 <- loadGene(name = "nonhost resistance", type = "immunity")
gene4 <- data.frame(rbind(gene1, gene2, gene3, gene4), stringsAsFactors = FALSE)
genes <-
genes#> geneName efficiency time_to_activ_exp time_to_activ_var
#> 1 MG 1 1.0 0 0
#> 2 Lr34 1.0 30 30
#> 3 gene 3 0.5 0 0
#> 4 nonhost resistance 1.0 0 0
#> mutation_prob Nlevels_aggressiveness fitness_cost tradeoff_strength
#> 1 1e-07 2 0.5 1
#> 2 1e-07 2 0.5 1
#> 3 1e-04 6 0.5 1
#> 4 0e+00 1 0.0 1
#> target_trait
#> 1 IR
#> 2 IR
#> 3 IR
#> 4 IR
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”:
$geneName == "MG 1", "mutation_prob"] <- 1e-3
genes[genes
genes#> geneName efficiency time_to_activ_exp time_to_activ_var
#> 1 MG 1 1.0 0 0
#> 2 Lr34 1.0 30 30
#> 3 gene 3 0.5 0 0
#> 4 nonhost resistance 1.0 0 0
#> mutation_prob Nlevels_aggressiveness fitness_cost tradeoff_strength
#> 1 1e-03 2 0.5 1
#> 2 1e-07 2 0.5 1
#> 3 1e-04 6 0.5 1
#> 4 0e+00 1 0.0 1
#> target_trait
#> 1 IR
#> 2 IR
#> 3 IR
#> 4 IR
Finally, the dataframe genes
can also be generated entirely from scratch:
data.frame(geneName = c("MG1", "MG2"),
genes_new <-efficiency = c(1.0 , 0.8 ),
time_to_activ_exp = c(0.0 , 0.0 ),
time_to_activ_var = c(0.0 , 0.0 ),
mutation_prob = c(1E-7 , 1E-4),
Nlevels_aggressiveness = c(2 , 2 ),
fitness_cost = c(0.50 , 0.75 ),
tradeoff_strength = c(1.0 , 1.0 ),
target_trait = c("IR" , "LAT"),
stringsAsFactors = FALSE)
genes_new#> geneName efficiency time_to_activ_exp time_to_activ_var mutation_prob
#> 1 MG1 1.0 0 0 1e-07
#> 2 MG2 0.8 0 0 1e-04
#> Nlevels_aggressiveness fitness_cost tradeoff_strength target_trait
#> 1 2 0.50 1 IR
#> 2 2 0.75 1 LAT
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.
The object simul_params
can be updated with setGenes()
and setCultivars()
:
setGenes(simul_params, dfGenes = genes)
simul_params <- setCultivars(simul_params, dfCultivars = cultivars)
simul_params <-@Genes
simul_params#> geneName efficiency time_to_activ_exp
#> MG 1 MG 1 1.0 0
#> Lr34 Lr34 1.0 30
#> gene 3 gene 3 0.5 0
#> nonhost resistance nonhost resistance 1.0 0
#> time_to_activ_var mutation_prob Nlevels_aggressiveness
#> MG 1 0 1e-03 2
#> Lr34 30 1e-07 2
#> gene 3 0 1e-04 6
#> nonhost resistance 0 0e+00 1
#> fitness_cost tradeoff_strength target_trait
#> MG 1 0.5 1 IR
#> Lr34 0.5 1 IR
#> gene 3 0.5 1 IR
#> nonhost resistance 0.0 1 IR
@Cultivars
simul_params#> cultivarName initial_density max_density growth_rate
#> Susceptible Susceptible 0.1 2 0.2
#> Resistant1 Resistant1 0.1 2 0.1
#> Resistant2 Resistant2 0.1 2 0.1
#> Resistant3 Resistant3 2.0 2 0.0
#> Forest Forest 0.0 2 0.0
#> reproduction_rate death_rate yield_H yield_L yield_I yield_R
#> Susceptible 0 0 2.5 0 0 0
#> Resistant1 0 0 2.5 0 0 0
#> Resistant2 0 0 2.5 0 0 0
#> Resistant3 0 0 2.5 0 0 0
#> Forest 0 0 0.0 0 0 0
#> planting_cost market_value
#> Susceptible 225 200
#> Resistant1 225 200
#> Resistant2 225 200
#> Resistant3 225 200
#> Forest 0 0
Then the function allocateCultivarGenes()
allows the attribution of resistance genes to cultivars:
allocateCultivarGenes(simul_params
simul_params <-cultivarName = "Resistant1"
, listGenesNames = c("MG 1"))
, allocateCultivarGenes(simul_params
simul_params <-cultivarName = "Resistant2"
, listGenesNames = c("Lr34", "gene 3"))
, allocateCultivarGenes(simul_params
simul_params <-cultivarName = "Resistant3"
, listGenesNames = c("nonhost resistance"))
, @Cultivars
simul_params#> cultivarName initial_density max_density growth_rate
#> Susceptible Susceptible 0.1 2 0.2
#> Resistant1 Resistant1 0.1 2 0.1
#> Resistant2 Resistant2 0.1 2 0.1
#> Resistant3 Resistant3 2.0 2 0.0
#> Forest Forest 0.0 2 0.0
#> reproduction_rate death_rate yield_H yield_L yield_I yield_R
#> Susceptible 0 0 2.5 0 0 0
#> Resistant1 0 0 2.5 0 0 0
#> Resistant2 0 0 2.5 0 0 0
#> Resistant3 0 0 2.5 0 0 0
#> Forest 0 0 0.0 0 0 0
#> planting_cost market_value
#> Susceptible 225 200
#> Resistant1 225 200
#> Resistant2 225 200
#> Resistant3 225 200
#> Forest 0 0
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.
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:
loadCroptypes(simul_params, names = c("Susceptible crop"
croptypes <-"Pure resistant crop"
, "Mixture"
, "Other"))
,
croptypes#> croptypeID croptypeName Susceptible Resistant1 Resistant2 Resistant3
#> 1 0 Susceptible crop 0 0 0 0
#> 2 1 Pure resistant crop 0 0 0 0
#> 3 2 Mixture 0 0 0 0
#> 4 3 Other 0 0 0 0
#> Forest
#> 1 0
#> 2 0
#> 3 0
#> 4 0
Then croptype
is updated by allocateCroptypeCultivars()
to specify the composition of every croptype with regard to cultivars (and proportion for mixtures):
allocateCroptypeCultivars(croptypes
croptypes <-croptypeName = "Susceptible crop"
, cultivarsInCroptype = "Susceptible")
, allocateCroptypeCultivars(croptypes
croptypes <-croptypeName = "Pure resistant crop"
, cultivarsInCroptype = "Resistant1")
, allocateCroptypeCultivars(croptypes
croptypes <-croptypeName = "Mixture"
, cultivarsInCroptype = c("Resistant2","Resistant3")
, prop = c(0.4, 0.6))
, allocateCroptypeCultivars(croptypes
croptypes <-croptypeName = "Other"
, cultivarsInCroptype = "Forest")
,
croptypes#> croptypeID croptypeName Susceptible Resistant1 Resistant2 Resistant3
#> 1 0 Susceptible crop 1 0 0.0 0.0
#> 2 1 Pure resistant crop 0 1 0.0 0.0
#> 3 2 Mixture 0 0 0.4 0.6
#> 4 3 Other 0 0 0.0 0.0
#> Forest
#> 1 0
#> 2 0
#> 3 0
#> 4 1
Finally the object simul_params
is updated using setCroptypes()
:
setCroptypes(simul_params, dfCroptypes = croptypes)
simul_params <-@Croptypes
simul_params#> croptypeID croptypeName Susceptible Resistant1 Resistant2 Resistant3
#> 0 0 Susceptible crop 1 0 0.0 0.0
#> 1 1 Pure resistant crop 0 1 0.0 0.0
#> 2 2 Mixture 0 0 0.4 0.6
#> 3 3 Other 0 0 0.0 0.0
#> Forest
#> 0 0
#> 1 0
#> 2 0
#> 3 1
Alternatively, the dataframe croptypes
can be generated from scratch:
data.frame(croptypeID = c(0, 1, 2, 3)
croptypes <-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)
, setCroptypes(simul_params, croptypes) simul_params <-
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:
list(c(0,1,3), c(0,2,3))
rotation_sequence <- 2 # number of years before rotation of the landscape
rotation_period <- list(rep(1/3, 3), rep(1/3, 3)) # proportion (in surface) of each croptype
prop <- 1 # level of spatial aggregation
aggreg <- allocateLandscapeCroptypes(simul_params
simul_params <-rotation_period = rotation_period
, rotation_sequence = rotation_sequence
, prop = prop
, aggreg = aggreg
, graphic = FALSE)
, # plot(simul_params@Landscape)
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 ?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_rel” : Relative Area Under Disease Progress Curve (average proportion of diseased hosts relative to the total number of existing hosts)
- “gla” : Green Leaf Area (average number of healthy hosts per square meter)
- “gla_rel” : Relative Green Leaf Area (average proportion of healthy hosts relative to the total number of existing hosts)
- “eco_yield” : total crop yield (in weight or volume units per ha)
- “eco_cost” : operational crop costs (in monetary units per ha)
- “eco_product” : total crop products (in monetary units per ha)
- “eco_margin” : Margin (products - operational costs, in monetary units per ha) - “contrib”: contribution of pathogen genotypes to LIR dynamics - “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 include:
- “evol_patho”: Dynamics of pathogen genotype frequencies
- “evol_aggr”: Evolution of pathogen aggressiveness
- “durability”: Durability of resistance genes
- “all”: compute all these outputs (default)
- "": none of these outputs will be generated.
A list of outputs can be generated using loadOutputs()
:
loadOutputs(epid_outputs = "all", evol_outputs = "all")
outputlist <-
outputlist#> $epid_outputs
#> [1] "all"
#>
#> $evol_outputs
#> [1] "all"
#>
#> $thres_breakdown
#> [1] 50000
#>
#> $GLAnoDis
#> [1] 1.48315
#>
#> $audpc100S
#> [1] 0.76
Then simul_params
can be updated via setOutputs()
:
setOutputs(simul_params, outputlist) simul_params <-
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
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)
saveDeploymentStrategy(simul_params) simul_params <-
Then, the function runSimul()
launches the simulation. Use ?runSimul
to get all available options.
runSimul(simul_params, graphic = TRUE, videoMP4 = FALSE)