# Relating Eulerian and Lagrangian spatial models for vector-host disease
dynamics through a fundamental matrix

Esteban Vargas Bernal1111Corresponding author. E-mail:
<EMAIL_ADDRESS>Omar Saucedo2, Joseph Hua Tien1
1Department of Mathematics, The Ohio State University.
2Department of Mathematics, Virginia Tech.

###### Abstract

We explore the relationship between Eulerian and Lagrangian approaches for
modeling movement in vector-borne diseases for discrete space. In the Eulerian
approach we account for the movement of hosts explicitly through movement
rates captured by a graph Laplacian matrix $L$. In the Lagrangian approach we
only account for the proportion of time that individuals spend in foreign
patches through a mixing matrix $P$. We establish a relationship between an
Eulerian model and a Lagrangian model for the hosts in terms of the matrices
$L$ and $P$. We say that the two modeling frameworks are consistent if for a
given matrix $P$, the matrix $L$ can be chosen so that the residence times of
the matrix $P$ and the matrix $L$ match. We find a sufficient condition for
consistency, and examine disease quantities such as the final outbreak size
and basic reproduction number in both the consistent and inconsistent cases.
In the special case of a two-patch model, we observe how similar values for
the basic reproduction number and final outbreak size can occur even in the
inconsistent case. However, there are scenarios where the final sizes in both
approaches can significantly differ by means of the relationship we propose.

Keywords— Vector-borne disease, Movement, Spatial models, Lagrangian approach,
Eulerian approach, Fundamental matrix.

##  1 Introduction

Understanding how connectivity between different spatial locations affects
vector-borne disease dynamics is a fundamental issue in disease ecology and
public health. In particular, wide variation in disease transmission between
locales is commonplace, reflecting heterogeneity in breeding site
availability, geography, climate, availability of bed nets and window screens,
demography, and many other factors. How this spatial heterogeneity interacts
with connectivity through host and vector movement to inform disease dynamics
is not obvious. For example, empirical studies have shown how a vector-borne
disease may persist in cities where mosquito abundance is low or zero (for
example for malaria in [38]). Some authors have shown that this persistence
can be explained by host movement. Stoddard et al. [49] used a conceptual
model to show that even when the vector density is low, the risk of acquiring
the disease may be high due to movement. Similarly, Cosner et al. [17]
constructed a spatial vector-borne disease model to study how human movement
can affect transmission, and discovered that human movement between
heterogeneous locations was sufficient to sustain disease persistence. In
other situations, movement can lead to disease extinction, even in areas with
high local transmission [50]. Movement patterns also affect the spatial spread
of vector-borne diseases such as Lyme disease [21], West Nile Virus [35],
Dengue [20], Zika [39, 61], and Malaria [50, 58]. Both connectivity and local
conditions for disease transmission are important considerations when
designing disease surveillance and control efforts [60].

Many different approaches have been taken for modeling spatial vector-borne
disease dynamics, including PDEs [34], ODEs [1, 40, 44], stochastic models
[30, 57], and agent-based models [12, 29]. A widely-used building block for
modeling mosquito-borne disease is the Ross-MacDonald model [41, 43]. See [48]
for derivation of the Ross-Macdonald system that is considered here. This
paper concerns extensions of the Ross-MacDonald framework to include multiple
discrete spatial locations. These locations might correspond to villages,
cities, health districts, or the like, with linkages between them through
movement of host and vector. Note that the connectivity patterns for host and
vector may be different, for example reflecting different movement scales of
each. Regarding mosquito-borne diseases, Service [45] provides a review of the
types of mosquitoes movement (long/short dispersal), which may vary
significantly among different species. Empirical studies using different
capture methods (bed net catches, exit trap catches, oviposition traps) allow
us to have an idea of the spatial scale of mosquito movement [25, 52]. Given
the importance that host and vector movement may have for the spread of a
vector-borne disease, in this paper we consider both movements.

We now encounter a dichotomy in the modeling approaches: whether to treat
individuals (host or vector) as residents of a particular patch and commuting
to the others, versus migration between patches without a fixed sense of home.
Following the terminology of [16, 17], we will refer to the former approach as
Lagrangian, and the latter as Eulerian. This terminology stems from
similarities to the Lagrangian and Eulerian modeling approaches for the
description of fluid motion in fluid mechanics [31]. For more on Eulerian and
Lagrangian approaches for modeling movement, see [22, 23]. Each of these
approaches has its strengths. Lagrangian models are often a natural choice on
small spatial scales and in settings where individuals have a sense of home.
Eulerian models may be suitable on large spatial scales, and in migratory
settings where the origin of the individuals is less important than their
current location. This may include situations where we are interested in
introduction, reintroduction, or global spread of disease, as discussed in
[49].

Eulerian and Lagrangian approaches are widely-used to model vector-borne
disease dynamics on discrete space. See [19, 24, 42, 44] for some studies that
have used the Lagrangian approach. For example, in [19, 24], each study used a
Lagrangian framework to examine how heterogeneity in the distribution of
_Anopheles_ mosquitoes could lead to a larger reproduction number compared to
the setting when the mixing is homogeneous. See also [3, 21, 26, 35] for
applications of the Eulerian approach. For example, in [35] an Eulerian model
was used to study how the long range dispersal of birds may explain
discontinuities in the spread of West Nile Virus. Combinations of the Eulerian
and Lagrangian approaches are taken in [4, 5, 6, 27].

What are the functional implications of using one approach versus the other?
Do the Eulerian and Lagrangian approaches yield similar results regarding
important disease quantities such as the basic reproduction number and
outbreak size? Eulerian models tend to be analytically tractable, for example
allowing establishment of asymptotic disease profiles [3], global stability
[46], and application of techniques from spectral graph theory to estimate the
basic reproduction number [53]. We would like to know whether one can safely
use an Eulerian approach to model vector-borne disease in a setting where
Lagrangian data are available (e.g., perhaps in terms of proportions of time
spent in different locations) or where spatial scales are small and commuting
(vs. migration) is typical. Comparing the Eulerian and Lagrangian approaches
to modeling movement, and the functional implications of using one approach
versus the other in vector-borne disease models, is the focus of this paper.

Consider vector-host disease dynamics on $n$ spatial locations. We define two
families of models for this setting, where both use an Eulerian framework for
modeling vector movement. In the first model we consider a Lagrangian approach
for host movement, and we refer to this as the Lagrangian model [17, 36]. In
the second model we consider an Eulerian framework for host movement, and we
refer to this model as the Eulerian model [17, 36]. The Lagrangian character
of the first model is captured by a mixing matrix $P=(p_{ij})_{i,j\leq n}$,
where $p_{ij}$ is the fraction of time that a resident of patch $j$ spends in
patch $i$. The Eulerian character of the second model is captured by an
adjacency matrix $A=(m_{ij})_{i,j\leq n}$, with $m_{ij}$ being the per capita
movement rate from patch $j$ to patch $i$. We will work extensively with the
unnormalized (‘combinatorial’) graph Laplacian $L=W-A$, where $W$ diagonal
with $W_{ii}=\sum_{k=1}^{n}A_{ki}$. The graph Laplacian $L$ is a basic object
in graph theory that conveys a great deal of structural information about the
network associated with $A$, including the number of connected components,
spanning trees, community structure, and more [15, 37, 56].

The main goal of this paper is studying the relationship between the
Lagrangian and Eulerian models. Specifically, we obtain a relationship between
the matrices $L$ and $P$ through a fundamental matrix that captures the
expected time that an individual from one location spends in another. We give
criteria for when the Eulerian and Lagrangian frameworks are consistent,
meaning that the two frameworks can exactly match in terms of this fundamental
matrix, and consider the functional implications for disease dynamics in both
the consistent and inconsistent settings. These results can serve as a guide
for when one framework can be substituted for the other.

The following is the distribution of the content of this paper. In Section 2
we define the Lagrangian and Eulerian systems that we study, and we give the
disease-free equilibria and the basic reproduction numbers
$\mathcal{R}_{0}^{\operatorname{lag}}$, $\mathcal{R}_{0}^{\operatorname{eul}}$
for both systems. In Section 3 we relate the Lagrangian and Eulerian systems
through a minimization problem involving the matrices $P$ and $L$. We say that
the two systems are consistent if this minimization problem can be solved
exactly. Following some preliminaries in Section 3.1, we formulate this
minimization problem in Section 3.2. In Section 3.3 we look at the consistent
scenario. We first provide a sufficient condition for the consistency of the
proposed relationship in Section 3.3.1, and then, we give an example of a
consistent relationship in Section 3.3.2. In Section 3.4 we give examples
where the relationship is inconsistent. In Section 4 we study the relationship
between the Lagrangian and Eulerian approaches for a simple network consisting
of two patches. In Section 4.1 we compare the final outbreak sizes and basic
reproduction numbers of both systems under an inconsistent example. In Section
4.2 we compare the basic reproduction numbers of consistent examples when we
vary the entries of the mixing matrix. In Section 5 we explore some examples
of matrices $P$ from empirical and hypothetical data. In Section 6 we present
the main conclusions of the paper. Finally, we give the details of some
results of the previous sections in an appendix in Section 7.

##  2 Modeling frameworks

The basic building block for the modeling frameworks considered in this paper
is the Ross-MacDonald vector-host model, as considered by [48]. For a single
spatial location, the model equations are:

$\begin{array}[]{ccl}\dot{S}&=&\Lambda-\beta\frac{S}{N}I_{v}-\mu S\\\
\dot{I}&=&\beta\frac{S}{N}I_{v}-(\gamma+\mu)I\\\ \dot{R}&=&\gamma I-\mu R\\\
\dot{S}_{v}&=&\Lambda_{v}-\beta_{v}\frac{I}{N}S_{v}-\mu_{v}S_{v}\\\
\dot{I}_{v}&=&\beta_{v}\frac{I}{N}S_{v}-\mu_{v}I_{v}\,.\end{array}$ (1)

System (1) takes a Susceptible-Infectious-Recovered (SIR) framework for host
and Susceptible-Infectious (SI) framework for vector. $S$ represents the
number of susceptible hosts, $I$ represents the number of infectious hosts,
and $R$ represents the number of recovered hosts. Similarly, $S_{v}$ denotes
the number of susceptible vectors and $I_{v}$ denotes the number of infected
vectors. The total host population is denoted by $N$. $\Lambda$ and
$\Lambda_{v}$ are the constant recruitment rates for host and vector,
respectively. The transmission rate from vector to host is $\beta_{i}$ and the
transmission rate from host to vector is $\beta_{v,i}$. The parameters $\mu$
and $\mu_{v}$ correspond to the mortality rates for host and vectors. Many
modifications to this framework are possible, including incubation periods,
seasonal forcing, and much more (see [14] for an example of seasonal forcing
and [47] for a more general review). We consider here the very simple system
(1) in order to focus on the impact of connectivity between different spatial
locations.

Consider $n$ distinct spatial locations, each with local Ross-MacDonald
dynamics as in (1) with patch-specific parameters. We will consider two
modeling frameworks that differ in how the spatial locations are coupled
through host movement: one using a Lagrangian approach, and the other an
Eulerian approach [17, 36].

In the Lagrangian approach, coupling is a mixing matrix $P=(p_{ij})_{i,j\leq
n}$, where $p_{ij}$ is the proportion of time a resident of patch $j$ spends
in patch $i$. By contrast, in the Eulerian approach, coupling is via the
adjacency matrices $M^{\operatorname{X}}=(m_{ij}^{\operatorname{X}})_{i,j\leq
n}$ of weighted, directed graphs, where $m_{ij}^{\operatorname{X}}$ is the per
capita migration rate of hosts in state
$\operatorname{X}\in\\{\operatorname{S},\operatorname{I},\operatorname{R}\\}$.

In both modeling frameworks, vector movement is modeled using an Eulerian
framework, with $M^{v}$ the weighted adjacency matrix for vector movement.
Thus the two frameworks considered are Lagrangian (host) / Eulerian (vector),
and Eulerian (host) / Eulerian (vector). For brevity we will refer to these
frameworks as simply Lagrangian in the former, and Eulerian in the latter.

In the ensuing analysis we will make extensive use of the (unnormalized) graph
Laplacian $L$ [2]. Let $A$ be the adjacency matrix for a (weighted, directed)
graph $G$, with $A_{ij}$ the weight of the edge from $j$ to $i$. Let $W$ be
the diagonal out-degree matrix with $W_{ii}=\sum_{i=1}^{n}A_{ij}$. Then the
graph Laplacian is defined as

$L=W-A\,.$ (2)

The graph Laplacian is a fundamental quantity in graph theory that conveys
structural information about $G$, including the number of connected
components, spanning trees, community structure, and more [2, 15, 37]. In the
context of infectious disease dynamics, the Laplacian arises naturally in the
calculation of ${\cal R}_{0}$ for migration models [53].

All model parameters throughout are assumed to be non-negative. The adjacency
matrices $M^{\operatorname{X}}$, for
$\operatorname{X}\in{\\{\operatorname{S},\operatorname{I},\operatorname{R}\\}}$,
and $M^{v}$ are assumed to have zero diagonal (the corresponding graphs for
Eulerian movement have no self-edges). By definition, the columns of the
mixing matrix $P$ sum to one.

###  2.1 Lagrangian for host, Eulerian for vector

The Lagrangian approach can be viewed as a multigroup model with groups
corresponding to hosts that are residents of the different spatial locations.
The number of susceptible, infectious and recovered individuals that are
residents of patch $i$ are denoted by $S_{i}$, $I_{i}$ and $R_{i}$
respectively, and the number of susceptible and infectious vectors in patch
$i$ are denoted by $S_{v,i}$ and $I_{v,i}$ respectively. The proportion
$p_{ij}$ is the ratio between the time spent by a resident of $j$ in patch $i$
to the whole time spent by a resident of $j$ in all the visited patches. All
of these proportions are collected by the mixing matrix $P=(p_{ij})_{i,j\leq
n}$. The movement rate of a vector from patch $j$ to patch $i$ is given by
$m_{ij}^{v}$, and all these rates are collected in the adjacency matrix
$M^{v}=(m_{ij}^{v})_{i,j\leq n}$. The transmission rates $\beta_{i}$ (from
vectors to hosts) and $\beta_{v,i}$ (from hosts to vectors) are intrinsic to
the corresponding patch and determine the transmission rates in (3) by
averaging according to the mixing rates. Namely, the transmission rate for
host residents of patch $i$ is
$\sum_{j=1}^{n}\beta_{j}p_{ji}\left(S_{i}/N_{i}\right)I_{v,j}$ and the
transmission rate for vectors in patch $i$ is
$\beta_{v,i}\left[\left(\sum_{j=1}^{n}p_{ij}I_{j}\right)\middle/\left(\sum_{j=1}^{n}p_{ij}N_{j}\right)\right]S_{v,i}$.
The transmission rate
$\sum_{j=1}^{n}\beta_{j}p_{ji}\left(S_{i}/N_{i}\right)I_{v,j}$ has been
considered in Lagrangian models such as in [17, 19, 42, 44]. There exist other
models that incorporate more complex transmission rates for the host. For
example, [11] considered $\sum_{k=1}^{n}p_{jk}N_{k}$ instead of $N_{i}$ in the
denominator of the transmission rate for host individuals under other
assumptions. However, the transmission rate considered here is justified by
its inclusion in other studies and its analytical tractability. The
transmission rate
$\beta_{v,i}\left[\left(\sum_{j=1}^{n}p_{ij}I_{j}\right)\middle/\left(\sum_{j=1}^{n}p_{ij}N_{j}\right)\right]S_{v,i}$
for vectors has been adopted in studies such as [44]. Moreover, the whole
Lagrangian system that we consider was also studied in [44].

The recruitment, mortality and recovery rates for host residents of patch $i$
are $\Lambda_{i}^{\operatorname{lag}},\mu_{i}^{\operatorname{lag}},$ and
$\gamma_{i}^{\operatorname{lag}}$, respectively. The recruitment and mortality
rates for vectors in patch $i$ are $\Lambda_{v,i}$ and $\mu_{v,i},$
respectively. Let
$\delta_{i}^{\operatorname{lag}}=\mu_{i}^{\operatorname{lag}}+\gamma_{i}^{\operatorname{lag}}$
and $\delta_{v,i}=\mu_{v,i}$ denote the host and vector removal rates,
respectively. Equations for the Lagrangian framework are given in system (3):

$\begin{array}[]{ccl}\dot{S_{i}}&=&\Lambda_{i}^{\operatorname{lag}}-\sum_{j=1}^{n}\beta_{j}p_{ji}\frac{S_{i}}{N_{i}}I_{v,j}-\mu_{i}^{\operatorname{lag}}S_{i}\\\
\dot{I_{i}}&=&\sum_{j=1}^{n}\beta_{j}p_{ji}\frac{S_{i}}{N_{i}}I_{v,j}-\left(\gamma_{i}^{\operatorname{lag}}+\mu_{i}^{\operatorname{lag}}\right)I_{i}\\\
\dot{R_{i}}&=&\gamma_{i}^{\operatorname{lag}}I_{i}-\mu_{i}^{\operatorname{lag}}R_{i}\\\
\dot{S}_{v,i}&=&\Lambda_{v,i}-\beta_{v,i}\frac{\sum_{j=1}^{n}p_{ij}I_{j}}{\sum_{j=1}^{n}p_{ij}N_{j}}S_{v,i}+\sum_{j=1}^{n}m_{ij}^{v}S_{v,j}-\sum_{j=1}^{n}m_{ji}^{v}S_{v,i}-\mu_{v,i}S_{v,i}\\\
\dot{I}_{v,i}&=&\beta_{v,i}\frac{\sum_{j=1}^{n}p_{ij}I_{j}}{\sum_{j=1}^{n}p_{ij}N_{j}}S_{v,i}+\sum_{j=1}^{n}m_{ij}^{v}I_{v,j}-\sum_{j=1}^{n}m_{ji}^{v}I_{v,i}-\mu_{v,i}I_{v,i}\,,\end{array}$
(3)

