This tutorial explains step-by-step the main features of dynamAedes package, a unified modelling framework for invasive Aedes mosquitoes. Users can apply the stochastic, time-discrete and spatially-explicit 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 intra-specific 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 1-year long temperature time series.

1. Input data

1.1. Simulate temperature data with seasonal trend

We first simulate a 1-year 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.UTF-8")

# [1] "en_GB.UTF-8"

ndays <- 365*1 #length of the time series in days
set.seed(123)
sim_temp <- create_sims(n_reps = 1, 
  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,
 sim_temp[[1]]$x,
 main="Simulated temperatures seasonal trend", 
 xlab="Date", ylab="Temperature (°C)"
 )

2. Modelling

2.1 Model settings

We define a few model parameters which need to be defined by the user.

## Define the day of introduction (July 1st is day 1)
str <- "2000-07-01"
## Define the end-day of life cycle (August 1st is the last day)
endr <- "2000-08-01"
## Define the number of eggs to be introduced
ie <- 1000
## Define the number of model iterations
it <- 1 # The higher the number of simulations the better
## Define the number of liters for the larval density-dependent mortality
habitat_liters <- 1
#Define latitude and longitude for the diapause process
myLat <- 42
myLon <- 7
## Define the number of parallel processes (for sequential iterations set nc=1)
cl <- 1

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).

df_temp <- data.frame("Date" = sim_temp[[1]]$date, "temp" = sim_temp[[1]]$x)
w <- t(as.integer(df_temp$temp*1000)[format(as.Date(str)+1,"%j"):format(as.Date(endr)+1,"%j")])

2.2 Run the model

Running the model with the settings specified in this example takes about 2 minutes.

simout <- dynamAedes.m(species="albopictus", 
 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%

3. Analyse the results

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:

simout@n_iterations

The simulation output is stored in:

simout@simulation

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:

simout@simulation[[1]][[1]]

#              [,1]
# egg           645
# juvenile      355
# adult           0
# diapause_egg    0

simout@simulation[[1]][[15]]

#              [,1]
# egg           128
# juvenile      108
# adult           8
# diapause_egg    0

We can use the auxiliary functions of the package to analyse the results.

3.1 Derive probability of a successfull introduction at the end of the simulated period

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

3.2 Derive abundance 95% CI for each life stage and in each day

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
dd <- max(simout)

# Compute the inter-quartile of abundances along the iterations
breaks=c(0.25,0.50,0.75)
ed=1:dd

outdf <- rbind(
  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.

outdf$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)

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"))

01. The punctual scale model