This tutorial explains stepbystep the main features of dynamAedes package, a unified modelling framework for invasive Aedes mosquitoes. Users can apply the stochastic, timediscrete and spatiallyexplicit population dynamical model initially developed in Da Re et al., (2021) for Aedes aegypti and then expanded for other three species: Ae. albopictus, Ae. japonicus and Ae. koreicus Da Re et al., (2022).
The model is driven by temperature, photoperiod and intraspecific larval competition and can be applied to three different spatial scales: punctual, local and regional. These spatial scales consider different degrees of spatial complexity and data availability, by accounting for both active and passive dispersal of the modelled mosquito species as well as for the heterogeneity of input temperature data.
We will describe the model applications for Ae. albopictus and for all spatial scales by using a simulated temperature dataset.
At the punctual scale, the model only requires a weather station temperature time series provided as a numerical matrix (in degree Celsius). For the purpose of this tutorial, we simulate a 1year long temperature time series.
We first simulate a 1year temperature time series with seasonal trend. For the time series we consider a mean value of 16°C and standard deviation of 2°C.
#Load packages
# simulate temperatures
library(eesim)
# modelling
library(dynamAedes)
# data manipulation and plotting
library(ggplot2)
Sys.setlocale("LC_TIME", "en_GB.UTF8")
# [1] "en_GB.UTF8"
365*1 #length of the time series in days
ndays <set.seed(123)
create_sims(n_reps = 1,
sim_temp <n = ndays,
central = 16,
sd = 2,
exposure_type = "continuous",
exposure_trend = "cos1", exposure_amp = 1.0,
average_outcome = 12,
outcome_trend = "cos1",
outcome_amp = 0.8,
rr = 1.0055)
A visualisation of the distribution of temperature values and temporal trend.
hist(sim_temp[[1]]$x,
xlab="Temperature (°C)",
main="Histogram of simulated temperatures")
plot(sim_temp[[1]]$date,
1]]$x,
sim_temp[[main="Simulated temperatures seasonal trend",
xlab="Date", ylab="Temperature (°C)"
)
We define a few model parameters which need to be defined by the user.
## Define the day of introduction (July 1st is day 1)
"20000701"
str <## Define the endday of life cycle (August 1st is the last day)
"20000801"
endr <## Define the number of eggs to be introduced
1000
ie <## Define the number of model iterations
1 # The higher the number of simulations the better
it <## Define the number of liters for the larval densitydependent mortality
1
habitat_liters <#Define latitude and longitude for the diapause process
42
myLat < 7
myLon <## Define the number of parallel processes (for sequential iterations set nc=1)
1 cl <
Float numbers in the temperature matrix would slow the computational speed, they must be multiplied by 1000 and then transformed in integer numbers. We also transpose the matrix from long to wide format, since we conceptualized the model structure considering the rows as the spatial component (e.g., observations; here = 1) and the columns as the temporal one (e.g., variables).
data.frame("Date" = sim_temp[[1]]$date, "temp" = sim_temp[[1]]$x)
df_temp < t(as.integer(df_temp$temp*1000)[format(as.Date(str)+1,"%j"):format(as.Date(endr)+1,"%j")]) w <
Running the model with the settings specified in this example takes about 2 minutes.
dynamAedes.m(species="albopictus",
simout <scale="ws",
jhwv=habitat_liters,
temps.matrix=w,
startd=str,
endd=endr,
n.clusters=cl,
iter=it,
intro.eggs=ie,
compressed.output=TRUE,
lat=myLat,
long=myLon,
verbose=FALSE,
seeding=TRUE)
#

  0%
A first summary of simulations can be obtained with:
summary(simout)
# Length Class Mode
# 1 dynamAedesClass S4
The simout object is a S4 object where simulation outputs and related details are saved in different slot:
For example, the number of model iterations is saved in:
@n_iterations simout
The simulation output is stored in:
@simulation simout
Which is a list where the the first level stores simulation of different iteration, while the secondcorresponds to the simulated days in the corresponding iteration. If we inspect the first iteration, we observe that the model has computed length(simout[[1]])
days, since we have started the simulation on the 1st of July and ended on the 1st of August.
length(simout@simulation[[1]])
# [1] 30
The third level corresponds to the quantity of individuals for each stage (rows) in each day for each iteration. If we inspect the 1st and the 15th day within the first iteration, we obtain a matrix having:
@simulation[[1]][[1]] simout
# [,1]
# egg 645
# juvenile 355
# adult 0
# diapause_egg 0
@simulation[[1]][[15]] simout
# [,1]
# egg 128
# juvenile 108
# adult 8
# diapause_egg 0
We can use the auxiliary functions of the package to analyse the results.
First, we can retrieve the “probability of successful introduction”, computed as the proportion of model iterations that resulted in a viable mosquito population at a given date. In this case the results is 1, since we have only one iteration and the population is still viable at the end of the simulation.
psi(input_sim = simout, eval_date = 30)
# Days_after_intro p_success stage
# 1 Day 30 1 Population
We now compute the interquantile range abundance for all the stages of the simulated population using the function adci.
# Retrieve the maximum number of simulated days
max(simout)
dd <
# Compute the interquartile of abundances along the iterations
c(0.25,0.50,0.75)
breaks=1:dd
ed=
rbind(
outdf <adci(simout, eval_date=ed, breaks=breaks, stage="Eggs"),
adci(simout, eval_date=ed, breaks=breaks, stage="Juvenile"),
adci(simout, eval_date=ed, breaks=breaks, stage="Adults"),
adci(simout, eval_date=ed, breaks=breaks, stage="Dia"))
Then we can look at the time series of the population dynamics stage by stage.
$stage < factor(outdf$stage, levels= c('Egg', 'DiapauseEgg', 'Juvenile', 'Adult'))
outdf$Date < rep(seq.Date(as.Date(str), as.Date(endr)  2, by="day"), 4)
outdf
ggplot(outdf, aes(x=Date, y=X50., group=factor(stage), col=factor(stage))) +
ggtitle("Ae. albopictus Interquantile range abundance") +
geom_ribbon(aes(ymin=X25., ymax=X75., fill=factor(stage)),
col="white",
alpha=0.2,
outline.type="full") +
geom_line(linewidth=0.8) +
labs(x="Date", y="Interquantile range abundance", col="Stage", fill="Stage") +
facet_wrap(~stage, scales = "free_y") +
theme_light() +
theme(legend.position="bottom",
text = element_text(size=16),
strip.text = element_text(face = "italic"))