\begin{matrix}c \\ ◯ \end{matrix} Revista de Matemática: Teoría y Aplicaciones 2020 27(1) : 23–48 cimpa – ucr issn: 1409-2433 (Print), 2215-3373 (Online)

doi: https:// doi. org/ 10.15517/ rmta. v27i1. 39946

a two-patch epidemic model with nonlinear reinfection

un modelo epidémico de dos poblaciones con reinfección no lineal

Juan G. Calvo^{∗}  Alberto Hernández^{†}
Mason A. Porter^{‡}  Fabio Sanchez^{§}

Received: 27/Jul/2019; Revised: 17/Sep/2019; Accepted: 31/Oct/2019

∗University of Costa Rica, CIMPA-School of Mathematics, San José, Costa Rica. E-Mail: juan.calvo@ucr.ac.cr

†Misma dirección que/Same address as: J.G. Calvo. E-Mail: albertojose.hernandez@ucr.ac.cr ‡University of California Los Angeles, Department of Mathematics, Los Angeles CA, United

States of America. E-Mail: mason@math.ucla.edu §Misma dirección que/Same address as: J.G. Calvo. E-Mail: fabio.sanchez@ucr.ac.cr

24 j.g. calvo — a. hernández — m.a. porter — f. sanchez

Abstract

The propagation of infectious diseases and its impact on individuals play a major role in disease dynamics, and it is important to incorporate population heterogeneity into efforts to study diseases. As a simplistic but illustrative example, we examine interactions between urban and rural populations on the dynamics of disease spreading. Using a compartmental framework of susceptible–infected–susceptible (SIS̃) dynamics with some level of immunity, we formulate a model that allows nonlinear reinfection. We investigate the effects of population movement in a simple scenario: a case with two patches, which allows us to model population movement between urban and rural areas. To study the dynamics of the system, we compute a basic reproduction number for each population (urban and rural). We also compute steady states, determine the local stability of the disease-free steady state, and identify conditions for the existence of endemic steady states. From our analysis and computational experiments, we illustrate that population movement plays an important role in disease dynamics. In some cases, it can be rather beneficial, as it can enlarge the region of stability of a disease-free steady state.

Keywords: dynamical systems; population dynamics; mathematical modeling; biological contagions; population movement.

Resumen

La propagación de enfermedades infecciosas y su impacto en individuos juega un gran rol en la dinámica de enfermedades, y es importante incorporar heterogeneidad en la población en los esfuerzos por estudiar enfermedades. De manera simplística pero ilustrativa, se examinan interacciones entre una población urbana y una rural en la dinámica de la propagación de una enfermedad. Utilizando un sistema compartimental de dinámicas entre susceptibles–infectados–susceptibles (SIS̃) con cierto nivel de inmunidad, se formula un modelo que permite reinfecciones no lineales. Se investiga los efectos de movimiento de poblaciones en un escenario simple: un caso con dos poblaciones, que permite modelar movimiento entre un área urbana y otra rural. Con el fin de estudiar la dinámica del sistema, se calcula el número básico reproductivo para cada comunidad (rural y urbana). Se calculan también puntos de equilibrio, la estabilidad local del estado libre de enfermedad, y se identifican condiciones para la existencia de estados de equilibrio endémicos. Del análisis y experimentos computacionales, se ilustra que el movimiento en la población juega un rol importante en la dinámica del sistema. En algunos casos, puede ser beneficioso, pues incrementa la región de estabilidad del punto de equilibrio del estado libre de infección.

Palabras clave: sistemas dinámicos; dinámica de poblaciones; modelado matemático; contagios biológicos; movimiento de poblaciones.

Mathematics Subject Classification: 92D25, 92D30.

1 Introduction

It is relatively easy for individuals to move between towns, cities, countries, and even continents; and incorporating movement between populations has become increasingly prevalent in the modeling and analysis of disease spreading [3, 22]. It is also important to consider movement in which individuals travel to a distinct location from their place of origin and then return to their original location in a relatively short period of time. Such movement can lead to rapid spreading of infectious diseases, and examination of connected environments can give clues about the types of strategies that are needed to control disease propagation [2, 20].

In 2003, the severe acute respiratory syndrome (SARS) epidemic was a major concern among public health officials worldwide [4]. This new infectious disease spread rapidly, and scientists and researchers scrambled to try to discern how to contain its spread (e.g., by reducing the spreading rate) and to seek treatments and a vaccine. The best control measure that was found at the time was to isolate individuals who had been in contact with infected individuals. The rapid spread of the disease was associated with the movement of a doctor who was identified as “patient 0” for SARS [7, 14]. Population movement has also played an important role in subsequent events, such as the spread of ebola to the western hemisphere [5], the spread of measles in some parts of the world by travelers [6], and the resurgence of malaria through the mass migration of Nicaraguans to the northern part of Costa Rica [1].

The use of compartmental models to describe the spreading of diseases has been explored thoroughly in numerous scenarios [3, 22]. For example, when there is nonlinear reinfection, an individual who was infected previously can become infected again through contact with an infectious individual after losing immunity [30]. Several models that allow individuals to lose immunity and become infectious again also exhibit a backward bifurcation, such that a stable endemic steady state coexists with a stable disease-free steady state when the associated basic reproduction number is smaller than 1 [12, 24, 27, 28, 26]. Moreover, in some models of social contagion processes (e.g., the spread of drug use, the adoption of products, and so on), after an initial “contagion”, a backward bifurcation can arise via social inputs [27].

In the present paper, we generalize the compartmental model from Sanchez et al. [24], who studied a continuous dynamical system (in the form of coupled ordinary differential equations) that describes interactions between susceptible and infected individuals with the possibility of reinfection after loss of immunity.

In our generalization, we incorporate population movement between urban and rural environments. Taking movement into account is important for studies of disease dynamics in practice, and it changes the qualitative dynamics of disease spreading. We explore a simple case in which a population has two patches, and we obtain insights that will be useful for subsequent explorations of disease spreading in a population that includes a network of patches.

Our paper proceeds as follows. In Section 2, we describe our compartmental model of disease spreading (including nonlinear reinfection) between urban and rural environments. In Section 3, we give a formula for the model’s basic reproduction number R_{0}, analyze the existence and local stability of the diseasefree steady state, and study the existence of endemic steady states. In Section 4, we illustrate several example scenarios with numerical computations. Finally, in Section 5, we conclude and discuss biological insights of our model.

2 A two-patch compartmental model

We present a two-patch model of disease spreading in humans that incorporates nonlinear reinfection and population movement. Specifically, we generalize the model of Sanchez et al. [24] by incorporating the idea of urban versus rural environments. We use the subscript u for urban variables and parameters and the subscript r for rural variables and parameters. For j ∈ {u, r}, let S_{j}, I_{j}, and S̃_{j} denote the numbers of susceptible, infected, and post-recovery susceptible individuals, respectively. New susceptible individuals enter the system in proportion to the total population N_{j}, where N_{j} = S_{j} + I_{j} + S̃_{j}. Let μ_{j} denote the rate of both births and deaths. (For simplicity, we assume that they are the same.) Susceptible individuals become infected at rate β_{j}, infected individuals transition to a state of post-recovery susceptibility at rate γ_{j}, and post-recovery susceptible individuals are reinfected at rate ρ_{j}. Such reinfection corresponds to infectious diseases (such as tuberculosis and malaria [15, 16]) in which subsequent infections are possible after loss of immunity. Typically, in such diseases, an initial infection tends to produce stronger symptoms than subsequent reinfections [9]. We model population movement between patches using the functions δ_{ij}(t), which denotes the fraction of individuals who travel from patch i∈{u, r} to patch j ∈{u, r} (with i≠j) at time t.