for $i=1,\dots,n$.

We assume the following throughout:

  * A1:

The adjacency matrix for vector movement has zeros on the diagonal (i.e.
$m^{v}_{ii}=0$ for all $i$).

  * A2:

The parameters $\beta_{i},\beta_{v,i}$ are non-negative.

  * A3:

The parameters
$\Lambda_{i}^{\operatorname{lag}},\Lambda_{v,i},\gamma_{i}^{\operatorname{lag}}$,
$\mu_{i}^{\operatorname{lag}}$ and $\mu_{v,i}$ are all positive.

  * A4:

The mixing matrix $P$ is non-negative, with $\sum_{i=1}^{n}p_{ij}=1$ for all
$j$.

Table 1 contains the parameters used for system (3).

Table 1: Parameters for systems (3) and (8).  Parameter | Meaning | Units
---|---|---
$\Lambda_{i}^{\operatorname{lag}},\Lambda_{i}^{\operatorname{eul}}$ | Recruitment rate of susceptible host in patch $i$ | Hosts $\times$ Days-1
$\Lambda_{v,i}$ | Recruitment rate of susceptible vectors in patch $i$ | Vectors $\times$ Days-1
$\gamma_{i}^{\operatorname{eul}}$ | Per capita recovery rate of hosts in patch $i$. | Days-1
$\mu_{i}^{\operatorname{eul}}$ | Per capita mortality rate of hosts in patch $i$ | Days-1
$\delta_{i}^{\operatorname{eul}}$ | Per capita removal rate of infectious hosts in patch $i$ | Days-1
$\gamma_{i}^{\operatorname{lag}}$ | Per capita recovery rate of hosts from patch $i$ | Days-1
$\mu_{i}^{\operatorname{lag}}$ | Per capita mortality rate of hosts from patch $i$ | Days-1
$\delta_{i}^{\operatorname{lag}}$ | Per capita removal rate of infectious hosts from patch $i$ | Days-1
$p_{ji}$ | Proportion of time that a host from patch $i$ spends in patch $j$ | Dimensionless
$\beta_{i}$ | Transmission rate to hosts per vector in patch $i$ | Hosts $\times$ Days${}^{-1}\times$ Vectors-1
$\beta_{v,i}$ | Per capita transmission rate to vectors in patch $i$ | Days-1
$m_{ji}^{\operatorname{X}}$ | Per capita movement rate of hosts in state $\operatorname{X}$ | Days-1
| from patch $i$ to $j$, for $\operatorname{X}\in\\{\operatorname{S},\operatorname{I},\operatorname{R}\\}$ |
$m_{ji}^{v}$ | Per capita movement rate of vectors from patch $i$ to $j$ | Days-1

Let $G_{v}\mathrel{\mathop{\mathchar 58\relax}}=L_{v}+D_{\delta_{v}}$, where
$L_{v}$ is the graph Laplacian corresponding to the adjacency matrix $M^{v}$
that captures the vector movement, $D_{\delta_{v}}\mathrel{\mathop{\mathchar
58\relax}}=\operatorname{diag}\left\\{\delta_{v,i}\right\\}$ and
$D_{\mu}^{\operatorname{lag}}\mathrel{\mathop{\mathchar
58\relax}}=\operatorname{diag}\left\\{\mu_{i}^{\operatorname{lag}}\right\\}$.
Note that $G_{v}$ has the Z-sign pattern [10], and under assumption A2 has
positive column sums. Thus $G_{v}$ is a non-singular $M$-matrix [10]. The
disease-free equilibrium (DFE) of the susceptible compartments of model (3) is
then

$\begin{array}[]{rcl}\left(S^{\operatorname{lag}}\right)^{*}=\left(N^{\operatorname{lag}}\right)^{*}&=&\left(D_{\mu}^{\operatorname{lag}}\right)^{-1}\Lambda^{\operatorname{lag}}\,,\\\
\left(S_{v}^{\operatorname{lag}}\right)^{*}=\left(N_{v}^{\operatorname{lag}}\right)^{*}&=&G_{v}^{-1}\Lambda_{v}\,,\end{array}$
(4)

where
$\left(S^{\operatorname{lag}}\right)^{*},\left(N^{\operatorname{lag}}\right)^{*}$,
$\left(S_{v}^{\operatorname{lag}}\right)^{*},\left(N_{v}^{\operatorname{lag}}\right)^{*}$,
$\Lambda_{v}$ and $\Lambda^{\operatorname{lag}}$ are column vectors with
components
$\left(S_{i}^{\operatorname{lag}}\right)^{*},\left(N_{i}^{\operatorname{lag}}\right)^{*}$,
$\left(S_{v,i}^{\operatorname{lag}}\right)^{*},\left(N_{v,i}^{\operatorname{lag}}\right)^{*}$,
$\Lambda_{v,i}$ and $\Lambda_{i}^{\operatorname{lag}},$ respectively. (The
superscript $*$ indicates evaluation at the DFE.)

Consider the basic reproduction number ${\cal R}_{0}^{\operatorname{lag}}$ for
system (3), computed using the next generation matrix approach [54]. Then
$\left({\cal
R}_{0}^{\operatorname{lag}}\right)^{2}=\rho\left(\left(F^{\operatorname{lag}}\right)\left(V^{\operatorname{lag}}\right)^{-1}\right)$,
where $\rho$ denotes the spectral radius, and $F^{\operatorname{lag}}$ and
$V^{\operatorname{lag}}$ denote the fecundity and transfer matrices for system
(3). Let $D_{\beta}\mathrel{\mathop{\mathchar
58\relax}}=\operatorname{diag}\\{\beta_{i}\\}$,
$\hat{N}^{*}\mathrel{\mathop{\mathchar
58\relax}}=P(N^{\operatorname{lag}})^{*}$,
$D_{\beta_{v}}^{\operatorname{lag}}\mathrel{\mathop{\mathchar
58\relax}}=\operatorname{diag}\left\\{\beta_{v,i}\left(N_{v,i}^{\operatorname{lag}}\right)^{*}/\hat{N_{i}}^{*}\right\\}$
and $D_{\delta}^{\operatorname{lag}}\mathrel{\mathop{\mathchar
58\relax}}=\operatorname{diag}\left\\{\delta_{i}^{\operatorname{lag}}\right\\}$.
The resulting fecundity and transfer matrices are

$F^{\operatorname{lag}}=\begin{pmatrix}0&P^{\operatorname{T}}D_{\beta}\\\
D_{\beta_{v}}^{\operatorname{lag}}P&0\end{pmatrix}$ (5)

and

$V^{\operatorname{lag}}=\begin{pmatrix}D_{\delta}^{\operatorname{lag}}&0\\\
0&G_{v}\end{pmatrix}\,,$ (6)

with basic reproduction number

$\left(\mathcal{R}_{0}^{\operatorname{lag}}\right)^{2}=\rho\left(\left(F^{\operatorname{lag}}\right)\left(V^{\operatorname{lag}}\right)^{-1}\right)=\rho\left(P^{\operatorname{T}}D_{\beta}G_{v}^{-1}D_{\beta_{v}}^{\operatorname{lag}}P\left(D_{\delta}^{\operatorname{lag}}\right)^{-1}\right)\,.$
(7)

Details are provided in Appendix 7.1. The next generation matrix is the
product of two terms: $P^{\operatorname{T}}D_{\beta}G_{v}^{-1}$ corresponding
to secondary host infections created by infectious vectors, and
$D_{\beta_{v}}^{\operatorname{lag}}P\left(D_{\delta}^{\operatorname{lag}}\right)^{-1}$
corresponding to secondary vector infections created by infectious hosts. Note
that for the Lagrangian model, host movement influences the next generation
matrix via the mixing matrix $P$ appearing in the fecundity matrix
$F^{\operatorname{lag}}$.

###  2.2 Eulerian for host, Eulerian for vector

Consider now the setting where spatial locations are coupled via migration of
both host and vector. The abundances of susceptible, infectious and recovered
hosts in patch $i$ are $S_{i}$, $I_{i}$ and $R_{i}$, and the number of
susceptible and infectious vectors are $S_{v,i}$ and $I_{v,i}$ respectively.
The per capita movement rate of hosts in state $\operatorname{X}$ from patch
$j$ to patch $i$ is $m_{ij}^{\operatorname{X}}$, for
$\operatorname{X}\in\\{\operatorname{S},\operatorname{I},\operatorname{R}\\}$,
and these rates are recorded in the host movement adjacency matrix
$M^{\operatorname{X}}=\left(m_{ij}^{\operatorname{X}}\right)_{i,j\leq n}$. The
assumption that the movement rate between two patches depends on the state
$\operatorname{X}\in\\{\operatorname{S},\operatorname{I},\operatorname{R}\\}$
has been considered in studies such as [26]. Similarly, the per capita
movement rate of vectors from patch $j$ to patch $i$ is $m_{ij}^{v}$ and these
rates are collected in the vector movement adjacency matrix
$M^{v}=\left(m_{ij}^{v}\right)_{i,j\leq n}$. As before, the host and vector
transmission rates for patch $i$ are $\beta_{i}$ and $\beta_{v,i}$. We treat
these rates as intrinsic to the patch, and thus take them to be the same as in
model (3). The recruitment, mortality and recovery rate for hosts in patch $i$
are $\Lambda_{i}^{\operatorname{eul}},\mu_{i}^{\operatorname{eul}}$ and
$\gamma_{i}^{\operatorname{eul}}$ respectively. The recruitment and mortality
rate for vectors in patch $i$ are $\Lambda_{v,i}$ and $\mu_{v,i}$
respectively, taken as the same as in model (3). We also define the removal
rates
$\delta_{i}^{\operatorname{eul}}=\mu_{i}^{\operatorname{eul}}+\gamma_{i}^{\operatorname{eul}}$
and $\delta_{v,i}=\mu_{v,i}$. This comprises a modeling framework where an
Eulerian approach is used to model both host and vector movement. The
corresponding equations are shown in (8):

$\begin{array}[]{ccl}\dot{S_{i}}&=&\Lambda_{i}^{\operatorname{eul}}-\beta_{i}\frac{S_{i}}{N_{i}}I_{v,i}+\sum_{j=1}^{n}m_{ij}^{\operatorname{S}}S_{j}-\sum_{j=1}^{n}m_{ji}^{\operatorname{S}}S_{i}-\mu_{i}^{\operatorname{eul}}S_{i}\\\
\dot{I_{i}}&=&\beta_{i}\frac{S_{i}}{N_{i}}I_{v,i}+\sum_{j=1}^{n}m_{ij}^{\operatorname{I}}I_{j}-\sum_{j=1}^{n}m_{ji}^{\operatorname{I}}I_{i}-\left(\gamma_{i}^{\operatorname{eul}}+\mu_{i}^{\operatorname{eul}}\right)I_{i}\\\
\dot{R_{i}}&=&\gamma_{i}^{\operatorname{eul}}I_{i}+\sum_{j=1}^{n}m_{ij}^{\operatorname{R}}R_{j}-\sum_{j=1}^{n}m_{ji}^{\operatorname{R}}R_{i}-\mu_{i}^{\operatorname{eul}}R_{i}\\\
\dot{S}_{v,i}&=&\Lambda_{v,i}-\beta_{v,i}\frac{I_{i}}{N_{i}}S_{v,i}+\sum_{j=1}^{n}m_{ij}^{v}S_{v,j}-\sum_{j=1}^{n}m_{ji}^{v}S_{v,i}-\mu_{v,i}S_{v,i}\\\
\dot{I}_{v,i}&=&\beta_{v,i}\frac{I_{i}}{N_{i}}S_{v,i}+\sum_{j=1}^{n}m_{ij}^{v}I_{v,j}-\sum_{j=1}^{n}m_{ji}^{v}I_{v,i}-\mu_{v,i}I_{v,i}\,,\end{array}$
(8)

for $i=1,\dots,n$.

A summary of the parameters of system (8) is given in Table 1. We assume A1-A2
hold for system (8). In addition, we assume A5-A7:

  * A5:

The adjacency matrix for host movement has zeros on the diagonal (i.e.
$m^{\operatorname{X}}_{ii}=0$ for all $i$ and
$\operatorname{X}\in\\{\operatorname{S},\operatorname{I},\operatorname{R}\\}$).

  * A6:

The parameters $\beta_{i},\beta_{v,i},\Lambda_{v,i},\mu_{v,i}$ are intrinsic
to the patch $i$, so they are considered to be the same as in system (3).

  * A7:

The parameters
$\Lambda_{i}^{\operatorname{eul}},\Lambda_{v,i},\gamma_{i}^{\operatorname{eul}}$,
$\mu_{i}^{\operatorname{eul}}$ and $\mu_{v,i}$ are all positive.

Note that in system (8), individuals assume the characteristics of the patch
they are currently located in. For example, a host individual that migrates
from $j$ to $i$ now recovers from infection at rate
$\gamma_{i}^{\operatorname{eul}}$. The force of infection for patch $i$ in
system (8) depends only upon population abundances in patch $i$ (in contrast
to the Lagrangian system (3), where the force of infection in $i$ involves
contributions from other patches weighted by the mixing matrix $P$).

Let $L^{\operatorname{X}},L_{v}$ denote the graph Laplacians corresponding to
the host ($M^{\operatorname{X}}$, for
$\operatorname{X}\in\\{\operatorname{S},\operatorname{I},\operatorname{R}\\}$)
and vector ($M^{v}$) movement matrices, respectively. Let
$G_{v}\mathrel{\mathop{\mathchar 58\relax}}=L_{v}+D_{\delta_{v}}$ as in
Section 2.1, and let
$D_{\mu}^{\operatorname{eul}}=\operatorname{diag}\left\\{\mu_{i}^{\operatorname{eul}}\right\\}$.
Then at the DFE, the susceptible compartments of system (8) are given by

$\begin{array}[]{rcl}\left(S^{\operatorname{eul}}\right)^{*}=\left(N^{\operatorname{eul}}\right)^{*}&=&\left(L^{\operatorname{S}}+D_{\mu}^{\operatorname{eul}}\right)^{-1}\Lambda^{\operatorname{eul}}\,,\\\
\left(S_{v}^{\operatorname{eul}}\right)^{*}=\left(N_{v}^{\operatorname{eul}}\right)^{*}&=&G_{v}^{-1}\Lambda_{v}\,,\end{array}$
(9)

where
$\left(S^{\operatorname{eul}}\right)^{*},\left(N^{\operatorname{eul}}\right)^{*},\left(S_{v}^{\operatorname{eul}}\right)^{*},\left(N_{v}^{\operatorname{eul}}\right)^{*},\Lambda^{\operatorname{eul}}$
and $\Lambda_{v}$ are column vectors as before. Notice that
$\left(N_{v}^{\operatorname{eul}}\right)^{*}=\left(N_{v}^{\operatorname{lag}}\right)^{*}=G_{v}^{-1}\Lambda_{v}$,
so we define the vector $N_{v}^{*}\mathrel{\mathop{\mathchar
58\relax}}=\left(N_{v}^{\operatorname{eul}}\right)^{*}=\left(N_{v}^{\operatorname{lag}}\right)^{*}$
with entries
$N_{v,i}^{*}=\left(N_{v,i}^{\operatorname{eul}}\right)^{*}=\left(N_{v,i}^{\operatorname{lag}}\right)^{*}$.

Consider the basic reproduction number ${\cal R}_{0}^{\operatorname{eul}}$ for
system (8), computed using the next generation matrix approach. Define
$D_{\beta}\mathrel{\mathop{\mathchar
58\relax}}=\operatorname{diag}\\{\beta_{i}\\}$ and
$D_{\beta_{v}}^{\operatorname{eul}}\mathrel{\mathop{\mathchar
58\relax}}=\operatorname{diag}\left\\{\beta_{v,i}\left(N_{v,i}^{\operatorname{eul}}\right)^{*}/\left(N_{i}^{\operatorname{eul}}\right)^{*}\right\\}$.
Let $D_{\delta}^{\operatorname{eul}}\mathrel{\mathop{\mathchar
58\relax}}=\operatorname{diag}\\{\delta_{i}^{\operatorname{eul}}\\}$, and
$G\mathrel{\mathop{\mathchar
58\relax}}=L^{\operatorname{I}}+D_{\delta}^{\operatorname{eul}}.$ Note that
assumption A2 implies that $G$ is a non-singular $M$-matrix. We show in
Appendix 7.1 that if

$F^{\operatorname{eul}}=\begin{pmatrix}0&D_{\beta}\\\
D_{\beta_{v}}^{\operatorname{eul}}&0\end{pmatrix}$ (10)

and

$V^{\operatorname{eul}}=\begin{pmatrix}G&0\\\ 0&G_{v}\end{pmatrix}\,,$ (11)

then

$\left(\mathcal{R}_{0}^{\operatorname{eul}}\right)^{2}=\rho\left(\left(F^{\operatorname{eul}}\right)\left(V^{\operatorname{eul}}\right)^{-1}\right)=\rho\left(D_{\beta}G_{v}^{-1}D_{\beta_{v}}^{\operatorname{eul}}G^{-1}\right)\,.$
(12)

