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.
#Load packages
require(spatstat)
require(gstat)
require(parallel)
require(eesim)
require(dplyr)
require(geosphere)
require(ggplot2)
require(terra)
require(rgdal)
require(raster)
require(dynamAedes)
Sys.setlocale("LC_TIME", "en_GB.UTF8")
# [1] "en_GB.UTF8"
At 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.
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)"
)
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).
We are now left with 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 itarations set nc=1)
1
cl =## Set output name for the *.RDS output will be saved
#outname= paste0("dynamAedes.m_albo_ws_dayintro_",str,"_end",endr,"_niters",it,"_neggs",ie)
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%
We first explore the model output structure: the simout object is a nested list.
The first level corresponds to the number of model iterations
print(it)
# [1] 1
print(length(simout))
# [1] 1
The second level corresponds to the simulated days. So if we inspect the first iteration, we observe that the model has computed rlength(simout[[1]])
days, since we have started the simulation on the 1st of July and ended on the 1st of August.
length(simout[[1]])
# [1] 30
The third level corresponds to the amount of individuals for each stage (rows). So if we inspect the 1st and the 15th day within the first iteration, we obtain a matrix having
dim(simout[[1]][[1]])
# [1] 4 1
1]][[1]] simout[[
# [,1]
# egg 645
# juvenile 355
# adult 0
# diapause_egg 0
1]][[15]] simout[[
# [,1]
# egg 128
# juvenile 108
# adult 8
# diapause_egg 0
We can now 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(sapply(simout, function(x) length(x)))
dd <
# Compute the abundance
rbind(data.frame(cbind(
outdf <data.frame(adci(simout, eval_date=1:dd, breaks=c(0.25,0.50,0.75), st=1)),
myStage = 'Egg')),
data.frame(cbind(
adci(simout, eval_date=1:dd, breaks=c(0.25,0.50,0.75), st=2),
myStage='Juvenile')),
data.frame(cbind(
adci(simout, eval_date=1:dd, breaks=c(0.25,0.50,0.75), st=3),
myStage='Adult')),
data.frame(cbind(
adci(simout, eval_date=1:dd, breaks=c(0.25,0.50,0.75), st=4),
myStage='Diapausing egg'))
)
# Plot
%>%
outdf mutate(myStage=factor(myStage, levels= c('Egg', 'Diapausing egg', 'Juvenile', 'Adult')),
X25.=as.numeric(X25.),
X50.=as.numeric(X50.),
X75.=as.numeric(X75.),
Date=rep(seq.Date(as.Date(str), as.Date(endr)2, by="day"),4)) %>%
ggplot( aes(x=Date, y=X50., group=factor(myStage),col=factor(myStage))) +
ggtitle("Ae. albopictus Interquantile range abundance")+
geom_line(linewidth=0.8)+
geom_ribbon(aes(ymin=X25.,ymax=X75.,fill=factor(myStage)),
col="white",
alpha=0.2,
outline.type="full")+
labs(x="Date", y="Interquantile range abundance", col="Stage", fill="Stage")+
facet_wrap(~myStage, scales = "free_y")+
theme_light()+
theme(legend.pos="bottom",
text = element_text(size=16) ,
strip.text = element_text(face = "italic"))