In this vignette, we show how to run **MGDrivE2** simulations using a variety of sampling techniques, including deterministic integration with ODE solvers from the `deSolve`

package. We demonstrate a logistic (carrying capacity) form of larval mortality as well as the Lotka-Volterra form. While the behavior of the two types of density-dependent mortality will be nearly identical, because they are both quadratic in the number of larvae, this is meant to demonstrate the relative ease at which new density-dependent (and in general, any hazard, see “MGDrivE2: Simulation of Time-inhomogeneous Stochastic Processes (Seasonality)” for a more complex example) functions can be swapped as necessary for the simulation task at hand.

We start by loading the **MGDrivE2** package, as well as the **MGDrivE** package for access to inheritance cubes and **ggplot2** for graphical analysis. We will use the basic cube to simulate Mendelian inheritance for this example.

```
# simulation functions
library(MGDrivE2)
# inheritance patterns
library(MGDrivE)
# plotting
library(ggplot2)
# basic inheritance pattern
MGDrivE::cubeMendelian() cube <-
```

Next, we set entomological parameters that dictate life-history characteristics, as well as the equilibrium number of adult female mosquitoes. The density dependent per-capita larval mortality is thus: \(\mu_{L}(1 + \frac{L}{K})\), where \(L\) is the current larval population size. Note that when the population is zero there is no effect of density dependence.

The parameters `qE`

, `qL`

, and `qP`

are the mean dwell times in the egg, larval, and pupae stages, respectively. The Erlang distributed aquatic stages allow users to select the mean and variance of the dwell time distribution for each stage. For example, \(\frac{1}{q_{E}}\) will give the mean dwell time in the egg stage, and variance \(\frac{1}{n_{E} \cdot q_{E}^{2}}\). A table of (case-sensitive) biological parameters the user needs to specify is given below (Note that all parameters must be specified as **rates per day**). The equilibrium constant for larval density dependence (`K`

) is returned when equilibrium is calculated.

The `nu`

parameter gives the rate at which unmated female mosquitoes mate when males become available again. Ordinarily, female mosquitoes are assumed to mate upon emergence if males are present. However if male mosquitoes become locally extinct and are later replenished either by currently developing egg batches or releases (or for metapopulation simulations, inbound male migration), then there will be some time where newly emerging adult females are unable to find mates. These unmated female mosquitoes experience mortality at the same rate as mated females, `muF`

. Because they are unmated, they do not produce eggs, nor do they blood-feed (which is done in order to provision developing egg batches). While blood-feeding behavior doesn’t matter for lifecycle simulations, this means that for the epidemiological simulations we don’t need to worry about unmated females becoming infected, because they will not be biologically primed for this behavior until after they have mated and need to produce egg batches. Once males become available, they will mate at rate `nu`

.

Parameter | Description |
---|---|

`qE` |
inverse of mean duration of egg stage |

`nE` |
shape parameter of Erlang-distributed egg stage |

`qL` |
inverse of mean duration of larval stage |

`nL` |
shape parameter of Erlang-distributed larval stage |

`qP` |
inverse of mean duration of pupal stage |

`nP` |
shape parameter of Erlang-distributed pupal stage |

`muE` |
density-independent egg mortality |

`muL` |
density-independent larval mortality |

`muP` |
density-independent pupal mortality |

`muF` |
density-independent adult female mortality |

`muM` |
density-independent adult male mortality |

`beta` |
rate of egg production in adult females |

`nu` |
mating rate of unmated females |

```
# number of adult female mosquitoes
500
NF <-
# entomological parameters
list(
theta <-qE = 1/4,
nE = 2,
qL = 1/3,
nL = 3,
qP = 1/6,
nP = 2,
muE = 0.05,
muL = 0.15,
muP = 0.05,
muF = 0.09,
muM = 0.09,
beta = 16,
nu = 1/(4/24)
)
```

