# Disease Persistence in Epidemiological Models:
The Interplay between Vaccination and Migration

Jackson Burton  Lora Billings<EMAIL_ADDRESS>Derek A. T.
Cummings  Ira B. Schwartz  Montclair State University, Department of
Mathematical Sciences, Montclair, NJ 07043  Johns Hopkins Bloomberg School of
Public Health, Department of International Health, Baltimore, MD, 21205  US
Naval Research Laboratory, Code 6792, Nonlinear System Dynamics Section,
Plasma Physics Division, Washington, DC 20375

###### Abstract

We consider the interplay of vaccination and migration rates on disease
persistence in epidemiological systems. We show that short-term and long-term
migration can inhibit disease persistence. As a result, we show how migration
changes how vaccination rates should be chosen to maintain herd immunity. In a
system of coupled SIR models, we analyze how disease eradication depends
explicitly on vaccine distribution and migration connectivity. The analysis
suggests potentially novel vaccination policies that underscore the importance
of optimal placement of finite resources.

###### keywords:

epidemics, migration, vaccination, herd immunity

††journal: Mathematical Biosciences

##  1 Introduction

Countries are increasingly connected by travel and economics. Due to economic
disparities and political turmoil, extreme heterogeneity exists in childhood
vaccination coverage across the two sides of multiple national boundaries. It
has been suggested that the immunization coverage of neighboring countries or
those countries well connected by travel can or should be used when crafting
national level immunization policy. In the case of hepatitis B, Gay and
Edmunds [1] argue that it would be four times more cost effective for the
United Kingdom to sponsor a vaccination program in Bangladesh than to
introduce its own universal program. When indigenous wild poliovirus was
eradicated in all but four endemic countries in 2005: India, Nigeria, Pakistan
and Afghanistan, it was exported from northern Nigeria and northern India and
subsequently caused $>50$ outbreaks and paralyzed $>1500$ children in
previously polio-free countries across Asia and Africa [2]. And in 2007, the
WHO estimated that there were 197,000 measles deaths, despite the 82%
worldwide vaccination coverage. In countries where measles has been largely
eliminated, cases imported from other countries remain an important source of
infection [3]. It is clear that a country needs to be concerned with the
vaccination rate of a neighboring country as well as its own.

On another scale, vaccination policies must also take into consideration the
subpopulation dynamics within a country. Wilson, et al. [4] models linked
urban and rural epidemics of HIV and discusses how to optimize a limited
treatment supply to minimize new infections. Cummings et al. [5] uses data to
identify a distinct pattern in the periodicity of measles outbreaks in
Cameroon before the widespread vaccination efforts of the Measles Initiative.
The southern part of Cameroon experienced a significant measles epidemic
approximately every three years. In contrast, the three northern provinces
contend with annual measles epidemics. In 2000 and 2001, these cyclic
outbreaks coincided, exacerbating the situation and causing a much more severe
epidemic [5].

Noting that a small contribution of infections from one population to another
could drive a new type of epidemic that would not normally occur, we study how
migration between populations could change dynamics and respective herd
immunity levels in metapopulation models. We analyze a model of a disease
imported between subpopulations of a region by short-term and long-term
migration with limited vaccination coverage. Our initial study is based on the
analysis of a system of canonical SIR compartmental models. The system allows
the rigorous proof of the qualitative affects migration has on herd immunity.
The model can be enhanced to include more compartments or seasonal forcing,
but most of these systems will require numerical exploration of trends in
spatial synchrony and bifurcation analysis, which will be explored in future
papers. In this article, we revisit the fundamental ways migration is modeled
in metapopulation models and how it fundamentally affects herd immunity.

Migration is often treated as a phenomenological input to maintain incidence
in a population that might experience local fade-out [6]. Long-term migration
has been analyzed by Liebovitch and Schwartz [7], with a thorough derivation
of the linear flux term coupling the patches. This approach also agrees with
the classes of models proposed by Sattenspiel and Dietz [8] and Lloyd and
Jansen [9]. Keeling and Rohani [10] investigated the spatial coupling of
dynamics exhibited in models using multiple formulations of migration
including mass-action coupling and linear flux terms. However, they did not
explore the impact of coupling in the presence of vaccination needed to
maintain disease free states. Additional analysis of mixed long-term and
short-term migration in transport-related disease spread can be found in [11,
12, 13]. These papers derive the global asymptotic stability of the disease
free state for a new disease. Because there is no vaccination, the papers
conclude that it is essential to strengthen restrictions of passenger travel
as soon as the infectious diseases appear.

Our paper considers how migration directly affects the vaccination levels
needed for herd immunity against a known disease and how that would impact
optimum usage of limited vaccination supplies. We investigate the dynamics of
models that include mass-action coupling, an assumption that assumes mixing
occurs at fast time scales, and linear migration, which is more consistent
with mixing occurring at long time scales. The organization of this paper is
as follows: We introduce a coupled compartmental model in Section 2 and
perform stability analysis of the disease free state as a function of the
migration and vaccination rates. We also consider normal forms of the
bifurcations created by the short-term and long-term migration dynamics.
Section 3 describes how vaccination rates should be adjusted with respect to
short-term and long-term migration levels to preserve herd immunity. Section 4
has a summary of our observations and conclusions.