Our model consists of the following coupled system of ordinary differential equations:

\begin{aligned}\frac{dS_{u}}{dt} & = & μ_{u}N_{u} − β_{u} \frac{I_{u}}{N_{u}} S_{u} − μ_{u}S_{u} + δ_{ru}S_{r} − δ_{ur}S_{u} , \\ \frac{dI_{u}}{dt} & = & β_{u} \frac{I_{u}}{N_{u}} S_{u} − (μ_{u} + γ_{u})I_{u} + ρ_{u} \frac{I_{u}}{N_{u}} S̃_{u} + δ_{ru}I_{r} − δ_{ur}I_{u} , \\ \frac{dS̃_{u}}{dt} & = & γ_{u}I_{u} − ρ_{u} \frac{I_{u}}{N_{u}} S̃_{u} − μ_{u}S̃_{u} + δ_{ru}S̃_{r} − δ_{ur}S̃_{u} , \\ \frac{dS_{r}}{dt} & = & μ_{r}N_{r} − β_{r} \frac{I_{r}}{N_{r}} S_{r} − μ_{r}S_{r} + δ_{ur}S_{u} − δ_{ru}S_{r} , \\ \frac{dI_{r}}{dt} & = & β_{r} \frac{I_{r}}{N_{r}} S_{r} − (μ_{r} + γ_{r})I_{r} + ρ_{r} \frac{I_{r}}{N_{r}} S̃_{r} + δ_{ur}I_{u} − δ_{ru}I_{r} , \\ \frac{dS̃_{r}}{dt} & = & γ_{r}I_{r} − ρ_{r} \frac{I_{r}}{N_{r}} S̃_{r} − μ_{r}S̃_{r} + δ_{ur}S̃_{u} − δ_{ru}S̃_{r} , \end{aligned} (1)

with initial conditions S_{j}(0) = S_{j0}, I_{j}(0) = I_{j0}, and S̃_{j}(0) = S̃_{j0} (for j ∈{u, r}). Additionally, N_{j0} = S_{j0} + I_{j0} + S̃_{j0} for j ∈{u, r}.

By adding the first three and last three equations in (1), we see that N_{u} and N_{r} satisfy the linear dynamical system

\frac{d}{dt} [ \begin{matrix}N_{u} \\ N_{r} \end{matrix} ] = [ \begin{matrix}−δ_{ur} & δ_{ru} \\ δ_{ur} & −δ_{ru} \end{matrix} ] [ \begin{matrix}N_{u} \\ N_{r} \end{matrix} ] .

Its solution is

N_{u}(t) = e^{− ∫_{0}^{t} δ(s)ds} [N_{u,0} + (N_{u,0} + N_{r,0}) ∫_{0}^{t} δ_{ru}(s)e^{∫_{0}^{s} δ(h)dh} ds] ,
N_{r}(t) = e^{− ∫_{0}^{t} δ(s)ds} [N_{r,0} + (N_{u,0} + N_{r,0}) ∫_{0}^{t} δ_{ur}(s)e^{∫_{0}^{s} δ(h)dh} ds] ,

where δ(t) = δ_{ur}(t) + δ_{ru}(t) is the net movement of individuals at time t. The initial conditions are N_{u,0} = N_{u}(0) and N_{r,0} = N_{r}(0). Note that N_{u} + N_{r} is constant for all t.

In our analysis, we assume at first that δ_{ru} and δ_{ur} are constant, so N_{u} and N_{r} simplify to

N_{u}(t) = N_{u,0}e^{−t(δ_{ur}+δ_{ru})} + (N_{u,0} + N_{r,0})\frac{δ_{ru}}{δ_{ru} + δ_{ur}} (1 − e^{−t(δ_{ur}+δ_{ru})}) ,
N_{r}(t) = N_{r,0}e^{−t(δ_{ur}+δ_{ru})} + (N_{u,0} + N_{r,0})\frac{δ_{ur}}{δ_{ru} + δ_{ur}} (1 − e^{−t(δ_{ur}+δ_{ru})}) .

There are three cases:

1. There is no population movement (i.e., δ_{ur} = δ_{ru} = 0), which corresponds to the case of independent patches that was studied in [24]. In this case, each patch has a backward bifurcation, and the steady state can depend on the number of initially infected individuals.

2. Population movement occurs exclusively from one patch to the other. For example, suppose that δ_{ur} = 0 and δ_{ru} > 0. In this case, lim_{t→∞} N_{r}(t) = 0 (i.e., eventually, the entire population is in u). There is a backward bifurcation in this case as well, but now there is a total population of N_{u,0}+N_{r,0} in the urban patch. We explore this case in Example 1 in Section 4.

3. There is population movement in both directions between the two patches (i.e., δ_{ur} > 0 and δ_{ru} > 0). This is the primary scenario (and the principal novel contribution) of the present paper.

3 Analysis of our model

3.1 Disease-free steady state and basic reproduction number

Because the total population is constant, we can eliminate S̃_{r} from the dynamical system (1). We then compute the Jacobian matrix of the reduced system and evaluate it at the disease-free steady state

(S_{u}^{∗}, I_{u}^{∗}, S̃_{u}^{∗}, S_{r}^{∗}, I_{r}^{∗}, S̃_{r}^{∗}) = (S_{u}^{∗}, 0, 0, S_{r}^{∗}, 0, 0)  (2)

to obtain the matrix

[ \begin{matrix}−δ_{ur} & μ_{u} − β_{u} & μ_{u} & δ_{ru} & 0 \\ 0 & η_{ur} & 0 & 0 & δ_{ru} \\ −δ_{ru} & γ_{u} − δ_{ru} & −δ_{ru} − δ_{ur} − μ_{u} & −δ_{ru} & −δ_{ru} \\ δ_{ur} − μ_{r} & −μ_{r} & −μ_{r} & −μ_{r} − δ_{ru} & −β_{r} \\ 0 & δ_{ur} & 0 & 0 & η_{ru} \end{matrix} ] , (3)

where η_{ur} = β_{u}−δ_{ur} −γ_{u}−μ_{u} and η_{ru} = β_{r} −δ_{ru}−γ_{r} −μ_{r}. For the dynamical system (1), steady states with I_{u} = I_{r} = 0 necessarily also satisfy S̃_{u} = S̃_{r} = 0.

The five eigenvalues of the matrix (3) are

λ_{1} = −(δ_{ru} + δ_{ur}) ,
λ_{2±} = \frac{1}{2} [−(δ_{ru} + δ_{ur} + μ_{r} + μ_{u}) ± \sqrt{(δ_{ru} − δ_{ur} + μ_{r} − μ_{u})^{2} + 4δ_{ru}δ_{ur}}] ,
λ_{3±} = \frac{1}{2} [η_{ur} + η_{ru} ± \sqrt{(η_{ur} − η_{ru})^{2} + 4δ_{ru}δ_{ur}}] .

All eigenvalues are real, and λ_{1} and λ_{2±} are negative. The two remaining eigenvalues λ_{3±} are negative as long as

δ_{ru}δ_{ur} < η_{ur}η_{ru} , η_{ur} + η_{ru} < 0 .  (4)

When there is only one population (i.e., one patch), the basic reproduction number is R_{0} = \frac{β}{μ+γ}. For our multiple-patch case, we define a “local basic reproduction number” for each patch:

R_{0u} = \frac{β_{u}}{γ_{u} + μ_{u}} , R_{0r} = \frac{β_{r}}{γ_{r} + μ_{r}} .

We can then express the conditions in (4) as