As for the Lagrangian model, the next generation matrix is the product of two
terms, one ($D_{\beta}G_{v}^{-1}$) corresponding to secondary host infections
created by infectious vectors, and the other
($D_{\beta_{v}}^{\operatorname{eul}}G^{-1}$) corresponding to secondary vector
infections created by infectious hosts. Note that for the Eulerian model, host
movement affects the next generation matrix via the transfer matrix
$V^{\operatorname{eul}}$. This is in contrast with the Lagrangian model, where
host movement appears in the fecundity matrix $F^{\operatorname{lag}}$.
Additionally, in the Lagrangian model the host transmission rates in the next
generation matrix are scaled by host movement [i.e. the
$P^{\operatorname{T}}D_{\beta}G_{v}^{-1}$ term in (7)], whereas in the
Eulerian approach they are not [i.e. the $D_{\beta}G_{v}^{-1}$ term in (12)].
As we will see in the next sections, these differences lead to different
values of $\mathcal{R}_{0}^{\operatorname{lag}}$ and
$\mathcal{R}_{0}^{\operatorname{eul}}$ (in particular, see the end of Appendix
7.2 for intuition on the difference between
$\mathcal{R}_{0}^{\operatorname{lag}}$ and
$\mathcal{R}_{0}^{\operatorname{eul}}$ in a two-patch example).

##  3 Model comparison through a fundamental matrix

###  3.1 Preliminaries

Our objective is to compare the Lagrangian (3) and Eulerian (8) frameworks. As
pointed out by Cosner et al. (see Section 2.2.3 of [17]), the frameworks are
in general distinct if we try to relate the number of individuals of both
systems. Specifically, Cosner et al. showed that the dynamical system
resulting from modeling the number of individuals currently located in each
patch under the Lagrangian framework does not correspond to an Eulerian model.
Here, we take an alternative approach: we compare systems (3) and (8) by
tuning the host mobility matrix $M$ so that the expected amount of time spent
in one location starting from another matches between the two frameworks as
closely as possible.

We begin with the analysis of Cosner et al. [17], who considered when the
equilibria between an Eulerian and Lagrangian model can be matched. Let $X$
denote a host population type in $\\{S,I,N\\}$. To match population sizes, we
should have that the population size $X_{i}^{\operatorname{eul}}$ in any patch
$i$ for the Eulerian model is equal to the total combined proportions of
number of residents $X_{j}^{\operatorname{lag}}$ of any other patch $j$ that
are in $i$ for the Lagrangian model, which is
$\sum_{j}p_{ij}X_{j}^{\operatorname{lag}}$. In other words, if
$X^{\operatorname{eul}}$ and $X^{\operatorname{lag}}$ are column vectors with
entries $X_{i}^{\operatorname{eul}}$ and $X_{i}^{\operatorname{lag}}$
respectively, we should have:

$X^{\operatorname{eul}}=PX^{\operatorname{lag}}\,.$ (13)

We will assume in the remainder of the paper that (13) holds at the DFE. Using
$\left(N^{\operatorname{eul}}\right)^{*}=\left(L^{\operatorname{S}}+D_{\mu}^{\operatorname{eul}}\right)^{-1}\Lambda^{\operatorname{eul}}$
and
$\left(N^{\operatorname{lag}}\right)^{*}=\left(D_{\mu}^{\operatorname{lag}}\right)^{-1}\Lambda^{\operatorname{lag}}$,
we set the following assumption:

  * A8:

We assume that the DFEs of systems (3) and (8) match in the sense of (13), so
$\left(N^{\operatorname{eul}}\right)^{*}=P\left(N^{\operatorname{lag}}\right)^{*}$.
Specifically, let $\Lambda^{\operatorname{eul}}$ and
$\Lambda^{\operatorname{lag}}$ be such that

$\left(L^{\operatorname{S}}+D_{\mu}^{\operatorname{eul}}\right)^{-1}\Lambda^{\operatorname{eul}}=P\left(D_{\mu}^{\operatorname{lag}}\right)^{-1}\Lambda^{\operatorname{lag}}\,.$
(14)

The condition in (14) states that at the DFE the amount of susceptible
individual in patch $i$ (Eulerian equilibrium
$\left(N_{i}^{\operatorname{eul}}\right)^{*}$) is the average (weighted by the
row $i$ of $P$) over all different patches $j$ of the number of susceptible
individuals that are residents of $j$ (Lagrangian equilibrium
$\left(N_{j}^{\operatorname{lag}}\right)^{*}$). As pointed out by Cosner et
al. [17], in general it may not be possible to satisfy equation (14). Namely,
for a given $\Lambda^{\operatorname{lag}}$ we may find that
$\left(L^{\operatorname{S}}+D_{\mu}^{\operatorname{eul}}\right)P\left(D_{\mu}^{\operatorname{lag}}\right)^{-1}\Lambda^{\operatorname{lag}}$
has some negative entries, and since the equation
$\Lambda^{eul}=\left(L^{\operatorname{S}}+D_{\mu}^{\operatorname{eul}}\right)P\left(D_{\mu}^{\operatorname{lag}}\right)^{-1}\Lambda^{\operatorname{lag}}$
is equivalent to (14), then (14) would not be true for
$\Lambda^{\operatorname{eul}}$ with positive entries. However, since
$L^{\operatorname{S}}+D_{\mu}^{\operatorname{eul}}$ is an $M$-matrix (see
[10], page 137), then $L^{\operatorname{S}}+D_{\mu}^{\operatorname{eul}}$ is
semi-positive, i.e., there exists $x>0$ such that
$\left(L^{\operatorname{S}}+D_{\mu}^{\operatorname{eul}}\right)x>0$ (see [10],
page 136). Therefore, equation (14) can be satisfied for appropriate
$\Lambda^{\operatorname{lag}},\Lambda^{\operatorname{eul}}$.

In addition, since
$D_{\beta_{v}}^{\operatorname{lag}}\mathrel{\mathop{\mathchar
58\relax}}=\operatorname{diag}\left\\{\beta_{v,i}\left(N_{v,i}^{\operatorname{lag}}\right)^{*}/\hat{N_{i}}^{*}\right\\}$,
$D_{\beta_{v}}^{\operatorname{eul}}\mathrel{\mathop{\mathchar
58\relax}}=\operatorname{diag}\left\\{\beta_{v,i}\left(N_{v,i}^{\operatorname{eul}}\right)^{*}/\left(N_{i}^{\operatorname{eul}}\right)^{*}\right\\}$,
$\left(N_{v,i}^{\operatorname{lag}}\right)^{*}=\left(N_{v,i}^{\operatorname{eul}}\right)^{*}$
and
$\hat{N}^{*}=P\left(N^{\operatorname{lag}}\right)^{*}=\left(N^{\operatorname{eul}}\right)^{*}$
(under (14)), then

$D_{\beta_{v}}\mathrel{\mathop{\mathchar
58\relax}}=D_{\beta_{v}}^{\operatorname{lag}}=D_{\beta_{v}}^{\operatorname{eul}}\,.$

Therefore, under A14 we have

$\displaystyle\left(\mathcal{R}_{0}^{\operatorname{lag}}\right)^{2}$
$\displaystyle=\rho\left(P^{\operatorname{T}}D_{\beta}G_{v}^{-1}D_{\beta_{v}}P\left(D_{\delta}^{\operatorname{lag}}\right)^{-1}\right)\,,$
$\displaystyle\left(\mathcal{R}_{0}^{\operatorname{eul}}\right)^{2}$
$\displaystyle=\rho\left(D_{\beta}G_{v}^{-1}D_{\beta_{v}}G^{-1}\right)\,.$
(15)

We can also compare the recovery and mortality rates of both models. In
general, these parameters are distinct between the two modeling frameworks.
Consider, for example, the mortality rate for patch $i$ in each of the
frameworks. In the Eulerian model, $\mu_{i}^{\operatorname{eul}}$ reflects
only the characteristics of location $i$. By contrast, the mortality rate
$\mu_{i}^{\operatorname{lag}}$ in the Lagrangian model reflects
characteristics of all the spatial locations, weighted according to the
proportion of time that a resident of $i$ spends in each location. Thus, we
model this relationship by averaging the rates $\mu_{j}^{\operatorname{eul}}$
and setting
$\mu_{i}^{\operatorname{lag}}=\sum_{j}p_{ji}\mu_{j}^{\operatorname{eul}}$.
These considerations lead us to the following assumption:

  * A9:

We assume that the parameters $\mu_{i}^{\operatorname{lag}}$,
$\gamma_{i}^{\operatorname{lag}}$, $\delta_{i}^{\operatorname{lag}}$ are
related to $\mu_{i}^{\operatorname{eul}}$, $\gamma_{i}^{\operatorname{eul}}$,
$\delta_{i}^{\operatorname{eul}}$ by

$\begin{array}[]{rcl}\mu_{i}^{\operatorname{lag}}&=&\sum_{j}p_{ji}\mu_{j}^{\operatorname{eul}}\,,\\\
\gamma_{i}^{\operatorname{lag}}&=&\sum_{j}p_{ji}\gamma_{j}^{\operatorname{eul}}\,,\\\
\delta_{i}^{\operatorname{lag}}&=&\sum_{j}p_{ji}\delta_{j}^{\operatorname{eul}}\,.\end{array}$
(16)

In matrix form, this means
$D_{\delta}^{\operatorname{lag}}=\operatorname{diag}\left\\{\mathbf{1}^{\operatorname{T}}D_{\delta}^{\operatorname{eul}}P\right\\}$,
where $\mathbf{1}=(1,1,\ldots,1)^{\operatorname{T}}$.

###  3.2 Problem formulation

We compare the Lagrangian and Eulerian modeling frameworks through a
relationship that can be interpreted in terms of a fundamental matrix for a
Markov process. Systems (3) and (8) are deterministic, not stochastic.
However, aspects of both systems can be interpreted in terms of a fundamental
matrix for a continuous time random walk [18]. We will use this fundamental
matrix to relate the two systems. Consider the next generation matrix
$F^{\operatorname{eul}}\left(V^{\operatorname{eul}}\right)^{-1}$ for system
(8) [54]. The transfer matrix $V^{\operatorname{eul}}$ includes a block
$G=L^{\operatorname{I}}+D_{\delta}^{\operatorname{eul}}$ that generates an
absorbing random walk on the host movement network. $G^{-1}$ corresponds to
the fundamental matrix of this random walk, with $(i,j)$ entry giving the
expected time that an infectious individual starting in $j$ spends in $i$
before being absorbed (removed) from the system [51, 54]. This interpretation
of the transfer matrix underlies intuition for $FV^{-1}$ as giving the number
of infectious individuals in the ‘next generation’, and corresponding
threshold of $\rho\left(FV^{-1}\right)>1$ for disease invasion. See [54].

Now consider the mixing matrix $P$ for the Lagrangian system (3). The entries
$p_{ij}$ give the probability that a resident of $j$ is in patch $i$, and
$1/\delta_{j}^{\operatorname{lag}}$ gives the expected time that a resident of
$j$ stays infectious, so $p_{ij}/\delta_{j}^{\operatorname{lag}}$ represents
the expected time that an infectious host from patch $j$ spends in patch $i$
according to the Lagrangian approach. Matching the expected ‘residence times’
(times spent in $i$, starting from $j$) for the Eulerian and Lagrangian
frameworks and applying A9, we have:

$\left(L^{\operatorname{I}}+D_{\delta}^{\operatorname{eul}}\right)^{-1}=P\left(D_{\delta}^{\operatorname{lag}}\right)^{-1}=(P)\left(\operatorname{diag}^{-1}\left\\{\mathbf{1}^{\operatorname{T}}D_{\delta}^{\operatorname{eul}}P\right\\}\right)\,.$
(17)

As we will see in Section 3.4, for a given $D_{\delta}^{\operatorname{eul}}$
and $P$ it may not always be possible to find a graph Laplacian matrix
$L^{\operatorname{I}}$ such that (17) holds. This fact leads us to the
following definition.

######  Definition 1.

(Consistency) Assume A1-A9 and suppose that $P$ and
$D_{\delta}^{\operatorname{eul}}$ are given. We say that systems (3) and (8)
are consistent if there exists a graph Laplacian matrix $L^{\operatorname{I}}$
satisfying (17), and inconsistent if such a matrix does not exist.

Note that as $L^{\operatorname{I}}+D_{\delta}^{\operatorname{eul}}$ and
$D_{\delta}^{\operatorname{lag}}$ are non-singular, if systems (3) and (8) are
consistent, then $P$ is non-singular. However, as we will see in Section 3.4,
the converse is not necessarily true. We will assume in the following sections
that $P$ is non-singular.

  * A10:

$P$ is non-singular.

Notice that under A10, relationship (17) is equivalent to

$L^{\operatorname{I}}=D_{\delta}^{\operatorname{lag}}P^{-1}-D_{\delta}^{\operatorname{eul}}\,.$
(18)

Let $p_{ij}^{\prime}$ denote the $(i,j)$ entry of the inverse of $P$, that is,
let $P^{-1}=(p_{ij}^{\prime})_{i,j\leq n}$. Then the $j^{th}$ column sum of
$D_{\delta}^{\operatorname{lag}}P^{-1}-D_{\delta}^{\operatorname{eul}}$ is

$\left(\sum_{i}p_{ij}^{\prime}\delta_{i}^{\operatorname{lag}}\right)-\delta_{j}^{\operatorname{eul}}=\left(\sum_{i}p_{ij}^{\prime}\sum_{k}p_{ki}\delta_{k}^{\operatorname{eul}}\right)-\delta_{j}^{\operatorname{eul}}=\left(\sum_{k}\delta_{k}^{\operatorname{eul}}\sum_{i}p_{ki}p_{ij}^{\prime}\right)-\delta_{j}^{\operatorname{eul}}=0\,.$
(19)

The above condition is consistent with the graph Laplacian having zero column
sums. However, we also want the off-diagonal elements of
$D_{\delta}^{\operatorname{lag}}P^{-1}-D_{\delta}^{\operatorname{eul}}$ to be
non-positive, which is not in general the case as we will show in the example
of Section 3.4. In consequence, we are interested in finding

$E=\inf_{L\in\mathbb{R}^{n\times n}}\left\\{\lVert
L-(D_{\delta}^{\operatorname{lag}}P^{-1}-D_{\delta}^{\operatorname{eul}})\rVert_{F}\mathrel{\mathop{\mathchar
58\relax}}L\text{ is a graph Laplacian matrix }\right\\},$ (20)

where we use the Frobenius norm defined by $\lVert
B\rVert_{F}=\sqrt{\operatorname{tr}(BB^{\operatorname{T}})}$ for a given
matrix $B$. The Frobenius norm allows us to treat (20) as a non-negative least
squares problem. For example, in the 2 by 2 case, we have that

$E=\inf_{m\geq 0}\lVert L^{*}m-\bar{x}\rVert_{2}\,,$ (21)

where

$L^{*}=\begin{pmatrix}1&0\\\ -1&0\\\ 0&1\\\ 0&-1\end{pmatrix},\quad
m=\begin{pmatrix}m_{21}\\\ m_{12}\end{pmatrix},\quad
D_{\delta}^{\operatorname{lag}}P^{-1}-D_{\delta}^{\operatorname{eul}}=\begin{pmatrix}a&c\\\
b&d\end{pmatrix},\quad\bar{x}=(a,b,c,d)^{\operatorname{T}}\,.$

Notice that $\lVert m\rVert_{2}\leq\lVert L\rVert_{F}=\lVert
L^{*}m\rVert_{2}\leq\lVert L^{*}m-\bar{x}\rVert_{2}+\lVert\bar{x}\rVert_{2}$.
Therefore, if $\\{m_{k}\\}_{k\geq 1}$ is a sequence such that
$\lim_{k\rightarrow\infty}\lVert L^{*}m_{k}-\bar{x}\rVert_{2}=E$, then
$\limsup_{k\geq 1}\lVert m_{k}\rVert_{2}\leq E+\lVert\bar{x}\rVert_{2}$. In
consequence, $\\{\lVert m_{k}\rVert_{2}\\}_{k\geq 1}$ is bounded and therefore
$\inf_{m\geq 0}\lVert L^{*}m-\bar{x}\rVert_{2}$ is attained, i.e.,
$E=\min_{m\geq 0}\lVert L^{*}m-\bar{x}\rVert_{2}.$

In general, $L^{*}$ is an $n^{2}\times n(n-1)$ matrix, and from the Karush-
Kuhn-Tucker conditions the optimum $\bar{m}$ satisfies
$(L^{*}\bar{m}-\bar{x})^{\operatorname{T}}L^{*}\bar{m}=0,\bar{m}\geq 0$ (see
Section 10.10 of [13]), giving

$\bar{x}^{\operatorname{T}}L^{*}\bar{m}=\lVert L^{*}\bar{m}\rVert_{2}^{2}\geq
0,E^{2}=(L^{*}\bar{m}-\bar{x})^{\operatorname{T}}(L^{*}\bar{m}-\bar{x})=\lVert
D_{\delta}^{\operatorname{lag}}P^{-1}-D_{\delta}^{\operatorname{eul}}\rVert_{F}^{2}-\bar{x}^{\operatorname{T}}L^{*}\bar{m}\,.$
(22)

Hence we have the following upper bound for $E$:

$E\leq\lVert
D_{\delta}^{\operatorname{lag}}P^{-1}-D_{\delta}^{\operatorname{eul}}\rVert_{F}\,.$
(23)

###  3.3 Consistency

####  3.3.1 Sufficient condition

By direct calculation we can show that the condition $p_{12}+p_{21}<1$
guarantees that systems (3) and (8) are consistent in the two-patch setting.
This suggests that the off-diagonal terms of the mixing matrix $P$ must be
sufficiently small for the Lagrangian and Eulerian systems to be consistent.
On the other hand, if

$P=\begin{pmatrix}9/10&0&0\\\ 0&1&1/10\\\ 1/10&0&9/10\end{pmatrix},\text{ then
}D_{\delta}^{\operatorname{lag}}P^{-1}-D_{\delta}^{\operatorname{eul}}=\begin{pmatrix}10/9\delta_{1}^{\operatorname{lag}}-\delta_{1}^{\operatorname{eul}}&0&0\\\
1/81\delta_{2}^{\operatorname{lag}}&\delta_{2}^{\operatorname{lag}}-\delta_{2}^{\operatorname{eul}}&-1/9\delta_{2}^{\operatorname{lag}}\\\
-10/81\delta_{3}^{\operatorname{lag}}&0&10/9\delta_{3}^{\operatorname{lag}}-\delta_{3}^{\operatorname{eul}}\end{pmatrix}.$
(24)