The next thing we need to do is set up the structural properties of the SPN model (technically just a Petri Net, PN, at this point). First, we set all of the “places” in the model. As this is a compartmental model, “places” are simply distinct life stages (and locations for metapopulation versions). Thus, the “places” include all aquatic stages (number of bins in each egg, larval, and pupal stages), adult males, and adult females (further separated by mate genotype). All of these stages are unique for each genotype as well, so the number of “places” (i.e., compartments) grows with the number of genotypes and Erlang-distributed stage length.

Once the places for the SPN are setup, we define the transitions between each place. Transitions are how objects move from one compartment to another. These are not the rates that objects change compartments, only whether or not such a move is possible. The rate of movement from one “place” to another are calculated later by the `spn_hazards()`

function.

Finally, as not all transitions apply to all “places”, we create a summary of possible transitions to and from each “place”. This is handled by the `spn_S()`

function.

```
# "Places"
# These are defined by Erlang-distributed life stages and genotypes
spn_P_lifecycle_node(params = theta,cube = cube)
SPN_P <-
# Transitions
# This is a list of viable transitions from one "place" to another
spn_T_lifecycle_node(spn_P = SPN_P,params = theta,cube = cube)
SPN_T <-
# Stoichiometry matrix
# A sparse matrix representing the effect of each transition on each place
spn_S(spn_P = SPN_P, spn_T = SPN_T) S <-
```

Once the “places” are defined, we use `equilibrium_lifeycle()`

to get initial conditions (\(M(0)\)), solved at the dynamic equilibrium for `NF`

number of adult females, as well as the density-dependent constants for larval competition, either `K`

(logistic dependence) or `gamma`

(Lotka-Volterra dependence).

```
# calculate equilibrium and setup initial conditions
# outputs required parameters in the named list "params"
# outputs intial equilibrium for adv users, "init
# outputs properly filled initial markings, "M0", vectorized over:
# all patches
# any genotypes( includes XX/XY inheritance patterns)
# allows initial ratios to vary
# properly mates males/females
equilibrium_lifeycle(params = theta, NF = NF, phi = 0.5,
initialCons <-log = TRUE, spn_P = SPN_P, cube = cube)
```

The transitions in our SPN define possible movement from one compartment to another, but they do not define the rate that these changes occur. The rate of change from one compartment to another is called a “hazard”, calculated by the `spn_hazards()`

function. This calculates the rate of change from one compartment to another. Below, we make both a set of exact hazards, for sampling algorithms with integer valued state spaces, and a set of approximate hazards, for algorithms using a continuous approximation.

```
# approximate hazards for continous approximation
spn_hazards(spn_P = SPN_P, spn_T = SPN_T, cube = cube,
approx_hazards <-params = initialCons$params, type = "life",
log = TRUE, exact = FALSE, tol = 1e-8,
verbose = FALSE)
# exact hazards for integer-valued state space
spn_hazards(spn_P = SPN_P, spn_T = SPN_T, cube = cube,
exact_hazards <-params = initialCons$params, type = "life",
log = TRUE, exact = TRUE, tol = NaN,
verbose = FALSE)
```

Finally, before we run simulations, we set up a release scheme to generate interesting dynamics. We will release 50 adult females with homozygous recessive alleles 5 times, every 10 days, starting at day 20. Releases are possible at every life-stage, so long as the event name matches the place in the simulation. It is critically important that **the event names match a place name** in the simulation. The simulation function checks this and will throw an error if the event name does not exist as a place in the simulation. This format is used in **MGDrivE2** for consistency with solvers in `deSolve`

.

```
# releases
seq(from = 20, length.out = 5, by = 10)
r_times <- 50
r_size <- data.frame("var" = paste0("F_", cube$releaseType, "_", cube$wildType),
events <-"time" = r_times,
"value" = r_size,
"method" = "add",
stringsAsFactors = FALSE)
```