R_{0u} < 1 + \frac{δ_{ur}}{μ_{u} + γ_{u}} ,  (5a)
R_{0r} < 1 + \frac{δ_{ru}}{μ_{r} + γ_{r}} ,  (5b)
\frac{δ_{ur}}{μ_{u} + γ_{u}} \frac{δ_{ru}}{μ_{r} + γ_{r}} < (R_{0u} − 1 − \frac{δ_{ur}}{μ_{u} + γ_{u}} ) (R_{0r} − 1 − \frac{δ_{ru}}{μ_{r} + γ_{r}} ) . (5c)

See Figure 1 for an illustration of a typical region in which all eigenvalues are negative. For progressively smaller δ_{ur} and δ_{ru}, the shaded region approaches the unit square.

Figure 1: A typical region of local asymptotic stability of the disease-free steady state (2) of (1) that satisfies the conditions in (5) as a function of the local basic reproduction numbers R_{0u} and R_{0r}.

We have thus established the following lemma.

Lemma 1 Assume that the inequalities (5) hold. It then follows that the diseasefree steady state (2) of the dynamical system (1) is locally asymptotically stable.

Remark 1 The conditions in (5) are satisfied when R_{0u} < 1 and R_{0r} < 1. We then have local stability in the rural and urban patches if we treat them as independent. Additionally, from Lemma 1, we see that it is possible to obtain local asymptotic stability for the disease-free steady state (2) even when one or both local basic reproduction numbers are larger than 1. In such a scenario, population movement is beneficial, as it leads to local asymptotic stability of the disease-free steady state in situations in which this is not the case for independent patches. We illustrate such a scenario in Example 2 in Section 4.

3.2 Technical tool: the Poincaré–Miranda theorem

As preparation for analyzing the existence of steady states in our model, we briefly recall some classical results by Poincaré and others. See [13, 17, 29, 31] for detailed accounts of the relevant theory. In 1817, Bolzano proved the following well-known theorem.

Theorem 1 Suppose that f : [a, b] ⟶ ℝ is a continuous function such that f(a)·f(b)<0. There then exists c∈(a, b) such that f(c)=0.

Definition 1 Let I^{n} = [0, 1]^{n}, and let ∂I^{n} denote its boundary. For each i∈{1, …, n}, let

I_{i}^{−} = {x∈I^{n}|x_{i} = 0} , I_{i}^{+} = {x∈I^{n}|x_{i} = 1}

be the opposite i-th faces of the boundary ∂I^{n}.

In 1883–1884, Poincaré announced a generalization of Bolzano’s theorem without providing a proof.

Theorem 2 (Poincaré) Let F : I^{n} ⟶ ℝ^{n}, with F = (F_{1}, …, F_{n}), be a continuous map such that

F_{i}(I_{i}^{−}) ⊆ (−∞, 0]

and

F_{i}(I_{i}^{+}) ⊆ [0, +∞)

for every i ∈ {1, …, n}. It then follows that there exists c ∈ I^{n} such that F (c) = 0.

In the 1940s, Miranda rediscovered Poincaré’s theorem and showed that it is logically equivalent to Brouwer’s fixed-point theorem. Since then, this result has often been called the Poincaré–Miranda theorem. We require a modified version of the Poincaré–Miranda theorem in two dimensions.

Proposition 1 Let F : [0, 1]^{2} ⟶ ℝ^{2} such that F (x, y) = (F_{1}(x, y), F_{2}(x, y))^{T} is continuous and F (0, 0) = 0. Assume that

F_{1}(x, 0) > 0 for all x ∈ (0, 1] , F_{1}(x, 1) < 0 for all x ∈ [0, 1] ,
F_{2}(0, y) > 0 for all y ∈ (0, 1] , F_{2}(1, y) < 0 for all y ∈ [0, 1] .

Assume additionally that ∂F_{1}/∂y and ∂F_{2}/∂x are both continuous from the right at (0, 0), with \frac{∂F_{1}}{∂y} (0, 0) > 0 and \frac{∂F_{2}}{∂x} (0, 0) > 0. It then follows that there exists (x_{0}, y_{0})∈(0, 1)^{2} such that F (x_{0}, y_{0}) = (0, 0)^{T}.

Proof. Because F is continuous, both F_{1} and F_{2} are also continuous. Using the facts that F_{1}(x, 0) > 0 for all x ∈ (0, 1] and that ∂F_{1}/∂x is continuous from the right, it follows that there exists ϵ_{0} > 0 such that \frac{∂F_{1}}{∂x} (x, 0) > 0 for all x ∈ (0, ϵ_{0}). By the implicit function theorem, for x ∈ (0, ϵ_{0}), one can write y = g_{1}(x), where the function g_{1} is differentiable. Therefore,

g′_{1}(0) = − \frac{\frac{∂F_{1}}{∂x} (0, 0)}{\frac{∂F_{1}}{∂y} (0, 0)} .  (6)

By the continuity of F_{1} in the compact set [0, 1]^{2} and using Equation (6), it follows that there exists ϵ_{1} > 0 such that F_{1}(x, ϵ_{1}) > 0 for all x ∈ [0, 1]. By the same argument, there also exists ϵ_{2} ∈ (0, 1) such that F_{2}(ϵ_{2}, y) > 0 for all y ∈ [0, 1]. By the Poincaré–Miranda theorem, there must exist (x_{0}, y_{0})∈[ϵ_{1}, 1] × [ϵ_{2}, 1] such that F (x_{0}, y_{0}) = (0, 0).

3.3 Existence of multiple-population endemic steady states

We now examine endemic steady states, in which there are infected individuals at steady state in both the urban and the rural patches.

We start by proving the following lemma.

Lemma 2 Suppose that R_{0u} > 1 + \frac{δ_{ur}}{μ_{u} + γ_{u}} and R_{0r} > 1 + \frac{δ_{ru}}{μ_{r} + γ_{r}}. It then

follows that there exists at least one endemic steady state, for which both I_{u}^{∗} > 0 and I_{r}^{∗} > 0. In other words, there are infected individuals at steady state in both the urban and the rural patches.

Proof. We deduce the existence of a solution of the following nonlinear system of algebraic equations:

\begin{aligned}0 & = & μ_{u}N_{u} − β_{u} \frac{I_{u}}{N_{u}} S_{u} − μ_{u}S_{u} + δ_{ru}S_{r} − δ_{ur}S_{u} , \\ 0 & = & β_{u} \frac{I_{u}}{N_{u}} S_{u} − (μ_{u} + γ_{u})I_{u} + ρ_{u} \frac{I_{u}}{N_{u}} S̃_{u} + δ_{ru}I_{r} − δ_{ur}I_{u} , \\ 0 & = & γ_{u}I_{u} − ρ_{u} \frac{I_{u}}{N_{u}} S̃_{u} − μ_{u}S̃_{u} + δ_{ru}S̃_{r} − δ_{ur}S̃_{u} , \\ 0 & = & μ_{r}N_{r} − β_{r} \frac{I_{r}}{N_{r}} S_{r} − μ_{r}S_{r} + δ_{ur}S_{u} − δ_{ru}S_{r} , \\ 0 & = & β_{r} \frac{I_{r}}{N_{r}} S_{r} − (μ_{r} + γ_{r})I_{r} + ρ_{r} \frac{I_{r}}{N_{r}} S̃_{r} + δ_{ur}I_{u} − δ_{ru}I_{r} , \\ 0 & = & γ_{r}I_{r} − ρ_{r} \frac{I_{r}}{N_{r}} S̃_{r} − μ_{r}S̃_{r} + δ_{ur}S̃_{u} − δ_{ru}S̃_{r} . \end{aligned}  (7)

We introduce a pair of parameters, (b_{u}, b_{r}) ∈ [0, 1]^{2}, and we consider the auxiliary linear system