From (24),
$D_{\delta}^{\operatorname{lag}}P^{-1}-D_{\delta}^{\operatorname{eul}}$ has
positive off-diagonal entries, and thus systems (3) and (8) are inconsistent.
In this example, some entries $p_{ij}$ of $P$ are small (for example
$p_{21}=0$), which suggests that non-diagonal entries of the mixing matrix
must also be sufficiently large in order for the Eulerian and Lagrangian
frameworks to be consistent. In Proposition 1 below, we prove that a
sufficient condition for the consistency of systems (3) and (8) is that the
off-diagonal entries of the mixing matrix belong to an intermediate range
$p_{*}<p_{ij}<p^{*}$.

######  Proposition 1.

Assume A1-A10 and let $P=(p_{ij})_{i,j\leq n}$ be the mixing matrix associated
with system (3), where $n\geq 2$. Let $p_{*}$ and $p^{*}$ be constant numbers
in the interval $(0,1)$ such that

$\frac{4(n-1)^{2}p^{*}}{1-p^{*}}<1\text{ and
}p_{*}=\frac{4(n-1)^{2}(p^{*})^{2}}{1-p^{*}}\,.$ (25)

In addition, suppose that

$p_{*}\leq p_{ij}\leq p^{*},\text{ for }i\neq j\,.$ (26)

Then, the systems (3) and (8) are consistent.

###### Proof.

From (19), it follows that the column sums of
$D_{\delta}^{\operatorname{lag}}P^{-1}-D_{\delta}^{\operatorname{eul}}$ are
all zero. Consequently, we just need to show that the off-diagonal entries of
$D_{\delta}^{\operatorname{lag}}P^{-1}-D_{\delta}^{\operatorname{eul}}$ are
non-positive in order to prove that systems (3) and (8) are consistent.
Observe that if $p^{*}$ goes to zero, then $4(n-1)^{2}p^{*}/(1-p^{*})$
approaches to zero, so the first part of (25) can be satisfied by small enough
$p^{*}$. The condition (25) then implies that

$\frac{p_{*}}{p^{*}}=\frac{4(n-1)^{2}p^{*}}{1-p^{*}}<1\,.$

Therefore, the interval $[p_{*},p^{*}]$ is non-empty and we can pick $p_{ij}$
such that (26) holds. In order to show that the off-diagonal entries of
$D_{\delta}^{\operatorname{lag}}P^{-1}-D_{\delta}^{\operatorname{eul}}$ are
non-positive, it suffices to prove that the off-diagonal entries of $P^{-1}$
are non-positive. Let $\Delta P$ be a matrix such that $P=I-\Delta P$, and
consider the matrix 1-norm defined by $\lVert
B\rVert_{1}=\max_{j}\sum_{i}|b_{ij}|$ for a given matrix $B$. We can write
$P^{-1}=(I-\Delta P)^{-1}=\mathcal{P}^{(0)}+\mathcal{P}^{(1)}$, where
$\mathcal{P}^{(0)}=I+\Delta P$ and $\mathcal{P}^{(1)}=(\Delta
P)^{2}\sum_{k\geq 0}(\Delta P)^{k}$. If $i\neq j$, then
$\mathcal{P}^{(0)}_{ij}=-p_{ij}$, so we want to show
$|\mathcal{P}^{(1)}_{ij}|\leq p_{ij}$ to get that the off-diagonal elements of
$P^{-1}$ are non-positive. Indeed, from (25) and (26) it follows that
$\lVert\Delta P\rVert_{1}\leq 2(n-1)p^{*}$ and

$\displaystyle|\mathcal{P}^{(1)}_{ij}|$
$\displaystyle\leq\lVert\mathcal{P}^{(1)}\rVert_{1}$
$\displaystyle\leq\lVert\Delta P\rVert_{1}^{2}\sum_{k\geq 0}\lVert\Delta
P\rVert_{1}^{k}$ $\displaystyle=\lVert\Delta P\rVert_{1}^{2}/(1-\lVert\Delta
P\rVert_{1})$ $\displaystyle\leq\left[2(n-1)p^{*}\right]^{2}/(1-p^{*})$
$\displaystyle=p_{*}\leq p_{ij}\,,$

as we desired. ∎

Proposition 1 establishes that systems (3) and (8) are consistent when the
off-diagonal entries of the mixing matrix $P$ lie in an interval
$[p_{*},p^{*}]\subset(0,1)$. The upper and lower bounds for this interval
satisfy (25) and (26). We note that these bounds are not unique, as more than
one $p^{*}$ can satisfy (25) and (26). Furthermore, the width of the resulting
interval $(p_{*},p^{*})$ may be small for large $n$.

For example, let $n=3$. From (25) it suffices to choose $p^{*}$ such that
$p^{*}<1/17$. Let us consider, for instance, $p^{*}=0.05$. Consequently, from
(25) we have $p_{*}=0.0421$ and then $p_{*}=0.0239\leq p_{ij}\leq
p^{*}=0.1324$ is a sufficient condition for consistency of systems (3) and
(8). We note that this interval may not be the widest interval among those
obtained using other values of $p^{*}$ satisfying (25).

Additionally, (25) may be improved. Namely, we can write
$P^{-1}=\mathcal{P}^{(0)}+\mathcal{P}^{(1)}$ where $\mathcal{P}^{(0)}=I+\Delta
P+(\Delta P)^{2}$ and $\mathcal{P}^{(1)}=(\Delta P)^{3}\sum_{k\geq 0}(\Delta
P)^{k}$. By imposing the condition $p_{*}\leq p_{ij}\leq p^{*}$, with $i\neq
j$, we have that $p_{*}$ and $p^{*}$ must satisfy

$\displaystyle-\mathcal{P}^{(0)}_{ij}$ $\displaystyle\geq
p_{*}+2(n-1)p_{*}^{2}-(n-2)(p^{*})^{2}$
$\displaystyle\geq\frac{\left[2(n-1)p^{*}\right]^{3}}{1-p^{*}}$
$\displaystyle\geq|\mathcal{P}^{(1)}_{ij}|\,.$ (27)

Figure 1 shows the region of pairs $(p_{*},p^{*})$ that satisfy the inequality
(3.3.1) for $n=5$. Point A of Figure 1 indicates that $p_{*}=0.0065\leq
p_{ij}\leq p^{*}=0.022$, for $i\neq j$, is a sufficient condition for
consistency of systems (3) and (8). This condition improves (25), where
$p^{*}=0.015$ is associated to $p_{*}=0.0065$ and the interval
$(0.0065,0.015)$, which is smaller than the interval $(0.0065,0.022)$
corresponding to (3.3.1).

Figure 1: Inequality (3.3.1) holds in the blue region for $n=5$. For every
point $(p_{*},p^{*})$ on the top boundary of the blue region we have a
consistency condition. Namely, under A1-A10, if $p_{ij}\in[p_{*},p^{*}]$ for
all $i\neq j$, then systems (3) and (8) are consistent. For example, at the
black point A we have $p_{*}=0.0065$ and $p^{*}=0.022$. Therefore, if
$0.0065\leq p_{ij}\leq 0.022$ for $i\neq j$, then the systems (3) and (8) are
consistent.

####  3.3.2 A consistent example: star graphs

In this section we give an example where the Eulerian and Lagrangian systems
are consistent. Specifically, we consider a setting where the mixing matrix
$P$ corresponds to a star graph, where the ‘hub’ node is the only location
that residents of other patches visit. This setting is motivated by empirical
networks where there exist nodes $k$ for which all the $p_{ki}$ are large. For
example, in the data analyzed in [44] on malaria in Namibia, non-residents are
much more likely to visit a few locations (e.g. the capital Windhoek) than
others. A schematic of the class of mixing matrices considered in this section
is shown in Figure 2. We will show that for such $P$, systems (3) and (8) are
consistent.

Figure 2: Star graph as example of consistency.

We will use the Sherman-Morrison formula

$\left(A+vu^{\operatorname{T}}\right)^{-1}=A^{-1}-\frac{A^{-1}vu^{\operatorname{T}}A^{-1}}{1+u^{\operatorname{T}}A^{-1}v}\,,$
(28)

where $A$ is non-singular and $v,u$ are column vectors such that
$1+u^{\operatorname{T}}A^{-1}v\neq 0$. Let us suppose that

$P=\operatorname{diag}\\{1,1-p_{12},\ldots,1-p_{1n}\\}+(1,0,\ldots,0)^{\operatorname{T}}(0,p_{12},\ldots,p_{1n})\,,$
(29)

$D_{\delta}^{\operatorname{lag}}=D_{\delta}^{\operatorname{eul}}=D_{\delta}=\delta
I$ and let $b_{ij}$ denote the $(i,j)$ entry of
$D_{\delta}^{\operatorname{lag}}P^{-1}-D_{\delta}^{\operatorname{eul}}$. Using
the Sherman-Morrison formula (28), we get

$P^{-1}=\operatorname{diag}\\{1,1/(1-p_{12}),\ldots,1/(1-p_{1n})\\}-(1,0,\ldots,0)^{\operatorname{T}}(0,p_{12}/(1-p_{12}),\ldots,p_{1n}/(1-p_{1n}))$

and

$b_{ij}=\frac{-\delta p_{ij}}{1-p_{ij}},\textrm{ for }i\neq j\,.$

Then, as $p_{ij}<1$ for $i\neq j$, systems (3) and (8) are consistent.

Moreover, let $f\mathrel{\mathop{\mathchar 58\relax}}\mathbb{R}^{n\times
n}\rightarrow\mathbb{R}^{n\times n}$ be defined by
$f(P)\mathrel{\mathop{\mathchar 58\relax}}=D_{\delta}(P^{-1}-I)$, and let $U$
denote the open set $U=\\{M\in\mathbb{R}^{n\times n}\mathrel{\mathop{\mathchar
58\relax}}\text{ the off-diagonal entries of }M\text{ are negative }\\}$. If
$P_{0}$ is of the form (29), then $f(P_{0})\in U$.

By continuity of $f$, for $P$ with small enough $p_{ij},i\neq j$ and $i\neq
1$, we have that $f(P)\in U$, i.e, systems (3) and (8) are consistent.

###  3.4 An inconsistent example

We now present an example where systems (3) and (8) are inconsistent, and in
fact the upper bound in (23) is attained.

Suppose that

$\begin{array}[]{ccl}p_{ij}&\mathrel{\mathop{\mathchar 58\relax}}=&p_{i},i\neq
j\,,\\\ p_{ii}&\mathrel{\mathop{\mathchar 58\relax}}=&1-\sum_{k\neq
i}p_{k}\,.\end{array}$ (30)

Let $\theta\mathrel{\mathop{\mathchar 58\relax}}=1-\sum_{k=1}^{n}p_{k}\neq 0$.
Then we have

$P=\theta I+(p_{1},\ldots,p_{n})^{\operatorname{T}}(1,\ldots,1)\,.$

As $A\mathrel{\mathop{\mathchar 58\relax}}=\theta I$ is non-singular and
$1+(1,\ldots,1)A^{-1}(p_{1},\ldots,p_{n})^{\operatorname{T}}=1+\frac{1}{\theta}\sum_{k=1}^{n}p_{k}=\frac{1}{\theta}\neq
0$, applying the Sherman-Morrison formula gives

$P^{-1}=\frac{1}{\theta}I-\frac{1}{\theta}(p_{1},\ldots,p_{n})^{\operatorname{T}}(1,\ldots,1)\,.$

In addition, suppose that
$\delta_{i}^{\operatorname{eul}}\mathrel{\mathop{\mathchar
58\relax}}=\delta,i=1,\ldots,n$, so that
$D_{\delta}^{\operatorname{lag}}=D_{\delta}^{\operatorname{eul}}=\delta I$.
Therefore, if
$D_{\delta}^{\operatorname{lag}}P^{-1}-D_{\delta}^{\operatorname{eul}}=\delta(P^{-1}-I)\mathrel{\mathop{\mathchar
58\relax}}=(b_{ij})_{i,j\leq n}$, then we have that

$b_{ij}=-\frac{\delta p_{i}}{\theta},\quad i\neq j\,.$ (31)

For sufficiently large $p_{k}$ (such that $\sum_{k}p_{k}>1$) we have
$\theta<0$, leading to positive $b_{ij}$ for $i\neq j$, which is inconsistent
with the off-diagonal entries of a Laplacian matrix. In this case the optimum
in (22) is $\bar{m}=0$ and the error is
$E=||D_{\delta}^{\operatorname{lag}}P^{-1}-D_{\delta}^{\operatorname{eul}}||_{F}$,
which is the largest possible error.

As a specific example, consider three patches with
$p_{1}=0.8,p_{2}=0.15,p_{3}=0.15$. Then $\theta=-0.1<0$ and

$P=\begin{pmatrix}0.7&0.8&0.8\\\ 0.15&0.05&0.15\\\
0.15&0.15&0.05\end{pmatrix},$ (32)

for which $L=0$ is the solution of (20). Thus it is possible for systems (3)
and (8) to not only be inconsistent, but in fact for the upper bound in (23)
to be attained. We note that the preceding example requires that some of the
off-diagonal entries of $P$ are large, which may not be realistic in
situations where host individuals spend the majority of their time in a
distinguished ‘home’ location.

##  4 Two-patch network

In this section we explore results obtained from consistency and inconsistency
of systems (3) and (8) for two-patch systems ($n=2$). In Section 4.1 we
compare the final outbreak size and the basic reproduction number obtained
from systems (3) and (8) for an example where the systems are inconsistent. In
Section 4.2 we state explicit bounds for the relative difference between the
basic reproduction numbers of systems (3) and (8) when the transmission and
removal rates are the same for both patches in Proposition 2. In Section 4.2
we also compare the basic reproduction number of both systems for a particular
example where the removal rates for the two patches are different.

###  4.1 Final outbreak size and basic reproduction number for an inconsistent
example

We now consider the functional implications of consistency / inconsistency of
the Eulerian and Lagrangian frameworks, in terms of important disease
quantities such as the basic reproduction number and final outbreak size.

Let us consider an example for a two-patch network where
$P=\begin{pmatrix}1-p_{21}&p_{12}\\\ p_{21}&1-p_{12}\end{pmatrix}$, with
$p_{12}+p_{21}>1$. Let us also assume that
$D_{\delta}^{\operatorname{lag}}=D_{\delta}^{\operatorname{eul}}=D_{\delta}=\delta
I$, $D_{\beta_{v}}=\beta_{v}I$, $D_{\beta}=\beta I$,
$D_{\delta_{v}}=\delta_{v}I$ and $L_{v}=m_{v}\begin{pmatrix}1&-1\\\
-1&1\end{pmatrix}$. Then the Eulerian and Lagrangian frameworks are
inconsistent by the argument in Section 3.4, and $L=0$ optimizes (20) in this
case, meaning that the solution to (20) corresponds to a completely
disconnected set of nodes in the Eulerian framework. Thus the networks in the
Lagrangian and Eulerian frameworks are wildly different. Here we examine how
this difference in connectivity corresponds to differences in ${\cal R}_{0}$
and final outbreak size.

In Figure 3 we use the parameters $\delta=1/150$, $\beta=0.3\times 0.1$,
$\delta_{v}=0.05$, $m_{v}=0.02$ [44] and define $p_{12}=0.95$, $p_{21}=0.1$
(therefore $p_{12}+p_{21}>1$, which implies inconsistency). In Figure 3(a) we
observe that the final outbreak size obtained from the Lagrangian system is
larger than the final outbreak size obtained from the Eulerian system.
Furthermore, for values of ${\cal R}_{0}^{\operatorname{eul}}$ around one, we
get a significant relative difference between the outbreak sizes for the two
systems. For example, when ${\cal R}_{0}^{\operatorname{eul}}=1.2$, the
outbreak size of Lagrangian system is $220\%$ larger than the outbreak size of
the Eulerian system. Thus inconsistency in terms of Definition 1 might lead to
significant differences in outbreak size between the Eulerian and Lagrangian
frameworks. Moreover, we can modify the example of Figure 3 so the systems are
consistent and there is still a significant difference in outbreak size for
intermediate values of $\mathcal{R}_{0}^{\operatorname{eul}}$. On the other
hand, when ${\cal R}_{0}^{\operatorname{eul}}=1.6$, the percentage change is
$25.2\%$, and when ${\cal R}_{0}^{\operatorname{eul}}=0.5$ the final sizes are
almost the same. In Figure 3(b) we observe that ${\cal
R}_{0}^{\operatorname{lag}}>{\cal R}_{0}^{\operatorname{eul}}$ if $0\leq{\cal
R}_{0}^{\operatorname{eul}}\leq 2$. Additionally, ${\cal
R}_{0}^{\operatorname{lag}}$ is linear with respect to ${\cal
R}_{0}^{\operatorname{eul}}$, where ${\cal R}_{0}^{\operatorname{lag}}$ is at
most $18.4\%$ larger than ${\cal R}_{0}^{\operatorname{eul}}$.

(a)  (b)