The first simulation we will run is a mean-field deterministic approximation to the full stochastic model. Formally, this is a moment-closure technique where we model the 1st moment and ignore contributions from all higher order moments, so we can expect that the ODEs will approximate the mean of stochastic simulations well only when nonlinear effects and influence of higher order moments are small. For this continuous-state approximation, we need to use the approximate hazards, created by setting the following parameters `exact = FALSE, tol = 1e-8`

(above) in the function `spn_hazards()`

. When used in ODE approximation, the interpretation of hazard functions simplifies to simple rate functions.

Internally, **MGDrivE2** uses the high quality numerical solvers in from `deSolve`

to integrate a mean-field approximation to the stochastic model. Effectively, it evaluates the ODEs over the interval between time points where the user requested model output, calculating the rate of that hazard, given the current simulation state, multiplied with the stoichiometry matrix. The default method used is `lsoda`

, though for highly variable rates and inhomogeneous systems, consider using `rk4`

.

Now that the deterministic model is fully specified, let’s run the simulation for 200 days, outputting values of the state variables every day.

```
# max simulation time
125
tmax <-# time-step for output return, not the time-step of the sampling algorithm
1
dt <-
# run deterministic simulation
sim_trajectory_R(x0 = initialCons$M0, tmax = tmax, dt = dt, S = S,
ODE_out <-hazards = approx_hazards, sampler = "ode", method = "lsoda",
events = events, verbose = FALSE)
```

**MGDrivE2** has several helper functions for summarizing simulation data in a nice format for plotting. Here, we will look at the adult male and female data, though there are many more options explored below.

```
# summarize females by genotype
summarize_females(out = ODE_out$state, spn_P = SPN_P)
ODE_out_f <-
# summarize males by genotype
summarize_males(out = ODE_out$state)
ODE_out_m <-
# add sex for plotting
$sex <- "Female"
ODE_out_f$sex <- "Male"
ODE_out_m
# plot
ggplot(data = rbind(ODE_out_f, ODE_out_m)) +
geom_line(aes(x = time, y = value, color = genotype)) +
facet_wrap(facets = vars(sex), scales = "fixed") +
theme_bw() +
ggtitle("SPN: ODE Solution")
```

We see the initial equilibrium at 500 females, as was set above, and 500 males, because we defined an equal sex breakdown in the population. At time `t = 25`

, we see the first release of `aa`

females into the population, and then the 4 subsequent releases. Notice that releases only appear in the females, nothing was released into the males. Since there are no fitness costs, we see the `aa`

individuals come to equilibrium predominantly as heterozygotes, i.e. `Aa`

, with a very small fraction of homozygotes, and the majority of the population remaining `AA`

.

For stochastic simulations, we generally rely on tau-leaping; this method retains the integer valued state space, thus necessitating the use of exact hazards. The theoretical basis of tau-leaping depends upon realizing that each transition can be represented as a Poisson process; tau-leaping is therefore an Euler-method for solving continuous-time Markov chains where the “tau-leaping” refers to a fixed time-step of size “tau”, where the number of each transition that fires over such a time-step (\(\Delta t\)) is approximated by a Poisson random variable with a linear approximation to the cumulative hazard. This is provided as the `dt_stoch`

an argument to the function `sim_trajectory_R()`

.

We set `tau`

as our sampler type, use the same release events, conditions, and stoichiometry matrix as above, set the exact hazards instead of the approximate, and finally run the simulation and plot the output.

```
# delta t
0.1
dt_stoch <-
# tau sampling
sim_trajectory_R(x0 = initialCons$M0, tmax = tmax, dt = dt,
PTS_out <-dt_stoch = dt_stoch, S = S, hazards = exact_hazards,
sampler = "tau", events = events, verbose = FALSE)
# summarize females/males
summarize_females(out = PTS_out$state, spn_P = SPN_P)
PTS_out_f <- summarize_males(out = PTS_out$state)
PTS_out_m <-
# add sex for plotting
$sex <- "Female"
PTS_out_f$sex <- "Male"
PTS_out_m
# plot adults
ggplot(data = rbind(PTS_out_f, PTS_out_m)) +
geom_line(aes(x = time, y = value, color = genotype)) +
facet_wrap(facets = vars(sex), scales = "fixed") +
theme_bw() +
ggtitle("SPN: Tau-leaping Approximation")
```

