\documentclass[11pt]{article}
\usepackage{fullpage}
\usepackage{setspace}
\usepackage{parskip}
\usepackage{titlesec}
\usepackage[section]{placeins}
\usepackage{xcolor}
\usepackage{breakcites}
\usepackage{lineno}
\usepackage{hyphenat}
\setlength\columnsep{25pt}
\PassOptionsToPackage{hyphens}{url}
\usepackage[colorlinks = true,
linkcolor = blue,
urlcolor = blue,
citecolor = blue,
anchorcolor = blue]{hyperref}
\usepackage{etoolbox}
\makeatletter
\patchcmd\@combinedblfloats{\box\@outputbox}{\unvbox\@outputbox}{}{%
\errmessage{\noexpand\@combinedblfloats could not be patched}%
}%
\makeatother
\usepackage{natbib}
\renewenvironment{abstract}
{{\bfseries\noindent{\abstractname}\par\nobreak}\footnotesize}
{\bigskip}
\titlespacing{\section}{0pt}{*3}{*1}
\titlespacing{\subsection}{0pt}{*2}{*0.5}
\titlespacing{\subsubsection}{0pt}{*1.5}{0pt}
\usepackage{authblk}
\usepackage{graphicx}
\usepackage[space]{grffile}
\usepackage{latexsym}
\usepackage{textcomp}
\usepackage{longtable}
\usepackage{tabulary}
\usepackage{booktabs,array,multirow}
\usepackage{amsfonts,amsmath,amssymb}
\providecommand\citet{\cite}
\providecommand\citep{\cite}
\providecommand\citealt{\cite}
% You can conditionalize code for latexml or normal latex using this.
\newif\iflatexml\latexmlfalse
\providecommand{\tightlist}{\setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}%
\AtBeginDocument{\DeclareGraphicsExtensions{.pdf,.PDF,.eps,.EPS,.png,.PNG,.tif,.TIF,.jpg,.JPG,.jpeg,.JPEG}}
\usepackage[utf8]{inputenc}
\usepackage[T2A]{fontenc}
\usepackage[ngerman,greek,polish,english]{babel}
\newcommand{\ie}{\emph{i.e.},~}
\newcommand{\eg}{\emph{e.g.},~}
\newcommand{\etc}{\emph{etc.}~}
\newcommand{\ddt}[1]{\frac{\textrm{d} #1}{\textrm{d}t}}
\renewcommand{\figurename}{Fig.}
\newcommand{\beginsupplement}{%
\setcounter{table}{0}
\renewcommand{\thetable}{S\arabic{table}}%
\setcounter{figure}{0}
\renewcommand{\thefigure}{S\arabic{figure}}%
\setcounter{section}{0}
\renewcommand{\thesection}{S\arabic{section}}%
}
\begin{document}
\title{Multi-strain disease dynamics on metapopulation networks}
\author[1,2]{Matthew Michalska-Smith}%
\author[1]{Kimberly VanderWaal}%
\author[1]{Montserrat Torremorell}%
\author[1]{Cesar Corzo}%
\author[1]{Meggan E Craft}%
\affil[1]{Department of Veterinary Population Medicine, University of Minnesota, St. Paul, MN USA}%
\affil[2]{Department of Plant Pathology, University of Minnesota, St. Paul MN, USA}%
\vspace{-1em}
\date{\today}
\begingroup
\let\center\flushleft
\let\endcenter\endflushleft
\maketitle
\endgroup
\linenumbers
\doublespacing
\newpage
\selectlanguage{english}
\begin{abstract}
Many pathogens have clusters of variation in their genotypes~that we
refer to as strain structure. Importantly, when considering related
pathogen strains, host immunity to one strain~is often neither
independent from nor equivalent to immunity to other strains. This
partial cross-reactive immunity can thus allow repeated infection with
(different strains of) the same pathogen and shapes disease dynamics
across a population, in turn influencing the effectiveness of
intervention strategies. To better understand the dynamics governing
multi-strain pathogens in complex landscapes, we combine two frameworks
well-studied in their own right: multi-strain disease dynamics and
metapopulation network structure. We simulate the dynamics of a
multi-strain disease on a network of populations connected by migration
and characterize the joint effects of disease model parametrization and
network structure on these dynamics. We find that the movement of
(partially) immune individuals tends to have a larger impact than the
movement of infectious individuals, dampening infection dynamics in
populations further along a chain. When disease parameters differ
between populations, we find that dynamics can propagate from one
population to another, alternatively stabilizing or destabilizing
destinations populations based on the dynamics of origin populations. In
addition to providing novel insights into the role of host movement on
disease dynamics, this work provides a framework for future predictive
modelling of multi-strain diseases across generalized~population
structures.%
\end{abstract}%
\sloppy
\section{Introduction}
{\label{776665}}
Many of the most impactful infectious diseases that affect humans
(influenza, malaria, human papillomavirus,~\emph{etc.}), livestock
(porcine reproductive and respiratory syndrome, foot-and-mouth
disease,~\emph{etc.}), and wildlife (anthrax, plague,~\emph{etc.}) have
clusters in their population-genetic variability that we classify as
strains. This variation in pathogen genotype is often associated with
differences in phenotype, which can have profound effects on the
efficacy of host immune defenses. While the human immune system is
usually capable of preventing re-infection---\emph{i.e.} infection with
something to which it has been previously exposed---sufficient,
divergent evolution among pathogen strains can reduce the ability of the
host to recognize, and thus mount an immunological response to,
subsequent exposures. In some cases, this change is not sufficient to
completely avoid recognition by the host's immune system, yielding an
immune response that is neither as strong as it would be in the case of
re-exposure to the same strain, nor as weak as in the case of exposure
to a novel pathogen. This partial cross-reactive immunity can likewise
lead to reduced transmissibility, affecting disease dynamics across the
population.
While the study of multi-strain diseases goes back
decades~\hyperref[csl:1]{(1}; \hyperref[csl:2]{2)}, the resulting modelling framework has not
yet been generalized to a collection of sub-populations connected
through host movement,~\emph{i.e.} a metapopulation (but
see~\selectlanguage{ngerman}\hyperref[csl:3]{(3)}). Initially introduced through the concepts of
island biogeography~\hyperref[csl:4]{(4)}, the network approach of
metapopulations can be applied to a variety of systems, including human
movement between cities, livestock transport between farms, and wildlife
living in fragmented natural habitats. In each case, there exist
relatively high-density areas which are connected to one another through
a network of individuals' movement. A metapopulation framework allows
the application of network analyses to characterize patterns of
connection within the larger system, and can provide unique insights
across scales.
Historically, metapopulation studies have been been divided into two
main camps: those that model within-population dynamics and ``cell
occupancy'' models. The latter of which, where only the presence or
absence of a given species within a population is recorded in a given
timepoint~\hyperref[csl:5]{(5)}, has received much more theoretical
attention. Importantly, cell occupancy models rest on an assumption of
temporal separation in which local dynamics occur on a timescale that
can be treated as instantaneous relative to that of the
between-population dynamics~\hyperref[csl:6]{(6)}. When considering
pathogens in systems with relatively high migration rates, however, this
assumption rarely holds, and the presence-absence approach can
significantly limit model accuracy~\hyperref[csl:7]{(7}; \hyperref[csl:8]{8}; \hyperref[csl:9]{9)}.
The presence of metapopulation structure has been repeatedly associated
with increased stability~\hyperref[csl:10]{(10}; \hyperref[csl:11]{11}; \hyperref[csl:12]{12)}. This is due in part to the
ability of migration between asynchronous populations to rescue
temporarily low density populations from extinction~\hyperref[csl:13]{(13)}.
This is particularly relevant when populations are undergoing cyclical
or chaotic dynamics, where repeated instances of low density are
generally considered to be at greater risk of extinction than a
population maintaining steady state dynamics~\hyperref[csl:14]{(14}; \hyperref[csl:15]{15)}.
Here, we build on the strain theory of host-pathogen systems proposed
by~\hyperref[csl:16]{(16)}, considering a scenario where a collection of
populations undergoing local dynamics are furthermore interconnected
through the movement of individuals between populations. We simulate
disease dynamics on this system, characterizing the effects of
parametrization and network structure on these dynamics. This work is
divided into three sections: first, we explore a case of interconnected
populations which have been parametrized to display identical dynamics
in the absence of host migration. Second, we consider cases where
parameters differ between populations. Finally, we explore the role of
network structure on disease dynamics in larger networks of connected
populations.
\section{Methods}
{\label{312439}}
\subsection{Model framework for one
population}
{\label{826836}}
We work from a system of ordinary differential equations which delineate
a population into classes based on current and past exposure to
different strains of a pathogen. Pathogens with strain structure can
differ in both the number of strains and the level of cross-reactive
immunity afforded by past exposure to similar strains. To model the
number of strains, we signify a strain i = \{x\textsubscript{1},
x\textsubscript{2}, \ldots{}, x\textsubscript{n}\} as a set of n loci,
each of which can take on a finite number of alleles. For instance, a
pathogen with two loci (a and b) and two alleles at each loci has a
total of four potential strains: \{a\textsubscript{1},
b\textsubscript{1}\}, \{a\textsubscript{1}, b\textsubscript{2}\},
\{a\textsubscript{2}, b\textsubscript{1}\}, \{a\textsubscript{2},
b\textsubscript{2}\}. For cross-reactive immunity, we use a parameter \selectlanguage{greek}γ
\selectlanguage{english}which indicates the degree of reduced susceptibility a host has to
strains that are similar to (\emph{i.e.} strains that share at least one
allele with one another) past exposures. Importantly, in this framework,
the number of strains is fixed and finite. While strains may go extinct
over time, there is no process for the generation of new strains or to
re-introduce strains that had previously gone extinct (but
see~\hyperref[csl:16]{(16)}).
The model consists of sets of three nested equations (one set for each
strain i): y, z, and w. See~\hyperref[csl:17]{(17)} for a more comprehensive
discussion of the model framework, including a graphical representation.
y\textsubscript{i} represents the proportion of the population currently
infectious with strain i. z\textsubscript{i} represents the proportion
of the population that has been exposed to strain i. These individuals
harbor complete immunity to future infections with strain i and include
those currently infected,~\emph{i.e.~}y\textsubscript{i}, those that
have recovered but were previously infectious, and those that were
exposed, but protected from becoming infectious due to partial
cross-protective immunity. Finally, w\textsubscript{i} represents the
proportion of the population which has been exposed to any strain j
which has at least one allele in common with strain i (including strain
i itself),~\selectlanguage{ngerman}\emph{i.e.~\(j\cap i\ne\text{Ø}\)}. These individuals have at
least partial immunity to strain i.~\emph{N.b.} these equations are
nested such that any individual in y\textsubscript{i} is also in
z\textsubscript{i} and any individual in z\textsubscript{i} is also in
w\textsubscript{i}, and~\(y_i\le z_i\le w_i\ \forall\ \) strain i. In traditional
Susceptible-Infected (SI), Susceptible-Infected-Recovered
(SIR),~\emph{etc.} single-strain mathematical frameworks: the y class is
analogous to the I class, while w and z are composed of combinations of
I and R classes. The susceptible population is not modelled explicitly
in this framework.
Explicitly, these three equations (for a given strain i) are:
\begin{equation}
\label{eqn:onepop}
\begin{aligned}
\ddt{y_i} &= \beta \left( (1 - w_i) + (1 - \gamma) (w_i - z_i) \right) y_i - \sigma y_i - \mu y_i\\
\ddt{z_i} &= \beta (1 - z_i) y_i - \mu z_i\\
\ddt{w_i} &= \beta (1 - w_i) \selectlanguage{ngerman}\sum_{j \ni j \cap i \neq \text{Ø}} y_j - \mu w_i
\end{aligned}
\end{equation}
As above, we denote strains as subscripts and, in the equation for
w\textsubscript{i}, we sum over all strains j which share at least one
allele with the focal strain i. \selectlanguage{greek}β, σ, \selectlanguage{english}and \selectlanguage{greek}μ \selectlanguage{english}are the infection, recovery,
and death rates, respectively.~\selectlanguage{greek}γ \selectlanguage{english}(as mentioned above) is an indicator of
the level of cross-reactive immunity gained by prior exposure to alleles
in the target strain. Note that while we depict only one value per
demographic parameter (\emph{i.e.} all strains are functionally
equivalent) for clarity of notation, these values could also be written
to vary by strain (\emph{e.g.} \selectlanguage{greek}β\selectlanguage{english}\textsubscript{i}).
Immunity in this framework is non-waning: exposure to a strain yields
consistent protection from future infection over the lifespan of the
individual. Moreover, this protection is trichotomous: an individual can
either have no protection from a given strain (it has not seen any of
the alleles before), complete protection (it has seen this exact
combination of alleles before), or a set point in-between according to
the parameter \selectlanguage{greek}γ \selectlanguage{english}(it has seen at least one, but not all alleles before).
Put another way, we do not distinguish between loci, assuming that
sharing an allele at one locus is functionally identical to sharing an
allele at any other locus, or indeed all other loci except one.
\subsection{Extensions to more than one
population}
{\label{917367}}
Following~\hyperref[csl:18]{(18)}, we model movement between populations
using a dispersal matrix~\selectlanguage{greek}\textbf{Δ} \selectlanguage{english}=~\textbf{A} -~\textbf{E},
where~\textbf{A} is the weighted adjacency matrix containing elements
A\textsubscript{kl} indicating the proportion of population k (row)
moving to population l (column) per unit time and~\textbf{E} is a
diagonal matrix representing emigration, where each
entry~\(\mathrm{E}_{kk}=\sum_{k=1}^n\mathrm{A}_{kl}\) where n is the number of populations. Thus, the
whole system can be depicted by a set of three equations per strain i
per population k:
\begin{equation}
\label{eqn:multipop}
\begin{aligned}
\ddt{y_{i,l}} &= \beta \left( (1 - w_{i,k}) + (1 - \gamma) (w_{i,k} - z_{i,k}) \right) y_{i,k} - \sigma y_{i,k} - \mu y_{i,k} + \sum_k \Delta_{kl} y_{i,k}\\
\ddt{z_{i,l}} &= \beta (1 - z_{i,k}) y_{i,k} - \mu z_{i,k} + \sum_k \Delta_{kl} z_{i,k}\\
\ddt{w_{i,l}} &= \beta (1 - w_{i,k}) \selectlanguage{ngerman}\sum_{j \ni j \cap i \neq \text{Ø}} y_{j,k} - \mu w_{i,k} + \sum_k \Delta_{kl} w_{i,k}
\end{aligned}
\end{equation}
Where~\selectlanguage{greek}\textbf{Δ}\selectlanguage{english}\textsuperscript{T} signifies the transpose of
\selectlanguage{greek}\textbf{Δ}, \selectlanguage{english}and each equation from
Section~{\ref{826836}} is now additionally indexed
according to population and has an additional term to account for
migration between populations. While in principle the elements
of~\selectlanguage{greek}\textbf{Δ} \selectlanguage{english}can take any value {[}0, 1{]}, signifying a (continuous)
movement of between 0 and 100\% of individuals, for simplicity we use a
constant value \selectlanguage{greek}δ \selectlanguage{english}for the strength of each movement,~\emph{i.e.} for each
non-zero off-diagonal element of \selectlanguage{greek}\textbf{Δ}. \selectlanguage{english}Sensitivity to this value
is explored in the supporting information
(Fig~{\ref{175706}}).
This framework can be applied to a metapopulation of arbitrary size and
complexity, with the number of equations being linearly related to the
number of populations. The dynamics of each population are governed by a
set of three equations per pathogen strain, and these equations are
interlinked within populations by partial, cross-reactive immunity, and
between populations through a movement network. The total number of
differential equations for any given system will be three times the
number of strains multiplied by the number of populations in the
metapopulation.
\subsection{Simulation Procedure}
{\label{898391}}
All simulations were carried out in Julia version
1.3.0~\hyperref[csl:19]{(19)}, with graphics produced using the ggplot2
package~\hyperref[csl:20]{(20)} in R version 3.6.1~\hyperref[csl:21]{(21)}. For
simplicity of presentation, we fix the values of all variables other
than~\selectlanguage{greek}β \selectlanguage{english}(the infection rate) and~\selectlanguage{greek}\textbf{Δ} \selectlanguage{english}(the network of movement
information) to be identical for all populations in the metapopulation.
To demonstrate the variety of dynamics obtainable in this modeling
framework, we vary~\[\tilde{R_0} = \frac{\beta}{\sigma + \mu + E_{kk}},\]where E\textsubscript{kk} =
-\selectlanguage{greek}Δ\selectlanguage{english}\textsubscript{kk}, as noted above, signifies the total outgoing
movement from the population of interest. We add a \textasciitilde{}
over R\textsubscript{0} to denote that this is an approximation of the
true reproductive number, the precise form of which would additionally
take into account the inflow of infectious individuals from other
populations. We additionally vary~\selectlanguage{greek}\textbf{Δ} \selectlanguage{english}according to the number
and~interconnectedness of the populations. For the figures of the main
text, we utilize a strain structure of two loci, each with two alleles.
Sensitivity to these parameter choices is explored in the supporting
information (Fig~{\ref{124825}}).
\subsubsection{Populations with identical
parametrizations}
{\label{454156}}
To assess the effect of migration on population dynamics, we first
consider the simplest case of a set of populations sharing the same
disease parametrization: \selectlanguage{greek}β \selectlanguage{english}such that~\(\tilde{R_0}\) = 2, \selectlanguage{greek}σ \selectlanguage{english}= 8, \selectlanguage{greek}μ \selectlanguage{english}=
0.1, and \selectlanguage{greek}γ \selectlanguage{english}= 0.66. We use a movement network described by a chain of
populations,~\emph{i.e.}
A\(\rightarrow\)B\(\rightarrow\)C\(\rightarrow\)D or
\begin{equation*}
\mathbf{\Delta} = \begin{bmatrix} -\delta & \delta & 0 & 0 \\ 0 & -\delta & \delta & 0 \\ 0 & 0 & -\delta & \delta \\ 0 & 0 & 0 & 0 \end{bmatrix},
\end{equation*}
where \selectlanguage{greek}δ \selectlanguage{english}= 0.05, and ask how the dynamics of populations further down the
chain (\emph{i.e.} B, C, D) differ from those of the origin population
(\emph{i.e.} A), recalling that, without migration, all populations
would have identical dynamics.
\subsubsection{Populations with varying
parametrizations}
{\label{720243}}
We next consider the case where parameters differ between connected
populations, we restrict our consideration to a system of two
populations, identical in all respects other than the parameter \selectlanguage{greek}β, \selectlanguage{english}which
is set to either induce a steady state of coexisting strains (\selectlanguage{greek}β \selectlanguage{english}such
that~\(\tilde{R_0}\) = 5 in population A) or cyclically coexisting
strains (\selectlanguage{greek}β \selectlanguage{english}such that~\(\tilde{R_0}\) = 2 in population B). We then
display three potential patterns of connection: no migration (left
column), A\(\rightarrow\)B (middle column), and B\(\rightarrow\)A
(right column). Specifically, we set
\begin{equation*} \mathbf{\Delta} = \begin{bmatrix} 0 & 0 \\ 0 & 0 \end{bmatrix}, \mathbf{\Delta} = \begin{bmatrix} -\delta & \delta \\ 0 & 0 \end{bmatrix}, \textrm{and~} , \mathbf{\Delta} = \begin{bmatrix} 0 & 0 \\ \delta & -\delta \end{bmatrix},\end{equation*}
respectively, again with~\selectlanguage{greek}δ \selectlanguage{english}= 0.05, \selectlanguage{greek}σ \selectlanguage{english}= 8, \selectlanguage{greek}μ \selectlanguage{english}= 0.1, and \selectlanguage{greek}γ \selectlanguage{english}= 0.66.
To address the case of multiple origin populations feeding into a single
destination population, we consider a system of three populations:
A\(\rightarrow\)C\(\leftarrow\)B, or
\begin{equation*} \mathbf{\Delta} = \begin{bmatrix} -\delta & 0 & \delta \\ 0 & -\delta & \delta \\ 0 & 0 & -\delta \end{bmatrix}, \end{equation*}
where populations A and C have~\selectlanguage{greek}β \selectlanguage{english}corresponding to an~\(\tilde{R_0}\)
= 5, but population B has~\selectlanguage{greek}β \selectlanguage{english}corresponding to an \(\tilde{R_0}\) = 2;
all other parameters as above.
\subsubsection{Larger network structure}
{\label{930618}}
Finally, we characterize the role of global network structure through
considering the impact of degree distribution on a few summary
statistics of system-wide disease burden: the mean proportion of
infectious individuals (area under the currently infectious (\emph{i.e.}
y) curve), the mean level of strain-specific immunity (average z value),
and the mean time between epidemic peaks (\emph{i.e.} between local
maxima in y) over the course of the final 33\% of the simulation. We
omit the initial period of the simulation to reduce the impact of
transient dynamics.
We perform 100 simulations for each of five generic network ensembles
each with 25 populations and a connectedness of approximately 0.15.
Specifically, we examine~Erd\selectlanguage{polish}ő\selectlanguage{english}s-R\selectlanguage{ngerman}ényi (links randomly assigned between
populations), stochastic block (a metapopulation consisting of two
groups of populations which have high migration within the group, but
low migration to populations in the other group), tree-like (where there
are many chains of populations and no potential for
cycles),~Barabasi-Albert (a scale-free network in which there tends to
be a few populations with very many links, and many populations with few
links), and Watts-Strogatz (a small-world network structure which is
produced by partially re-wiring a spatially connected grid of
populations) network structures. To generate these networks, we utilize
functions from the tidygraph R package~\hyperref[csl:22]{(22)}, except in the
case of the tree and Watts-Strogatz configuration for which we use
custom algorithms. Note that we use a parameter of attachment of 4 for
the Barabasi-Albert random graphs. This allows for comparable
connectences to the other random graphs as well as distinguishing these
randomizations from trees (as would result from a default parameter of
attachment of 1). In all cases, each migration strength is set to a
constant~\selectlanguage{greek}δ \selectlanguage{english}= 0.01, only the pattern of connections varies. Each
population is assigned a random \selectlanguage{greek}β \selectlanguage{english}value corresponding to
a~\(\tilde{R_0}\) between {[}1, 6{]}. These results are qualitatively
similar if instead every population is assigned the same value of \selectlanguage{greek}β.
\selectlanguage{english}All code is made available on GitHub: \url{https://git.io/JeqMc}.
\section{Results}
{\label{482316}}
In the following sections, we provide figures to demonstrate the effect
of metapopulation structure on disease dynamics. In these figures, we
plot a time series for each of three subsets of the population: those
currently infected with a particular strain of the pathogen, those
having (complete) specific immunity against the focal strain, and those
who have at least partial cross-reactive immunity to the focal strain,
due to past exposure to a similar strain (see Methods). Populations
differ in their \selectlanguage{greek}β \selectlanguage{english}(and thus R\textsubscript{0}) value. This can be
considered, for example, as differences in population density, which
affects the probability of disease transmission. We only depict one
representative strain in each plot for visual clarity and parametrize
the model such that all strains are functionally equivalent
(\emph{i.e.~}they all have the same transmission and recovery rates
within any given population).
\subsection{Cyclical dynamics are dampened along chains in the
metapopulation
network}
{\label{623830}}
We found that even when all populations share the same parametrizations
and initial conditions, that populations further along network chains
have reduced proportions of currently infectious individuals and
dampened oscillatory dynamics compared to those they would exhibit in
isolation (Fig~{\ref{621939}}). This is due to the
movement of (partially) immune individuals between the populations,
increasing the proportion of individuals with specific and
cross-reactively immunity in populations further along the chain. While
infectious individuals move at an equal rate, the proportion of the
population that is currently infectious at any given time is much
smaller than the proportion with immunity.\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=500]{figures/dampening-cycles/dampening-cycles}
\caption{{Connecting multiple populations with the same parameters results in
reduced pathogen prevalence and dampened cycles in populations further
down the chain. Here, populations are connected such
that~\(A\rightarrow B\rightarrow C\rightarrow D\). Each column indicates a population, while each
row is one of the three population classes laid out above and in the
Methods,~\emph{i.e.} those currently infectious with the given strain,
those with (complete) specific immunity to the given strain, and those
with partial cross-protective immunity to the given strain. The mean
level of immunity (both specific (middle row) and cross-reactive (bottom
row)) increases in each sequential population, while the mean level of
currently infectious individuals (top row) decreases. All populations
have parameters \selectlanguage{greek}σ \selectlanguage{english}= 8,~\selectlanguage{greek}μ \selectlanguage{english}= 0.1,~\selectlanguage{greek}δ \selectlanguage{english}= 0.05, \selectlanguage{greek}γ \selectlanguage{english}= 0.66 and a \selectlanguage{greek}β \selectlanguage{english}chosen to
make \(\tilde{R_0}\) = 2 for all populations. We use a two-loci,
two-allele strain structure, but show only one strain for clarity (but
see supporting information Fig~{\ref{134222}}).
{\label{621939}}%
}}
\end{center}
\end{figure}
\subsection{Dynamics propagate through metapopulation
networks}
{\label{615462}}
We found that in the case of a simple chain of populations, the dynamics
of destination populations can be overridden by the dynamics of origin
populations (Fig~{\ref{438608}}). Interestingly, this
is true both of cyclical dynamics overruling stable dynamics
and~\emph{vice versa}, though the required amount of migration differs
according to the origin and destination dynamics (see supporting
information Fig~{\ref{175706}}).~This migration can
also allow for strain coexistence even in populations where the local
parameters would suggest extinction of one or more strains.\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=375]{figures/two-pop-comparison/two-pop-comparison}
\caption{{Destination populations tend to inherit origin population dynamics when
linking populations with different model parametrizations. As in
Fig~{\ref{621939}}, rows correspond to population
classes, but here, columns indicate network structure. While in
isolation (left column), population A has cyclical dynamics and
population B has steady-state dynamics, when the two populations are
linked by migration, the destination population inherits the dynamics of
the origin population (center and right columns). This is true
regardless of the direction of the movement (depending on the level of
migration; see supporting information
Fig~{\ref{175706}}). Populations A and B have
parameters \selectlanguage{greek}σ \selectlanguage{english}= 8,~\selectlanguage{greek}μ \selectlanguage{english}= 0.1,~\selectlanguage{greek}δ \selectlanguage{english}= 0.05, and \selectlanguage{greek}γ \selectlanguage{english}= 0.66 in common and \selectlanguage{greek}β \selectlanguage{english}chosen
to reflect~\(\tilde{R_0}\) = 2 and 5, respectively. We use a
two-loci, two-allele strain structure, but show only one strain for
clarity.
{\label{438608}}%
}}
\end{center}
\end{figure}
The issue of dynamics propagation gets more complicated when there are
multiple, varying origin populations for a given destination population.
We found that there is a hierarchy of dynamics in their propagation
through the network: when there are origin populations with both
cyclical and steady state dynamics, the destination population inherits
the cyclical dynamics (Fig~{\ref{441914}}), albeit
dampened from what they would have been without migration from a steady
state population. This asymetrical inheritance is robust to imbalance in
the relative contributions of the origins. Put another way, if just one
of many origin populations (or a small proportion of the total movement)
has cyclical dynamics, the destination population will also have
cyclical dynamics.\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=375]{figures/dynamics-hierarchy/dynamics-hierarchy}
\caption{{When multiple origin populations differ in their dynamics, the
destination population inherits cycles over steady states. As in
Fig~{\ref{621939}}, rows indicate population classes,
and columns the component populations. Here, we have populations A and B
feeding into population C at the same rate of~\selectlanguage{greek}δ \selectlanguage{english}= 0.05. Populations A
and C are parametrized to produce steady state dynamics in the absence
of migration, with~\selectlanguage{greek}σ \selectlanguage{english}= 8,~\selectlanguage{greek}μ \selectlanguage{english}= 0.1, \selectlanguage{greek}γ \selectlanguage{english}= 0.66, and \selectlanguage{greek}β \selectlanguage{english}corresponding to
a~\(\tilde{R_0} = 2\). Population B shows cyclical dynamics with~\selectlanguage{greek}β
\selectlanguage{english}corresponding to a~\(\tilde{R_0} = 5\) and all other parameters the same.
Note that, even though the parameters of population C would lead to a
steady state in the absence of migration, we see cyclical dynamics being
inherited from population B. We use a two-loci, two-allele strain
structure, but show only one strain for clarity.
{\label{441914}}%
}}
\end{center}
\end{figure}
\subsection{Degree distribution affects pathogen prevalence and
immunity}
{\label{311338}}
These simple patterns in the effects of origin population dynamics on
those in the destination population have clear implications when pieced
together into larger network structures. For instance, the propagation
of immune individuals through the metapopulation suggests that
populations further ``up the chain'' will tend to have higher on-average
disease incidence and also greater variability. The inheritance of
dynamical regimes combined with a hierarchy of dynamics in that
inheritence suggests that chaos and cycles should be more common,
especially in populations further ``down the chain.'' That is, except in
cases where the ultimate origin populations are all disposed toward
steady states, in which case the stabilizing effect could overrule
downstream local parametrizations, leading to an overall stable system.
In Fig~{\ref{948197}}, we report the effect of various
network structures on three summary statistics of pathogen prevalence
(and levels of immunity) using five common network ensembles. Depending
on the system being explored, empirical network structures might have
elements in common with one or more of these ensembles, for instance,
many social networks are considered to be ``small-world'' in structure
like Watts-Strogatz random graphs, while ecological networks are often
commented on for their formation of ``modules'' or clusters of more
densely interacting species as in stochastic block random graphs.
Networks were parametrized to have approximately equal connectance and
size in order to reduce uninformative variation (see
Section~{\ref{930618}}). This is because metapopulation
size and connectance have known effects on pathogen persistence,
independent of further network structure~\selectlanguage{ngerman}\hyperref[csl:23]{(23}; \hyperref[csl:24]{24}; \hyperref[csl:25]{25)}.~
We found that the network configurations with higher variation in
indegree (\emph{i.e.} the number of other populations each population
receives migration from) distributions (supporting information
Fig~{\ref{403698}}), such as those found in the tree
and Barabasi-Albert networks, tend to have higher levels of infection
over time, despite similar levels of immunity as the other three network
types (Fig {\ref{948197}}). We saw similar patterns to
those in infectious individuals when looking at time between epidemic
peaks across network types. While fewer of the populations ended up with
cyclical dynamics in the Barabasi-Abert graphs, the mean period of the
cycles tended to be slightly higher and have higher variance, but this
was not robust to alternative parametrizations (supporting information
Fig~{\ref{376031}}).\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=500]{figures/large-network-comparison/large-network-comparison}
\caption{{The effect of network structure on pathogen prevalence and levels of
immunity through time. In the top row, we depict a representative
network from each of the five ensembles. The second row shows the
distributions of each of three response variables for prevalence of the
pathogen and specific immunity over the course of the simulation, with
the units for the horizontal axes given by the panel headings. We depict
one point for each randomized network structure and box-plots indicating
the median and inter-quartile range of each network-type's distribution.
Network generating algorithms were tuned to produce networks of the same
size and approximate connectance and model parameters were either the
same for all populations and across simulations (\selectlanguage{greek}σ \selectlanguage{english}= 8, \selectlanguage{greek}μ \selectlanguage{english}= 0.1, \selectlanguage{greek}δ \selectlanguage{english}=
0.01, and \selectlanguage{greek}γ \selectlanguage{english}= 0.66) or randomized for each population in each simulation
(initial densities of infectious and immune individuals {[}0, 1{]} and~\selectlanguage{greek}β
\selectlanguage{english}value corresponding to a~\(\tilde{R_0}\) within {[}1, 6{]} for each
population).
{\label{948197}}%
}}
\end{center}
\end{figure}
\section{Discussion}
{\label{945890}}
Both metapopulation~\hyperref[csl:26]{(26)} and strain~\hyperref[csl:16]{(16)}
structure have long been known to be important to disease dynamics and
are increasingly being recognized as ubiquitous. Yet the combination of
these two areas of theory has been underexplored. We show here that this
lacuna can have real consequences for our understanding of disease
dynamics in empirical systems.
In probing the relationship between origin and destination dynamics in
simple metapopulations, we have demonstrated several patterns that
expand our understanding of disease dynamics in these systems. By
directly incorporating a movement network into our model framework, we
have constructed a very general approach that lends itself to
arbitrarily large and complex systems. This is noteworthy, as more and
more natural systems are being thought of in terms of networks of
interacting components (\emph{e.g.} separate species in ecological
communities~\hyperref[csl:27]{(27)} or host individuals exchanging
parasites~\hyperref[csl:28]{(28)}). By adjusting the scale of our
metapopulation, we can ask and answer different questions about the
forces influencing disease dynamics.
We found that the dynamics of prevalence and immunity among
migrationally connected populations are not independent, and that even
very small rates of population movement can have profound effects on a
population's disease dynamics: from reducing pathogen prevalence to
changing the dynamical regime of destination populations entirely. Our
findings regarding the reduction in cycle amplitude
(Section~{\ref{623830}}) echo results in dispersal
networks in ecology, where population dynamics were dampened following
the introduction of migration~\hyperref[csl:29]{(29)}.~
Contrary to prior focus in the literature on the role of migrating
infectious individuals~\hyperref[csl:30]{(30}; \hyperref[csl:31]{31}; \hyperref[csl:26]{26)}, we found that the migration
of immune individuals can be equally (or even more) important. This is
noteworthy, as the few previous studies relating multi-strain diseases
and metapopulation structure only allow pathogen transmission between
populations, not the movement of individuals
explicitly~\selectlanguage{ngerman}\hyperref[csl:3]{(3}; \hyperref[csl:32]{32)}---an approach that is more mathematically
tractable, but omits the potentially influential transmission of immune
individuals.
Finally, we show that larger network structure also has a part to play
in disease dynamics, resulting in significant differences in pathogen
prevalence across network types. Our results are in agreement with
previous results suggesting increased epidemic size in scale-free
network structures (such as those found in Barabasi-Albert random
graphs) when the spreading rate is sufficiently slow \hyperref[csl:33]{(33}; \hyperref[csl:34]{34)}
due to the high-degree nodes serving as
``super-spreaders''~\hyperref[csl:35]{(35}; \hyperref[csl:25]{25)}. Along these lines, there has
been some previous research indicating that node degree (the number of
other populations a given population is connected to) is directly
related to pathogen prevalence in that focal population
(\hyperref[csl:36]{(36)}, but see~\hyperref[csl:37]{(37)}), however a complete
investigation into network structure at the node-level is beyond the
scope of this work. A comprehensive investigation of the role of more
complex network structures in disease dynamics, however, remains a topic
for further investigation.
In this work, we have utilized a relatively simple model for disease
dynamics in an effort to maximize interpretability. Such simplification
inevitably comes with a cost, and several of our assumptions can be
critiqued as unrealistic. Perhaps foremost is the assumption of
continuous movement. While continuous movement might be appropriate for
very large populations with frequent, relatively small migrations
between them, when any of these three components is not present, we
would expect deviation from these predictions. Future work could explore
the importance of discrete movement regimes on these patterns. Likewise,
in this work we omit strain mutation and
recombination~\hyperref[csl:17]{(17)} (yet the latter is included in the
original framework of~\hyperref[csl:16]{(16)}). The generation of novel
strains is likely important to the global persistence of diseases in
humans~\hyperref[csl:38]{(38)} and animals~~\hyperref[csl:39]{(39)}. Finally, in
representing movement by adding a proportion of the origin population to
the destination, we introduce an assumption that the two populations in
question are of approximately the same size. In a metapopulation where
populations vary widely in size, the proportion leaving one population
would not correspond to the proportion entering another.
This work should not be seen as an attempt at comprehensive
categorization of the role of meta-population structure on the dynamics
of multi-strain diseases, but rather as an initial step in exploring the
complex interplay between the population structure of hosts and strain
structure of pathogens. Our results suggest there may be simple rules
underlying this relationship, at least for a wide range of parameter
values, but it remains to be seen where networks based on empirical data
fall in these parameter regimes, as well as how such systems might
deviate from theoretical expectations.~
\section{Acknowledgments}
{\label{519417}}
We would like to thank José Lourenço for helpful discussions regarding
the mathematical framework of this work. This work was supported by the
USDA National Institute of Food and Agriculture, Animal Health project
\#MINV-62-057.
\beginsupplement
\section{Supporting Information}
{\label{586469}}\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=500]{figures/dampening-cycles-all-strains/dampening-cycles-all-strains}
\caption{{Considering all the dynamics of all four strains from population A in
Fig {\ref{621939}}. Note that lines are colored
according to strain rather than population. Strains can be divided into
two discordant sets of non-overlapping alleles: \{1, 1\} and \{2, 2\},
and \{1, 2\} and \{2, 1\}. Each strain of a discordant set~ behaves
identically due to identical parametrization and no interaction between
strains that do not share at least one allele, but discordant sets
interact with one another due to partial cross-reactive immunity. Thus,
when one set is abundant, the other is rare and~\emph{vice versa}. We
highlight the maximum value of each discordant set's cycle with a
vertical line in order to facilitate comparisons between strains and
sets.
{\label{134222}}%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=500]{figures/variable-movement-rate/variable-movement-rate}
\caption{{The effect of variable migration rate on transference of dynamical
regime from origin to destination in a simple metapopulation of two
populations linked by unidirectional movement. In the top panel are
cases where the origin population has cyclical dynamics and the
destination has steady state dynamics when in isolation, while in the
bottom panel the opposite is true. As the movement rate (\selectlanguage{greek}δ; \selectlanguage{english}horizontal
axis) increases, there is a phase-transition at which point the
destination population's dynamics (indicated by the vertical axis)
switch to match those of the origin. We fit a binomial spline to
highlight this transition point.~We see that even with very small rates
of migration, a stable population can be converted to a cyclical one
(top panel). Yet, it is more difficult to convert a cyclical population
to one with steady state dynamics (bottom panel). Three parametrizations
are recorded here (color; see legend), with additional parameter~\selectlanguage{greek}μ \selectlanguage{english}=
0.15 being the same for both populations. Finally, note that the two
values of \(\tilde{R_0}\) listed in the legend correspond to the two
populations, with the larger value corresponding to the steady state
population and the smaller value the cyclical population (when in
isolation) .
{\label{175706}}%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=500]{figures/state-space-exploration/state-space-exploration}
\caption{{Effect of parametrization on dynamic regime for a population in
isolation. Here, columns indicate values of \selectlanguage{greek}γ \selectlanguage{english}and rows indicate values
of \selectlanguage{greek}σ. \selectlanguage{english}Depending on the combination of \selectlanguage{greek}β, σ, μ, \selectlanguage{english}and \selectlanguage{greek}γ, \selectlanguage{english}a population can
exhibit a range of dynamics including steady states for all strains
(pink), cyclical or chaotic dynamics for all strains (orange), or
partial extinction of some strains (blue). Green cells indicate a
numerical failure in integration. All simulations here utilize a
two-loci, two-allele strain structure. See \protect\hyperref[csl:16]{(16)} for a
similar figure for alternative parameterizations.
{\label{124825}}%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=500]{figures/large-network-degree-distributions/large-network-degree-distributions}
\caption{\selectlanguage{polish}{Summary statistics for the degree distributions of each randomized
network used for Fig~{\ref{948197}} in the main text.
Networks were constructed to have the same size and approximate
connectance, but with the network structure (which populations are
connected to which other popuations) otherwise generated according to
one of five algorithms: Erdős-Rényi, Barabasi-Albert, and
Watts-Strogatz, stochastic block, and tree (see
Section~{\ref{930618}} of the main text). Some
algorithms allowed perfect matching of connectance (Erdős-Rényi,
Barabasi-Albert, and Watts-Strogatz), while others necessitated some
minor variation (stochastic block and tree).
{\label{403698}}%
}\selectlanguage{ngerman}}
\end{center}
\end{figure}\selectlanguage{ngerman}\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=500]{figures/large-network-comparison-sigma=40/large-network-comparison-sigma=20}
\caption{{Similar to the lower row of Fig~{\ref{948197}},but with
\selectlanguage{greek}γ \selectlanguage{english}and~\selectlanguage{greek}σ \selectlanguage{english}equal to 0.55 and 32, respectively. All other parameters are
equal to or set randomly as in Fig~{\ref{948197}}. All
other parameters are equal to or set randomly as in
Fig~{\ref{948197}}. While the observed differences in
total infected are robust, note that here the mean time between epidemic
peaks is approximately equal across randomizations.
{\label{376031}}%
}}
\end{center}
\end{figure}
\selectlanguage{english}
\FloatBarrier
\section*{References}\sloppy
\phantomsection
\label{csl:1}1. Couch RB, Kasel JA. {Immunity to influenza in man}. Annual review of microbiology. 1983;37(1):529–49.
\phantomsection
\label{csl:2}2. Castillo-Chavez C, Hethcote HW, Andreasen V, Levin SA, Liu WM. {Epidemiological models with age structure, proportionate mixing, and cross-immunity.}. J Math Biol. 1989;27:233–58.
\phantomsection
\label{csl:3}3. Lourenço J, Recker M. {Natural, persistent oscillations in a spatial multi-strain disease system with application to dengue.}. PLoS Comput Biol. 2013;9:e1003308.
\phantomsection
\label{csl:4}4. Wilson EO, MacArthur RH. {The theory of island biogeography}. Princeton University Press; 1967.
\phantomsection
\label{csl:5}5. Taylor AD. {Large-scale spatial structure and population dynamics in arthropod predator-prey systems}. Annales Zoologici Fennici. 1988;:63–74.
\phantomsection
\label{csl:6}6. Hanski I. {A practical model of metapopulation dynamics}. The Journal of Animal Ecology. January 1994;63(1):151.
\phantomsection
\label{csl:7}7. Hess G. {Disease in metapopulation models: implications for conservation}. Ecology. July 1996;77(5):1617–32.
\phantomsection
\label{csl:8}8. Cross PC, Lloyd-Smith JO, Johnson PLF, Getz WM. {Duelling timescales of host movement and disease recovery determine invasion of disease in structured populations}. Ecology Letters. June 2005;8(6):587–95.
\phantomsection
\label{csl:9}9. Lloyd AL, Jansen VAA. {Spatiotemporal dynamics of epidemics: synchrony in metapopulation models}. Mathematical Biosciences. March 2004;188(1-2):1–16.
\phantomsection
\label{csl:10}10. Ruxton GD. {Low levels of immigration between chaotic populations can reduce system extinctions by inducing asynchronous regular cycles}. Proceedings of the Royal Society of London Series B: Biological Sciences. May 1994;256(1346):189–93.
\phantomsection
\label{csl:11}11. Earn DJD, Rohani P, Grenfell BT. {Persistence chaos and synchrony in ecology and epidemiology}. Proceedings of the Royal Society of London Series B: Biological Sciences. January 1998;265(1390):7–10.
\phantomsection
\label{csl:12}12. Hanski I. {Metapopulation dynamics}. In: Metapopulation Biology. Elsevier; 1997. p. 69–91.
\phantomsection
\label{csl:13}13. Brown JH, Kodric-Brown A. {Turnover Rates in Insular Biogeography: Effect of Immigration on Extinction}. Ecology. March 1977;58(2):445–9.
\phantomsection
\label{csl:14}14. Rosenzweig ML. {Paradox of Enrichment: Destabilization of Exploitation Ecosystems in Ecological Time}. Science. January 1971;171(3969):385–7.
\phantomsection
\label{csl:15}15. Hilker FM, Schmitz K. {Disease-induced stabilization of predator{\textendash}prey oscillations}. Journal of Theoretical Biology. December 2008;255(3):299–306.
\phantomsection
\label{csl:16}16. Gupta S, Ferguson N, Anderson R. {Chaos persistence, and evolution of strain structure in antigenically diverse infectious agents}. Science. May 1998;280(5365):912–5.
\phantomsection
\label{csl:17}17. Louren{\c{c}}o J, Wikramaratna PS, Gupta S. {{MANTIS}: an R package that simulates multilocus models of pathogen evolution}. {BMC} Bioinformatics. May 2015;16(1).
\phantomsection
\label{csl:18}18. Xiao Y, Zhou Y, Tang S. {Modelling disease spread in dispersal networks at two levels}. Mathematical Medicine and Biology. September 2011;28(3):227–44.
\phantomsection
\label{csl:19}19. Bezanson J, Edelman A, Karpinski S, Shah VB. {Julia: a fresh approach to numerical computing}. {SIAM} Review. January 2017;59(1):65–98.
\phantomsection
\label{csl:20}20. Wickham H. {ggplot2: Elegant Graphics for Data Analysis}. Springer-Verlag New York; 2016.
\phantomsection
\label{csl:21}21. {R Core Team}. {R: a language and environment for statistical computing}. Vienna, Austria: R Foundation for Statistical Computing; 2019.
\phantomsection
\label{csl:22}22. Pedersen TL. {tidygraph: a tidy {API} for graph manipulation}. 2019.
\phantomsection
\label{csl:23}23. Ganesh A, Massouli{\'e} L, Towsley D. {The effect of network topology on the spread of epidemics}. In: Proceedings IEEE 24th Annual Joint Conference of the IEEE Computer and Communications Societies. IEEE; 2005. p. 1455–66.
\phantomsection
\label{csl:24}24. Salathé M, Jones JH. {Dynamics and control of diseases in networks with community structure.}. PLoS Comput Biol. 2010;6:e1000736.
\phantomsection
\label{csl:25}25. Keeling MJ, Eames KTD. {Networks and epidemic models}. Journal of The Royal Society Interface. June 2005;2(4):295–307.
\phantomsection
\label{csl:26}26. Grenfell B, Harwood J. {(Meta)population dynamics of infectious diseases}. Trends in Ecology {\&} Evolution. October 1997;12(10):395–9.
\phantomsection
\label{csl:27}27. Proulx SR, Promislow DE, Phillips PC. {Network thinking in ecology and evolution.}. Trends Ecol Evol. 2005;20:345–53.
\phantomsection
\label{csl:28}28. Craft ME, Caillaud D. {Network models: an underutilized tool in wildlife epidemiology?}. Interdiscip Perspect Infect Dis. 2011;2011:676949.
\phantomsection
\label{csl:29}29. Holland MD, Hastings A. {Strong effect of dispersal network structure on ecological dynamics.}. Nature. 2008;456:792–4.
\phantomsection
\label{csl:30}30. Jesse M, Ezanno P, Davis S, Heesterbeek JAP. {A fully coupled mechanistic model for infectious disease dynamics in a metapopulation: Movement and epidemic duration}. Journal of Theoretical Biology. September 2008;254(2):331–8.
\phantomsection
\label{csl:31}31. North AR, Godfray HCJ. {The dynamics of disease in a metapopulation: The role of dispersal range}. Journal of Theoretical Biology. April 2017;418:57–65.
\phantomsection
\label{csl:32}32. Wikramaratna PS, Pybus OG, Gupta S. {Contact between bird species of different lifespans can promote the emergence of highly pathogenic avian influenza strains.}. Proc Natl Acad Sci U S A. 2014;111:10767–72.
\phantomsection
\label{csl:33}33. Pastor-Satorras R, Vespignani A. {Epidemics and immunization in scale-free networks}. In: Handbook of Graphs and Networks [Internet]. Wiley-{VCH} Verlag {GmbH} {\&} Co. {KGaA}; 2004. p. 111–30. Available at: \url{https://doi.org/10.1002\%2F3527602755.ch5}
\phantomsection
\label{csl:34}34. Lloyd-Smith JO, Schreiber SJ, Kopp PE, Getz WM. {Superspreading and the effect of individual variation on disease emergence.}. Nature. 2005;438:355–9.
\phantomsection
\label{csl:35}35. Shirley MDF, Rushton SP. {The impacts of network topology on disease spread}. Ecological Complexity. September 2005;2(3):287–99.
\phantomsection
\label{csl:36}36. Godfrey SS, Bull CM, James R, Murray K. {Network structure and parasite transmission in a group living lizard the gidgee skink, Egernia stokesii}. Behavioral Ecology and Sociobiology. April 2009;63(7):1045–56.
\phantomsection
\label{csl:37}37. VanderWaal KL, Atwill ER, Hooper S, Buckle K, McCowan B. {Network structure and prevalence of Cryptosporidium in Belding's ground squirrels}. Behavioral Ecology and Sociobiology. July 2013;67(12):1951–9.
\phantomsection
\label{csl:38}38. Antia R, Regoes RR, Koella JC, Bergstrom CT. {The role of evolution in the emergence of infectious diseases}. Nature. December 2003;426(6967):658–61.
\phantomsection
\label{csl:39}39. Tompkins DM, Carver S, Jones ME, Krkošek M, Skerratt LF. {Emerging infectious diseases of wildlife: a critical perspective.}. Trends Parasitol. 2015;31:149–59.
\end{document}