Figure 3: In (a) we have the final outbreak size for the inconsistent example
given in Section 4.1 (the recovered individuals $R_{1}(T)+R_{2}(T)$ for large
enough $T$). In (b) we have the comparison between the basic reproduction
number of both models. The used parameters are
$a=0.3,b=0.1,\mu=0,\delta=\gamma=1/150,\delta_{v}=0.05,\beta=ab,m_{v}=0.02$
[44], and the units are as in Table 1. We also assume
$G_{v}=m_{v}(2I-\mathbf{1}\mathbf{1}^{\operatorname{T}}),N_{v}=80,\Lambda_{v}$
such that
$G_{v}^{-1}\Lambda_{v}=N_{v}\mathbf{1},\Lambda^{\operatorname{lag}}=\Lambda^{\operatorname{eul}}=(0,0)^{\operatorname{T}},N_{0}=100.$
For a given value of $\mathcal{R}_{0}^{\operatorname{eul}}$ we choose
$\beta_{v,1}^{\operatorname{eul}}=\beta_{v,2}^{\operatorname{eul}}$ such that
$\beta_{v}\mathrel{\mathop{\mathchar
58\relax}}=\beta_{v,i}^{\operatorname{eul}}=\left(\mathcal{R}_{0}^{\operatorname{eul}}\right)^{2}\delta_{v}\delta
N_{0}/(\beta N_{v})$. We also define $\beta_{v,1}^{\operatorname{lag}}$ and
$\beta_{v,2}^{\operatorname{lag}}$ such that
$D_{\beta_{v}}^{\operatorname{lag}}=D_{\beta_{v}}^{\operatorname{eul}}=\beta_{v}N_{v}/N_{0}I$.
The initial conditions are
$S_{1}(0)=N_{0},S_{2}(0)=N_{0}-1,I_{1}(0)=0,I_{2}(0)=1,R_{1}(0)=R_{2}(0)=0,S_{v,1}(0)=S_{v,2}(0)=50,I_{v,1}=I_{v,2}=0$
and the final outbreak is taken at time $T=3000$.

###  4.2 Homogeneous infection

In this section, we consider the special case where parameter values for
transmission and removal are the same across both locations in the two-patch
network. We will show that in this case we can bound the difference between
the basic reproduction numbers for the Eulerian and Lagrangian frameworks.
Furthermore, we show that in this ‘homogeneous infection’ setting, the basic
reproduction number for the Lagrangian model is greater than or equal to the
reproduction number for the Eulerian model.

Consider the case where $0<p_{12},p_{21}<1/2$. Under these conditions, the
systems (3) and (8) are consistent. We assume A1–A10, and pick the same
transmission and recovery parameters for both patches. We refer to this
setting as the homogeneous infection scenario. The following proposition
describes the behavior of the quantities
$\mathcal{R}_{0}^{\operatorname{eul}}$ and
$\mathcal{R}_{0}^{\operatorname{lag}}$ for $0<p_{12}<1/2$ in the homogeneous
infection scenario. In this proposition, we can interchange $p_{12}$ by
$p_{21}$ and obtain symmetric results. A proof and additional details are
given in Appendix 7.2.

######  Proposition 2.

Assume A1–A10 and suppose that $D_{\beta}=\beta I$,
$D_{\beta_{v}}=\beta_{v}I$,
$D_{\delta}^{\operatorname{eul}}=D_{\delta}^{\operatorname{lag}}=D_{\delta}\mathrel{\mathop{\mathchar
58\relax}}=\delta I$, $D_{\delta_{v}}=\delta_{v}I$ and
$m_{12}^{v}=m_{21}^{v}=m_{v}$. Define $\mathcal{M}\mathrel{\mathop{\mathchar
58\relax}}=D_{\beta}G_{v}^{-1}D_{\beta_{v}}P\left(D_{\delta}^{\operatorname{lag}}\right)^{-1}$
and fix $p_{21}$ in the interval $(0,1/2)$. Then
$\rho(P^{\operatorname{T}}\mathcal{M})$ is a function of $p_{12}$ where
$0<p_{12}<1/2$, and we have that:

  1. a)

$\left({\cal
R}_{0}^{\operatorname{eul}}\right)^{2}=\rho(\mathcal{M})=(\beta\beta_{v})/(\delta\delta_{v})$
is constant on $0<p_{12}<1/2$.

  2. b)

$\left({\cal R}_{0}^{\operatorname{lag}}\right)^{2}=\rho(P^{t}\mathcal{M})$ is
decreasing on $(0,p_{21})$, increasing on $(p_{21},1/2)$ and attains its
absolute minimum over $[0,1/2]$ with value $\left({\cal
R}_{0}^{\operatorname{eul}}\right)^{2}=(\beta\beta_{v})/(\delta\delta_{v})$ at
$p_{12}=p_{21}$.

  3. c)

In addition, we have the inequality

$\begin{split}\frac{\left({\cal
R}_{0}^{\operatorname{lag}}\right)^{2}-\left({\cal
R}_{0}^{\operatorname{eul}}\right)^{2}}{\left({\cal
R}_{0}^{\operatorname{eul}}\right)^{2}}=\frac{\rho(P^{\operatorname{T}}\mathcal{M})-\rho(\mathcal{M})}{\rho(\mathcal{M})}\leq\frac{1}{4}\frac{\delta_{v}}{(2m_{v}+\delta_{v})}\,.\end{split}$
(33)

From (33) and using that ${\cal R}_{0}^{\operatorname{lag}}>{\cal
R}_{0}^{\operatorname{eul}}$, we get the following bound for the relative
difference $\mathcal{R}_{0}^{\operatorname{lag}}$ with respect to
$\mathcal{R}_{0}^{\operatorname{eul}}$:

$\frac{\mathcal{R}_{0}^{\operatorname{lag}}-\mathcal{R}_{0}^{\operatorname{eul}}}{\mathcal{R}_{0}^{\operatorname{eul}}}\leq\frac{1}{4}\frac{{\cal
R}_{0}^{\operatorname{eul}}}{{\cal R}_{0}^{\operatorname{lag}}+{\cal
R}_{0}^{\operatorname{eul}}}\leq\frac{1}{8}\,.$

In consequence, the percentage difference between the basic reproduction
numbers for systems (3) and (8) is at most $12.5\%$ under homogeneous
infection if $0<p_{12},p_{21}<1/2$. In addition, the larger $|p_{12}-p_{21}|$
is, the larger the difference between $\mathcal{R}_{0}^{\operatorname{eul}}$
and $\mathcal{R}_{0}^{\operatorname{lag}}$ is as well, as Figure 4 shows.

Figure 4: Basic reproduction number for system (3) in the two-patch case as a
function of $p_{12}$ with
$a=0.3,b=0.1,c=0.214,\delta=r=1/150,\delta_{v}=0.1,m\mathrel{\mathop{\mathchar
58\relax}}=\left(N_{v,i}^{\operatorname{lag}}\right)^{*}/\hat{N_{i}}^{*}=\left(N_{v,i}^{\operatorname{eul}}\right)^{*}/\left(N_{i}^{\operatorname{eul}}\right)^{*}=1,\beta=ab,\beta_{v}=acm$
[44], and units as in Table 1. In this case, the largest percentage difference
has value
$100\left(\mathcal{R}_{0}^{\operatorname{lag}}-\mathcal{R}_{0}^{\operatorname{eul}}\right)/\mathcal{R}_{0}^{\operatorname{eul}}=3.36\%$
and is attained at $p_{12}=1/2$.

In conclusion, under homogeneous conditions for both patches, the introduction
of infectious individuals creates more secondary infections according to the
Lagrangian dynamics for any matrix $P$. By contrast, if we suppose
$\delta_{v,1}>\delta_{v,2}$, it is then possible that ${\cal
R}_{0}^{\operatorname{lag}}<{\cal R}_{0}^{\operatorname{eul}}$ when
$p_{12}>p_{21}$ (see Figure 5).

Figure 5:  Comparison reproduction numbers for $\delta_{v,1}>\delta_{v,2}$.
The used parameters are
$a=0.3,b=0.1,c=0.214,\delta_{1}=\delta_{2}=1/150,\delta_{v,1}=0.1,\beta_{1}=\beta_{2}=ab,\beta_{v,1}=\beta_{v,2}=ac,m_{v}=0.02$
[44], and the corresponding units are as in Table 1.

##  5 Examples using data

We now turn to applying our analytical results and definition of consistency
to empirical data. There is an abundance of empirical data on mobility and
connectivity between locations [7, 33, 58, 59], and these data are being
increasingly incorporated into mathematical and computational models of
vector-borne disease dynamics [27, 44]. Many factors are involved for deciding
whether to use a Lagrangian or Eulerian modeling framework, including the
spatial scale involved, mathematical tractability, and type of data available.
Here we consider two empirical data sets (Sections 5.1, 5.2) and two
hypothetical data sets (Sections 5.3, 5.4) on host movement in the context of
vector-borne disease on spatial scales ranging from within village [55] to
country-wide [44], with data sources including long-term GPS trackers [55] and
mobile phones [44]. Studies in [27, 44] incorporate the corresponding data
into a Lagrangian modeling framework. Here we examine whether the data are
consistent with an Eulerian framework in the sense of Definition 1, and
discuss possible implications of using an Eulerian approach for these specific
settings in terms of the basic reproduction number.

###  5.1 Malaria in Namibia

Ruktanonchai et al. [44] use mobile phone records to examine movement between
health districts in Namibia, in the context of malaria control efforts.
Specifically, the authors identified mobility sources / sinks from mobile
phone records, together with transmission hot spots from malaria parasite
maps. The mobility data and local transmission parameters are combined in a
Lagrangian framework for vector-host dynamics [44] .

Anonymized mobile phone records were collected over a year from 1.2 million
phones, corresponding to approximately 85% of the adult population in Namibia.
Call and SMS data were used to identify locations at the health district
level. Home health districts and location changes were estimated, and
aggregated to produce a Lagrangian mixing matrix $P$ [44]. For most of the
locations in these data we have that the quantities
$F_{i}^{(\textrm{in})}\mathrel{\mathop{\mathchar 58\relax}}=\sum_{j\neq
i}p_{ij}$ and $F_{i}^{(\textrm{out})}\mathrel{\mathop{\mathchar
58\relax}}=\sum_{j\neq i}p_{ji}$ are small.

Using the data considered in [44], we take the ten health districts for which
the quantity
$\left(F_{i}^{(\textrm{in})}\right)^{2}+\left(F_{i}^{(\textrm{out})}\right)^{2}$
is the largest and define the $10\times 10$ mixing matrix $P$ for these
locales. Parameter values used in this example are as in [44]: $a=0.3\hskip
5.69046pt\text{Hosts}\times\text{Days}^{-1}$, $b=0.1\hskip
5.69046pt\text{Vectors}^{-1}$, $c=0.214\hskip 5.69046pt\text{Hosts}^{-1}$,
$\delta=r=1/150\hskip 5.69046pt\text{Days}^{-1}$, $\delta_{v}=0.1\hskip
5.69046pt\text{Days}^{-1}$, $m\mathrel{\mathop{\mathchar
58\relax}}=\left(N_{v,i}^{\operatorname{lag}}\right)^{*}/\hat{N_{i}}^{*}=\left(N_{v,i}^{\operatorname{eul}}\right)^{*}/\left(N_{i}^{\operatorname{eul}}\right)^{*}=1$,
$\beta=ab\hskip 5.69046pt\text{Days}^{-1}$, $\beta_{v}=acm\hskip
5.69046pt\text{Days}^{-1}$. In this case, systems (3) and (8) are inconsistent
with small relative error $\lVert
L-\left(D_{\delta}^{\operatorname{lag}}P^{-1}-D_{\delta}^{\operatorname{eul}}\right)\rVert_{F}/\lVert
D_{\delta}^{\operatorname{lag}}P^{-1}-D_{\delta}^{\operatorname{eul}}\rVert_{F}=0.018$,
where $L$ is the solution of (20). In addition, the Lagrangian basic
reproduction number $\mathcal{R}_{0}^{\operatorname{lag}}=1.7049$ is slightly
larger than the Eulerian basic reproduction number
$\mathcal{R}_{0}^{\operatorname{eul}}=1.6997$. In conclusion, the basic
reproduction numbers for both systems are similar, even though the systems are
not consistent (with small relative error).

###  5.2 Dengue in Brazil

Iggidr et al [27] use a Lagrangian framework to model dengue in eight
locations forming the metropolitan area of Rio de Janeiro. The host density in
each location was determined from the national census, and a 20% host-vector
ratio was assumed in their simulations. The mixing matrix was estimated from
data provided by the transport authority of Rio de Janerio (see Appendix B of
[27]).

Consider the parameter values $\delta_{i}=1/10.5\hskip
5.69046pt\text{Days}^{-1}$, $\delta_{v,i}=1/5.5\hskip
5.69046pt\text{Days}^{-1}$, $m\mathrel{\mathop{\mathchar
58\relax}}=\left(N_{v,i}^{\operatorname{eul}}\right)^{*}/\left(N_{i}^{\operatorname{eul}}\right)^{*}=1$,
$\beta_{i}=0.3750\hskip 5.69046pt\text{Days}^{-1}$, $\beta_{v,i}=0.3750\hskip
5.69046pt\text{Days}^{-1}$ used in [32]. In this case, systems (3) and (8) are
inconsistent albeit with relative error $\lVert
L-\left(D_{\delta}^{\operatorname{lag}}P^{-1}-D_{\delta}^{\operatorname{eul}}\right)\rVert_{F}/\lVert
D_{\delta}^{\operatorname{lag}}P^{-1}-D_{\delta}^{\operatorname{eul}}\rVert_{F}=0.0277$,
where $L$ is the solution of (20). Moreover, the Lagrangian basic reproduction
number $\mathcal{R}_{0}^{\operatorname{lag}}=2.8586$ is nearly equal to the
Eulerian basic reproduction number
$\mathcal{R}_{0}^{\operatorname{eul}}=2.8498$. Thus despite the systems being
inconsistent, similar reproduction numbers are obtained for the Eulerian and
Lagrangian modeling frameworks, suggesting flexibility in using either
framework in terms of the domain basic reproduction number.

###  5.3 Dengue in Iquitos, hypothetical mixing matrix

In [55], neck trap GPS devices were used to register movements of a carpenter
and mototaxi driver in Iquitos, Peru for approximately two weeks. The spatial
scale of movement here is small, suggesting a Lagrangian framework. The second
column of Table 2 (‘Time’) corresponds to the time spent by the carpenter in
$4$ houses (P1, P2, P4, P5) during 14 days according to the GPS data. The
third column of Table 2 (‘Proportions’) corresponds to the ratios between the
times in the second column and the total hours in 14 days.

House | Time (hours) | Proportion
---|---|---
P1 | 12.1 | 0.036
P2 | 3 | 0.008
P4 | 5.4 | 0.016
P5 | 18 | 0.053
Table 2: Time and proportion of time spent in the most frequented four visited
houses (other than home) by a carpenter during two weeks in Iquitos [55].

The mixing matrix $P$ is a hypothetical arrangement based on the proportions
in Table 2.

$P=\begin{pmatrix}0.887&0.036&0.036&0.036&0.036\\\
0.036&0.887&0.008&0.008&0.008\\\ 0.008&0.008&0.887&0.016&0.016\\\
0.016&0.016&0.016&0.887&0.053\\\ 0.053&0.053&0.053&0.053&0.887\end{pmatrix},$
(34)

$L=\begin{pmatrix}0.0241&-0.0077&-0.0077&-0.0077&-0.0077\\\
-0.0082&0.0236&-0.0014&-0.0014&-0.0014\\\
-0.0015&-0.0015&0.0235&-0.0034&-0.0034\\\
-0.0028&-0.0028&-0.0028&0.0241&-0.0121\\\
-0.0116&-0.0116&-0.0116&-0.0116&0.0246\end{pmatrix}.$

In this example we use the parameters $\delta_{i}=1/10.5\hskip
5.69046pt\text{Days}^{-1}$, $\delta_{v,i}=1/5.5\hskip
5.69046pt\text{Days}^{-1}$, $m\mathrel{\mathop{\mathchar
58\relax}}=\left(N_{v,i}^{\operatorname{lag}}\right)^{*}/\hat{N_{i}}^{*}=\left(N_{v,i}^{\operatorname{eul}}\right)^{*}/\left(N_{i}^{\operatorname{eul}}\right)^{*}=1$,
$\beta_{i}=0.3750\hskip 5.69046pt\text{Days}^{-1}$, $\beta_{v,i}=0.3750\hskip
5.69046pt\text{Days}^{-1}$ given in [32]. Here we obtain that systems (3) and
(8) are consistent, i.e, $\lVert
L-\left(D_{\delta}^{\operatorname{lag}}P^{-1}-D_{\delta}^{\operatorname{eul}}\right)\rVert_{F}/\lVert
D_{\delta}^{\operatorname{lag}}P^{-1}-D_{\delta}^{\operatorname{eul}}\rVert_{F}=0$.
Furthermore, the Lagrangian basic reproduction number
$\mathcal{R}_{0}^{\operatorname{lag}}=2.8536$ is approximately the Eulerian
basic reproduction number $\mathcal{R}_{0}^{\operatorname{eul}}=2.8498$.

###  5.4 Migratory hosts, hypothetical mixing matrix

We consider the hypothetical $3\times 3$ mixing matrix (32) from Section 3.4,
where the hosts in all the three patches spend most of their time in patch
$1$. A mixing matrix such as (32) could correspond to the movement of
migratory hosts that do not have a sense of home. In this situation an
Eulerian modeling framework is natural. Specifically, consider systems (3) and
(8) with $P$ and $L$ from (32) and patch parameters corresponding to West Nile
Virus in migratory birds, as used in [9]. These parameters are
$\delta_{i}=0.2222\hskip 5.69046pt\text{Days}^{-1}$,
$\delta_{v,i}=0.0666\hskip 5.69046pt\text{Days}^{-1}$,
$m\mathrel{\mathop{\mathchar
58\relax}}=\left(N_{v,i}^{\operatorname{lag}}\right)^{*}/\hat{N_{i}}^{*}=\left(N_{v,i}^{\operatorname{eul}}\right)^{*}/\left(N_{i}^{\operatorname{eul}}\right)^{*}=1$,
$\beta_{i}=0.2479\hskip 5.69046pt\text{Days}^{-1}$, $\beta_{v,i}=0.2479\hskip
5.69046pt\text{Days}^{-1}$. For this example, systems (3) and (8) with the
largest possible relative error, ($L=0$ the solution to the minimization
problem (20)). The Lagrangian basic reproduction number
$\mathcal{R}_{0}^{\operatorname{lag}}=1.5796$ is approximately $21\%$ larger
than the Eulerian basic reproduction number
$\mathcal{R}_{0}^{\operatorname{eul}}=1.3135$. In conclusion, in this example
the systems are inconsistent and the connected movement network for the
Lagrangian system is not reflected in the disconnected movement network for
the Eulerian system. In addition, the difference in the values of the basic
reproduction numbers from both systems may be significant.