##  2 The Model

We start with the classic Susceptible, Infected, Recovered (SIR) model. Let
$S$, $I$, and $R$ denote the number of people in each of the disease classes
for a population of size $N$. Let the parameters $\beta>0$ denote the contact
rate, $\mu>0$ denote the birth/death rate, and $\kappa>0$ denote the recovery
rate. The vaccination rate, $0\leq v\leq 1$, represents the removal of a
percentage of the incoming newborn population to recovered. The standard form
for this system is

$\displaystyle\frac{dS}{dt}$ $\displaystyle=$ $\displaystyle(1-v)\,\mu
N-\frac{\beta SI}{N}-\mu S,$ $\displaystyle\frac{dI}{dt}$ $\displaystyle=$
$\displaystyle\frac{\beta SI}{N}-\kappa I-\mu I,$ (1)
$\displaystyle\frac{dR}{dt}$ $\displaystyle=$ $\displaystyle v\mu N+\kappa
I-\mu R.$

The death rates in the classes balance the births so that the population size
$N>0$ is constant. For a detailed analysis of the single patch formulation of
this system, see Hethcote [14].

We now consider two coupled subpopulations where the disease dynamics of each
population are described by the SIR model. Let $S_{k}$, $I_{k}$, and $R_{k}$
denote the number of people in each of the disease classes, $\mu_{k}>0$ denote
the birth/death rates, and $v_{k}$ denote the vaccination rates of
subpopulations $N_{k}$ for $k=1,2$. To model long-term movement (linear
mixing), let $c_{1}\geq 0$ denote the rate of migration from population two to
population one and vice versa for the rate $c_{2}\geq 0$. To model short-term
movement (mass action mixing), let $0\leq c_{3}\leq 1$ be a scaling of the
number of infectives from one population who move into the other population
for a short time and mix with the susceptibles to produce additional
infections. Because $\beta$ is proportional to the average number of contacts
a person can make per unit time, we distribute the contacts for the
susceptibles between the infected people by mass action within and outside the
population by using the prefactors $c_{3}$ and $(1-c_{3})$ respectively as in
Keeling and Rohani [10]. The coupled two population model is as follows:

$\displaystyle\frac{dS_{1}}{dt}$ $\displaystyle=$
$\displaystyle(1-v_{1})\,\mu_{1}N_{1}-\frac{(1-c_{3})\,\beta\,S_{1}I_{1}}{N_{1}}-\frac{c_{3}\,\beta\,S_{1}I_{2}}{N_{1}}-\mu_{1}S_{1}+c_{1}S_{2}-c_{2}S_{1},$
$\displaystyle\frac{dI_{1}}{dt}$ $\displaystyle=$
$\displaystyle\frac{(1-c_{3})\,\beta\,S_{1}I_{1}}{N_{1}}+\frac{c_{3}\,\beta\,S_{1}I_{2}}{N_{1}}-\kappa
I_{1}-\mu_{1}I_{1}+c_{1}I_{2}-c_{2}I_{1},$ $\displaystyle\frac{dR_{1}}{dt}$
$\displaystyle=$ $\displaystyle v_{1}\mu_{1}N_{1}+\kappa
I_{1}-\mu_{1}R_{1}+c_{1}R_{2}-c_{2}R_{1},$ (2)
$\displaystyle\frac{dS_{2}}{dt}$ $\displaystyle=$
$\displaystyle(1-v_{2})\,\mu_{2}N_{2}-\frac{(1-c_{3})\,\beta\,S_{2}I_{2}}{N_{2}}-\frac{c_{3}\,\beta\,S_{2}I_{1}}{N_{2}}-\mu_{2}S_{2}+c_{2}S_{1}-c_{1}S_{2},$
$\displaystyle\frac{dI_{2}}{dt}$ $\displaystyle=$
$\displaystyle\frac{(1-c_{3})\,\beta\,S_{2}I_{2}}{N_{2}}+\frac{c_{3}\,\beta\,S_{2}I_{1}}{N_{2}}-\kappa
I_{2}-\mu_{2}I_{2}+c_{2}I_{1}-c_{1}I_{2},$ $\displaystyle\frac{dR_{2}}{dt}$
$\displaystyle=$ $\displaystyle v_{2}\mu_{2}N_{2}+\kappa
I_{2}-\mu_{2}R_{2}+c_{2}R_{1}-c_{1}R_{2}.$

We keep the number of people in the subpopulations constant by letting
$\rho=N_{2}/N_{1}$ and setting the constraint $c_{2}=c_{1}\rho$. This system
is overdetermined by the subpopulation constraints, $S_{k}+I_{k}+R_{k}=N_{k}$
for $k=1,2$, and therefore the analysis omits the variables $R_{k}$ for
$k=1,2$.

