This vignette provides a comprehensive guide to the dynamAedes package, a unified modelling framework for invasive Aedes mosquitoes. Users can employ this spatially-explicit, time-discrete, stochastic population dynamical model, which was initially formulated for *Aedes aegypti* in Da Re et al., (2021) and subsequently adapted for Ae. albopictus, Ae. japonicus, and Ae. koreicus Da Re et al., (2022).

This model integrates factors like temperature, photoperiod, and intra-specific larval competition and is versatile across three spatial scales: punctual, local, and regional. The various scales provide flexibility in addressing spatial complexity and data availability, accommodating both active and passive mosquito dispersal and varying input temperature data.

For illustrative purposes, this guide focuses on model applications for *Ae. albopictus* across all spatial scales using a simulated temperature dataset.

At the regional scale, the model operates similarly to the “punctual” scale, processing each grid cell separately without considering active or passive dispersal. This makes each cell an isolated mosquito population. To execute the model in this context, two datasets are essential:

- A temperature matrix (in °C) spatially and temporally defined.
- A two-column matrix showcasing the centroid coordinates (in meters) for every cell.

For the purpose of this guide, we will employ simulated datasets:

- A grid with a 5 km lattice and 250 m cell size.
- A spatially and temporally correlated temperature time series spanning one year.

To begin, we will establish a squared lattice arena, which will serve as the site for our mosquito introduction. The chosen dimensions are 5 km on each side with a resolution of 250 m (comprising 20 columns and 20 rows, totalling 400 cells).

```
# Libraries
library(gstat)
library(terra)
library(eesim)
library(dynamAedes)
library(ggplot2)
Sys.setlocale("LC_TIME", "en_GB.UTF-8")
```

`# [1] "en_GB.UTF-8"`

```
20 # 5000m/250 m = 20 columns and rows
gridDim <- expand.grid(x=1:gridDim, y=1:gridDim) xy <-
```

To infuse a spatial element into the lattice area, we will overlay a spatial pattern, which will subsequently be used to impart spatial correlation to the temperature time series. This spatially correlated pattern is derived using a semivariogram model with a specified sill and range. We then predict this model over the grid using unconditional Gaussian simulation.

```
vgm(psill=0.5, range=100, model='Exp') # psill = partial sill = (sill-nugget)
varioMod <-# Set up an additional variable from simple kriging
gstat(formula=z~1,
zDummy <-locations = ~x+y,
dummy=TRUE,
beta=1,
model=varioMod,
nmax=1)
# Generate a randomly autocorrelated predictor data field
set.seed(123)
predict(zDummy, newdata=xy, nsim=1) xyz <-
```

`# [using unconditional Gaussian simulation]`

Creating a spatially correlated raster with the SA variable enhances our understanding of environmental features like urban vegetation distribution.

```
"+proj=utm +zone=32 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
utm32N <- terra::rast(nrow=gridDim, ncol=gridDim, crs=utm32N, ext=terra::ext(1220000,1225000, 5700000,5705000))
r <-::values(r) <- xyz$sim1
terraplot(r, main="SAC landscape")
```

```
# convert to a data.frame
data.frame("id"=1:nrow(xyz), terra::crds(r))
df <- terra::as.polygons(terra::ext(r), crs=utm32N)
bbox <-
# Store Parameters for autocorrelation
terra::values(r) autocorr_factor <-
```

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

```
365
ndays =set.seed(123)
eesim::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 Temperature Seasonal Trend",
xlab="Date", ylab="Temperature (°C)"
)
```

We can then “expand onto space” the temperature time series by multiplying it with the autocorrelated surface simulated above.

```
do.call(rbind, lapply(1:ncell(r), function(x) {
mat <- sim_temp[[1]]$x*autocorr_factor[[x]]
d_t <-return(d_t)
}))
```

` par(mfrow = c(1,2)) oldpar <-`

A comparison between the distribution of the initial temperature time series and autocorrelated temperature surface