\begin{aligned}0 & = & μ_{u}N_{u} − β_{u}b_{u}S_{u} − μ_{u}S_{u} + δ_{ru}S_{r} − δ_{ur}S_{u} , \\ 0 & = & β_{u}b_{u}S_{u} − (μ_{u} + γ_{u})I_{u} + ρ_{u}b_{u}S̃_{u} + δ_{ru}I_{r} − δ_{ur}I_{u} , \\ 0 & = & γ_{u}I_{u} − ρ_{u}b_{u}S̃_{u} − μ_{u}S̃_{u} + δ_{ru}S̃_{r} − δ_{ur}S̃_{u} , \\ 0 & = & μ_{r}N_{r} − β_{r}b_{r}S_{r} − μ_{r}S_{r} + δ_{ur}S_{u} − δ_{ru}S_{r} , \\ 0 & = & β_{r}b_{r}S_{r} − (μ_{r} + γ_{r})I_{r} + ρ_{r}b_{r}S̃_{r} + δ_{ur}I_{u} − δ_{ru}I_{r} , \\ 0 & = & γ_{r}I_{r} − ρ_{r}b_{r}S̃_{r} − μ_{r}S̃_{r} + δ_{ur}S̃_{u} − δ_{ru}S̃_{r} . \end{aligned}  (8)

Suppose that (S_{u}, I_{u}, S̃_{u}, S_{r}, I_{r}, S̃_{r}) is a solution of the system (8), which is linear in its variables S_{i}, I_{i}, and S̃_{i} (with i ∈ {u, r}). This solution also satisfies the nonlinear system (7) if I_{u}/N_{u} = b_{u} and I_{r}/N_{r} = b_{r}. By adding all of the equations in (8), we see that δ_{ur}N_{u} = δ_{ru}N_{r}. We then rescale variables in (8) by substituting

s_{u} = \frac{S_{u}}{N_{u}} , s_{r} = \frac{S_{r}}{N_{r}} ,
i_{u} = \frac{I_{u}}{N_{u}} , i_{r} = \frac{I_{r}}{N_{r}} ,
s̃_{u} = \frac{S̃_{u}}{N_{u}} , s̃_{r} = \frac{S̃_{r}}{N_{r}}

a two-patch epidemic model with nonlinear reinfection 33 to obtain a linear system of algebraic equations with parameters b_{u} and b_{r}. This system is

μ_{u} − (β_{u}b_{u} + μ_{u} + δ_{ur}) s_{u} + δ_{ur}s_{r} = 0 ,  (9a)
β_{u}b_{u}s_{u} − (μ_{u} + γ_{u} + δ_{ur})i_{u} + ρ_{u}b_{u}s̃_{u} + δ_{ur}i_{r} = 0 ,  (9b)
γ_{u}i_{u} − ρ_{u}b_{u}s̃_{u} − (μ_{u} + δ_{ur}) s̃_{u} + δ_{ur}s̃_{r} = 0 ,  (9c)
μ_{r} − (β_{r}b_{r} + μ_{r} + δ_{ru}) s_{r} + δ_{ru}s_{u} = 0 ,  (9d)
β_{r}b_{r}s_{r} − (μ_{r} + γ_{r} + δ_{ru})i_{r} + ρ_{r}b_{r}s̃_{r} + δ_{ru}i_{u} = 0 ,  (9e)
γ_{r}i_{r} − ρ_{r}b_{r}s̃_{r} − (μ_{r} + δ_{ru}) s̃_{r} + δ_{ru}s̃_{u} = 0 .  (9f)

We now solve Equations (9a) and (9d) to obtain

s_{u}(b_{u}, b_{r}) = \frac{β_{r}μ_{u}b_{r} + δ_{ur}μ_{r} + δ_{ru}μ_{u} + μ_{r}μ_{u}}{(β_{u}b_{u} + μ_{u} + δ_{ur}) (β_{r}b_{r} + μ_{r} + δ_{ru}) − δ_{ru}δ_{ur}} ,
s_{r}(b_{u}, b_{r}) = \frac{β_{u}μ_{r}b_{u} + δ_{ur}μ_{r} + δ_{ru}μ_{u} + μ_{r}μ_{u}}{(β_{u}b_{u} + μ_{u} + δ_{ur}) (β_{r}b_{r} + μ_{r} + δ_{ru}) − δ_{ru}δ_{ur}} .

Using Equations (9b), (9c), (9e), and (9f), we solve for i_{u} and i_{r} to obtain

i_{u}(b_{u}, b_{r}) = \frac{1}{Δ} (s_{r}β_{r}δ_{ur}b_{r}Δ_{i_{u}}^{(1)} + s_{u}β_{u}b_{u}Δ_{i_{u}}^{(2)}) ,
i_{r}(b_{u}, b_{r}) = \frac{1}{Δ} (s_{u}β_{u}δ_{ru}b_{u}Δ_{i_{r}}^{(1)} + s_{r}β_{r}b_{r}Δ_{i_{r}}^{(2)}) ,

where

Δ_{i_{u}}^{(1)} = δ_{ur}μ_{r} + δ_{ru}μ_{u} + μ_{u}μ_{r} + ρ_{r}(μ_{u} + δ_{ur})b_{r} + ρ_{u}(γ_{r} + μ_{r} + δ_{ru})b_{u}
+ ρ_{u}ρ_{r}b_{u}b_{r} ,
Δ_{i_{u}}^{(2)} = (δ_{ru} + γ_{r} + μ_{r})(δ_{ur}μ_{r} + δ_{ru}μ_{u} + μ_{r}μ_{u}) + (δ_{ru} + μ_{r})(δ_{ur} + μ_{u})ρ_{r}b_{r}
+ (δ_{ru} + μ_{r})(δ_{ru} + γ_{r} + μ_{r})ρ_{u}b_{u} + (δ_{ru} + μ_{r})ρ_{u}ρ_{r}b_{u}b_{r} ,
Δ_{i_{r}}^{(1)} = δ_{ur}μ_{r} + δ_{ru}μ_{u} + μ_{u}μ_{r} + ρ_{u}(μ_{r} + δ_{ru})b_{u} + ρ_{r}(γ_{u} + μ_{u} + δ_{ur})b_{r}
+ ρ_{u}ρ_{r}b_{u}b_{r} ,
Δ_{i_{r}}^{(2)} = (δ_{ur} + γ_{u} + μ_{u})(δ_{ur}μ_{r} + δ_{ru}μ_{u} + μ_{r}μ_{u}) + (δ_{ur} + μ_{u})(δ_{ru} + μ_{r})ρ_{u}b_{u}
+ (δ_{ur} + μ_{u})(δ_{ur} + γ_{u} + μ_{u})ρ_{r}b_{r} + (δ_{ur} + μ_{u})ρ_{u}ρ_{r}b_{u}b_{r} ,
Δ = (δ_{ur}μ_{r} + δ_{ru}μ_{u} + μ_{r}μ_{u})[(δ_{ru} + γ_{r} + μ_{r} + b_{r}ρ_{r})
× (δ_{ur} + γ_{u} + μ_{u} + b_{u}ρ_{u}) − δ_{ur}δ_{ru}] .

Rev.Mate.Teor.Aplic. (ISSN print: 1409-2433; online: 2215-3373) Vol. 27(1): 23–48, Jan–Jun 2020

One can write similar expressions for s̃_{u}(b_{u}, b_{r}) and s̃_{r}(b_{u}, b_{r}). The solution (s_{u}, i_{u}, s̃_{u}, s_{r}, i_{r}, s̃_{r}) of (9) satisfies the nonlinear algebraic system (7) if and only if i_{u}(b_{u}, b_{r}) = b_{u} and i_{r}(b_{u}, b_{r}) = b_{r}. We then define