Motivated by the distinct subpopulation dynamics of Cameroon described in
Cummings et al. [5], numerical simulations will use parameters based on
Cameroon demographics. The values are listed in Table 1. The subpopulation
sizes are totals for the northern and southern regions based on data in [5].
The birth/death rates are averages over the northern and southern regions
based on data in [5]. The recovery rate is a parameter that is derived from
the biological characteristics of measles. The contact rate was estimated
using the average age of incident measles cases over the period 1998-2006 from
passive surveillance data [5]. The specific results here are fairly
insensitive to changes to $\beta$. An SIR model is used here without the
exposed class but we expect the inclusion of an exposed class would not
substantively change our qualitative results.

Table 1: Parameter Values for Model Based on Cameroon Data Parameter | Value | Unit | Description
---|---|---|---
$N_{1}$ | 4,451,000 | people | Northern subpopulation size
$N_{2}$ | 10,212,000 | people | Southern subpopulation size
$\rho$ | $2.2943$ | none | Ratio of $N_{2}/N_{1}$
$\beta$ | $700$ | year-1 | Contact rate
$\kappa$ | $100$ | year-1 | Measles recovery rate
$\mu_{1}$ | $.0428$ | year-1 | Birth and death rate for $N_{1}$
$\mu_{2}$ | $.0329$ | year-1 | Birth and death rate for $N_{2}$

###  2.1 General system analysis

We start with a general analysis of the system to determine the conditions
necessary for the populations to be disease free.

######  Proposition 1.

System (2) has a disease free equilibrium (DFE) and is given by

$(S_{1},I_{1},S_{2},I_{2})=(N_{1}\hat{S}_{1},0,N_{2}\hat{S}_{2},0),$ (3)

for

$\displaystyle\hat{S}_{1}=\frac{(1-v_{1})(\mu_{1}c_{1}+\mu_{1}\mu_{2})+(1-v_{2})\mu_{2}c_{1}\rho}{\mu_{1}c_{1}+\mu_{1}\mu_{2}+\mu_{2}c_{1}\rho},$
(4)
$\displaystyle\hat{S}_{2}=\frac{(1-v_{1})\mu_{1}c_{1}+(1-v_{2})(\mu_{1}\mu_{2}+\mu_{2}c_{1}\rho)}{\mu_{1}c_{1}+\mu_{1}\mu_{2}+\mu_{2}c_{1}\rho}.$
(5)

Note that the DFE does not depend on the short-term migration parameter
$c_{3}$. Without long-term migration ($c_{1}=0$), the DFE simplifies to
$(S_{1},I_{1},S_{2},I_{2})=\left(N_{1}(1-v_{1}),0,N_{2}(1-v_{2}),0\right)$,
which is the steady state for the uncoupled system. Also, there is no steady
state for which the disease dies out in only one of the two subpopulations if
$c_{1}>0$.

The local stability of the DFE can be determined by the eigenvalues of the
Jacobian of the system evaluated at the DFE. The resulting eigenvalues are

$\displaystyle\lambda_{1}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\left(c_{1}+c_{1}\rho+\mu_{1}+\mu_{2}+\sqrt{(c_{1}+c_{1}\rho+\mu_{1}-\mu_{2})^{2}-4c_{1}(\mu_{1}-\mu_{2})}\right),$
(6) $\displaystyle\lambda_{2}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\left(c_{1}+c_{1}\rho+\mu_{1}+\mu_{2}-\sqrt{(c_{1}+c_{1}\rho+\mu_{1}-\mu_{2})^{2}-4c_{1}(\mu_{1}-\mu_{2})}\right),$
(7) $\displaystyle\lambda_{3}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\left(c_{1}+c_{1}\rho+\mu_{1}+\mu_{2}+2\kappa-(1-c_{3})\,\beta\,(\hat{S}_{1}+\hat{S}_{2})\,+\sqrt{W}\right),$
(8) $\displaystyle\lambda_{4}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\left(c_{1}+c_{1}\rho+\mu_{1}+\mu_{2}+2\kappa-(1-c_{3})\,\beta\,(\hat{S}_{1}+\hat{S}_{2})\,-\sqrt{W}\right),$
(9)

for

$W=4\,\beta\,c_{3}\,c_{1}\left(\hat{S}_{2}+\rho\,\hat{S}_{1}\right)+4\left(\beta^{2}c_{3}^{2}\hat{S}_{1}\hat{S}_{2}+c_{1}^{2}\rho\right)+\left((1-c_{3})\,\beta\left(\hat{S}_{1}-\hat{S}_{2}\right)+c_{1}-c_{1}\rho-\mu_{1}+\mu_{2}\right)^{2}.$
(10)

The DFE is locally stable if the maximum value of the real parts of this set
of eigenvalues is negative.

######  Proposition 2.