Table 3 summarizes the comparison between the estimated reproduction numbers
in Sections 5.1, 5.2, 5.3 and 5.4, showing that in the inconsistent cases the
Lagrangian model gives a larger basic reproduction number, and the differences
${\cal R}_{0}^{\operatorname{lag}}-{\cal R}_{0}^{\operatorname{eul}}$
increases as the error increases for these examples. For each of these
examples, the basic reproduction numbers for the Eulerian and Lagrangian
frameworks are similar to one another. By contrast, in Section 5.4, we have an
example where host movement in Eulerian system that we obtain from (20) is
totally disconnected (i.e., $L=0$), and the difference between the basic
reproduction numbers may be significant.

Example | $E_{r}$ | $\mathcal{R}_{0}^{\operatorname{lag}}$ | $\mathcal{R}_{0}^{\operatorname{eul}}$
---|---|---|---
Malaria data in Namibia in [44] | 0.018 | 1.7049 | 1.6997
Brazil transportation data in [27] | 0.0277 | 2.8586 | 2.8498
Dengue in Iquitos, hypothetical $P$ | 0 | 2.8536 | 2.8498
Migratory hosts, hypothetical $P$ | 1 | 1.5796 | 1.3135
Table 3:  Relative error $E_{r}=\lVert
L-\left(D_{\delta}^{\operatorname{lag}}P^{-1}-D_{\delta}^{\operatorname{eul}}\right)\rVert_{F}/\lVert
D_{\delta}^{\operatorname{lag}}P^{-1}-D_{\delta}^{\operatorname{eul}}\rVert_{F}$
and values of the basic reproduction numbers
$\mathcal{R}_{0}^{\operatorname{lag}}$ and
$\mathcal{R}_{0}^{\operatorname{eul}}$ in the examples of Sections 5.1, 5.2,
5.3 and 5.4.

##  6 Discussion

Lagrangian and Eulerian approaches are important modeling tools for studying
the effects of heterogeneity and movement in disease dynamics [16]. We have
presented an approach for relating the Eulerian and Lagrangian systems through
a fundamental matrix by matching the time that infectious individuals reside
in other patches. We define the Eulerian and Lagrangian systems to be
consistent when the fundamental matrices match, in the sense that the minimum
value of the optimization problem (20) is zero.

As the star graph example in Section 3.3 and mixing matrix in Section 5.4
show, both consistency and inconsistency between the two frameworks is
possible. While we do not have a complete characterization of when the
Eulerian and Lagrangian frameworks are consistent, Proposition 1 gives a
sufficient condition. Specifically, Proposition 1 gives intervals
$[p_{*},p^{*}]$ such that if all the off-diagonal elements of the mixing
matrix are in $[p_{*},p^{*}]$, then the systems are consistent. The upper
bound in the sufficiency criterion can be interpreted in terms of Lagrangian
models being suitable for situations where individuals commute from a
distinguished home location. This setting often corresponds to individuals
spending the majority of their time in the home location, meaning that the
off-diagonal entries of the mixing matrix are small [16]. Inconsistent
examples with large off-diagonal elements of the mixing matrix such as in
Section 5.4 thus conflict with the sense of home that Lagrangian models try to
capture. Inconsistency is also possible when the off-diagonal entries are
small, as for the mixing matrix in (24). Identifying additional necessary
criteria for consistency is an area for future work.

In [27] it is discussed how to go from an Eulerian to a Lagrangian framework
where the movement rates are relatively larger than the removal rates. The
Eulerian framework considered in [27] is different from the Eulerian framework
considered in this paper. Specifically, [27] use an Eulerian framework with
$n^{2}$ variables corresponding to residents of patch $i$ that are currently
located in patch $j$. For example, $S^{h}_{ij}(t)$ represents the number of
susceptible hosts whose home is patch $i$ and are in patch $j$ at time $t$.
The movement rates corresponding to $S^{h}_{ij}$ are $m_{kj}^{i}$ for $j\neq
k$ (movement from $j$ to $k$), from which we can define a graph Laplacian
$L_{i}$ for each home patch $i$. A mixing matrix $P$ can be obtained from
$L_{1},\ldots,L_{n}$ under the assumption that the movement rates are large
compared to the removal rates. This is different than the Eulerian framework
we consider, where the movement rates are captured by a single graph Laplacian
$L^{\operatorname{X}}$ for every state
$\operatorname{X}\in\\{\operatorname{S},\operatorname{I},\operatorname{R}\\}$.
In addition, the consistency definition here requires the off-diagonal entries
of $D_{\delta}^{\operatorname{lag}}P^{-1}-D_{\delta}^{\operatorname{eul}}$ to
be non-positive. Therefore, for given $P$ and
$D_{\delta}^{\operatorname{eul}}$, consistency depends only upon the sign of
the off-diagonal elements of $P^{-1}$ and does not depend on the removal rates
in $D_{\delta}^{\operatorname{eul}}$. In consequence, in contrast to the time
scales assumption in [27], the concept of consistency that we present does not
depend on the relative timescales of movement to removal.

The functional implications of using a Lagrangian versus an Eulerian approach
are important to consider. We find that the domain ${\cal R}_{0}$ values are
similar under various scenarios when the two approaches are consistent. In the
homogeneous consistent case (Section 4.2), we obtain explicit bounds
(Proposition 2) for the difference in ${\cal R}_{0}$ between the Eulerian and
Lagrangian frameworks. Furthermore, the behavior of the Lagrangian basic
reproduction number in Proposition 2 of Section 4.2 is consistent with studies
such as [32] (in the sense of attaining a minimum value when $p_{12}=p_{21}$,
see Fig. 4 of [32]).

Although there is inconsistency in the examples of Sections 3.4, 5.1, 5.4, we
observe that the obtained values of basic reproduction number are still alike.
This suggests using (20) to relate mixing matrices such as those given in [27,
44] to Eulerian systems, and then studying the reproduction number of the
resulting Eulerian systems using techniques such as those in [28, 53], can be
an effective approach. The graph Laplacian matrix $L$ obtained from the
optimization problem (20) (as in Sections 5.1 and 5.2) allows series
expansions for the basic reproduction number, and derivation of important
quantities such as the absorption inverse $L^{d}$ [28] for analyzing the
mobility network. For example, $L^{d}$ captures the effect of absorption (that
are the removal rates in this case) on the movement network and also has
applications in community detection and node centrality [8, 28], leading to
extensions of the analysis of the mixing matrix $P$.

In contrast to ${\cal R}_{0}$, we observe significant differences in outbreak
size between the Eulerian and Lagrangian approaches. Indeed, differences in
outbreak size can be large not only when the systems are inconsistent (e.g.
Figure 3), but for the consistent case as well. Therefore, care must be taken
when going from one approach to the other. It would be useful to have a bound
for the differences in outbreak size between the two approaches. Analytical
results quantifying how different the final sizes are is an area for future
work.

## Declaration of Competing Interest

None.

## Acknowledgements

The authors are grateful to Chris Cosner for feedback on an early version of
this manuscript. This work was supported by the National Science Foundation
(DMS 1814737, DMS 1440386), and the Fullbright International Fellow Program.

##  7 Appendix

###  7.1 Basic reproduction number

In this section of the Appendix we use the next generation matrix approach
[54] to derive the expressions for the basic reproduction numbers (7) and (12)
of systems (3) and (8) respectively.

We first compute the next generation matrix
$\left(F^{\operatorname{lag}}\right)\left(V^{\operatorname{lag}}\right)^{-1}$
of the Lagrangian system. Let us consider the equations

$\displaystyle\dot{I_{i}}$
$\displaystyle=\sum_{j=1}^{n}\beta_{j}p_{ji}\frac{S_{i}}{N_{i}}I_{v,j}-\left(\gamma_{i}^{\operatorname{lag}}+\mu_{i}^{\operatorname{lag}}\right)I_{i}$
$\displaystyle\dot{I}_{v,i}$
$\displaystyle=\beta_{v,i}\frac{\sum_{j=1}^{n}p_{ij}I_{j}}{\sum_{j=1}^{n}p_{ij}N_{j}}S_{v,i}+\sum_{j=1}^{n}m_{ij}^{v}I_{v,j}-\sum_{j=1}^{n}m_{ji}^{v}I_{v,i}-\mu_{v,i}I_{v,i}$
(35)

corresponding to the infectious compartments of system (3). From (7.1) we
define the function
$\mathcal{F}^{\operatorname{lag}}\mathrel{\mathop{\mathchar
58\relax}}\mathbb{R}^{2n}\rightarrow\mathbb{R}^{2n}$ by

$\mathcal{F}_{i}^{\operatorname{lag}}(I_{1},\ldots,I_{n},I_{v,1},\ldots,I_{v,n})\mathrel{\mathop{\mathchar
58\relax}}=\sum_{j=1}^{n}\beta_{j}p_{ji}\frac{S_{i}}{N_{i}}I_{v,j}\,,$

for $i=1,\ldots,n$, and

$\mathcal{F}_{i}^{\operatorname{lag}}(I_{1},\ldots,I_{n},I_{v,1},\ldots,I_{v,n})\mathrel{\mathop{\mathchar
58\relax}}=\beta_{v,i}\frac{\sum_{j=1}^{n}p_{ij}I_{j}}{\sum_{j=1}^{n}p_{ij}N_{j}}S_{v,i}\,,$

for $i=n+1,\ldots,2n$. The DFE of the Lagrangian systems is determined by
$\left(S_{i}^{\operatorname{lag}}\right)^{*}=\left(N_{i}^{\operatorname{lag}}\right)^{*}$
and
$\left(S_{v,i}^{\operatorname{lag}}\right)^{*}=\left(N_{v,i}^{\operatorname{lag}}\right)^{*}$
defined by (4). We then define the Jacobian matrix

$F^{\operatorname{lag}}\mathrel{\mathop{\mathchar
58\relax}}=\left.\partial\mathcal{F}^{\operatorname{lag}}/\partial(I_{1},\ldots,I_{n},I_{v,1},\ldots,I_{v,n})\right|_{DFE}\,.$

We have that
$\left.\left(\partial\mathcal{F}_{i}^{\operatorname{lag}}/\partial
I_{v,j}\right)\right|_{DFE}=\beta_{j}p_{ji}$, so the upper-right block of
$F^{\operatorname{lag}}$ is $P^{\operatorname{T}}D_{\beta}$, where
$D_{\beta}\mathrel{\mathop{\mathchar
58\relax}}=\operatorname{diag}\\{\beta_{i}\\}$. We also have that
$\left.\left(\partial\mathcal{F}_{n+i}^{\operatorname{lag}}/\partial
I_{j}\right)\right|_{DFE}=\beta_{v,i}\left(N_{v,i}^{*}/\hat{N_{i}}^{*}\right)p_{ij}$,
so the lower-left block of $F^{\operatorname{lag}}$ is
$D_{\beta_{v}}^{\operatorname{lag}}P$, where
$\hat{N_{i}}^{*}\mathrel{\mathop{\mathchar
58\relax}}=\sum_{j=1}^{n}p_{ij}\left(N_{j}^{\operatorname{lag}}\right)^{*}$
and $D_{\beta_{v}}^{\operatorname{lag}}\mathrel{\mathop{\mathchar
58\relax}}=\operatorname{diag}\left\\{\beta_{v,i}N_{v,i}^{*}/\hat{N_{i}}^{*}\right\\}$.
Therefore, we have

$F^{\operatorname{lag}}=\begin{pmatrix}0&P^{\operatorname{T}}D_{\beta}\\\
D_{\beta_{v}}^{\operatorname{lag}}P&0\end{pmatrix}.$

Similarly, we define the function
$\mathcal{V}^{\operatorname{lag}}\mathrel{\mathop{\mathchar
58\relax}}\mathbb{R}^{2n}\rightarrow\mathbb{R}^{2n}$ by

$\mathcal{V}_{i}^{\operatorname{lag}}\mathrel{\mathop{\mathchar
58\relax}}=\left(\gamma_{i}^{\operatorname{lag}}+\mu_{i}^{\operatorname{lag}}\right)I_{i}\,,$

for $i=1,\ldots,n$, and

$\mathcal{V}_{i}^{\operatorname{lag}}\mathrel{\mathop{\mathchar
58\relax}}=-\sum_{j=1}^{n}m_{ij}^{v}I_{v,j}+\sum_{j=1}^{n}m_{ji}^{v}I_{v,i}+\mu_{v,i}I_{v,i}\,,$

for $i=n+1,\ldots,2n$. We also define the Jacobian matrix

$V^{\operatorname{lag}}\mathrel{\mathop{\mathchar
58\relax}}=\left.\partial\mathcal{V}^{\operatorname{lag}}/\partial(I_{1}.\ldots,I_{n},I_{v,1},\ldots,I_{v,n})\right|_{DFE}\,.$

If $L_{v}$ is the graph Laplacian of the vector movement (with adjacency
matrix $M^{v}=\left(m_{ij}^{v}\right)_{i,j\leq n}$, see (2)) and
$D_{\delta}^{\operatorname{lag}}\mathrel{\mathop{\mathchar
58\relax}}=\operatorname{diag}\left\\{\gamma_{i}^{\operatorname{lag}}+\mu_{i}^{\operatorname{lag}}\right\\}$,
then the upper-left block of $V^{\operatorname{lag}}$ is
$D_{\delta}^{\operatorname{lag}}$ and the lower-right block of
$V^{\operatorname{lag}}$ is $G_{v}=L_{v}+D_{\delta_{v}}$. Therefore,

$V^{\operatorname{lag}}=\begin{pmatrix}D_{\delta}^{\operatorname{lag}}&0\\\
0&G_{v}\end{pmatrix}.$

Consequently,

$\left(F^{\operatorname{lag}}\right)\left(V^{\operatorname{lag}}\right)^{-1}=\begin{pmatrix}0&P^{\operatorname{T}}D_{\beta}G_{v}^{-1}\\\
D_{\beta_{v}}^{\operatorname{lag}}P\left(D_{\delta}^{\operatorname{lag}}\right)^{-1}&0\end{pmatrix}$

and

$(\mathcal{R}_{0}^{\operatorname{lag}})^{2}=\rho\left(P^{\operatorname{T}}D_{\beta}G_{v}^{-1}D_{\beta_{v}}^{\operatorname{lag}}P\left(D_{\delta}^{\operatorname{lag}}\right)^{-1}\right)\,.$

We now compute the next generation matrix
$\left(F^{\operatorname{eul}}\right)\left(V^{\operatorname{eul}}\right)^{-1}$
of system (8). The equations of the infectious compartments of system (8) are

$\displaystyle\dot{I_{i}}$
$\displaystyle=\beta_{i}\frac{S_{i}}{N_{i}}I_{v,i}+\sum_{j=1}^{n}m_{ij}^{\operatorname{I}}I_{j}-\sum_{j=1}^{n}m_{ji}^{\operatorname{I}}I_{i}-\left(\gamma_{i}^{\operatorname{eul}}+\mu_{i}^{\operatorname{eul}}\right)I_{i}$
$\displaystyle\dot{I_{v,i}}$
$\displaystyle=\beta_{v,i}\frac{I_{i}}{N_{i}}S_{v,i}+\sum_{j=1}^{n}m_{ij}^{v}I_{v,j}-\sum_{j=1}^{n}m_{ji}^{v}I_{v,i}-\mu_{v,i}I_{v,i}\,.$
(36)

From the equations in (7.1) we define the function

$\mathcal{F}^{\operatorname{eul}}(I_{1},\ldots,I_{n},I_{v,1},\ldots,I_{v,n})\mathrel{\mathop{\mathchar
58\relax}}=\beta_{i}\frac{S_{i}}{N_{i}}I_{v,i}\,,$

for $i=1,\ldots,n$, and

$\mathcal{F}^{\operatorname{eul}}(I_{1},\ldots,I_{n},I_{v,1},\ldots,I_{v,n})=\beta_{v,i}\frac{I_{i}}{N_{i}}S_{v,i}\,,$

for $i=n+1,\ldots,2n$.

Using the DFE (9) of system (8), we have that
$\left.\left(\partial\mathcal{F}_{i}^{\operatorname{eul}}/\partial
I_{v,j}\right)\right|_{DFE}=\beta_{i}$ and
$\left.\left(\partial\mathcal{F}_{n+i}^{\operatorname{eul}}/\partial\mathcal{F}_{j}^{\operatorname{eul}}\right)\right|_{DFE}=\beta_{v,i}N_{v,i}^{*}/\left(N_{i}^{\operatorname{eul}}\right)^{*}\,.$
Therefore, if

$F^{\operatorname{eul}}\mathrel{\mathop{\mathchar
58\relax}}=\left.\partial\mathcal{F}^{\operatorname{eul}}/\partial(I_{1},\ldots,I_{n},I_{v,1},\ldots,I_{v,n})\right|_{DFE}\,,$

we then have

$F^{\operatorname{eul}}=\begin{pmatrix}0&D_{\beta}\\\
D_{\beta_{v}}^{\operatorname{eul}}&0\end{pmatrix}\,,$

where $D_{\beta}\mathrel{\mathop{\mathchar
58\relax}}=\operatorname{diag}\\{\beta_{i}\\}$ and
$D_{\beta_{v}}^{\operatorname{eul}}\mathrel{\mathop{\mathchar
58\relax}}=\operatorname{diag}\left\\{\beta_{v,i}N_{v,i}^{*}/\left(N_{i}^{\operatorname{eul}}\right)^{*}\right\\}$.
Similarly, we define the function