G_{1}(b_{u}, b_{r}) = i_{u}(b_{u}, b_{r}) − b_{u} ,
G_{2}(b_{u}, b_{r}) = i_{r}(b_{u}, b_{r}) − b_{r} .

We will show that the system

G_{1}(b_{u}, b_{r}) = G_{2}(b_{u}, b_{r}) = 0  (10)

has at least one solution (b_{u}, b_{r})∈(0, 1)^{2}. This implies that there exists a solution of the nonlinear algebraic system (7) with i_{u} > 0 and i_{r} > 0. This solution corresponds to a steady state with a nonzero infected population in both the urban and the rural patches.

Consider the surface z_{1} =G_{1}(b_{u}, b_{r}). We have that G_{1}(0, b_{r})=i_{u}(0, b_{r})≥0 for b_{r} ∈ [0, 1], with equality when b_{r} = 0. We seek a positive solution of Equation (10). A straightforward calculation yields

\frac{∂G_{1}}{∂b_{u}} (0, 0) = \frac{−δ_{ur}(γ_{r} + μ_{r}) + (β_{u} − γ_{u} − μ_{u})(δ_{ru} + γ_{r} + μ_{r})}{δ_{ur}(γ_{r} + μ_{r}) + (δ_{ru} + γ_{r} + μ_{r})(γ_{u} + μ_{u})}

> \frac{δ_{ur}δ_{ru}}{δ_{ur}(γ_{r} + μ_{r}) + (δ_{ru} + γ_{r} + μ_{r})(γ_{u} + μ_{u})} > 0 ,

because β_{u} > γ_{u} + μ_{u} + δ_{ur} (by hypothesis). Analogously, we compute that

\frac{∂G_{2}}{∂b_{r}} (0, 0) > 0 .  (11)

We also compute that G_{1}(1, b_{r}) = i_{u}(1, b_{r}) − 1 < 0 for b_{r} ∈ [0, 1] and that G_{2}(b_{u}, 1) < 0 for b_{u} ∈ [0, 1]. By Proposition 1, we conclude that there exists (b_{u}, b_{r})∈(ϵ_{1}, 1) × (ϵ_{2}, 1) such that G_{1}(b_{u}, b_{r}) = G_{2}(b_{u}, b_{r}) = 0. Therefore, the nonlinear algebraic system (7) has a solution that corresponds to a steady state with positive values for both I_{u}^{∗} and I_{r}^{∗}.

(a) The curve G_{1}(b_{u}, b_{r}) = 0 has a tangent with slope m_{u} <0 at the origin.

(b) The curves G_{1}(b_{u}, b_{r}) = 0 and G_{2}(b_{u}, b_{r}) = 0 intersect at two points with positive coordinates.

Figure 2: Endemic steady states in urban (vertical axis) and rural (horizontal axis) environments. In these steady states, an infected population persists in both the rural and the urban patches.

Remark 2 There are at most 9 steady states, because Equation (10) has at most 9 solutions. (It is equivalent to a polynomial equation of degree 9.) Using numerical computations, we observe the existence of two or more distinct endemic steady states. However, we need to explore them further to characterize them; see Figure 2b and Examples 3 and 7 in Section 4.

4 Examples

We now present some numerical simulations of the dynamical system (1) for a variety of parameter values. Using these examples, we illustrate that population movement strongly influences how a disease can propagate. For our simulations, we use Matlab and its ordinary differential equation solver ode45 (see the appendix).

Example 1 We first explore the behavior of our model (1) when δ_{ur} = 0 and δ_{ru} ≠ 0, which describes population movement in one direction (specifically, from the rural patch to the urban one). In some countries, it is common for individuals in rural areas to travel to urban areas for work [8]. This is a type of short-term mobility. We consider the following parameter values:

\begin{matrix}μ_{u} = 1/(365 · 80) , & ρ_{u} = 0.08 , & γ_{u} = 0.01 , & β_{u} = 0.03 , \\ μ_{r} = 1/(365 · 70) , & ρ_{r} = 0.04 , & γ_{r} = 0.01 , & β_{r} = 0.02 , \end{matrix}

with initial conditions

\begin{matrix}S_{u0} = 999 , & I_{u0} = 1 , & S̃_{u0} = 0 , \\ S_{r0} = 300 , & I_{r0} = 0 , & S̃_{r0} = 0 . \end{matrix}

In this case,

R_{0u} > 1 and R_{0r} ≈ 2 .

In Figure 3, we show I_{u}(t) for different values of δ_{ru}. We show the solutions for both the urban and the rural patches when δ_{ru} = 0 (i.e., when there is no population movement) and δ_{ru} = 0.01 in Figure 4. We observe that population movement in one direction increases the number of infected individuals and that increasing the value of δ_{ru} affects only the speed of approach to a steady state.

Figure 3: The number I_{u} of infected individuals in the urban patch as a function of time (in days) for δ_{ur} = 0 and different values of δ_{ru}; see Example 1. In this example, there is movement only from the rural patch to the urban one.

(a) Urban patch with δ_{ru} = 0  (b) Rural patch with δ_{ru} = 0

(c) Urban patch with δ_{ru} = 0.01  (d) Rural patch with δ_{ru} = 0.01

Figure 4: Effect of population movement in just one direction. In this scenario, δ_{ru} ≠ 0 and δ_{ur} = 0, so there is movement only from the rural patch to the urban one; see Example 1. In these plots, we show the number of individuals in each state as a function of time (in days).

Example 2 We now compare the effects of conditions (5a) and (5b) to the standard condition R_{0} < 1 for local asymptotic stability of the disease-free steady state (2). We use the parameter values

\begin{matrix}μ_{u} = 1/(365 · 80) , & ρ_{u} = 0.08 , & γ_{u} = 0.01 , & β_{u} = 3 · 10^{−2} , \\ μ_{r} = 1/(365 · 70) , & ρ_{r} = 0.40 , & γ_{r} = 0.10 , & β_{r} = 2 · 10^{−5} \end{matrix}

and initial conditions

\begin{matrix}S_{u0} = 999 , & I_{u0} = 1 , & S̃_{u0} = 0 , \\ S_{r0} = 300 , & I_{r0} = 0 , & S̃_{r0} = 0 . \end{matrix}

We calculate that R_{0u} ≈ 2.9 and R_{0r} ≈ 2 · 10^{−4}. Therefore, in the absence of population movement, disease persists in the urban patch but dies out in the rural patch; see Figures 5a and 5b. For δ_{ur} = δ_{ru} = 0.05, we calculate that

\frac{δ_{ur}}{γ_{u} + μ} ≈ 5.98, so the conditions in (5) are satisfied. Therefore, the disease-free

steady state is locally asymptotically stable; see Figures 5c and 5d.

In this example, there exists one endemic steady state, for which

(\frac{I_{u}^{∗}}{N_{u}^{∗}} , \frac{I_{r}^{∗}}{N_{r}^{∗}})≈(0.82, 0.76). We study two variations:

(1) When the number of initially infected individuals is sufficiently large, we reach the endemic steady state. See Figure 6afor an illustration with initial conditions S_{u0} = 900 and I_{u0} = 100.

(2) When we increase β_{u} slightly to β_{u} = 5.3 · 10^{−2}, condition (5c) is no longer satisfied, and the system reaches an endemic steady state when we start with only a single infected individual in the urban patch; see Figure 6b.

(a) Urban patch  (b) Rural patch

(c) Urban patch  (d) Rural patch

Figure 5: Illustration of the influence of population movement on disease dynamics in the urban and rural patches. (Top) In the absence of population movement (i.e., δ_{ur} = δ_{ru} = 0), the urban patch reaches an endemic state. (Bottom) For δ_{ur} = δ_{ru} = 0.05, both populations reach a disease-free steady state. In this sense, population movement allows the disease to die out; see Example 2. In these plots, we show the number of individuals in each state as a function of time (in days).