The eigenvalue $\lambda_{4}$ determines the local stability of the DFE.

###### Proof.

We can show $\lambda_{1}$ and $\lambda_{2}$ are always negative. First, let

$\begin{array}[]{rl}\theta_{1}&=c_{1}+c_{1}\rho+\mu_{1}+\mu_{2},\\\\[7.22743pt]
\theta_{2}&=(c_{1}+c_{1}\rho+\mu_{1}-\mu_{2})^{2}-4c_{1}(\mu_{1}-\mu_{2}).\\\
\end{array}$ (11)

If we consider $\theta_{2}$ as a quadratic expression in $c_{1}$ with leading
coefficient $(\rho+1)^{2}$, it attains an absolute minimum at
$d\theta_{2}/dc_{1}=0$. Solving this equation gives
$c_{1}=(\mu_{1}-\mu_{2})/(\rho+1)^{2}$. Substituting this expression into
$\theta_{2}$ to find the absolute minimum gives
$4\rho(\mu_{1}-\mu_{2})^{2}/(\rho+1)^{2}>0$. Therefore $\theta_{2}>0$ for all
$c_{1}$, which implies $\lambda_{1}$ and $\lambda_{2}$ are real valued.

Upon inspection we see that $\theta_{1}>0$, $\theta_{1}+\sqrt{\theta_{2}}>0$,
and therefore $\lambda_{1}<0$. For $\lambda_{2}<0$, we require
$\theta_{1}>\sqrt{\theta_{2}}$. This is equivalent to
$c_{1}>-\mu_{1}\mu_{2}/(\mu_{1}+\mu_{2}\rho)$. Since we assume $c_{1}>0$, this
is always true. Therefore $\theta_{1}-\sqrt{\theta_{2}}>0$, which implies
$\lambda_{2}<0$.

For our parameter assumptions, we see that $W>0$ by inspection. This implies
$\lambda_{3}$ and $\lambda_{4}$ are real valued and $\lambda_{4}>\lambda_{3}$.
The only way for the DFE to be unstable is for $\lambda_{3}>0$ or
$\lambda_{4}>0$. Because of the ordering, a sign change would have to happen
for $\lambda_{4}$ first. Therefore, $\lambda_{4}$ determines the stability of
the DFE. ∎

Figure 1:  Contour plot of $\lambda_{4}$ values as a function of $c_{1}$ and
$c_{3}$. Parameters are given by the values in Table 1, with $v_{1}=0.82$ and
$v_{2}=0.90$. The DFE is unstable for $\lambda_{4}>0$, which occurs for
smaller values of $c_{1}$ and $c_{3}$.

To quantify how each migration type effects the stability of the DFE, we can
monitor the sign of $\lambda_{4}$ as we vary $c_{1}$ and $c_{3}$. As an
example, we show a contour plot of $\lambda_{4}$ in Fig. 1 using the
parameters in Table 1, with $v_{1}=0.82$ and $v_{2}=0.90$. When
$\lambda_{4}>0$, the DFE is unstable. From the figure, you can see that as
$c_{1}$ and $c_{3}$ decrease, $\lambda_{4}$ increases and the DFE becomes
unstable. We now explore the underlying conditions necessary in each
subpopulation for which migration can cause the die out or invasion of a
disease.

In the absence of migration, we recover the basic reproductive numbers scaled
by vaccination for each subpopulation for the uncoupled system [14].
Specifically when $c_{1}=c_{3}=0$,

$\hat{R}_{1}(v_{1})=\frac{\beta(1-v_{1})}{\kappa+\mu_{1}}~{}~{}\mbox{and}~{}~{}\hat{R}_{2}(v_{2})=\frac{\beta(1-v_{2})}{\kappa+\mu_{2}}$
(12)

We omit the arguments for $\hat{R_{k}}$ for $k=1,2$ from here on, unless
otherwise specified. The basic reproductive number is the quantity that
defines the threshold between disease absence and persistence, and for the
canonical SIR model without vaccination $R_{0}=\frac{\beta}{\kappa+\mu}$. We
write these expressions as a function of the vaccination rate in the
subpopulation noting that the inequalities $\hat{R}_{k}<1$ for $k=1,2$ implies
that the DFE in each subpopulation is locally stable. At $\hat{R}_{k}=1$ for
$k=1,2$, these two expressions also represent transcritical bifurcations for
the uncoupled system. As an example, for the parameters used in Fig. 1,
$\hat{R}_{1}(0.82)>1$ and the disease would be endemic in $N_{1}$. Conversely,
$\hat{R}_{2}(0.90)<1$ and the disease would die out in $N_{2}$. This motivates
us to examine the effect migration has on a simple system with a transcritical
bifurcation in each component.

###  2.2 Normal Form for Linear Mixing

To understand how long-term migration directly affects the stability of the
DFE, we rewrite the system as the normal form of a transcritical bifurcation
with linear coupling. Since the SIS model has the same topology near the DFE
as the SIR model, we consider the standard SIS model with births and deaths
[15]