$\mathcal{V}^{\operatorname{eul}}(I_{1},\ldots,I_{n},I_{v,1},\ldots,I_{v,n})\mathrel{\mathop{\mathchar
58\relax}}=-\sum_{j=1}^{n}m_{ij}^{\operatorname{I}}I_{j}+\sum_{j=1}^{n}m_{ji}^{\operatorname{I}}I_{i}+\left(\gamma_{i}^{\operatorname{eul}}+\mu_{i}^{\operatorname{eul}}\right)I_{i}\,,$

for $i=1,\ldots,n$, and

$\mathcal{V}^{\operatorname{eul}}(I_{1},\ldots,I_{n},I_{v,1},\ldots,I_{v,n})\mathrel{\mathop{\mathchar
58\relax}}=-\sum_{j=1}^{n}m_{ij}^{v}I_{v,j}+\sum_{j=1}^{n}m_{ji}^{v}I_{v,i}+\mu_{v,i}I_{v,i}\,,$

for $i=n+1,\ldots,2n$. Therefore, if

$V^{\operatorname{eul}}\mathrel{\mathop{\mathchar
58\relax}}=\left.\partial\mathcal{V}^{\operatorname{eul}}/\partial(I_{1},\ldots,I_{n},I_{v,1}\ldots,I_{v,n})\right|_{DFE}\,,$

we then have

$V^{\operatorname{eul}}=\begin{pmatrix}G&0\\\ 0&G_{v}\end{pmatrix}\,,$ (37)

where $G\mathrel{\mathop{\mathchar
58\relax}}=L^{\operatorname{I}}+D_{\delta}^{\operatorname{eul}}$,
$L^{\operatorname{I}}$ is the graph Laplacian of the host movement [with
adjacency matrix $M^{\operatorname{I}}=(m_{ij}^{\operatorname{I}})_{i,j\leq
n}$, see (2)], $D_{\delta}^{\operatorname{eul}}\mathrel{\mathop{\mathchar
58\relax}}=\operatorname{diag}\\{\gamma_{i}^{\operatorname{eul}}+\mu_{i}^{\operatorname{eul}}\\}$,
and $G_{v}=L_{v}+D_{\delta_{v}}$. In consequence,

$\left(F^{\operatorname{eul}}\right)\left(V^{\operatorname{eul}}\right)^{-1}=\begin{pmatrix}0&D_{\beta}G_{v}^{-1}\\\
D_{\beta_{v}}^{\operatorname{eul}}G^{-1}&0\end{pmatrix}$

and

$\left(\mathcal{R}_{0}^{\operatorname{eul}}\right)^{2}=\rho\left(D_{\beta}G_{v}^{-1}D_{\beta_{v}}^{\operatorname{eul}}G^{-1}\right).$

###  7.2 Comparison of basic reproduction numbers

In this section we prove Proposition 2 of Section 4.2. Let
$\beta_{v,1}=\beta_{v,2}$, $N_{v}^{*}\mathrel{\mathop{\mathchar
58\relax}}=N_{v,1}^{*}=N_{v,2}^{*}$,
$\left(N^{\operatorname{eul}}\right)^{*}\mathrel{\mathop{\mathchar
58\relax}}=\left(N_{1}^{\operatorname{eul}}\right)^{*}=\left(N_{2}^{\operatorname{eul}}\right)^{*}$,
and
$\beta_{v}=\beta_{v,1}N_{v,1}^{*}/\left(N_{1}^{\operatorname{eul}}\right)^{*}=\beta_{v,2}N_{v,2}^{*}/\left(N_{2}^{\operatorname{eul}}\right)^{*}$.
Define $\left(N_{1}^{\operatorname{lag}}\right)^{*}$ and
$\left(N_{2}^{\operatorname{lag}}\right)^{*}$ such that
$\left(N_{1}^{\operatorname{eul}}\right)^{*}=p_{11}\left(N_{1}^{\operatorname{lag}}\right)^{*}+p_{12}\left(N_{2}^{\operatorname{lag}}\right)^{*}$
and
$\left(N_{2}^{\operatorname{eul}}\right)^{*}=p_{21}\left(N_{1}^{\operatorname{lag}}\right)^{*}+p_{22}\left(N_{2}^{\operatorname{lag}}\right)^{*}$.
Hence A14 holds, and
$\delta=\delta_{1}^{\operatorname{eul}}=\delta_{2}^{\operatorname{eul}}$, so
we also get
$\delta=\delta_{1}^{\operatorname{lag}}=\delta_{2}^{\operatorname{lag}}$ by
A9. Let $D_{\beta}=\beta I$, $D_{\beta_{v}}=\beta_{v}I$,
$D_{\delta}^{\operatorname{eul}}=D_{\delta}^{\operatorname{lag}}=D_{\delta}\mathrel{\mathop{\mathchar
58\relax}}=\delta I$, $D_{\delta_{v}}=\delta_{v}I$,
$L_{v}=m_{v}\begin{pmatrix}1&-1\\\
-1&1\end{pmatrix},P=\begin{pmatrix}1-p_{21}&p_{12}\\\
p_{21}&1-p_{12}\end{pmatrix}$.

We fix all the parameters except $p_{12}$ (we obtain analogous and symmetric
results if we fix $p_{21}$). From (7) and (12), we obtain that
$\left(\mathcal{R}_{0}^{\operatorname{lag}}\right)^{2}=\rho\left(P^{\operatorname{T}}D_{\beta}G_{v}^{-1}D_{\beta_{v}}P\left(D_{\delta}\right)^{-1}\right)$
and
$\left(\mathcal{R}_{0}^{\operatorname{eul}}\right)^{2}=\rho\left(D_{\beta}G_{v}^{-1}D_{\beta_{v}}P(D_{\delta}\right)^{-1})$.
Therefore, if we define $\mathcal{M}\mathrel{\mathop{\mathchar
58\relax}}=D_{\beta}G_{v}^{-1}D_{\beta_{v}}P\left(D_{\delta}^{\operatorname{lag}}\right)^{-1}$,
we get

$\left(\mathcal{R}_{0}^{\operatorname{lag}}\right)^{2}=\rho(P^{\operatorname{T}}\mathcal{M})\quad\textrm{and}\quad\left(\mathcal{R}_{0}^{\operatorname{eul}}\right)^{2}=\rho(\mathcal{M})\,.$
(38)

The following is the statement of the proposition.

######  Proposition 2.

Assume A1–A10 and suppose that $D_{\beta}=\beta I$,
$D_{\beta_{v}}=\beta_{v}I$,
$D_{\delta}^{\operatorname{eul}}=D_{\delta}^{\operatorname{lag}}=D_{\delta}\mathrel{\mathop{\mathchar
58\relax}}=\delta I$, $D_{\delta_{v}}=\delta_{v}I$ and
$m_{12}^{v}=m_{21}^{v}=m_{v}$. Define $\mathcal{M}\mathrel{\mathop{\mathchar
58\relax}}=D_{\beta}G_{v}^{-1}D_{\beta_{v}}P\left(D_{\delta}^{\operatorname{lag}}\right)^{-1}$
and fix $p_{21}$ in the interval $(0,1/2)$. Then $\rho(P^{\operatorname{T}}M)$
is a function of $p_{12}$ where $0<p_{12}<1/2$ and we have that:

  1. a)

$\rho(\mathcal{M})=(\beta\beta_{v})/(\delta\delta_{v})$ is constant on
$0<p_{12}<1/2$.

  2. b)

$\rho(P^{t}\mathcal{M})$ is decreasing on $(0,p_{21})$, increasing on
$(p_{21},1/2)$ and attains one absolute minimum over $(0,1/2)$ with value
$(\beta\beta_{v})/(\delta\delta_{v})$ at $p_{12}=p_{21}$.

  3. c)

In addition, we have the inequality

$\displaystyle\rho(P^{\operatorname{T}}\mathcal{M})-\rho(\mathcal{M})$
$\displaystyle\leq\frac{\beta\beta_{v}}{\delta\delta_{v}(2m_{v}+\delta_{v})}|p_{12}-p_{21}|(1-p_{12}-p_{21})\delta_{v}$
$\displaystyle\leq\frac{\beta\beta_{v}}{\delta\delta_{v}(2m_{v}+\delta_{v})}\max\left\\{p_{21}(1-p_{21})\delta_{v},(1/2-p_{21})^{2}\delta_{v}\right\\}$
$\displaystyle\leq\frac{\beta\beta_{v}}{\delta\delta_{v}(2m_{v}+\delta_{v})}\frac{\delta_{v}}{4}\,.$
(39)

###### Proof.

We have that

$\mathcal{M}=D_{\beta}G_{v}^{-1}D_{\beta_{v}}P\left(D_{\delta}^{\operatorname{lag}}\right)^{-1}=\frac{\beta\beta_{v}}{\delta\delta_{v}(2m_{v}+\delta_{v})}\begin{pmatrix}m_{v}+(1-p_{21})\delta_{v}&m_{v}+p_{12}\delta_{v}\\\
m_{v}+p_{21}\delta_{v}&m_{v}+(1-p_{12})\delta_{v}\end{pmatrix}\,,$

so the eigenvalues of $\mathcal{M}$ in this case are

$\frac{\beta\beta_{v}}{\delta\delta_{v}}\quad\textrm{and}\quad\frac{\beta\beta_{v}}{\delta\delta_{v}}\frac{(1-p_{12}-p_{21})\delta_{v}}{2m_{v}+\delta_{v}}.$

In consequence, we get
$\rho(\mathcal{M})=(\beta\beta_{v})/(\delta\delta_{v})$. We can also get
$\rho(P^{\operatorname{T}}\mathcal{M})$ explicitly by

$\displaystyle\rho(P^{\operatorname{T}}\mathcal{M})=\frac{\beta\beta_{v}}{\delta\delta_{v}(2m_{v}+\delta_{v})}\left[m_{v}+\delta_{v}-\left(p_{12}(1-p_{12})+p_{21}(1-p_{21})\right)\delta_{v}+\right.$
$\displaystyle\left.\sqrt{\left[m_{v}+\left(p_{12}(1-p_{12})+p_{21}(1-p_{21})\right)\delta_{v}\right]^{2}+\left[(p_{12}-p_{21})(1-p_{12}-p_{21})\delta_{v}\right]^{2}}\right]\,.$
(40)

From this equation, it follows that when $p_{12}=p_{21}$, we get

$\rho(P^{\operatorname{T}}\mathcal{M})=\frac{\beta\beta_{v}}{\delta\delta_{v}}=\rho(\mathcal{M})\,.$

Moreover, $\partial\rho(P^{\operatorname{T}}\mathcal{M})/\partial
p_{12}=\kappa\beta\beta_{v}/\left[\delta\delta_{v}(2m_{v}+\delta_{v})\right]$,
where

$\displaystyle\kappa=2p_{12}-1+\frac{m_{v}(1-2p_{21})+\delta_{v}\left[2p_{12}(1-p_{12})(1-2p_{12})+(p_{12}-p_{21})\left(2p_{12}(1-p_{12}-p_{21})+p_{12}+p_{21}\right)\right]}{\sqrt{\left[m_{v}+\left(p_{12}(1-p_{12})+p_{21}(1-p_{21})\right)\delta_{v}\right]^{2}+\left[(p_{12}-p_{21})(1-p_{12}-p_{21})\delta_{v}\right]^{2}}}\,.$
(41)

In particular, if $p_{12}=p_{21}$, then

$\frac{\partial\rho(P^{\operatorname{T}}\mathcal{M})}{\partial p_{12}}=0\,.$

Assume that $p_{12}>p_{21}$. Define

$\displaystyle\alpha$ $\displaystyle\mathrel{\mathop{\mathchar
58\relax}}=\sqrt{\left[m_{v}+\left(p_{12}(1-p_{12})+p_{21}(1-p_{21})\right)\delta_{v}\right]^{2}+\left[(p_{12}-p_{21})(1-p_{12}-p_{21})\delta_{v}\right]^{2}}$
$\displaystyle\leq\left[m_{v}+\left(p_{12}(1-p_{12})+p_{21}(1-p_{21})\right)\delta_{v}\right]+\left[(p_{12}-p_{21})(1-p_{12}-p_{21})\delta_{v}\right]\,,$

and

$\eta\mathrel{\mathop{\mathchar
58\relax}}=m_{v}(1-2p_{21})+\delta_{v}\left[2p_{12}(1-p_{12})(1-2p_{12})+(p_{12}-p_{21})\left(2p_{12}(1-p_{12}-p_{21})+p_{12}+p_{21}\right)\right]\,.$

We then have that

$\frac{\partial\rho(P^{\operatorname{T}}\mathcal{M})}{\partial
p_{12}}=\frac{\beta\beta_{v}}{\delta\delta_{v}(2m_{v}+\delta_{v})}\left[2p_{12}-1+\frac{\eta}{\alpha}\right]$

and

$\displaystyle(2p_{12}-1)\alpha+\eta$
$\displaystyle\geq(2p_{12}-1)\left[m_{v}+\left(p_{12}(1-p_{12})+p_{21}(1-p_{21})\right)\delta_{v}+(p_{12}-p_{21})(1-p_{12}-p_{21})\delta_{v}\right]+\eta$
$\displaystyle=2(p_{12}-p_{21})\left[m_{v}+\left((p_{12}-1)^{2}+p_{12}^{2}\right)\delta_{v}\right]>0\,.$

Therefore, $p_{12}>p_{21}$ implies that

$\frac{\partial\rho(P^{\operatorname{T}}\mathcal{M})}{\partial
p_{12}}=\frac{\beta\beta_{v}}{\delta\delta_{v}(2m_{v}+\delta_{v})}\frac{(2p_{12}-1)\alpha+\eta}{\alpha}>0\,.$

Now, let us assume that $p_{12}<p_{21}$. We then have that

$\alpha\geq m_{v}+\left(p_{12}(1-p_{12})+p_{21}(1-p_{21})\right)\delta_{v}$

and

$\displaystyle(2p_{12}-1)\alpha+\eta$
$\displaystyle\leq(2p_{12}-1)\left[m_{v}+\left(p_{12}(1-p_{12})+p_{21}(1-p_{21})\right)\delta_{v}\right]+\eta$
$\displaystyle=-(p_{21}-p_{12})(m_{v}+\delta_{v}(1+(p_{21}-p_{12})(1-2p_{12})))<0\,.$

Therefore, $p_{12}<p_{21}$ implies that

$\frac{\partial\rho(P^{\operatorname{T}}\mathcal{M})}{\partial
p_{12}}=\frac{\beta\beta_{v}}{\delta\delta_{v}(2m_{v}+\delta_{v})}\frac{(2p_{12}-1)\alpha+\eta}{\alpha}<0\,.$

Using that
$\alpha\leq\left[m_{v}+\left(p_{12}(1-p_{12})+p_{21}(1-p_{21})\right)\delta_{v}\right]+|p_{12}-p_{21}|(1-p_{12}-p_{21})\delta_{v}$,
we get

$\displaystyle\rho(P^{\operatorname{T}}\mathcal{M})-\rho(\mathcal{M})$
$\displaystyle\leq\frac{\beta\beta_{v}}{\delta\delta_{v}(2m_{v}+\delta_{v})}|p_{12}-p_{21}|(1-p_{12}-p_{21})\delta_{v}$
$\displaystyle\leq\frac{\beta\beta_{v}}{\delta\delta_{v}(2m_{v}+\delta_{v})}\max\left\\{p_{21}(1-p_{21})\delta_{v},(1/2-p_{21})^{2}\delta_{v}\right\\}$
$\displaystyle\leq\frac{\beta\beta_{v}}{\delta\delta_{v}(2m_{v}+\delta_{v})}\frac{\delta_{v}}{4}\,.$

∎

Let us try to get some intuition for the inequality ${\cal
R}_{0}^{lag}\geq{\cal R}_{0}^{eul}$ in the previous proposition. Suppose that
$p_{21}<p_{12}$ and all the other parameters are assumed to be as in
Proposition 2. Assume that the systems (3) and (8) are at the DFE and suppose
we introduce the same number of infectious hosts in both patches, say
$I_{h}=I_{h,1}=I_{h,2}$. From the last equation in (3) of the Lagrangian
system and A14, the rates at which vectors get infected in patches $1$ and $2$
are

$\beta_{v,1}\frac{p_{11}I_{h,1}+p_{12}I_{h,2}}{p_{11}\left(N_{1}^{\operatorname{lag}}\right)^{*}+p_{12}\left(N_{2}^{\operatorname{lag}}\right)^{*}}N_{v,1}^{*}=\beta_{v}I_{h}(1+p_{12}-p_{21})$

and

$\beta_{v,2}\frac{p_{21}I_{h,1}+p_{22}I_{h,2}}{p_{21}\left(N_{1}^{\operatorname{lag}}\right)^{*}+p_{22}\left(N_{2}^{\operatorname{lag}}\right)^{*}}N_{v,2}^{*}=\beta_{v}I_{h}(1+p_{21}-p_{12})$

respectively. Since $p_{12}>p_{21}$, the vector infection rate in patch $1$ is
greater than in patch $2$ (because $1+p_{12}-p_{21}>1>1+p_{21}-p_{12}$).
Therefore, over a period $\Delta t$, the number of infected vectors at patch
$1$, which is approximately $\Delta
I_{v,1}^{\operatorname{lag}}\mathrel{\mathop{\mathchar
58\relax}}=\beta_{v}I_{h}(1+p_{12}-p_{21})\Delta t$, would be larger than the
amount of infected vectors at patch $2$, which is approximately $\Delta
I_{v,2}^{\operatorname{lag}}\mathrel{\mathop{\mathchar
58\relax}}=\beta_{v}I_{h}(1+p_{21}-p_{12})\Delta t$. From the second equation
in (3), the amount of new host infections in patches $1$ and $2$ caused by the
new infected vectors in DFE would be

