In a previous vignette, we demonstrated how to use the tdsa package to perform time-dependent sensitivity analysis (TDSA), with a simple example that only involved one state variable and relatively few parameters. In this vignette, we will present another example that is more complex. Again, our focus will be on the code syntax (input arguments and output structure) rather than on the biological interpretation of the results. Since this vignette builds upon the previous vignette, we recommend going through the previous vignette first.
Some pathogens are known to have multiple host species. For example, common bee pathogens typically associated with honey bees or bumble bees have been found to infect species from other bee families (Ngor et al, 2020). However, due to variations in host competency and behaviour, pathogen persistence in a natural community may rely on only one or a small subset of species known as a maintenance subcommunity, which can then act as a reservoir for pathogen spillover to the rest of the community (Roberts & Heesterbeek, 2020).
It is useful to regard the collection of species outside the maintenance subcommunity as a sink subcommunity (in the context of source-sink dynamics), since the pathogen cannot persist in this subcommunity and is only maintained via spillover. Within this subcommunity, the overall risk to a species of concern depends on both the transmission network structure, as well as the risk of spillover from the reservoir to each node (Ng et al, 2023). If spillover risk varies seasonally, the most effective management intervention (in terms of which node or transmission link to target) may depend on the timing within the season. In Ng et al. (in press), we used a model to demonstrate how TDSA can help identify the most effective course of action at different times. We present this model below.
The model comprises a sink subcommunity with five species. Interspecific transmission only occurs between “nearest-neighbour” species, leading to a linear network as shown in the figure below. Such a linear network might occur if species interactions are strongly structured by trait-matching (Truitt et al., 2019), so each species only interacts with species that are adjacent along a one-dimensional trait space. We assume that exogenous spillover can only occur during the active season, so the disease is re-introduced each year at the start of the season, and dies out after the end of the season.
Species 5 is a species of management concern because it provides an important ecosystem service but has a vulnerable population, due to its low birth rate and high disease-induced mortality. To reduce disease transmission to Species 5, we may want to cull an intermediate species or dilute a transmission link, so we will use TDSA to help us identify what and when to target.
The state variables are Sj and Ij, the number of susceptible and infected individuals in species j, so Nj ≡ Sj + Ij gives the species population size. The dynamic equations are first-order ordinary differential equations $$\begin{aligned} \frac{dS_j}{dt} &=B_j N_j (1 - a_j N_j) - S_j \left(\sigma_j + \sum_{k=1}^m b_{j,k} I_k \right) - \mu_j S_j + \gamma_j I_j,\\ \frac{dI_j}{dt} &= S_j \left(\sigma_j + \sum_{k=1}^m b_{j,k} I_k \right) - (\mu_j + \nu_j + \gamma_j) I_j, \end{aligned}$$ where Bj is the unregulated per-capita birth rate, and aj the coefficient for intraspecific reproductive competition. For a susceptible individual, μj is the mortality rate, while for an infected individual, νj is the additional disease-induced mortality rate and γj the recovery rate. bj, k is the transmission coefficient from species k to j, and σj the per-capita exogeneous spillover rate.
Each year, Species 5 provides the ecosystem service over an active season from t0 = 0 to t1 = 5 in units of lifespan. (We assume that all five species are relatively short-lived, so five lifespans span one active season.) The reward J represents the total value of this service, and is given by $$J = \int_{t_0}^{t_1} \underbrace{\sum_{j=1}^5 W_j N_j}_{\substack{\text{reward}\\\text{integrand}}}\, dt+ \underbrace{\sum_{j=1}^5 V_j N_j(t_1)}_{\text{terminal payoff}},$$ where Wj is the per-capita rate of contribution to the service. We also include a terminal payoff to ascribe a value to maintaining a large population at the end of the season, since this will affect the population size next season.
The tdsa package is already on CRAN, so it can be
easily installed using the command
install.packages("tdsa")
. For the rest of this vignette, we
assume that the package has already been installed.
Before proceeding further, we load the tdsa package.
Again, we need to prepare the input arguments for the function
state_sens
. However, now because there are more parameters,
we need to think about how we want the parameter values to be assigned
in dynamic_fn
. The following options are available.
Using the argument parms
as was done previously.
Only parameters assigned this way are recognised by the function
parm_sens
used to calculate parameter
sensitivities.
Using a different argument of our choice. This is permitted,
since the ...
in function(t, y, parms, ...)
means that we are free to introduce any additional arguments we
want.
Within the environment of the function itself.
In the global environment.
Here, we will use the first two options simultaneously, by splitting
the model parameters into two lists parms
and
parms2
. parms
contains parameters whose
sensitivities we are interested in, while parms2
is a new
input argument of dynamic_fn
that contains all other
parameters. Note that this split is not really necessary—it is perfectly
fine to assign all parameters using parms
alone. However,
keep in mind that the function parm_sens
will calculate the
sensitivities for all parameters assigned using
parms
, so we waste a lot of computation time if there are
many model parameters that we are not interested in.
parms = list(
sigma = c(0.2, 0, 0, 0, 0), # Per-capita exogenous spillover rate; only nonzero for Species 1.
mu = c(1, 1, 1, 1, 1), # Mortality rate of a susceptible individual.
b = matrix(c(1.976537, 1.976537, 0.000000, 0.000000, 0.000000,
1.976537, 1.976537, 1.976537, 0.000000, 0.000000,
0.000000, 1.976537, 1.976537, 1.976537, 0.000000,
0.000000, 0.000000, 1.976537, 1.976537, 1.976537,
0.000000, 0.000000, 0.000000, 1.976537, 1.976537),
nrow=5, byrow=TRUE) # Matrix of transmission coefficients; values chosen so that disease R_0 = 0.9.
)
parms2 = list(
B = c(5, 5, 5, 5, 1.02), # Unregulated per-capita birth rate; low for Species 5.
a = c(0.8, 0.8, 0.8, 0.8, 0.02), # Intraspecific competition; chosen so all species have disease-free carrying capacity = 1.
nu = c(0, 0, 0, 0, 5), # Disease-induced mortality of an infected individual; only nonzero for Species 5.
gamma = c(5, 5, 5, 5, 0) # Recovery rate of an infected individual; zero for species 5.
)
Now we specify dynamic_fn
. Since there are five species,
each with two compartments (susceptible and infected), the input
argument y
is a numeric vector of length 10, corresponding
to the state vector $$\vec{y}=\begin{pmatrix}\vec{S}\\\vec{I}\end{pmatrix},
\quad \vec{S}=(S_1, S_2, S_3, S_4, S_5)^{\mathsf{T}}, \quad \vec{I}=
(I_1, I_2, I_3, I_4, I_5)^{\mathsf{T}}.$$ dynamic_fn
must return a list, the first element a numeric vector of length 10
corresponding to the RHS of the dynamic equations (in vector form) $$\frac{d\vec{y}}{dt}
=\begin{pmatrix}\frac{d\vec{S}}{dt}\\ \frac{d\vec{I}}{dt}\end{pmatrix}
=\begin{pmatrix} \vec{B}
\circ\vec{N}\circ(1-\vec{a}\circ\vec{N})-\vec{S}\circ(\vec{\sigma}+\mathbf{b}\vec{I})
-\vec{\mu}\circ\vec{S}+\vec{\gamma}\circ\vec{I} \\
\vec{S}\circ(\vec{\sigma}+\mathbf{b}\vec{I})
-(\vec{\mu}+\vec{\nu}+\vec{\gamma})\circ\vec{I} \end{pmatrix}.$$
Here, ∘ corresponds to the Hadamard
(element-wise) product, while all the model parameters have been grouped
into vectors and matrices such as B⃗ or b.
dynamic_fn = function(t, y, parms, parms2){
# To make the lines below easier to read, we "extract" each coefficient from the parameter lists.
mu = parms[["mu"]]
b = parms[["b"]]
sigma = parms[["sigma"]]
B = parms2[["B"]]
a = parms2[["a"]]
nu = parms2[["nu"]]
gamma = parms2[["gamma"]]
# To make the lines below easier to read, we "extract" the susceptible and infected parts from the state vector.
SS = y[1:5 ]
II = y[6:10]
# Calculate the species population size.
NN = SS + II
# RHS of the dynamic equations.
dSS = B * NN * (1 - a*NN) - SS * (sigma + b%*%II) - mu * SS + gamma * II
dII = SS * (sigma + b%*%II) - (mu + nu + gamma) * II
return( list( c(dSS, dII) ) )
}
Next, we need to specify reward_fn
and
terminal_fn
for the reward integrands and terminal
payoffs.
reward_fn = function(t, y){
# Parameters.
W = c(0, 0, 0, 0, 1) # Per-capita rate of contribution to ecosystem service; only nonzero for Species 5.
# Split the state vector.
SS = y[1:5]
II = y[6:10]
# Return the reward integrand.
NN = SS + II
return( sum(W * NN) )
}
terminal_fn = function(y){
# Parameters.
V = c(0, 0, 0, 0, 1) # Per-capita terminal payoff; only nonzero for Species 5.
# Split the state vector.
SS = y[1:5]
II = y[6:10]
# Return the terminal payoff.
NN = SS + II
return( sum(V*NN) )
}
Next, we need to specify y_0
, the initial conditions
y⃗(t0). We
assume that the subcommunity starts off disease-free at the carrying
capacity Kj = (1 − μj/Bj)/aj.
SS_0 = (1-parms[["mu"]]/parms2[["B"]])/parms2[["a"]] # At carrying capacity.
II_0 = c(0, 0, 0, 0, 0) # Disease-free.
y_0 = c(SS_0, II_0)
Finally, since this is a continuous-time model, we need to discretise
the continuous time interval between the initial and final times t0 = 0 and t1 = 5. We choose 1001
time steps, so the step size is 0.005. When using the function
seq
, it is more “fail-safe” to specify the argument
length.out
(number of time steps) instead of
by
(step size), to ensure that t1 is always the final
time step.
We calculate time-dependent state sensitivities using the function
state_sens
. Since parms2
is a new argument of
dynamic_fn
, we use the argument
dynamic_fn_arglist
to specify it in
state_sens
.
state_sens_out = state_sens(
model_type = "continuous",
dynamic_fn = dynamic_fn,
parms = parms,
reward_fn = reward_fn,
terminal_fn = terminal_fn,
y_0 = y_0,
times = times,
dynamic_fn_arglist = list(parms2 = parms2),
verbose = FALSE
)
As before, the output is a list that contains some input arguments,
as well as the matrices state
and tdss
. These
matrices now have 10 columns, since there are now 10 state variables,
and also the following 10 state sensitivities $$\vec{\lambda}=\begin{pmatrix}\vec{\lambda}_S\\\vec{\lambda}_I
\end{pmatrix}, \quad
\vec{\lambda}_S=(\lambda_{S_1},\lambda_{S_2},\lambda_{S_3},\lambda_{S_4},\lambda_{S_5})^{\mathsf{T}},\quad
\vec{\lambda}_I=(\lambda_{I_1},\lambda_{I_2},\lambda_{I_3},\lambda_{I_4},\lambda_{I_5})^{\mathsf{T}}.$$
## List of 7
## $ model_type : chr "continuous"
## $ dynamic_fn :function (t, y, parms, parms2)
## $ parms :List of 3
## ..$ sigma: num [1:5] 0.2 0 0 0 0
## ..$ mu : num [1:5] 1 1 1 1 1
## ..$ b : num [1:5, 1:5] 1.98 1.98 0 0 0 ...
## $ dynamic_fn_arglist:List of 1
## ..$ parms2:List of 4
## .. ..$ B : num [1:5] 5 5 5 5 1.02
## .. ..$ a : num [1:5] 0.8 0.8 0.8 0.8 0.02
## .. ..$ nu : num [1:5] 0 0 0 0 5
## .. ..$ gamma: num [1:5] 5 5 5 5 0
## $ times : num [1:1001] 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 ...
## $ state : num [1:1001, 1:10] 1 0.999 0.998 0.997 0.996 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : NULL
## .. ..$ : chr [1:10] "1" "2" "3" "4" ...
## $ tdss : num [1:1001, 1:10] -0.0396 -0.0397 -0.0397 -0.0398 -0.0399 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : NULL
## .. ..$ : chr [1:10] "1" "2" "3" "4" ...
In the left panel below, we plot the disease prevalence Pj ≡ Ij/Nj
in each species, calculated from the state
matrix. Since
the disease has to be relayed down the transmission chain, not
surprisingly the prevalence is a lot lower in Species 5, and also the
rise is more delayed. In the right panel, we plot −λNj,
defined as λNj = (1 − Pj)λSj + PjλIj.
This represents the sensitivity of the reward J to random culling of Species j, regardless of infection status.
(It makes no sense to cull Species 5 so we have excluded it from the
plot.) We find that early in the season, it is more effective to cull a
species that is upstream along the transmission chain, whereas the
reverse is true late in the season. This is indeed intuitively what we
expect—culling a downstream species too early is ineffective since the
population would have recovered somewhat by the time the disease reaches
the species, while culling an upstream species too late is also
ineffective since the disease is already being transmitted among the
downstream species.
# Calculate the derived quantities to be plotted.
# Disease prevalence:
SS = state_sens_out[["state"]][, 1:5] # Number of susceptible individuals.
II = state_sens_out[["state"]][, 6:10] # Number of infected individuals.
NN = SS + II # Species population size.
PP = II / NN # Disease prevalence.
# State sensitivities:
lambda_SS = state_sens_out[["tdss"]][, 1:5] # State sensitivities for S_j.
lambda_II = state_sens_out[["tdss"]][, 6:10] # State sensitivities for I_j.
lambda_NN = (1 - PP) * lambda_SS + PP * lambda_II
# Generate the plots.
palette = c("#42049EFF", "#900DA4FF", "#CC4678FF", "#F1844BFF", "#FCCE25FF") # Colour palette.
par(mar=c(5,6.7,4,1.7)+0.1, mfrow=c(1,2)) # Set graphical parameters.
# Plot PP.
plot(NA, xlim=c(0, 5), ylim=c(0.00067, 1), log="y", xaxs="i", yaxs="i", yaxt="n",
main="Infection dynamics",
xlab="Time (in units of lifespan)",
ylab="Fraction of population infected (log scale)",
cex.main=2, cex.lab=1.5, cex.axis=1.5)
axis(side=2, cex.axis=1.5, at=10^((-4):(0)), labels=c(0.0001, 0.001, 0.01, 0.1, 1))
for(i in 1:5){
lines(times, PP[,i], lwd=3, col=palette[i])
}
legend("topright", legend=paste("Species",1:5), lwd=3, col=palette[1:5], bty="n")
# Plot lambda_NN.
plot(NA, xlim=c(0, 5), ylim=c(0, 0.11), xaxs="i", yaxs="i",
main="Time-dep. state sensitivities",
xlab="Time (in units of lifespan)",
ylab=bquote(atop("Sensitivity of J to sudden decrease", "in species population "~(-lambda[N["j"]]))),
cex.main=2, cex.lab=1.5, cex.axis=1.5)
for(i in (1:4)){
lines(times, -lambda_NN[,i], lwd=3, col=palette[i])
}
legend("topright", legend=paste("Species",1:4), lwd=3, col=palette[1:4], bty="n")
Again, calculating the parameter sensitivities is easy, since it only
involves using the output of state_sens
as the input
argument of the function parm_sens
. To speed up the
calculation of numerical derivatives, we will also use the optional
argument numDeriv_arglist
to reduce the number of
Richardson iterations to two (see the numDeriv package
help for details).
parm_sens_out = parm_sens(state_sens_out = state_sens_out,
numDeriv_arglist = list(method.args=list(r=2)),
verbose = FALSE)
As before, the output is a list containing times
, and an
object called tdps
that gives the time-dependent parameter
sensitivity. In this example, parms
is a list containing
two vectors sigma
and mu
and a matrix
b
, so tdps
is a list containing two matrices
and a 3-index array that give the parameter sensitivities of
sigma
, mu
and b
respectively.
## List of 2
## $ times: num [1:1001] 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 ...
## $ tdps :List of 3
## ..$ sigma: num [1:1001, 1:5] -0.695 -0.692 -0.689 -0.687 -0.684 ...
## ..$ mu : num [1:1001, 1:5] 0.0396 0.0403 0.0411 0.0418 0.0426 ...
## ..$ b : num [1:1001, 1:5, 1:5] 0 -0.000685 -0.00135 -0.001996 -0.002624 ...
## .. ..- attr(*, "dimnames")=List of 3
## .. .. ..$ : NULL
## .. .. ..$ : NULL
## .. .. ..$ : NULL
In the left panel below, we plot (excluding Species 5) κμj, the parameter sensitivities of μj. We note that the shape is similar to the negative of the state sensitivities −λNj, which is not surprising because a brief spike in mortality leads to a sudden decrease in the population size. In the right panel, we plot −κbj + 1, j, the negative of the parameter sensitivities of the forward transmission coefficients b2, 1, b3, 2, etc. We find that early in the season, it is more effective to target an upstream transmission link, whereas the reverse is true late in the season. This is again an intuitive result.
# Extract the list containing the parameter sensitivities.
tdps = parm_sens_out[["tdps"]]
# Parameter sensitivities for the mortality rate of a susceptible individual.
kappa_mu = tdps[["mu"]]
# Parameter sensitivities for the forward transmission rate b_{2,1}, b_{3,2}, etc.
# These are given by tdps[["b"]][,2,1], tdps[["b"]][,3,2], etc.
# A more systematic way to extract these elements is to use mapply.
kappa_b = mapply(function(i,j){tdps[["b"]][,i,j]}, 2:5, 1:4)
# Generate the plots.
palette = c("#42049EFF", "#900DA4FF", "#CC4678FF", "#F1844BFF", "#FCCE25FF") # Colour palette.
par(mar=c(5,6.7,4,1.7)+0.1, mfrow=c(1,2)) # Set graphical parameters.
# Plot kappa_mu.
plot(NA, xlim=c(0, 5), ylim=c(0, 0.11), xaxs="i", yaxs="i",
main="Time-dep. parm. sensitivities",
xlab="Time (in units of lifespan)",
ylab=bquote(atop("Sensitivity of J to brief increase", "in mortality of susceptibles "~(kappa[mu["j"]]))),
cex.main=2, cex.lab=1.5, cex.axis=1.5)
for(i in 1:4){
lines(times, kappa_mu[,i], lwd=3, col=palette[i])
}
legend("topright", legend=paste("Species",1:4), lwd=3, col=palette[1:4], bty="n")
# Plot kappa_b
plot(NA, xlim=c(0, 5), ylim=c(0, 0.077), xaxs="i", yaxs="i",
main="Time-dep. parm. sensitivities",
xlab="Time (in units of lifespan)",
ylab=bquote(atop("Sensitivity of J to brief decrease", "in forward transmission "~(-kappa[b["j+1, j"]]))),
cex.main=2, cex.lab=1.5, cex.axis=1.5)
for(i in 1:4){
lines(times, -kappa_b[,i], lwd=3, col=palette[i])
}
legend("topright", legend=paste("Species", 1:4, "to", 2:5), lwd=3, col=palette[1:4], bty="n")
Finally, we plot −κσ1, the negative of the parameter sensitivity of σ1, the exogenous spillover rate to Species 1. (We omit the other σj ≠ 1 since the model assumes no exogenous spillover to other species.) We see that it is most effective to target the exogenous disease source and reduce the spillover rate σ1 as early as possible. Later in the season the disease is already being transmitted among the “downstream” species, so reducing the exogenous spillover to Species 1 (which is upstream) becomes a lot less effective.
# Parameter sensitivity for the exogenous spillover rate to Species 1.
kappa_sigma_1 = tdps[["sigma"]][,1]
# Generate the plot.
par(mar=c(5,6.7,4,1.7)+0.1) # Set graphical parameters.
plot(times, -kappa_sigma_1, xlim=c(0, 5), ylim=c(0, 0.72), xaxs="i", yaxs="i",
type="l", lwd=3,
main="Time-dep. parm. sensitivity",
xlab="Time (in units of lifespan)",
ylab=bquote(atop("Sensitivity of J to brief decrease", "in exogeneous spillover to Species 1 "~(kappa[sigma["1"]]))),
cex.main=2, cex.lab=1.5, cex.axis=1.5)
Ng, W. H., Myers, C. R., McArt, S., & Ellner, S. P. (2023). Predicting and Controlling Spillover in Multispecies Disease Transmission Networks: Steady-State Analysis. American Naturalist, 201, 880–894. doi: 10.1086/724009.
Ng, W. H., Myers, C. R., McArt, S., & Ellner, S. P. (in press). A time for every purpose: using time-dependent sensitivity analysis to help understand and manage dynamic ecological systems. American Naturalist. doi: 10.1086/726143. eprint doi: 10.1101/2023.04.13.536769.
Ng, W. H., Myers, C. R., McArt, S., & Ellner, S. P. (in review). tdsa: An R package to perform time-dependent sensitivity analysis. Methods in Ecology and Evolution.
Ngor, L., Palmer-Young, E. C., Burciaga Nevarez, R., Russell, K. A., Leger, L., Giacomini, S. J., Pinilla-Gallego, M. S., Irwin, R. E., & McFrederick, Q. S. (2020). Cross-infectivity of honey and bumble bee-associated parasites across three bee families. Parasitology, 12, 1290–1304. doi: 10.1017/S0031182020001018
Roberts, M. G., & Heesterbeek, J. A. P. (2020). Characterizing reservoirs of infection and the maintenance of pathogens in ecosystems. Journal of The Royal Society Interface, 17, 2019054. doi: 10.1098/rsif.2019.0540.
Truitt, L. L., McArt, S. H., Vaughn, A. H., & Ellner, S. P. (2019). Trait-based modeling of multi-host pathogen transmission: Plant-pollinator networks. American Naturalist, 193, E149–E167. doi: 10.1086/702959.