$\displaystyle\frac{ds}{dt}$ $\displaystyle=\mu-\beta si+\kappa i-\mu s,$
(13a) $\displaystyle\frac{di}{dt}$ $\displaystyle=\beta si-\kappa i-\mu i,$
(13b)

with nondimensional variables representing percentages of the population.
Since $s+i=1$, the system is overdetermined and we need only to solve $di/dt$.

In the first subpopulation, let $x=i$ and $1-x=s$. Substitute these variables
into Eq. (13b) and rescale time by $\beta$. We repeat the process using the
variable $y$ to represent the second population. By adding the linear
migration terms to these equations, we find

$\displaystyle\dot{x}$ $\displaystyle=$ $\displaystyle r_{1}x-x^{2}-\alpha
x+\alpha y,$ (14) $\displaystyle\dot{y}$ $\displaystyle=$ $\displaystyle
r_{2}y-y^{2}-\alpha y+\alpha x.$

Here, the bifurcation parameter $r_{k}=(\hat{R}_{k}-1)/\hat{R}_{k}$ for
$k=1,2$ from Eq. (12). In addition, we rescaled the long-term migration rate
as the parameter $\alpha=c_{1}/\beta$.

The steady state $(x,y)=(0,0)$ is equivalent to the DFE in the full system in
Eq. (2). In the absence of coupling ($\alpha=0$), a transcritical bifurcation
will occur in the $x$ system, transferring the stability from $x=0$ to
$x=r_{1}$ at $r_{1}=0$. The dynamics are similar for $y$, respectively.
Linearizing about the steady state $(x,y)=(0,0)$ yields two eigenvalues,

$\displaystyle\Lambda_{1}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(r_{1}+r_{2}-2\alpha+\sqrt{4\alpha^{2}+(r_{1}-r_{2})^{2}}\right),$
(15) $\displaystyle\Lambda_{2}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(r_{1}+r_{2}-2\alpha-\sqrt{4\alpha^{2}+(r_{1}-r_{2})^{2}}\right).$

The following analysis uses this linearization approach to conclude when long-
term migration can change the stability of this steady state.

We start by considering the case of two isolated endemic populations, which
have basic reproduction numbers greater than one. We ask if it is possible to
stabilize the die out state through the coupling parameter, $\alpha$.

######  Proposition 3.

If $r_{1},r_{2}>0$, then the fixed point $(x,y)=(0,0)$ is unstable for all
$\alpha\in[0,\infty)$.

###### Proof.

Upon inspection, $\Lambda_{1}$ is the dominant eigenvalue. Assuming
$r_{1},r_{2}>0$, then $\Lambda_{1}>0$ implies

$(r_{1}+r_{2})+\left(\sqrt{4\alpha^{2}+(r_{1}-r_{2})^{2}}\right)>2\alpha.$
(16)

Squaring both sides and simplifying, we find

$r_{1}^{2}+r_{2}^{2}+(r_{1}+r_{2})\sqrt{4\alpha^{2}+(r_{1}-r_{2})^{2}}>0,$
(17)

which is always true. Therefore, $(x,y)=(0,0)$ is unstable for all
$\alpha\in[0,\infty)$. ∎

We can interpret this abstract result in the original system by concluding
that for two isolated endemic populations, the amount of long-term migration
is irrelevant to the persistence of the disease. The stability of the DFE
cannot be changed by migration and intervention by vaccination is necessary
for disease die out.

Next, consider the case where we have sufficient vaccination so that one of
the basic reproductive numbers is less than one, while the other is not.
Again, we ask under what conditions the coupling parameter can stabilize the
die out state.

######  Proposition 4.

Without loss of generality, we assume $r_{1}>0$ and $r_{2}<0$. Case 1: If
$-r_{2}<r_{1}$, then the fixed point $(x,y)=(0,0)$ is unstable for all
$\alpha\in[0,\infty)$. Case 2: If $-r_{2}>r_{1}$, then there exists some
$\alpha^{*}\in[0,\infty)$ such that the fixed point $(0,0)$ is stable for all
$\alpha>\alpha^{*}$.

###### Proof.

In both cases, assume $r_{1}>0$ and $r_{2}<0$.

Case 1: For $-r_{2}<r_{1}$, $\Lambda_{1}>0$ reduces to the relationship in Eq.
(17). This is always true for $0<r_{1}+r_{2}$ and we conclude $(x,y)=(0,0)$ is
unstable for all $\alpha\in[0,\infty)$.

Case 2: For $-r_{2}>r_{1}$, $\Lambda_{1}<0$ reduces to the relationship

$\left(\sqrt{4\alpha^{2}+(r_{1}-r_{2})^{2}}\right)<2\alpha-(r_{1}+r_{2}).$
(18)

Because both sides are positive, we can square both sides to find

$\alpha>\frac{r_{1}r_{2}}{r_{1}+r_{2}}>0.$ (19)

