Simulation of Fully Specified SPN Model
Deterministic: Decoupled ODE model
To first demonstrate that the system is indeed at equilibrium, we show the deterministic solution using decoupled ODE sampling without events.
dt <- 1
ode_out <- sim_trajectory_R_decoupled(
x0 = initialCons$M0,
h0 = initialCons$H,
SPN_P = SPN_P,
theta = theta,
tmax = 100,
inf_labels = SPN_T$inf_labels,
dt = dt,
S = S,
hazards = approx_hazards,
sampler = "ode-decoupled",
events = NULL,
verbose = FALSE,
human_ode = "Imperial",
cube = cube
)
# summarize females/humans by genotype
ode_female <- summarize_females_epi(out = ode_out$state, spn_P = SPN_P)
ode_humans <- summarize_humans_epiImperial(out = ode_out$state, index=1)
# plot
ggplot(data = ode_female) +
geom_line(aes(x = time, y = value, color = inf)) +
facet_wrap(~ genotype, scales = "free_y") +
theme_bw() +
ggtitle("SPN: ODE Decoupled Approximation - Mosquito")
ggplot(data = ode_humans) +
geom_line(aes(x = time, y = value, color = inf)) +
facet_wrap(~ genotype, scales = "free_y") +
theme_bw() +
ggtitle("SPN: ODE Decoupled Approximation - Human")
### Stochastic: Tau Leaping Solutions {#soln_tau}
As a further example, we run a single stochastic realization of the
same simulation, using the tau sampler with \(\Delta t = 1\), approximating 10 jumps per
day. As the adult male mosquitoes do not contribute to infection
dynamics, we will only view the adult female mosquito and human dynamics
here. We do not show events in this simulation, as genotype-specific
dynamics have not yet been implemented.
# delta t - one day
dt_stoch <- 0.1
dt <- 1
# run ode-decoupled simulation
tau_out <- sim_trajectory_R_decoupled(
x0 = initialCons$M0,
h0 = initialCons$H,
SPN_P = SPN_P,
theta = theta,
tmax = 100,
inf_labels = SPN_T$inf_labels,
dt = dt,
dt_stoch = dt_stoch,
S = S,
hazards = approx_hazards,
sampler = "tau-decoupled",
events = NULL,
verbose = FALSE,
human_ode = "Imperial",
cube = cube,
maxhaz = 1e12
)
# summarize females/humans by genotype
tau_female <- summarize_females_epi(out = tau_out$state, spn_P = SPN_P)
tau_humans <- summarize_humans_epiImperial(out = tau_out$state, index=1)
# plot
ggplot(data = tau_female) +
geom_line(aes(x = time, y = value, color = inf)) +
facet_wrap(~ genotype, scales = "free_y") +
theme_bw() +
ggtitle("SPN: Tau-leaping Decoupled Approximation - Mosquito")
ggplot(data = tau_humans) +
geom_line(aes(x = time, y = value, color = inf)) +
facet_wrap(~ genotype, scales = "free_y") +
theme_bw() +
ggtitle("SPN: Tau-leaping Decoupled Approximation - Human")
In the stochastic simulation, we see that both values drift around
equilibrium. Varying delta_t can provide a visual
understanding of how the underlying stochastic dynamics converge to the
ODE solution as delta_t -> 0.
Finally, we show a stochastic simulation with events.
Similar to previous simulations, we will release 50 adult females
with homozygous recessive alleles 5 times, every 10 days, but starting
at day 20. Remember, 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. To see the effect on transmission of
genotype-specific releases, vary the c and b0
parameters in the inheritance cube (see line 187 of this vignette).
r_times <- seq(from = 20, length.out = 5, by = 10)
r_size <- 50
events <- data.frame("var" = paste0("F_", cube$releaseType, "_", cube$wildType, "_S"),
"time" = r_times,
"value" = r_size,
"method" = "add",
stringsAsFactors = FALSE)
tau_out <- sim_trajectory_R_decoupled(
x0 = initialCons$M0,
h0 = initialCons$H,
SPN_P = SPN_P,
theta = theta,
tmax = 100,
inf_labels = SPN_T$inf_labels,
dt = dt,
dt_stoch = dt_stoch,
S = S,
hazards = approx_hazards,
sampler = "tau-decoupled",
events = events,
verbose = FALSE,
human_ode = "Imperial",
cube = cube,
maxhaz = 1e12
)
# summarize females/humans by genotype
tau_female <- summarize_females_epi(out = tau_out$state, spn_P = SPN_P)
tau_humans <- summarize_humans_epiImperial(out = tau_out$state, index=1)
# plot
ggplot(data = tau_female) +
geom_line(aes(x = time, y = value, color = inf)) +
facet_wrap(~ genotype, scales = "free_y") +
theme_bw() +
ggtitle("SPN: Tau-leaping Decoupled Approximation - Mosquito")
ggplot(data = tau_humans) +
geom_line(aes(x = time, y = value, color = inf)) +
facet_wrap(~ genotype, scales = "free_y") +
theme_bw() +
ggtitle("SPN: Tau-leaping Decoupled Approximation - Human")