```
par(mfrow=c(2,1))
hist(mat, xlab="Temperature (°C)", main="Histogram of simulated spatial autocorreled temperature")
hist(sim_temp[[1]]$x, xlab="Temperature (°C)", main="Histogram of simulated temperatures", col="red")
```

`par(mfrow=c(1,1))`

`par(oldpar) `

```
names(mat) <- paste0("d_", 1:ndays)
cbind(df, mat) df_temp <-
```

Float numbers in the temperature matrix would slow the computational speed, thus we first multiply them by 1000 and then transform them in integer numbers. **w** will be subset to match the simulated time period below.

` sapply(df_temp[,-c(1:3)], function(x) as.integer(x*1000)) w <-`

We can now define a two-column matrix of coordinates to identify each cell in the lattice grid.

` df_temp[,c("x","y")] cc <-`

We are now left with a few model variables which need to be defined.

```
## Define the day of introduction (May 1st is day 1)
"2000-06-01"
str <-## Define the end-day of life cycle (July 2nd is the last day)
"2000-07-02"
endr <-## Define the number of eggs to be introduced
100
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 density-dependent mortality
1
habitat_liters <-## Define proj4 string for input coordinates
"+proj=utm +zone=32 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
utm32N <-## Define the number of parallel processes (for sequential iterations set nc=1)
1 cl <-
```

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

```
dynamAedes.m(species="albopictus",
simout <-scale="rg",
jhwv=habitat_liters,
temps.matrix=w[,as.numeric(format(as.Date(str),"%j")):as.numeric(format(as.Date(endr),"%j"))],
coords.proj4=utm32N,
cells.coords=as.matrix(cc),
startd=str,
endd=endr,
n.clusters=cl,
iter=it,
intro.eggs=ie,
compressed.output=TRUE,
seeding=TRUE,
verbose=FALSE)
```

#3. Analyse the results A first summary of simulations can be obtained with:

`summary(simout)`

```
# Summary of dynamAedes simulations:
# ----------------------------------
# Species: Aedes albopictus
# Scale: REGIONAL
# Start Date: 2000-06-01
# End Date: 2000-07-02
# Number of Iterations: 1
# Introduced Stage: egg
# Number Introduced: 100
# Is Output Compressed?: Yes
# Water in the System: 1 L
# Min days with population: 30
# Max days with population: 30
```

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 **second**corresponds 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 within each grid cell of the landscape (columns). If we inspect the first day within the first iteration, we obtain a matrix having:

`dim(simout@simulation[[1]][[1]])`

`# [1] 4 400`

We can now use the auxiliary functions of the package to Analyse the results.

First, we can retrieve the “probability of a successful introduction”, computed as the proportion of model iterations that resulted in a viable mosquito population (in any cells of the grid) at a given date.

`psi(input_sim = simout, eval_date = 30)`

```
# Days_after_intro p_success stage
# 1 Day 30 1 Population
```

We can also get a “spatial output”, using the function *psi_sp*, which requires as additional input only the matrix of the pixels coordinates

`plot(psi_sp(input_sim = simout, eval_date = 30))`

We can now compute the interquantile range abundance of the simulated population using the function *adci* over the whole landscape.

```
max(simout) #retrieve the maximum number of simulated days
dd <-
# Compute the inter-quartile of abundances along the iterations
c(0.25,0.50,0.75)
breaks=1:dd
ed=
# type "O" derives a non-spatial time series
rbind(
outdf <-adci(simout, eval_date=ed, breaks=breaks, stage="Eggs", type="O"),
adci(simout, eval_date=ed, breaks=breaks, stage="Juvenile", type="O"),
adci(simout, eval_date=ed, breaks=breaks, stage="Adults", type="O"),
adci(simout, eval_date=ed, breaks=breaks, stage="Dia", type="O")
)
```

Then we can look at the time series of the population dynamics stage by stage at the whole landscape level.

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