(a) Total number of infected individuals for I_{0u} =100 and β_{u} =3 · 10^{−2}.

(b) Total number of infected individuals for I_{0u} =1 and β_{u} =5.3 · 10^{−2}.

Figure 6: Effect of the conditions in (5) on infected classes of individuals; see Example 2. In these plots, the horizontal axis gives the time (in days).

Example 3 We now give an example with two endemic steady states. We use the parameter values

\begin{matrix}μ_{u} = 1/(365 · 70) , & ρ_{u} = 0.80 , & γ_{u} = 0.15 , & β_{u} = 3 · 10^{−4} , & δ_{ur} = 0.04 , \\ μ_{r} = 1/(365 · 70) , & ρ_{r} = 0.40 , & γ_{r} = 0.10 , & β_{r} = 2 · 10^{−5} , & δ_{ru} = 0.05 . \end{matrix}

We fix the initial conditions

\begin{matrix}N_{u0} = 1000 , & S̃_{u0} = 0 , \\ N_{r0} = 300 , & I_{r0} = 0 , & S̃_{r0} = 0 \end{matrix}

and vary the initial number of infected individuals in the urban patch from 0 to N_{u}. In this case, there are two endemic steady states, with (I_{u}^{∗}/N_{u}^{∗}, I_{r}^{∗}/N_{r}^{∗}) ≈ (0.09, 0.07) and (I_{u}^{∗}/N_{u}^{∗}, I_{r}^{∗}/N_{r}^{∗}) ≈ (0.49, 0.45), which we illustrate in Figure 2b. We show I_{u}^{∗} and I_{r}^{∗} as a function of I_{0u} in Figure 7. The steady state (I_{u}^{∗}/N_{u}^{∗}, I_{r}^{∗}/N_{r}^{∗}) ≈ (0.09, 0.07) is unstable.

Figure 7: The number of infected individuals in the urban and rural patches at an endemic steady state as a function of the initial condition I_{u}(0); see Example 3.

Example 4 In this example, we explore the dependence on R_{0u} and R_{0r} in our model to illustrate the region of local asymptotic stability from (5). We consider the parameter values

\begin{matrix}μ_{u} = 1/(365 · 70) , & ρ_{u} = 0.80 , & γ_{u} = 0.01 , & δ_{ur} = 0.005 , \\ μ_{r} = 1/(365 · 70) , & ρ_{r} = 0.40 , & γ_{r} = 0.05 , & δ_{ru} = 0.010 . \end{matrix}

We fix the initial conditions

\begin{matrix}N_{u0} = 999 , & I_{u0} = 1 , & S̃_{u0} = 0 , \\ N_{r0} = 300 , & I_{r0} = 0 , & S̃_{r0} = 0 \end{matrix}

and vary β_{u} and β_{r}; see our results in Figure 8. We observe that the conditions in (5) precisely describe the disease-free region (in gray). Therefore, R_{0u} or R_{0r} can be slightly larger than 1 in situations with a disease-free steady state. This is not the case when we consider just one patch, so population movement can be beneficial.

Figure 8: Population as a function of the local basic reproduction numbers in the urban R_{0u} (horizontal axis) and rural R_{0r} (vertical axis) environments; see Example 4. The gray area represents the region of local asymptotic stability of both (urban and rural) populations from Lemma 1.

Example 5 In this example, we study how a disease spreads through the two populations as we vary the movement parameters δ_{ur} and δ_{ru} when we start with a single infected individual in one patch and no infected individuals in the other patch. We consider the parameter values

\begin{matrix}μ_{u} = 1/(365 · 70) , & ρ_{u} = 0.80 , & γ_{u} = 0.15 , & β_{u} = 3 · 10^{−1} , \\ μ_{r} = 1/(365 · 70) , & ρ_{r} = 0.40 , & γ_{r} = 0.10 , & β_{r} = 2 · 10^{−3} . \end{matrix}

We fix the initial conditions

\begin{matrix}N_{u0} = 99999 , & I_{u0} = 1 , & S̃_{u0} = 0 , \\ N_{r0} = 30000 , & I_{r0} = 0 , & S̃_{r0} = 0 \end{matrix}

and vary δ_{ur} and δ_{ru}. We show our results in Figure 9.

Example 6 To model different population movement between patches on weekdays and weekends, we now take δ_{ur} and δ_{ru} to be piecewise-constant and periodic. Specifically, we use the functions δ_{ur}(t) and δ_{ru}(t) that we show in Figure 10. In Figure 11, we show the results of our numerical computations using the parameter values

\begin{matrix}μ_{u} = 1/(365 · 70) , & ρ_{u} = 0.80 , & γ_{u} = 0.15 , & β_{u} = 3 · 10^{−4} , \\ μ_{r} = 1/(365 · 70) , & ρ_{r} = 0.40 , & γ_{r} = 0.10 , & β_{r} = 2 · 10^{−5} \end{matrix}

(a) Steady-state fraction (I_{u}^{∗}) of infected individuals in the urban patch.

(b) Steady-state fraction (I_{r}^{∗}) of infected individuals in the rural patch.

Figure 9: The effect of a single infected individual in the urban patch as a function of the movement parameters δ_{ur} (horizontal axis) and δ_{ru} (vertical axis); see Example 5.

and initial conditions

\begin{matrix}N_{u0} = 999 , & I_{u0} = 1 , & S̃_{u0} = 0 , \\ N_{r0} = 300 , & I_{r0} = 0 , & S̃_{r0} = 0 . \end{matrix}

Figure 10: Periodic, piecewise-constant values of the movement parameters (left) δ_{ur} and (right) δ_{ru} as a function of time (in days). We use these functions in Example 6.

Example 7 We now revisit Remark 2, where we stated that we observe the existence of several distinct endemic steady states in our numerical computations. Equation (10) is equivalent to a polynomial equation (in (b_{u}, b_{r})) of degree 9, for which there can exist at most 9 real solutions for (i_{u}^{∗}, i_{r}^{∗}). For the parameter values

Figure 11: The number of infected individuals in the (left) urban (i.e.,I_{u}(t)) and (right) rural (i.e.,I_{r}(t)) populations as a function of time (in days) for the periodic, piecewise-constant movements δ_{ur} and δ_{ru} from Figure 10; see Example 6.

μ_{u} = 1/(365 · 70), ρ_{u} = 1, γ_{u} = 0.15, β_{u} = 3 · 10^{−4}, δ_{ur} = 4 · 10^{−6},
μ_{r} = 1/(365 · 70), ρ_{r} = 2, γ_{r} = 1.00, β_{r} = 2 · 10^{−3}, δ_{ru} = 5 · 10^{−5},

we see that all 9 solutions are in [0, 1]^{2}. They correspond to feasible steady states (i_{u}^{∗}, i_{r}^{∗}); see Figure 12. We see numerically that 4 of these points are locally stable.

Figure 12: We illustrate the existence of multiple stable and unstable steady states of (1) that we obtain by solving (9). We plot the curves G_{1}(i_{u}, i_{r}) = 0 (solid curves) and G_{2}(i_{u}, i_{r}) = 0 (dashed curves). We observe numerically that there are 4 locally asymptotically stable steady states, illustrating that population movement between patches can introduce multiple steady states. We show the locally stable steady states with solid disks and the other steady states with open disks.

5 Conclusions and discussion

The dynamics of spreading diseases are influenced significantly by spatial heterogeneity and population movement. In this paper, we illustrated the importance of incorporating movement into models of disease dynamics using a simple but biologically meaningful model. Specifically, we constructed a two-patch compartmental model that incorporates (1) population movement between urban and rural patches and (2) the possibility of reinfection after recovery.