$\Delta I_{h,1}^{\operatorname{lag}}\mathrel{\mathop{\mathchar
58\relax}}=\beta_{1}p_{11}\left(\Delta
I_{v,1}^{\operatorname{lag}}\right)+\beta_{2}p_{21}\left(\Delta
I_{v,2}^{\operatorname{lag}}\right)=\beta\left(p_{11}\left(\Delta
I_{v,1}^{\operatorname{lag}}\right)+p_{21}\left(\Delta
I_{v,2}^{\operatorname{lag}}\right)\right)$

and

$\Delta I_{h,2}^{\operatorname{lag}}\mathrel{\mathop{\mathchar
58\relax}}=\beta_{1}p_{12}\left(\Delta
I_{v,1}^{\operatorname{lag}}\right)+\beta_{2}p_{22}\left(\Delta
I_{v,2}^{\operatorname{lag}}\right)=\beta\left(p_{12}\left(\Delta
I_{v,1}^{\operatorname{lag}}\right)+p_{22}\left(\Delta
I_{v,2}^{\operatorname{lag}}\right)\right)$

respectively. Therefore, the total amount of new infected hosts would be

$\Delta I_{h}^{\operatorname{lag}}\mathrel{\mathop{\mathchar 58\relax}}=\Delta
I_{h,1}^{\operatorname{lag}}+\Delta
I_{h,2}^{\operatorname{lag}}=\beta\left(\Delta
I_{v,1}^{\operatorname{lag}}+\Delta
I_{v,2}^{\operatorname{lag}}\right)+\beta(p_{12}-p_{21})\left(\Delta
I_{v,1}^{\operatorname{lag}}-\Delta I_{v,2}^{\operatorname{lag}}\right)\,.$
(42)

On the other hand, from the last equation in (8) of the Eulerian system, the
rates at which vectors get infected in patches $1$ and $2$ are

$\beta_{v,1}\frac{I_{h,1}}{\left(N_{1}^{\operatorname{eul}}\right)^{*}}N_{v,1}^{*}=\beta_{v}I_{h}\quad\text{and}\quad\beta_{v,2}\frac{I_{h,2}}{\left(N_{2}^{\operatorname{eul}}\right)^{*}}N_{v,2}^{*}=\beta_{v}I_{h}$

respectively. Therefore, over a period $\Delta t$, the amounts of new infected
vectors are $\Delta I_{v,1}^{\operatorname{eul}}=\Delta
I_{v,2}^{\operatorname{eul}}=\beta_{v}I_{h}\Delta t$. Notice that

$\Delta I_{v,1}^{\operatorname{lag}}+\Delta
I_{v,2}^{\operatorname{lag}}=2\beta_{v}I_{h}\Delta t=\Delta
I_{v,1}^{\operatorname{eul}}+\Delta I_{v,2}^{\operatorname{eul}}\,.$

From the second equation in (8) of the Eulerian system, the total amount of
new infected hosts is

$\Delta I_{h}^{\operatorname{eul}}=\beta(\Delta
I_{v,1}^{\operatorname{eul}}+\Delta I_{v,2}^{\operatorname{eul}})=\beta(\Delta
I_{v,1}^{\operatorname{lag}}+\Delta I_{v,2}^{\operatorname{lag}})\,.$ (43)

Since $p_{12}>p_{21}$ and $\Delta I_{v,1}^{\operatorname{lag}}>\Delta
I_{v,2}^{\operatorname{lag}}$, we get that the term
$\beta(p_{12}-p_{21})\left(\Delta I_{v,1}^{\operatorname{lag}}-\Delta
I_{v,2}^{\operatorname{lag}}\right)$ in (42) is positive, and then using (43),
we have that $\Delta I_{h}^{\operatorname{lag}}>\Delta
I_{h}^{\operatorname{eul}}$, which corresponds to ${\cal
R}_{0}^{\operatorname{lag}}\geq{\cal R}_{0}^{\operatorname{eul}}$. Note that
this implicitly relies upon the fact that the removal rates, the rates
$\beta_{v,i}$ and the rates $\beta_{i}$ are the same across patches. As we
observed in Figure 5, imposing different removal rates, for instance, can lead
to ${\cal R}_{0}^{\operatorname{eul}}>{\cal R}_{0}^{\operatorname{lag}}$.

## References

  * [1] M. A. Acevedo, O. Prosper, N. Ruktanonchai, T. Caughlin, M. Martcheva, C. W. Osenberg, and D. L. Smith.  Spatial Heterogeneity, Host Movement and Mosquito-Borne Disease Transmission.  PLoS One, 6(10), 2015.
  * [2] C. C. Aggarwal.  Data classification: algorithms and applications.  CRC press, 2014.
  * [3] L. J. S. Allen, B. M. Bolker, Y. Lou, and A. L. Nevai.  Asymptotic profiles of the steady states for an SIS epidemic patch model.  SIAM Journal on Applied Mathematics, 67(5):1283–1309, 2007.
  * [4] J. Arino, J. R. Davis, D. Hartley, R. Jordan, J. M. Miller, and P. van Den Driessche.  A multi-species epidemic model with spatial dynamics.  Mathematical Medicine and Biology, 22(2):129–142, 2005.
  * [5] J. Arino and P. van den Driessche.  A multi-city epidemic model.  Mathematical Population Studies, 10(3):175–193, 2003.
  * [6] J. Arino and P. van den Driessche.  Disease Spread in Metapopulations.  Fields Institute Communications, 48(2006):1–12, 2006.
  * [7] L. Bengtsson, X. Lu, A. Thorson, R. Garfield, and J. Von Schreeb.  Improved response to disasters and outbreaks by tracking population movements with mobile phone network data: a post-earthquake geospatial study in Haiti.  PLoS Medicine, 8(8), 2011.
  * [8] M. Benzi, P. Fika, and M. Mitrouli.  Graphs with absorption: Numerical methods for the absorption inverse and the computation of centrality measures.  Linear Algebra and its Applications, 574:123–152, 2019.
  * [9] L. D. Bergsman, J. M. Hyman, and C. A. Manore.  A mathematical model for the spread of West Nile virus in migratory and resident birds.  Mathematical Biosciences & Engineering, 13(2):401–424, 2016.
  * [10] A. Berman and R. J. Plemmons.  Nonnegative matrices in the mathematical sciences, volume 9.  SIAM, 1994.
  * [11] D. Bichara and C. Castillo-Chavez.  Vector-borne diseases models with residence times–A Lagrangian perspective.  Mathematical Biosciences, 281:128–138, 2016.
  * [12] A. Bomblies.  Agent-based modeling of malaria vectors: the importance of spatial simulation.  Parasites and Vectors, 308(7):426–435, 2014.
  * [13] C. L. Byrne.  A First Course in Optimization.  Chapman and Hall/CRC, 2014.
  * [14] D. Z. Childs and M. Boots.  The interaction of seasonal forcing and immunity and the resonance dynamics of malaria.  Journal of the Royal Society Interface, 7(43):309–319, 2009.
  * [15] F. R. K. Chung and F. C. Graham.  Spectral Graph Theory.  American Mathematical Soc., 1997.
  * [16] C. Cosner.  Models for the effects of host movement in vector-borne disease systems.  Mathematical Biosciences, 270:192–197, 2015.
  * [17] C. Cosner, J. C. Beier, R. S. Cantrell, D. Impoinvil, L. Kapitanski, M. D. Potts, A. Troyo, and S. Ruan.  The effects of human movement on the persistence of vector-borne diseases.  Journal of Theoretical Biology, 258(4):550–560, 2009.
  * [18] R. P. Dobrow.  Introduction to Stochastic Processes with R.  John Wiley & Sons, 2016.
  * [19] C. Dye and G. Hasibeder.  Population dynamics of mosquito-borne disease: effects of flies which bite some people more frequently than others.  Transactions of the Royal Society of Tropical Medicine and Hygiene, 80(1):69–77, 1986.
  * [20] G. Espana, C. Hogea, A. Guignard, A. Guignard, Q. Bosch, A. Morrison, Smith D. L., T. W. Scott, A. Schmidt, and A. Perkins.  Biased efficacy estimates in phase-iii dengue vaccine trials due to heterogeneous exposure and differential detectability of primary infections across trial arms.  PLoS One, 14(1), 2019.
  * [21] H. D. Gaff and L. J. Gross.  Modeling tick-borne disease: a metapopulation model.  Bulletin of Mathematical Biology, 69(1):265–288, 2007.
  * [22] D. Grunbaum and A. Okubo.  Modelling social animal aggregations.  Lecture Notes in Biomathematics, 100:296–325, 1994.
  * [23] S. Gueron, S. A. Levin, and D. I. Rubenstein.  The dynamics of herds: from individuals to aggregations.  Journal of Theoretical Biology, 182(1):85–98, 1996.
  * [24] G. Hasibeder and C. Dye.  Population dynamics of mosquito-borne disease: persistence in a completely heterogeneous environment.  Theoretical Population Biology, 33(1):31–53, 1988.
  * [25] N. A. Honório, W. C. Silva, P. J. Leite, J. M. Gonçalves, L. P. Lounibos, and R. Lourenço-de Oliveira.  Dispersal of Aedes aegypti and Aedes albopictus (Diptera: Culicidae) in an urban endemic dengue area in the State of Rio de Janeiro, Brazil.  Memórias do Instituto Oswaldo Cruz, 98(2):191–198, 2003.
  * [26] Y. H. Hsieh, P. van den Driessche, and L. Wang.  Impact of travel between patches for spatial spread of disease.  Bulletin of Mathematical Biology, 69(4):1355–1375, 2007.
  * [27] A. Iggidr, J. Koiller, M. L. F. Penna, G. Sallet, M. A. Silva, and M. O. Souza.  Vector borne diseases on an urban environment: the effects of heterogeneity and human circulation.  Ecological Complexity, 30:76–90, 2017.
  * [28] K. A. Jacobsen and J. H. Tien.  A generalized inverse for graphs with absorption.  Linear Algebra and its Applications, 537:118–147, 2018.
  * [29] A. Jindal and S. Rao.  Agent-Based Modeling and Simulation Mosquito-Borne Disease Transmission.  AAMAS ’17 Proceedings of the 16th Conference on Autonomous Agents and MultiAgent Systems, pages 426–435, 2017.
  * [30] M. Jovanovic and M. Krstic.  Stochastically perturbed vector-borne disease models with direct transmission.  Applied Mathematical Modeling, 36(11):5214–5228, 2012.
  * [31] E. Krause.  Fluid Mechanics I.  Springer, 2005.
  * [32] S. Lee and C. Castillo-Chavez.  The role of residence times in two-patch dengue transmission dynamics and optimal strategies.  Journal of Theoretical Biology, 374:152–164, 2015.
  * [33] J. Lessler, W. J. Edmunds, M. E. Halloran, T. D. Hollingsworth, and A. L. Lloyd.  Seven challenges for model-driven data collection in experimental and observational studies.  Epidemics, 10:78–82, 2015.
  * [34] M. Lewis, J. Rencławowic, and P. van den Driessche.  Traveling Waves and Spread Rates for a West Nile Virus Model.  Bulletin of Mathematical Biology, 68(1):3–23, 2006.
  * [35] R. Liu, J. Shuai, J. Wu, and H. Zhu.  Modeling spatial spread of West Nile virus and impact of directional dispersal of birds.  Mathematical Biosciences and Engineering, 3(1):145, 2006.
  * [36] M. Martcheva.  An introduction to mathematical epidemiology, volume 61.  Springer, 2015.
  * [37] A. Y. Ng, M. I. Jordan, and Y. Weiss.  On Spectral Clustering: Analysis and an algorithm.  In Advances in Neural Information Processing Systems, pages 849–856, 2002.
  * [38] J. Ochoa and L. Osorio.  Epidemiology of urban malaria in Quibdo, Choco.  Biomedica, 26(2):278–285, 2006.
  * [39] K. M. O’Reilly, R. Lowe, W. J. Edmunds, P. Mayaud, A. Kucharski, R. M. Eggo, S. Funk, D. Bhatia, K. Khan, M. U. G. Kraemer, et al.  Projecting the end of the Zika virus epidemic in Latin America: a modelling analysis.  BMC Medicine, 16(1):180, 2018.
  * [40] D. K. Pindolia, A. J. Garcia, A. Wesoloski, D. L. Smith, C. O. Buckee, A. M. Noor, R. W. Snow, and A. J. Tatem.  Human Movement Data for malaria control and elimination strategic planning.  Malaria Journal, 205(11), 2012.
  * [41] R. C. Reiner Jr., T. A. Perkins, C. M. Barker, T. Niu, L. F. Chaves, A. M. Ellis, D. B. George, A. Le Menach, J. R. C. Pulliam, D. Bisanzio, et al.  A systematic review of mathematical models of mosquito-borne pathogen transmission: 1970–2010.  Journal of The Royal Society Interface, 10(81):20120921, 2013.
  * [42] D. J. Rodríguez and L. Torres-Sorando.  Models of infectious diseases in spatially heterogeneous environments.  Bulletin of Mathematical Biology, 63(3):547–571, 2001.
  * [43] R. Ross.  An Application of the Theory of Probabilities to the Study of a priori Pathometry.Part I.  Proceedings of the Royal Society of London. Series A, Containing papers of a mathematical and physical character, 92(638):204–230, 1916.
  * [44] N. W. Ruktanonchai, P. DeLeenheer, A. J. Tatem, V. A. Alegana, T. T. Caughlin, E. zu Erbach-Schoenberg, C. Lourenço, C. W. Ruktanonchai, and D. L. Smith.  Identifying Malaria Transmission Foci for Elimination Using Human Mobility Data.  PLoS Computational Biology, 12(4):e1004846, 2016.
  * [45] M. W. Service.  Mosquito (Diptera: Culicidae) dispersal—the long and short of it.  Journal of Medical Entomology, 34(6):579–588, 1997.
  * [46] Z. Shuai and P. van den Driessche.  Global Stability of Infectious Disease Models Using Lyapunov Functions.  SIAM Journal on Applied Mathematics, 73:1513–1532, 2013.
  * [47] D. L. Smith, K. E. Battle, S. I. Hay, C. M. Barker, T. W. Scott, and F. E. McKenzie.  Ross, Macdonald, and a theory for the dynamics and control of mosquito-transmitted pathogens.  PLoS Pathogens, 8(4):e1002588, 2012.
  * [48] D. L. Smith and F. E. McKenzie.  Statics and dynamics of malaria infection in Anopheles mosquitoes.  Malaria Journal, 3(1):13, 2004.
  * [49] S. T. Stoddard, A. C. Morrison, G. M. Vazquez-Prokopec, V. P. Soldan, T. J. Kochel, U. Kitron, J. P. Elder, and T. W. Scott.  The role of human movement in the transmission of vector-borne pathogens.  PLoS Neglected Tropical Diseases, 3(7):e481, 2009.
  * [50] A. Tatem and D. L. Smith.  International population movements and regional Plasmodium falciparum malaria elimination strategies.  PNAS, 27(107):12222–12227, 2010.
  * [51] S. Tavare.  A note on finite homogeneous continuous-time Markov chains.  Biometrics, pages 831–834, 1979.
  * [52] M. C. Thomson, S. J. Connor, M. L. Quiñones, M. Jawara, J. Todd, and B. M. Greenwood.  Movement of Anopheles gambiae s.l. malaria vectors between villages in The Gambia.  Medical and Veterinary Entomology, 9(4):413–419, 1995.
  * [53] J. H. Tien, Z. Shuai, M. C. Eisenberg, and P. van den Driessche.  Disease invasion on community networks with environmental pathogen movement.  Journal of Mathematical Biology, 70(5):1065–1092, 2015.
  * [54] P. van den Driessche and J. Watmough.  Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission.  Mathematical Biosciences, 180(1-2):29–48, 2002.
  * [55] G. M. Vazquez-Prokopec, S. T. Stoddard, V. Paz-Soldan, A. C. Morrison, J. P. Elder, T. J. Kochel, T. W. Scott, and U. Kitron.  Usefulness of commercially available GPS data-loggers for tracking human movement and exposure to dengue virus.  International Journal of Health Geographics, 8(1):68, 2009.
  * [56] U. Von Luxburg.  A Tutorial on Spectral Clustering.  Statistics and Computing, 17(4):395–416, 2007.
  * [57] D. Wanduku and G. S. Ladde.  Global properties of a two-scale network stochastic delayed human epidemic dynamic model.  Nonlinear Analysis: Real World Applications, 13(2):794–816, 2012\.
  * [58] A. Wesolowski, N. Eagle, A. Tatem, D. Smith, A. Noor, R. W. Snow, and C. O. Buckee.  Quantifying the impact of human mobility on malaria.  Science, 338(6104):267–270, 2012.
  * [59] A. Wesolowski, T. Qureshi, M. F. Boni, P. R. Sundsøy, M. A. Johansson, S. B. Rasheed, K. Engø-Monsen, and C. O. Buckee.  Impact of human mobility on the emergence of dengue epidemics in Pakistan.  Proceedings of the National Academy of Sciences, 112(38):11887–11892, 2015.
  * [60] M. E. J. Woolhouse, C. Dye, J. F. Etard, T. Smith, J. D. Charlwood, G. P. Garnett, P. Hagan, J. L. K. Hii, P. D. Ndhlovu, R. J. Quinnell, et al.  Heterogeneities in the transmission of infectious agents: implications for the design of control programs.  Proceedings of the National Academy of Sciences, 94(1):338–342, 1997.
  * [61] Q. Zhang, K. Sun, M. Chinazzi, A. P. Piontti, N. E. Dean, D. P. Rojas, S. Merler, D. Mistry, P. Poletti, L. Rossi, et al.  Spread of Zika virus in the Americas.  Proceedings of the National Academy of Sciences, 114(22):E4334–E4343, 2017.