There exists an $\alpha^{*}=r_{1}r_{2}/(r_{1}+r_{2})$ for
$\alpha^{*}\in[0,\infty)$. Therefore, for $\alpha>\alpha^{*}$, it follows that
$\Lambda_{1}<0$ and $(x,y)=(0,0)$ is stable. ∎

This result implies that in a system with one population supporting an endemic
state, there is a minimum amount of migration necessary for the system to
achieve stability of the DFE. In fact, we can interpret the requirement of
$-r_{2}>r_{1}$ as $y=0$ in the uncoupled system is more stable than $x=0$.
Therefore, $y$ is sharing its extra stability with $x$.

For completeness, we can also show that long-term migration cannot change the
stability of a stable die out state for two isolated populations that have
basic reproduction numbers less than one.

######  Proposition 5.

If $r_{1}$, $r_{2}<0$, then the fixed point $(x,y)=(0,0)$ is stable for all
$\alpha\in[0,\infty)$.

###### Proof.

For $r_{1}$, $r_{2}<0$, $\Lambda_{1}<0$ reduces to Eq. (18). Because both
sides are positive, we can square both sides to find

$\alpha>\frac{r_{1}r_{2}}{r_{1}+r_{2}}.$ (20)

Since $\frac{r_{1}r_{2}}{r_{1}+r_{2}}<0$, $\Lambda_{1}<0$ for all
$\alpha\in[0,\infty)$ and $(x,y)=(0,0)$ is stable. ∎

We conclude that long-term migration has a positive effect on the stability of
the DFE. The mixing in all classes diffuses the force of infection, making it
harder for the disease to persist. In applications where migration is common,
this effect might be significant.

###  2.3 Normal Form for Mass Action Mixing

To capture the effect of short-term migration for each subpopulation in Eq.
(2), we follow the construction of model for linear mixing by substituting $x$
and $y$ for $i$ in Eq. (13b). In this system, we use mass action coupling with
the parameter $\sigma=c_{3}$ controlling the mass action mixing strength.
Specifically, the term $\sigma(1-x)y$ represents the infectious person from
$y$ coming into contact with a susceptible from $x$. The system take the form

$\displaystyle\dot{x}$ $\displaystyle=$
$\displaystyle\left(\frac{r_{1}-\sigma}{1-\sigma}\right)x-x^{2}-\frac{\sigma}{1-\sigma}(1-x)y,$
(21) $\displaystyle\dot{y}$ $\displaystyle=$
$\displaystyle\left(\frac{r_{2}-\sigma}{1-\sigma}\right)y-y^{2}-\frac{\sigma}{1-\sigma}(1-y)x.$

Again, the bifurcation parameter $r_{k}=(\hat{R}_{k}-1)/\hat{R}_{k}$ for
$k=1,2$ from Eq. (12) and time has been rescaled by $\beta(1-\sigma)$.

Performing the linearization about the steady state $(x,y)=(0,0)$ yields two
eigenvalues,

$\displaystyle\Lambda_{1}$ $\displaystyle=$
$\displaystyle\frac{1}{2(1-\sigma)}\,\left(r_{1}+r_{2}-2\sigma+\sqrt{4\,{\sigma}^{2}+(r_{1}-r_{2})^{2}}\right),$
(22) $\displaystyle\Lambda_{2}$ $\displaystyle=$
$\displaystyle\frac{1}{2(1-\sigma)}\,\left(r_{1}+r_{2}-2\sigma-\sqrt{4\,{\sigma}^{2}+(r_{1}-r_{2})^{2}}\right).$
(23)

The eigenvalues for this system are of a similar form as those for linear
migration in Eq. (15), but multiplied by $1/(1-\sigma)$. Therefore, the mass
action mixing can change the stability of the DFE as a function of basic
reproduction numbers in the same settings as linear mixing.

##  3 Vaccination responses

This section directly considers how the migration rates change the vaccination
levels necessary to keep the DFE stable, which implies the occurrence of herd
immunity. It explores whether neglecting the short- and long-term migration
rates overestimates or underestimates the minimum vaccination rates necessary
for disease fade-out.

We first restrict our attention to long-term migration only; i.e., letting
$c_{3}=0$ in Eq. (2). For $c_{1}>0$, this is equivalent to identifying the
bifurcation points in $(v_{1},v_{2})$ when $\lambda_{4}=0$ in Eq. (9). Notice
that if $c_{1}=0$, the subpopulations are isolated. The vaccination levels
needed in each for the disease to die out is equivalent to solving for $v_{1}$
and $v_{2}$ in $\hat{R}_{1,2}\leq 1$ from Eq. (12). The constant solutions for
$\hat{R}_{1,2}=1$ using parameter values in Table 1 are shown in Figure 2 as
solid black lines, and the disease will die out in both populations in the top
right quadrant.