When there is sufficient movement, there are regions of local asymptotic stability of the disease-free steady-state even when the basic reproduction number R_{0j} > 1 for j ∈ {u, r}. This arises predominantly from the numerous individuals who move between patches.

The exploration of interacting populations plays an important role in the understanding of disease dynamics. Many models of disease spreading focus on a single population [3], but populations do not exist in isolation. Using our two-patch model with urban and rural environments, we illustrated several examples of plausible real-world scenarios in which incorporating population movement yields insightful information about disease spreading and epidemics. We expect that such dynamics will be relevant for studies of disease spreading on networks, such as when many people commute daily between their homes in rural areas and work in urban centers (as is the case in many countries in South and Central America).

Acknowledgements

We thank the Research Center in Pure and Applied Mathematics and the Mathematics Department at Universidad de Costa Rica for their support during the preparation of this manuscript. The authors gratefully acknowledge institutional support for project B8747 from an UCREA grant from the Vice Rectory for Research at Universidad de Costa Rica. We also acknowledge helpful discussions with Profs. Luis Barboza, Carlos Castillo-Chavez, and Esteban Segura.

A Matlab Code for Example 4.1

1 % f i n a l time
2 T = 4∗365;
3 % ODE i n i t i a l c o n d i t i o n s
4 Su0 = 999; Iu0 = 1 ; Ru0 = 0 ;
5 Sr0 = 300; I r 0 = 0 ; Rr0 = 0 ;

6 i n i t _c o n d = [Su0 ; Iu0 ;Ru0; Sr0 ; I r 0 ; Rr0 ] ;
7 % ODE c o e f f i c i e n t s
8 dur = @ ( t ) 0 . 0 5 ;
9 dru = @ ( t ) 0 . 0 5 ;
10 mu = @ ( t ) 1 / ( 3 6 5 ∗ 8 0 );
11 mr = @ ( t ) 1 / ( 3 6 5 ∗ 7 0 );
12 pu = @ ( t ) . 0 8 ;
13 pr = @ ( t ) . 0 4 ;
14 gu = @ ( t ) . 0 1 ;
15 gr = @ ( t ) . 0 1 ;
16 bu = @ ( t ) 0 . 0 3 ;
17 br = @ ( t ) 0 . 0 2 ;

18

19 % d e f i n e ODE
20 f1 = @ ( t ,x) mu( t ) . ∗ ( x ( 1 )+x ( 2 )+x ( 3 ) )−bu ( t ) . ∗x ( 1 ) . ∗x ( 2 )
. / (x ( 1 )+x ( 2 )+x ( 3 ) )−mu( t ) . ∗x ( 1 )+dru ( t ) . ∗x ( 4 )−dur ( t )
. ∗x ( 1 ) ;
21 f2 = @ ( t ,x) bu ( t ) . ∗x ( 1 ) . ∗x ( 2 ) . / (x ( 1 )+x ( 2 )+x ( 3 ) )−(mu( t
)+gu ( t ) ) . ∗x ( 2 ) + pu ( t ) . ∗x ( 2 ) . ∗x ( 3 ) . / (x ( 1 )+x ( 2 )+x
( 3 ) )+dru ( t ) . ∗x ( 5 )−dur ( t ) . ∗x ( 2 ) ;
22 f3 = @ ( t ,x) gu ( t ) . ∗x ( 2 )−pu ( t ) . ∗x ( 2 ) . ∗x ( 3 ) . / (x ( 1 )+x ( 2 )
+x ( 3 ) )−mu( t ) . ∗x ( 3 )+dru ( t ) . ∗x ( 6 )−dur ( t ) . ∗x ( 3 ) ;
23 f4 = @ ( t ,x) mr( t ) . ∗ ( x ( 4 )+x ( 5 )+x ( 6 ) )−br ( t ) . ∗x ( 4 ) . ∗x ( 5 )
. / (x ( 4 )+x ( 5 )+x ( 6 ) )−mr( t ) . ∗x ( 4 )+dur ( t ) . ∗x ( 1 )−dru ( t )
. ∗x ( 4 ) ;
24 f5 = @ ( t ,x) br ( t ) . ∗x ( 4 ) . ∗x ( 5 ) . / (x ( 4 )+x ( 5 )+x ( 6 ) )−(mr( t
)+ gr ( t ) ) . ∗x ( 5 ) + pr ( t ) . ∗x ( 5 ) . ∗x ( 6 ) . / (x ( 4 )+x ( 5 )+x
( 6 ) )+dur ( t ) . ∗x ( 2 )−dru ( t ) . ∗x ( 5 ) ;
25 f6 = @ ( t ,x) gr ( t ) . ∗x ( 5 )−pr ( t ) . ∗x ( 5 ) . ∗x ( 6 ) . / (x ( 4 )+x ( 5 )
+x ( 6 ) )−mr( t ) . ∗x ( 6 )+dur ( t ) . ∗x ( 3 )−dru ( t ) . ∗x ( 6 ) ;
26 f f = @ ( t ,x) [ f1 ( t ,x) ; f2 ( t ,x) ; f3 ( t ,x) ; f4 ( t ,x) ; f5 ( t ,x) ;
f6 ( t ,x) ] ;
27 % o p t i o n s f o r ode s o l v e r
28 o p t i o n s = o d e s e t ( ’ RelTol ’ ,1e−10, ’ S t a t s ’ , ’ o f f ’ , ’ AbsTol
’ ,1e−10) ;
29 % s o l v e ODE
30 [ t o u t , xout ] = ode45 ( f f , [ 0 ,T] , i n i t _c o n d , o p t i o n s ) ;

31

32 f i g u r e

33 % p l o t s o l u t i o n f o r urban p a t c h
34 s u b p l o t ( 3 , 1 , 1 )
35 p l o t ( t o u t , xout ( : , 1 ) , ’ : k ’ ) , hold on
36 p l o t ( t o u t , xout ( : , 2 ) , ’−k ’ )
37 p l o t ( t o u t , xout ( : , 3 ) , ’−k ’ )
38 l e g e n d ({ ’ $S_u$ ’ , ’ $I_u$ ’ , ’$ \w i d e t i l d e {S}_u$ ’ } , ’
I n t e r p r e t e r ’ , ’ Latex ’ )
39 t i t l e ( ’ Urban ’ )
40 a x i s ( [ 0 T 0 1 0 0 0 ] )
41 g r i d on
42 % p l o t s o l u t i o n f o r r u r a l p a t c h
43 s u b p l o t ( 3 , 1 , 2 )
44 p l o t ( t o u t , xout ( : , 4 ) , ’ : k ’ ) , hold on
45 p l o t ( t o u t , xout ( : , 5 ) , ’−k ’ )
46 p l o t ( t o u t , xout ( : , 6 ) , ’−k ’ )
47 l e g e n d ({ ’ $S_r$ ’ , ’ $ I _r $ ’ , ’$ \w i d e t i l d e {S}_r$ ’ } , ’
I n t e r p r e t e r ’ , ’ Latex ’ )
48 t i t l e ( ’ Rural ’ )
49 a x i s ( [ 0 T 0 1 0 0 0 ] )
50 g r i d on

References

[1] L. Alvarado, Costa Rica once again under malaria alert, The Costa Rica Star, 2018. Available at https:// news. co.cr/

       costa-rica-once-again-under-malaria-alert/      73681 ,

accessed 25/04/2019.

[2] D. Bichara, C. Castillo-Chavez, Vector-borne diseases models with residence times — A Lagrangian perspective, Math. Biosci. 281(2016),

       128–138. doi: 10.1016/ j.mbs. 2016. 09.006

[3] F. Brauer, C. Castillo-Chavez, Mathematical Models in Population Biology and Epidemiology, 2nd edition, Springer-Verlag, Providence RI, USA, 2012.