As Tau-leaping is a stochastic approximation of our system, we should see some variation around the solution provided by the ODE simulation above. These simulations begin at the proper equilibrium amount, 500 individuals of each sex, with releases performed in females and an equilibrium shortly thereafter. The plots are heuristically similar to the ODE solution, which is good, but notice the slow trend of `Aa`

individuals towards extinction. Since there is such a small amount of `a`

alleles in the population, they are steadily being lost through stochastic drift. If we carried the simulations out longer, it is possible that the `a`

allele will be completely lost from the population.

Now we will show an example with the Lotka-Volterra style density dependence in larval stages. Please note that we do not need to recreate the structural pieces of the Petri Net (places and transitions), or update any of the parameters, because they remain the same; we simply need to calculate a new equilibrium point and generate new hazards.

To use Lotka-Volterra density dependence, we need to recalculate the equilibrium and then rebuild the hazard functions. The per-capital larval mortality is \(\mu_{L} + \gamma L\), noted by `log = FALSE`

in the options to `equilibrium_lifeycle()`

.

```
# using the same parameters as above, along with the SPN_P object already created
# calculate equilibrium for lotka-volterra dynamics
equilibrium_lifeycle(params = theta, NF = NF, phi = 0.5,
initialCons <-log = FALSE, spn_P = SPN_P,cube=cube)
```

We now follow the same steps as before to make the vector of hazards (setting `log = FALSE`

).

```
# approximate hazards for continous approximation
spn_hazards(spn_P = SPN_P, spn_T = SPN_T, cube = cube,
approx_hazards <-params = initialCons$params, type = "life",
exact = FALSE, tol = 1e-8, verbose = FALSE,
log = FALSE)
# exact hazards for integer-valued state space
spn_hazards(spn_P = SPN_P, spn_T = SPN_T, cube = cube,
exact_hazards <-params = initialCons$params, type = "life",
exact = TRUE, tol = NaN, verbose = FALSE,
log = FALSE)
```

We will use the releases defined above, and again begin with an ODE solution to our system, before testing a different stochastic sampler.

As before, we run an ODE simulation, but this time we will plot the aquatic stages in two ways: looking at their genotypes over time, and looking at the progress through Erlang-distributed stages.

```
# deterministic simulation
sim_trajectory_R(x0 = initialCons$M0, tmax = tmax, dt = dt,
ODE_out <-S = S, hazards = approx_hazards, sampler = "ode",
events = events, verbose = FALSE)
# summarize aquatic stages by genotype
summarize_eggs_geno(out = ODE_out$state, spn_P = SPN_P)
ODE_out_e <- summarize_larvae_geno(out = ODE_out$state, spn_P = SPN_P)
ODE_out_l <- summarize_pupae_geno(out = ODE_out$state, spn_P = SPN_P)
ODE_out_p <-
# add stage name
$stage <- "Egg"
ODE_out_e$stage <- "Larvae"
ODE_out_l$stage <- "Pupae"
ODE_out_p
# plot by genotype
ggplot(data = rbind(ODE_out_e, ODE_out_l,ODE_out_p)) +
geom_line(aes(x = time, y = value, color = genotype)) +
facet_wrap(facets = vars(stage), scales = "free_y") +
theme_bw() +
ggtitle("SPN: ODE Solution - Genotypes")
```

Notice the vastly different equilibrium number for the aquatic stages. These values are parameterized from the adult female value, and the death rates at each stage. Notice that the “releases” here appear to be `Aa`

individuals - this is because our released females (`aa`

) mated with wild-type males (`AA`

), and produced heterozygous (`Aa`

) offspring. Also notice how the releases are smoothed out from Egg to Larvae and Larvae to Pupae stages, due to the death in each stage. We can also see a small dip in the pupal stage, stemming from the over-compensation of density-dependent mortality in the larval stage, before returning to equilibrium.