Figure 2:  Boundary of the region of stability for the DFE as we vary $c_{1}$.
The curves represent $\lambda_{4}=0$ using the parameter values in Table 1 and
$c_{3}=0$.The vertical and horizontal black lines represent the vaccination
rates necessary for a stable DFE in isolated populations ($c_{1}=0$). The
dashed line represents the limiting curve as we increase $c_{1}$.

As we increase $c_{1}$, the boundary for the region of stability for the DFE
spreads away from the $c_{1}=0$ case, increasing the die out region. The limit
is a line, shown by the dotted black line in Figure 2.

######  Proposition 6.

As $c_{1}\rightarrow\infty$, the bifurcation curve bounding the stable region
approaches the line

$v_{2}=-\left(\frac{\mu_{1}}{\mu_{2}\rho}\right)v_{1}+\frac{(\mu_{1}+\mu_{2}\rho)(\beta-\kappa)}{\beta\mu_{2}\rho}-\frac{(\mu_{1}+\mu_{2}\rho)^{2}}{\beta\mu_{2}\rho(1+\rho)}.$
(24)

This line is decreasing in $v_{1}$, with the slope depending on a ratio of the
birth rates and subpopulation sizes. Note that if we use the basic
reproductive numbers for the isolated subpopulations, Eq. (24) is equivalent
to

$v_{2}=-\left(\frac{\mu_{1}}{\mu_{2}\rho}\right)v_{1}+\frac{\left(1+\frac{\mu_{1}}{\mu_{2}\rho}\right)}{(1+\rho)}\left(\left(\frac{\hat{R}_{1}(0)-1}{\hat{R}_{1}(0)}\right)+\left(\frac{\hat{R}_{2}(0)-1}{\hat{R}_{2}(0)}\right)\right).$
(25)

As $\hat{R}_{1}(0)$ and/or $\hat{R}_{2}(0)$ increase, the $v_{2}$ intercept
increases and shifts the line up vertically. Therefore, the attainable stable
DFE region in $(v_{1},v_{2})$ space decreases, as expected.

Next, we consider only short-term migration, i.e. letting $c_{1}=0$.
Similarly, for $0\leq c_{3}\leq 1$, this is equivalent to identifying the
bifurcation points in $(v_{1},v_{2})$ when $\lambda_{4}=0$. We fix the
parameters to the values in Table 1 and vary $c_{3}$ to see the changes to the
boundary of the stable region. The solution for $c_{3}=0$ is shown as solid
black lines in Figure 3. Again the disease will die out in both populations in
the top right quadrant. As we increase $c_{3}$, the boundary for the region of
stability for the DFE smoothly pulls away from top right quadrant, increasing
the die out region.

######  Proposition 7.

The limit as we increase $c_{3}\rightarrow 1$ is

$v_{2}=1-{\frac{\left(\kappa+\mu_{{2}}\right)\left(\kappa+\mu_{{1}}\right)}{\left(1-v_{{1}}\right){\beta}^{2}}}.$

Therefore, in both cases, underestimating the migration between populations
causes an overestimation of the vaccination levels needed for herd immunity.

Figure 3:  Boundary of the region of stability for the DFE as we vary $c_{3}$.
The curves represent $\lambda_{4}=0$ using the parameter values in Table 1 and
$c_{1}=0$. The vertical and horizontal black lines represent the vaccination
rates necessary for a stable DFE in isolated populations ($c_{3}=0$).

##  4 Conclusions

In this paper, we consider the effects of short- and long-term migration in
coupled population models in the presence of vaccination. We study the
interplay between the independent vaccination and migration rates across
different populations. We conclude that neglecting migration effects
overestimates the vaccination levels necessary to achieve herd immunity.

We have proven that if two isolated populations support an endemic state
simultaneously, migration cannot change the stability of those endemic states.
Analogously, this is also true for two populations with stable disease free
equilibria. In contrast, migration can lead to disease die out in the mixed
case. If a single population has a vaccination rate sufficient for herd
immunity in isolation, low levels of migration from a population that is
endemic will not necessarily make the disease endemic in both. In fact,
increased levels of migration can lead to disease die out in both populations.
However, migration rates are only physically realistic when they are small.

Our results suggest more efficient vaccination strategies may be identified
for groups of countries with significant migration between them. For example,
instead of increasing the vaccination levels in a population that has already
achieved herd immunity, sending vaccine to the less vaccinated neighboring
country could have a greater impact on outbreak levels. The most efficient
control algorithm would be to target the stable die out region as shown in
Figure 2 or Figure 3.

Conversely, populations for which vaccine delivery is difficult may benefit to
a degree by vaccination of neighboring countries. More specifically, consider
decreasing the migration rates to a country with a lower vaccination rate. We
show in Figure 4 a policy where $N_{2}$, which has a vaccination rate
$v_{2}=0.9$, decreases the long-term migration rate $c_{1}$ with $N_{1}$,
which has a vaccination rate $v_{1}=0.7$. The short term migration rate is
held constant at $c_{3}=0.1$. The decrease in number of new infections for
$N_{2}$ is a small percentage of the increase in infections in $N_{1}$, and we
conclude that the policy meant to help $N_{2}$ has unintended negative
consequences for $N_{1}$.