[4] Center for Disease Control and Prevention, Severe acute respiratory system

       (SARS), 2019. Available at https:// www. cdc. gov/ sars/ index.
       html

[5] Center for Disease Control and Prevention, Ebola (Ebola virus dis-

ease), 2019. Available at https:// www. cdc. gov/ vhf/ ebola/

history/ 2014-2016-outbreak/   index. html

[6] Center for Disease Control and Prevention, Measles (Rubeola), 2019.

Available at https:// www. cdc. gov/ measles/ index. html

[7] G. Chowell, P.W. Fenimore, M.A. Castillo-Garsow C. Castillo-Chavez, SARS outbreaks in Ontario, Hong Kong and Singapore: The role of diagnosis and isolation as a control mechanism, J. Theor. Biol. 224(2003),

no. 1, 1–8. doi: 10.1016/ S0022-5193(03)00228-5

[8] M.P. Coffee, G.P. Garnett, M. Mlilo, H.A.C.M. Voeten, S. Chandiwana, S. Gregson, Patterns of movement and risk of HIV Infection in rural Zimbabwe, J. Infect. Dis. 191(2005), no. 1, S159–S167. doi: 10.1086/

425270

[9] J.M. Crutcher, S.L. Hoffman, Malaria, in: S. Baron (Ed.) Medical Microbiology, 4th edition, University of Texas Medical Branch at Galveston, Galveston TX, 1996, ch. 83. Available in: https:// www. ncbi. nlm.

nih. gov/ books/ NBK8584/

[10] R. DeVore, A. Ron, Approximation of functions, Proc. Sympos. Appl. Math. 36(1986), 34–56.

[11] Z. Feng, J. X. Velasco-Hernández, Competitive exclusion in a vector–host model for the dengue fever, J. Math. Biol. 35(1997), no. 5, 523–544. doi:

10.1007/ s002850050064

[12] Z. Feng, C. Castillo-Chavez, A. F. Capurro, A model for tuberculosis with exogenous reinfection, Theor. Popul. Biol. 57(2000), no. 3, 235–247. doi:

10.1006/ tpbi. 2000. 1451

[13] H. Frankowska, The Poincaré–Miranda theorem and viability condition, J. Math. Anal. Appl. 463(2018), no. 2, 832–837. doi: 10.1016/ j.jmaa.

2018. 03.047

[14] J. Gjorgjieva, K. Smith, G. Chowell, F. Sanchez, J. Snyder, C. Castillo-Chavez, The role of vaccination in the control of SARS, Math. Biosci. Eng. 2(2005), no. 4, 753–769. doi: 10.1006/ tpbi. 2000. 1451 [15] J.R. Glynn, J. Murray, A. Bester, G. Nelson, S. Shearer, P. Sonnenberg, Ef-

fects of duration of HIV infection and secondary tuberculosis transmission on tuberculosis incidence in the South African gold mines, AIDS 22(2008),

no. 14, 1859–1867. doi: 10.1097/ QAD. 0b013e3283097cfa

[16] J.L. Grun, W.P. Weidanz, Antibody-independent immunity to reinfection malaria in B-cell-deficient mice, Infect. and Immun. 41(1983), no. 3, 1197– 1204.

[17] W. Kulpa, The Poincaré–Miranda theorem, Amer. Math. Month. 104(1997), no. 6, 545–550. doi: 10.2307/ 2975081

[18] S. Lee, C. Castillo-Chavez, The role of residence times in two-patch dengue transmission dynamics and optimal strategies, J. Theor. Biol. 374(2015),

no. 7, 152–164. doi: 10.1016/ j.jtbi. 2015. 03.005

[19] C.A. Manore, K.S. Hickmann, S. Xu, H.J. Wearing, J.M. Hyman, Comparing dengue and chikungunya emergence and endemic transmission in A. aegypti and A. albopictus, J. Theor. Biol. 356(2014), 174–191. doi: 10.

1016/ j.jtbi. 2014. 04.033

[20] P. Martens, L. Hall, Malaria on the move: Human population movement and malaria transmission, Emerg. Infect. Diseas. 6(2000), no. 2, 103–109.

doi: 10.3201/ eid0602. 000202

[21] D. Murillo, S. Holechek, A. Murillo, F. Sanchez, C. Castillo-Chavez, Vertical transmission in a two-strain model of dengue fever, Lett. Biomath.

1(2014), no. 2, 249–271. doi: 10.1080/ 23737867. 2014. 11414484

[22] R. Pastor-Satorras, C. Castellano, P. van Mieghem, A. Vespignani, Epidemic processes in complex networks, Rev. Mod. Phys. 87(2015), no. 3, 925–979.

[23] F. Sanchez, M. Engman, L.C. Harrington, C. Castillo-Chavez, Models for dengue transmission and control, in: A.B. Gumel, C. Castillo-Chavez, R.E. Mickens & D.P. Clemence (Eds.) Mathematical Studies on Human Disease Dynamics. Emerging Paradigms and Challenges, Contemp. Math. 410, Amer. Math. Soc., Providence RI, USA, 2006, pp. 311–326. doi: 10.

1090/ conm/ 410/ 07734

[24] F. Sanchez, X. Wang, C. Castillo-Chavez, D. Gorman, P.J. Gruenewald, Drinking as an epidemic — A simple mathematical model with recovery and relapse, in: K.A. Witkiewitz & G.A. Marlatt (Eds.) Therapist’s Guide to Evidence-Based Relapse Prevention, Academic Press, Cambridge, MA, USA, 2007, pp. 353–368. doi: 10.1016/ B978-012369429-4/

50046-X

[25] F. Sanchez, D. Murillo, C. Castillo-Chavez, Change in host behavior and its impact on the transmission dynamics of dengue, in: R.P. Mondaini (Ed.) International Symposium on Mathematical and Computational Biology, BIOMAT 2011 (Santiago, Chile), 2012, pp. 191–203. doi: https://

doi. org/ 10.1142/ 9789814397711_0013

[26] F. Sanchez, J. G. Calvo, E. Segura, Z. Feng, A partial differential equation model with age-structure and nonlinear recidivism: Conditions for a backward bifurcation and a general numerical implementation, Computers and Mathematics with Applications 78(2018), no. 12, 3916–3930. doi: 10.

1016/ j.camwa. 2019. 06.021

[27] B. Song, M. Castillo-Garsow, K.R. Rios-Soto, M. Mejran, L. Henso, C. Castillo-Chavez, Raves, clubs and ecstasy: The impact of peer pressure, Math. Biosci. Eng. 3(2006), no. 1, 249–266. doi: 10.3934/ mbe. 2006.

3.249

[28] B. Song, W. Du, J. Lou, Different types of backward bifurcations due to density-dependent treatments, Math. Biosci. Eng. 10(2013), no. 5–6,

1651–1668. doi: 10.3934/ mbe. 2013. 10.1651

[29] K. Szymańska-Dȩbowska, On a generalization of the Miranda Theorem and its application to boundary value problems, J. Diff. Equ. 258(2015),

no. 8, 2686–2700. doi: 10.1016/ j.jde. 2014. 12.022

[30] A.J. Treno, P.J. Gruenewald, L.G. Remer, F. Johnson, E.A. LaScala, Examining multi-level relationships between bars, hostility and aggression: Social selection and social influence, Addiction 103(2007), no. 1, 66–77.

doi: 10.1111/ j.1360-0443. 2007. 02039. x

[31] M. Turzánski, The Bolzano–Poicaré–Miranda theorem — Discrete version, Topol. Appl. 159(2012), no. 13, 3130–3135. doi: 10.1016/ j.topol.

2012. 05.026