We can also explore the aquatic stages by dwell time instead of genotype.

```
# summarize aquatic stages by Erlang stage
summarize_eggs_stage(out = ODE_out$state, spn_P = SPN_P)
ODE_out_e <- summarize_larvae_stage(out = ODE_out$state, spn_P = SPN_P)
ODE_out_l <- summarize_pupae_stage(out = ODE_out$state, spn_P = SPN_P)
ODE_out_p <-
# add stage name
$stage <- "Egg"
ODE_out_e$stage <- "Larvae"
ODE_out_l$stage <- "Pupae"
ODE_out_p
# plot by Erlang stage
ggplot(data = rbind(ODE_out_e, ODE_out_l,ODE_out_p)) +
geom_line(aes(x = time, y = value, color = `Erlang-stage`)) +
facet_wrap(facets = vars(stage), scales = "free_y") +
theme_bw() +
ggtitle("SPN: ODE Solution - Erlang Dwell Stage")
```

This plot shows how individuals are moving through each aquatic stage. There are significantly more individuals in the first Erlang-stage, compared to the second Erlang-stage, and we clearly see how the large increase in eggs is attenuated by death through the larval stage, and the effect of density-dependence actually lowering the amount of larvae that make it to the pupal stage.

Finally, for comparison with the logistic density-dependent dynamics above, and a different stochastic sampler below, we plot the adult stages.

```
# summarize females/males
summarize_females(out = ODE_out$state, spn_P = SPN_P)
ODE_out_f <- summarize_males(out = ODE_out$state)
ODE_out_m <-
# add sex for plotting
$sex <- "Female"
ODE_out_f$sex <- "Male"
ODE_out_m
# plot adults
ggplot(data = rbind(ODE_out_f, ODE_out_m)) +
geom_line(aes(x = time, y = value, color = genotype)) +
facet_wrap(facets = vars(sex), scales = "fixed") +
theme_bw() +
ggtitle("SPN: ODE Solution - Adult Stages")
```

Notice the similarity with the logistic density-dependent solution. Since there are no fitness effects or complex dynamics happening here, the solutions are effectively the same.

As an alternative to tau-leaping, which is a discrete stochastic approximation, we investigate a continuous stochastic approximation (Chemical Langevin equation) to the full stochastic model. Theoretically, the CLE approximation is a second-order approximation using continuous state; it is the Fokker-Planck approximation to the Master equation (Kolmogorov Forwards Equations) of the integer valued stochastic process. Because we rely on a relatively simple Euler-Maruyama scheme to solve the system of stochastic differential equations, much like tau-leaping, the user must provide a \(\Delta t\) to the `sim_trajectory_R()`

function, and specify the `cle`

sampler.

Note the change in the sampler, to `cle`

, and the use of approximate hazards.

```
# chemical langevin sampler
sim_trajectory_R(x0 = initialCons$M0, tmax = tmax, dt = dt,
CLE_out <-dt_stoch = dt_stoch, S = S, hazards = approx_hazards,
sampler = "cle", events = events, verbose = FALSE)
# summarize females/males
summarize_females(out = CLE_out$state, spn_P = SPN_P)
CLE_out_f <- summarize_males(out = CLE_out$state)
CLE_out_m <-
# add sex for plotting
$sex <- "Female"
CLE_out_f$sex <- "Male"
CLE_out_m
# plot adults
ggplot(data = rbind(CLE_out_f, CLE_out_m)) +
geom_line(aes(x = time, y = value, color = genotype)) +
facet_wrap(facets = vars(sex), scales = "fixed") +
theme_bw() +
ggtitle("SPN: Chemical Langevin Approximation")
```

Looking at this final set of plots, we immediately see the similarity with the ODE solution above. This is an indication that the stochastic time-step is sufficient for an accurate approximation.