Figure 4:  Time series of infectives in both populations using the parameter
values in Table 1, with $v_{1}=0.7$, $v_{2}=0.9$, and $c_{3}=0.1$. The value
of $c_{1}$ is decreased to the constant noted in each window.

In future directions, the model can be extended to include the effects of
seasonality. A similar analysis of stable periodic behavior can reveal the
sensitivity of synchronization to short-term and long-term migration. For
example, the work of Schwartz [16] predicts new periodic orbits that can be
excited by the mass action coupling in models with seasonal forcing.
Specifically, these orbits exhibit long period outbreaks in small populations
due to mass action coupling. When applying time dependent vaccination
schedules, other parameters must be considered in addition to the average
vaccination rates, such as pulse frequency and phase with respect to periodic
application. Other techniques can be extended to migration models with vaccine
control, such as prediction of future outbreaks as reported in Schwartz, et
al. [17].

## Acknowledgments

We gratefully acknowledge support from the Office of Naval Research. The
authors were also supported by the National Institute of General Medical
Sciences (Award No. R01GM090204). DATC holds a Career Award at the Scientific
Interface from the Burroughs Wellcome Fund and received funding from the Bill
and Melinda Gates Foundation Vaccine Modeling Initiative. The content is
solely the responsibility of the authors and does not necessarily represent
the official views of the National Institute of General Medical Sciences or
the National Institutes of Health. We would also like to thank Leah Shaw and
Luis Mier-Y-Teran for their useful discussions.

## References

  * [1] N. J. Gay, W. J. Edmunds, Developed countries could pay for hepatitis B vaccination in developing countries, British Medical Journal 316 (1998) 7142.
  * [2] B. Aylward, R. Tangermann, The global polio eradication initiative: Lessons learned and prospects for success, Vaccine Suppl 4 (2011) D80–5.
  * [3] World Health Organization (WHO), Measles Fact Sheet No. 286, December 2009, http://www.who.int/mediacentre/factsheets/fs286/en//.
  * [4] D. P. Wilson, J. Kahn, S. M. Blower, Predicting the epidemiological impact of antiretroviral allocation strategies in KwaZulu-Natal: The effect of the urban-rural divide, PNAS 103 (2006) 14228–14233.
  * [5] D. A. T. Cummings, W. J. Moss, K. Long, C. S. Wiysonge, T. J. Muluh, B. Kollo, E. Nomo, N. D. Wolfe, D. S. Burke, Improved measles surveillance in Cameroon reveals two major dynamic patterns of incidence, International Journal of Infectious Diseases 10 (2) (2006) 148–155.
  * [6] N. M. Ferguson, C. A. Donnelly, R. M. Anderson, Transmission dynamics and epidemiology of dengue: Insights from age-stratified sero-prevalence surveys, Phil. Trans. R. Soc. London, Ser. B 354 (1999) 757–768.
  * [7] L. S. Liebovitch, I. B. Schwartz, Migration induced epidemics: Dynamics of flux-based multipatch models, Physics Letters A 332 (2004) 256–267.
  * [8] L. Sattenspiel, K. Dietz, A structured epidemic model incorporating geographic mobility among regions, Mathematical Biosciences 128 (1995) 71–91.
  * [9] A. L. Lloyd, V. A. A. Jansen, Spatiotemporal dynamics of epidemics: Synchrony in metapopulation models, Mathematical Biosciences 188 (2004) 1–16.
  * [10] M. J. Keeling, P. Rohani, Estimating spatial coupling in epidemiological systems: A mechanistic approach, Ecology Letters 5 (1) (2002) 20–29.
  * [11] J. Cui, Y. Takeuchi, Y. Saito, Spreading disease with transport-related infection, Journal of Theoretical Biology 239 (3) (2006) 376 – 390.
  * [12] Y. Takeuchi, X. Liu, J. Cui, Global dynamics of SIS models with transport-related infection, Journal of Mathematical Analysis and Applications 329 (2) (2007) 1460 – 1471.
  * [13] J. Liu, Y. Zhou, Global stability of an SIRS epidemic model with transport-related infection, Chaos, Solitons & Fractals 40 (1) (2009) 145 – 158\.
  * [14] H. W. Hethcote, The mathematics of infectious diseases, SIAM Review 42 (2000) 599–653.
  * [15] J. D. Murray, Mathematical Biology, Springer, Berlin, 1989.
  * [16] I. B. Schwartz, Small amplitude, long period outbreaks in seasonally driven epidemics, J. Math. Biology 30 (1992) 473–491.
  * [17] I. B. Schwartz, L. Billings, E. M. Bollt, Dynamical epidemic suppression using stochastic prediction and control, Physical Review E 70 (2004) 046220.