# Analysis and Applications of Population Flows
in a Networked SEIRS Epidemic Process

Brooks A. Butler
Elmore Family School of Electrical and Computer Engineering
Purdue University
<EMAIL_ADDRESS>
& Raphael Stern
Department of Civil Engineering
University of Minnesota
<EMAIL_ADDRESS>
& Philip E. Paré
Elmore Family School of Electrical and Computer Engineering
Purdue University
<EMAIL_ADDRESS>
This work was funded in part by the National Science Foundation, grants NSF-
CNS #2028738 (P.E.P.), NSF-CNS #2028946 (R.S.), and NSF-ECCS #2032258 and
#2238388 (P.E.P.).

###### Abstract

Transportation networks play a critical part in the spread of infectious
diseases between populations. In this work, we define a networked susceptible-
exposed-infected-recovered epidemic process with loss of immunity over time
(SEIRS) that explicitly models the flow of individuals between sub-
populations, which serves as the propagating mechanism for infection. We
provide sufficient conditions for local stability and instability of the
healthy state of the system and show that no perturbation of population flows
can change the local stability of any healthy state. We also provide
sufficient conditions for the existence and uniqueness of an endemic state. We
then develop tools and methods for applying our model to real-world data,
including spreading parameter estimation and disease arrival time prediction,
and apply them in a case study using both travel and infection data from
counties in Minnesota during the first year of the COVID-19 pandemic.

_Keywords_ Epidemic Spread $\cdot$ Transportation $\cdot$ Network Modeling

##  1 Introduction

With the prevalence of transportation networks that enable rapid travel
between locations globally, understanding how these transportation networks
facilitate the spread of infectious diseases only continues to grow in
importance. While certainly not the first globally transmitted disease, the
recent COVID-19 pandemic has demonstrated the negative effect that a highly
infectious disease can have on general health and safety of any population
connected via travel Ruan et al. (2006); Tatem et al. (2006); Lowe et al.
(2014); Findlater and Bogoch (2018); Hollingsworth et al. (2006); Tang et al.
(2009), as well as long-lasting impacts on both local and global economies
Maital and Barzani (2020); Kaye et al. (2021); Ahmad et al. (2020); Kumar et
al. (2020). Given the growing impact and scale of human mobility as global
interconnectivity continues to increase, as well as the near certainty of
other pandemic-level diseases emerging in the future, it is critical that we
develop a rigorous and sophisticated understanding of such phenomena in the
networked sense. Thus, in this work, we aim to contribute to a deeper
understanding of how transportation networks may be used to model and predict
the transmission of infectious diseases over a network of connected
populations.

This work builds off the deep body of literature dedicated to modeling
epidemic processes using compartmental models Kermack and McKendrick (1927);
Brauer (2008); Prasse and Van Mieghem (2020); Fagnani and Zino (2017); Chen et
al. (2019). Such models separate a population, or individuals in a population,
into distinct compartments which represent their infection state. Simple
compartmental models may have as few as two states, such as being either
susceptible to a given disease or infected. For diseases with more complex
behavior, more compartments may be needed to better describe the progression
of a population through the epidemic process of the disease. In this work, we
focus on a four-compartment model which separates individuals into one of four
states, namely: susceptible, exposed, infected, and recovered. Further, we
allow for a rate at which a population or individual may lose immunity to the
given disease, and therefore may become susceptible again. This model,
abbreviated to the SEIRS model, allows for the modeling of a disease that has
a period in which an individual may be exposed and will become infectious, but
is not immediately apparent. Additionally, the recovered state and eventual
loss of immunity model a disease that can evolve over time, but may allow for
temporary immunity after an initial recovery. This model is chosen, in part,
to better emulate what is now known of the behavior of COVID-19 and its
variants Xiang et al. (2021); Rahmandad et al. (2021); Wang (2020); Chen et
al. (2020). Compartmental models alone, however, do not account for the
transmission of a disease in a networked environment where a general
population can be divided into subpopulations connected via various
transportation networks Nowzari et al. (2016); Qian and Ukkusuri (2021);
Salehi et al. (2015); Qu and Wang (2017); Kang et al. (2018); Nadini et al.
(2018); Li et al. (2021); Yang et al. (2012); Hackl and Dubernet (2019); Paré
et al. (2020). Thus, an increase in scope is needed to account for this
important behavioral complexity in disease transmission.

Previous work that incorporates networked population flows in an epidemic
process model includes analysis of both networked SIS and SIR models with
flows Ye et al. (2020); Yin et al. (2020) as well as using a networked SIR
model with flows to predict arrival times for various epidemics using global
flight data Brockmann and Helbing (2013), where both models are developed in
continuous time. Additionally, a data-driven approach to modeling the effect
of population flow on epidemic spreading is performed in Di Giamberardino et
al. (2021), which uses data from lockdowns in Italy during the initial
spreading of COVID-19. Additional work on epidemic spreading with population
flow has been performed in modeling multi-city epidemics Arino and van den
Driessche (2003); Liu and Stechlinski (2013); Knipl (2016); Pujari and
Shekatkar (2020) and in the modeling of disease spread in metapopulations
Arino and Van den Driessche (2006); Wang and Li (2014); Shao and Han (2022);
Zhu et al. (2021). We contribute to the development of such networked models
by including the exposed state in our model formulation as well as a loss of
immunity. In our previous work, we consider capturing the effect of
transportation on the spread of COVID-19 using a networked discrete-time SEIR
model Vrabac et al. (2022); Butler et al. (2021); however, the key distinction
in this work is that, in addition to being a continuous-time model, infection
propagation over the network is modeled by the relocation of infected
individuals to other sub-populations rather than assuming direct contact and
intermingling between sub-populations. Further, we allow recovered individuals
to experience a loss of immunity over time, which creates the possibility of
an endemic state for the epidemic process.

In this paper, we make the following contributions:

  1. 1.

We define a networked SEIRS model that incorporates travel flows as a means of
disease propagation

  2. 2.

We provide conditions for local stability and instability of the healthy state
and prove that perturbations to population flows cannot change the stability
of the healthy state

  3. 3.

We provide conditions under which a unique endemic state will exist

  4. 4.

We present methods for estimating the spreading parameters of our model using
both infection and travel data as well as a method for predicting the arrival
time of a disease using travel data to compute the effective distance between
nodes

  5. 5.

We demonstrate how these tools can be applied to our model using real-world
mobility data from Minnesota during the COVID-19 pandemic.

We now provide the necessary definitions of notation and previous analytical
results used in our model construction and analysis sections, respectively.

##  2 Preliminaries

We denote the set of real numbers and positive real numbers as $\mathbb{R}$
and $\mathbb{R}_{>0}$, respectively. For any positive integer $n$, we have
$[n]=\\{1,2,...,n\\}$. A diagonal matrix is denoted as diag$(\cdot)$. The
transpose of a vector $x\in\mathbb{R}^{n}$ is $x^{\top}$. We use $\mathbf{0}$
and $\mathbf{1}$ to denote the vectors of the appropriate size whose entries
all equal 0 and 1, respectively. We let
$\mathcal{G}=(\mathbb{V},\mathbb{E},\mathbb{W})$ denote a weighted directed
graph where $\mathbb{V}=\\{v_{1},v_{2},...,v_{n}\\}$ is the set of nodes,
$\mathbb{E}\subseteq\mathbb{V}\times\mathbb{V}$ is the set of edges, and
$\mathbb{W}:\mathbb{E}\rightarrow\mathbb{R}_{>0}$ maps to the real valued edge
weights on each edge. We denote the configuration of edges in a directed graph
at time $t$ as $\mathcal{G}(t)=(\mathbb{V},\mathbb{E}(t),\mathbb{W})$, where
$\mathbb{E}(t)$ denotes the set of edges at time $t$. Furthermore, we denote
$\cup_{t\geq 0}\mathbb{E}(t)$ as the union of all non-zero edge configurations
on a graph for all $t\geq 0$. We define a graph $\mathcal{G}$ as being
strongly connected if there is a path consisting of nonzero edge weights from
every node to every other node in the graph.

For a complex number $x$ we use $|x|$ and $\text{Re}(x)$ to denote its
magnitude and real part, respectively. For a real square matrix $M$, we use
$\rho(M)$ to denote its spectral radius and $s(M)$ to denote the largest real
part amongst its eigenvalues, i.e.,

$\displaystyle\rho(M)$
$\displaystyle=\max\\{|\lambda|:\lambda\in\sigma(M)\\},$ $\displaystyle s(M)$
$\displaystyle=\max\\{\text{Re}(\lambda):\lambda\in\sigma(M)\\},$

where $\sigma(M)$ denotes the spectrum of $M$.

A real square matrix is called Metzler if its off-diagonal entries are all
non-negative. Thus, any non-negative matrix is Metzler.

######  Lemma 1.

For any matrix $M$ and any real number $\phi$, if $A:=M-\phi I$, then
$\sigma(M)=\sigma(A)+\phi$.

The following results from Chapter 2 of Varga (1962) for non-negative
matrices, which also hold for Metzler matrices by Lemma 1, with
$\phi=\min\\{0,m_{11},\dots,m_{nn}\\}$, will be used in the subsequent
analysis.

######  Lemma 2.

(Lemma 2.3 in Varga (1962)) Suppose that $M$ is an irreducible Metzler matrix.
Then, $s(M)$ is a simple eigenvalue of $M$ and there exists a unique (up to
scalar multiple) vector $x\gg\mathbf{0}$ such that $Mx=s(M)x$.

######  Proposition 1.

Suppose that $\Lambda$ is a negative diagonal matrix in $\mathbb{R}^{n\times
n}$ and $N$ is an irreducible non-negative matrix in $\mathbb{R}^{n\times n}$.
Let $M=\Lambda+N$. Then, $s(M)<0$ if and only if $\rho(-\Lambda^{-1}N)<1$,
$s(M)=0$ if and only if $\rho(-\Lambda^{-1}N)=1$, and $s(M)>0$ if and only if
$\rho(-\Lambda^{-1}N)>1$.

######  Lemma 3.

(Levy–Desplanques Theorem) A strictly diagonally dominant matrix is non-
singular. In other words, let $A\in C^{n\times n}$ be a matrix satisfying the
property

$|a_{ii}|>\sum_{j\neq i}|a_{ij}|,\,\,\forall i;$

then $\det(A)\neq 0$.

######  Lemma 4.

(Observation 6.3.1 in Horn and Johnson (2012)) Let $M\in{\rm I\\!R}^{n\times
n}$ be diagonalizable with $M=S\Lambda S^{-1}$ and
$\Lambda=diag(\lambda_{1},\dots,\lambda_{n})$ and let $E\in{\rm
I\\!R}^{n\times n}$. If $\hat{\lambda}$ is an eigenvalue of $M+E$, then there
exists some eigenvalue $\lambda_{i}$ of $M$ for which

$\displaystyle|\hat{\lambda}-\lambda_{i}|\leq{||S||}_{\infty}{||S^{-1}||}_{\infty}{||E||}_{\infty}=\kappa_{\infty}(S){||E||}_{\infty},$

where $\kappa_{\infty}(S)$ denotes the condition number with respect to the
infinity matrix norm ${||\cdot||}_{\infty}$.

##  3 Network Flows Model Definition

We now present a networked SEIRS model incorporating the population flow of
individuals between subpopulations. First, consider a group of $n$ sub-
populations in a graph, where each sub-population $i\in[n]$ is represented by
a node in the graph $\mathcal{G}$. We use the SEIRS model to describe how
susceptible individuals in sub-population $i$ become exposed, infected,
recover, and gradually lose immunity as the result of an infectious disease
Shu et al. (2012). We begin with defining the SEIRS model behavior without
graph connections for each sub-population $i\in[n]$. Let $S_{i}$, $E_{i}$,
$I_{i}$, and $R_{i}$ represent the number of susceptible, exposed, infected,
and recovered individuals in sub-population $i$, respectively, and their
dynamics evolve as

$\displaystyle\dot{S}_{i}(t)$
$\displaystyle=\alpha_{i}R_{i}(t)-\beta_{i}\frac{I_{i}(t)}{N_{i}}S_{i}(t)$
(1a) $\displaystyle\dot{E}_{i}(t)$
$\displaystyle=\beta_{i}\frac{I_{i}(t)}{N_{i}}S_{i}(t)-\sigma_{i}E_{i}(t)$
(1b) $\displaystyle\dot{I}_{i}(t)$
$\displaystyle=\sigma_{i}E_{i}(t)-\delta_{i}I_{i}(t)$ (1c)
$\displaystyle\dot{R}_{i}(t)$
$\displaystyle=\delta_{i}I_{i}(t)-\alpha_{i}R_{i}(t),$ (1d)

where $\beta_{i}$ is the infection rate, $\sigma_{i}$ is the transition rate
from exposed to infected, $\delta_{i}$ is the healing rate, and $\alpha_{i}$
is the rate of immunity loss. Although the state variables, except population,
will continue to vary with time, we remove the time-dependence notation for
convenience and ease of reading from this point forward.

To account for the flow of individuals between sub-populations we expand the
model in (1):

$\displaystyle\dot{S}_{i}$
$\displaystyle=-\beta_{i}\frac{I_{i}}{N_{i}}S_{i}+{\color[rgb]{0,0,0}\alpha_{i}}R_{i}+\sum_{j\neq
i}\left(F_{ij}\frac{S_{j}}{N_{j}}-F_{ji}\frac{S_{i}}{N_{i}}\right)$ (2a)
$\displaystyle\dot{E}_{i}$
$\displaystyle=\beta_{i}\frac{I_{i}}{N_{i}}S_{i}-\sigma_{i}E_{i}+\sum_{j\neq
i}\left(F_{ij}\frac{E_{j}}{N_{j}}-F_{ji}\frac{E_{i}}{N_{i}}\right)$ (2b)
$\displaystyle\dot{I}_{i}$
$\displaystyle=\sigma_{i}E_{i}-\delta_{i}I_{i}+\sum_{j\neq
i}\left(F_{ij}\frac{I_{j}}{N_{j}}-F_{ji}\frac{I_{i}}{N_{i}}\right)$ (2c)
$\displaystyle\dot{R}_{i}$
$\displaystyle=\delta_{i}I_{i}-{\color[rgb]{0,0,0}\alpha_{i}}R_{i}+\sum_{j\neq
i}\left(F_{ij}\frac{R_{j}}{N_{j}}-F_{ji}\frac{R_{i}}{N_{i}}\right),$ (2d)

where $F_{ij}$ represents the number of individuals flowing from sub-
population $j$ to $i$, with $F_{ii}=0$. By making a substitution of variables
where
$s_{i}=S_{i}/N_{i},e_{i}=E_{i}/N_{i},x_{i}=I_{i}/N_{i},r_{i}=R_{i}/N_{i}$ we
can model the proportion of individuals as follows

$\displaystyle\dot{s}_{i}$
$\displaystyle=\alpha_{i}r_{i}-\beta_{i}x_{i}s_{i}+\frac{1}{N_{i}}\sum_{j\neq
i}\left(F_{ij}s_{j}-F_{ji}s_{i}\right)$ (3a) $\displaystyle\dot{e}_{i}$
$\displaystyle=\beta_{i}x_{i}s_{i}-\sigma_{i}e_{i}+\frac{1}{N_{i}}\sum_{j\neq
i}\left(F_{ij}e_{j}-F_{ji}e_{i}\right)$ (3b) $\displaystyle\dot{x}_{i}$
$\displaystyle=\sigma_{i}e_{i}-\delta_{i}x_{i}+\frac{1}{N_{i}}\sum_{j\neq
i}\left(F_{ij}x_{j}-F_{ji}x_{i}\right)$ (3c) $\displaystyle\dot{r}_{i}$
$\displaystyle=\delta_{i}x_{i}-\alpha_{i}r_{i}+\frac{1}{N_{i}}\sum_{j\neq
i}\left(F_{ij}r_{j}-F_{ji}r_{i}\right),$ (3d)

where $s_{i}+e_{i}+x_{i}+r_{i}=1$. Note that both (2) and (3) assume the
subpopulations are well mixed and that the likelihood of an individual
traveling is independent of their infectious state, that is, whether they are
susceptible, exposed, infected, or recovered. We can compute the number of
individuals flowing from sub-population $j$ to $i$ as
$F_{ij}=\gamma_{j}w_{ij}N_{j}$, where $\gamma_{j}$ is the proportion of the
total population flowing out of node $j$ computed as
$\gamma_{j}=\frac{\sum_{i\neq j}F_{ij}}{N_{j}}$, and $w_{ij}$ is the
proportion of traveling individuals flowing from sub-population $j$ to $i$
computed as $w_{ij}=\frac{F_{ij}}{\sum_{l\neq j}F_{lj}}$, with $w_{ii}=0$.
Thus, we can further derive the dynamics for the susceptible proportion at
sub-population $i$ as

$\displaystyle\dot{s}_{i}$
$\displaystyle=\alpha_{i}r_{i}-\beta_{i}x_{i}s_{i}+\frac{1}{N_{i}}\sum_{j\neq
i}\left(F_{ij}s_{j}-F_{ji}s_{i}\right)$
$\displaystyle=\alpha_{i}r_{i}-\beta_{i}x_{i}s_{i}+\frac{1}{N_{i}}\sum_{j\neq
i}\left(\gamma_{j}w_{ij}N_{j}s_{j}-\gamma_{i}w_{ji}N_{i}s_{i}\right)$
$\displaystyle=\alpha_{i}r_{i}-\beta_{i}x_{i}s_{i}+\sum_{j\neq
i}\left(\frac{N_{j}}{N_{i}}w_{ij}\gamma_{j}s_{j}-w_{ji}\gamma_{i}s_{i}\right).$

Using the fact that $\sum_{j\neq i}w_{ji}=1$, we have that

$\dot{s}_{i}=\alpha_{i}r_{i}-(\beta_{i}x_{i}+\gamma_{i})s_{i}+\sum_{j\neq
i}\frac{N_{j}}{N_{i}}w_{ij}\gamma_{j}s_{j}.$

By similar derivations we can rewrite (3) as

$\displaystyle\dot{s}_{i}$
$\displaystyle=\alpha_{i}r_{i}-(\beta_{i}x_{i}+\gamma_{i})s_{i}+\sum_{j\neq
i}\frac{N_{j}}{N_{i}}w_{ij}\gamma_{j}s_{j}$ (4a) $\displaystyle\dot{e}_{i}$
$\displaystyle=\beta_{i}x_{i}s_{i}-(\sigma_{i}+\gamma_{i})e_{i}+\sum_{j\neq
i}\frac{N_{j}}{N_{i}}w_{ij}\gamma_{j}e_{j}$ (4b) $\displaystyle\dot{x}_{i}$
$\displaystyle=\sigma_{i}e_{i}-(\delta_{i}+\gamma_{i})x_{i}+\sum_{j\neq
i}\frac{N_{j}}{N_{i}}w_{ij}\gamma_{j}x_{j}$ (4c) $\displaystyle\dot{r}_{i}$
$\displaystyle=\delta_{i}x_{i}-(\gamma_{i}+\alpha_{i})r_{i}+\sum_{j\neq
i}\frac{N_{j}}{N_{i}}w_{ij}\gamma_{j}r_{j}$ (4d)

\begin{overpic}[width=260.17464pt]{Images/SEIRS_flows_1_2.png}
\put(40.0,43.5){$\gamma_{1}$} \put(18.0,35.0){$\beta_{1}$}
\put(29.0,32.0){$\sigma_{1}$} \put(40.0,29.0){$\delta_{1}$}
\put(22.0,5.0){$\alpha_{1}$} \put(65.0,43.0){$\beta_{2}$}
\put(74.0,41.0){$\sigma_{2}$} \put(85.0,38.0){$\delta_{2}$}
\put(70.5,15.0){$\alpha_{2}$} \put(63.0,0.0){$\gamma_{2}$} \end{overpic}
Figure 1: An illustration of a two-node SIERS network with travel flows as
defined by (4).

where an illustration of this model for two nodes is shown in Figure 1. Note
that our model does not allow for state transitions to occur during travel,
i.e., we do not model individuals becoming infected as a result of their
travel. Rather, individuals will change their state in the epidemic process at
the rate defined at each node. Of course, this may not always be the case in
reality, as some methods of public transportation may facilitate infection
spreading for those sharing a confined space with inadequate ventilation.
Therefore, this study implicitly focuses on more long-range modes of
transportation that transfer individuals between more distinctly separated
populations (such as travel via flights or long-rage commuting) where
traveling individuals may not be interacting with other travelers, or where
there is a sufficient amount of ventilation such that the spread of the
disease is significantly reduced during transit Hu et al. (2021). The
inclusion of infections occurring on the edges of transportation networks is
left as a direction for future work.

For the model in (4) to be well-defined, we required the following
assumptions.

######  Assumption 1.

Let $\sum_{i\neq j}F_{ji}=\sum_{i\neq j}F_{ij}$ for all $i,j\in[n]$.

This assumption requires that the total flow of individuals into a given sub-
population must be equal to the total flow out, which naturally follows from
the assumption that the population at each node remains constant. Furthermore,
we impose the following assumption on the model parameters and initial
conditions.

######  Assumption 2.

Let $\alpha_{i},\beta_{i},\sigma_{i},\delta_{i}\in\mathbb{R}_{>0}$,
$s_{i}(t_{0}),e_{i}(t_{0}),x_{i}(t_{0}),r_{i}(t_{0})\in[0,1]$, and
$s_{i}(t_{0})+e_{i}(t_{0})+x_{i}(t_{0})+r_{i}(t_{0})=1$ for all $i\in[n]$ and
$t_{0}\in\mathbb{R}_{\geq 0}$.

Under these assumptions, we can show that the model will always remain well-
defined.

######  Lemma 5.

Let Assumptions 1-2 hold, then $s_{i}(t),e_{i}(t),x_{i}(t),r_{i}(t)\in[0,1]$,
and $s_{i}(t)+e_{i}(t)+x_{i}(t)+r_{i}(t)=1$ for all $i\in[n]$ and for all
$t\in\mathbb{R}_{\geq 0}$.

###### Proof.

First, by taking the sum of the system dynamics as defined in (3) we have

$\displaystyle\dot{s}_{i}(t)+\dot{e}_{i}(t)+\dot{x}_{i}(t)+\dot{r}_{i}(t)=0,$
(5)

by Assumption 1 at every time $t$. Thus, by integrating (5) from $t_{0}$ to
$t$ we have

$\displaystyle\int_{t_{0}}^{t}\dot{s}_{i}(\tau)+\dot{e}_{i}(\tau)+\dot{x}_{i}(\tau)+\dot{r}_{i}(\tau)d\tau=0,$
$\displaystyle\implies
s_{i}(t)+e_{i}(t)+x_{i}(t)+r_{i}(t)=1,t\in\mathbb{R}_{\geq 0}$

by Assumption 2. To show that $s_{i}(t),e_{i}(t),x_{i}(t),r_{i}(t)\in[0,1]$,
we consider the system as each state approaches zero. Note that
$\lim_{q_{i}(t)\rightarrow 0}\dot{q}_{i}(t)\geq 0$ where
$q_{i}\in\\{s_{i},e_{i},x_{i},r_{i}\\}$. Thus, by the assumption that the
initial states are non-negative, we have that $q_{i}(t)\geq 0,\forall t\geq
0$. Further, since we have shown that $s_{i}(t)+e_{i}(t)+x_{i}(t)+r_{i}(t)=1$
it follows directly that $s_{i}(t),e_{i}(t),x_{i}(t),r_{i}(t)\in[0,1],\forall
t\geq 0$ and $\forall i\in[n]$. ∎

For the purpose of the forthcoming analysis, we also express the system in (4)
in matrix form. Let $\Phi=N^{-1}W\Gamma N$, where $N=diag(N_{i})$,
$W\in\mathbb{R}^{n\times n}$ is the matrix comprised of the entries $w_{ij}$,
and $\Gamma=diag(\gamma_{i})$. The vectorized equations then become

$\displaystyle\dot{s}$ $\displaystyle=Ar-(BX(t)+\Gamma)s+\Phi s$ (6a)
$\displaystyle\dot{e}$ $\displaystyle=BX(t)s-(\Sigma+\Gamma)e+\Phi e$ (6b)
$\displaystyle\dot{x}$ $\displaystyle=\Sigma e-(D+\Gamma)x+\Phi x$ (6c)
$\displaystyle\dot{r}$ $\displaystyle=Dx-(A+\Gamma)r+\Phi r,$ (6d)

where $A=diag(\alpha_{i})$, $B=diag(\beta_{i})$, $\Sigma=diag(\sigma_{i})$,
$D=diag(\delta_{i})$, $X(t)=diag(x_{i})$, $s=[s_{1},\dots,s_{n}]^{\top}$,
$e=[e_{1},\dots,e_{n}]^{\top}$, $x=[x_{1},\dots,x_{n}]^{\top}$, and
$r=[r_{1},\dots,r_{n}]^{\top}$. We further construct the matrix

$\underbrace{\begin{bmatrix}-BX(t)-\Gamma+\Phi&0&0&A\\\
BX(t)&-\Sigma-\Gamma+\Phi&0&0\\\ 0&\Sigma&-D-\Gamma+\Phi&0\\\
0&0&D&-A-\Gamma+\Phi\end{bmatrix}}_{H(t)},$ (7)

where

$H(t)=-Q(t)+M(t),$ (8)

with

$Q(t)=\begin{bmatrix}{\color[rgb]{0,0,0}BX(t)}+\Gamma&&&\\\
&\Sigma+\Gamma&&\\\ &&D+\Gamma&\\\ &&&A+\Gamma\end{bmatrix},$ (9)
$\displaystyle M(t)$ $\displaystyle=\begin{bmatrix}\Phi&0&0&A\\\
BX(t)&\Phi&0&0\\\ 0&\Sigma&\Phi&0\\\ 0&0&D&\Phi\end{bmatrix},$ (10)

and let the complete state vector be defined as

$z=\begin{bmatrix}s^{\top}&e^{\top}&x^{\top}&r^{\top}\end{bmatrix}^{\top}.$
(11)

Thus, we can describe the state dynamics as

$\dot{z}=\big{(}-Q(t)+M(t)\big{)}z,$ (12)

where $Q(t)$ is a positive diagonal matrix and $M(t)$ is a non-negative matrix
$\forall t\geq 0$. Note that the time-dependence of $Q(t)$ and $M(t)$ comes
from the fact that they are functions of $X(t)$, where the indication of time
dependence is left for emphasis of this fact.

######  Definition 1.

A graph $\mathcal{G}(t)=(\mathbb{V},\cup_{t\geq 0}\mathbb{E}(t),\mathbb{W})$
for $t\in\mathbb{R}_{\geq 0}$ is _$K$ -strongly connected_ if there exist some
bound $K$ such that $(\mathbb{V},\cup_{j=t}^{t+K-1}\mathbb{E}(j),\mathbb{W})$
is strongly connected, for all $t\in\mathbb{R}_{\geq 0}$.

######  Assumption 3.

Let the graph $\mathcal{G}(t)=(\mathbb{V},\cup_{t\geq
0}\mathbb{E}(t),\mathbb{W})$, where
$\mathbb{W}:\mathbb{E}(t)\rightarrow\mathbb{R}_{>0}$ is defined by
$w_{ij}(t)$, be $K$-strongly connected.

######  Lemma 6.

Given Assumptions 2-3, if $\exists t_{0}\geq 0$ and $\exists i\in[n]$ such
that $x_{i}(t_{0})>0$, the matrix $M(t_{0})$ is irreducible.

###### Proof.

Assume, by contradiction, that $M(t_{0})$ is reducible. Then
$M(t_{0})-Q(t_{0})$ will also be reducible, as $Q(t_{0})$ is a diagonal
matrix. Thus, the graph described by the adjacency matrix in (12) cannot be
strongly connected. However, we can partition the graph into $n$ strongly
connected components where each component describes the evolution of the
epidemic in each subpopulation. Consider the case where $n=1$, which reduces
to a non-networked SEIRS model. By Assumption 2, this system is strongly
connected in the sense of the epidemic states given $x_{1}(t_{0})>0$ where

$\displaystyle\begin{bmatrix}\dot{s}_{1}\\\ \dot{e}_{1}\\\ \dot{x}_{1}\\\
\dot{r}_{1}\end{bmatrix}=\begin{bmatrix}-\beta_{1}x_{1}(t_{0})&0&0&\alpha_{1}\\\
\beta_{1}x_{1}(t_{0})&-\sigma_{1}&0&0\\\ 0&\sigma_{1}&-\delta_{1}&0\\\
0&0&\delta_{1}&-\alpha_{1}\end{bmatrix}\begin{bmatrix}s_{1}\\\ e_{1}\\\
x_{1}\\\ r_{1}\end{bmatrix}.$

For $n>1$, the network connections are defined by the total flows between each
sub-population. Thus, for the system to be reducible, there must exist a sub-
population that is unreachable via travel from at least one other sub-
population in the graph, which contradicts Assumption 3. Therefore, $M(t_{0})$
must be irreducible. ∎

Note that the irreducibulity of $M(t)$ will play an important role in proving
the existence of an endemic equilibrium, which will be proven in Section 4.

##  4 Network Flows Model Analysis

In this section, we analyze the fundamental behaviors of our model, which is
divided into the analysis of the healthy and the endemic states of the system
in (4).

###  4.1 Healthy State Analysis

We first consider the existence of a disease-free, or healthy state of the
system and when the stability of the healthy state can be ensured. We then
examine the effect of manipulating travel flows on the stability of the
healthy state and prove that there is no permissible perturbation of flows
that will change the local stability of the healthy state for any system.

######  Proposition 2.

The healthy state,
$(s^{*},e^{*},x^{*},r^{*})=(\mathbf{1},\mathbf{0},\mathbf{0},\mathbf{0})$, is
always an equilibrium of the system.

###### Proof.

Substituting $(\mathbf{1},\mathbf{0},\mathbf{0},\mathbf{0})$ into (3) yields

$\displaystyle\dot{s}_{i}=\frac{1}{N_{i}}\sum_{j\neq
i}\left(F_{ij}-F_{ji}\right),\,\,\dot{e}_{i}=0,\,\,\dot{x}_{i}=0,\,\,\dot{r}_{i}=0$

for all $i\in[n]$, where $\sum_{j\neq i}\left(F_{ij}-F_{ji}\right)=0$ by
Assumption 3. Thus,
$(s^{*},e^{*},x^{*},r^{*})=(\mathbf{1},\mathbf{0},\mathbf{0},\mathbf{0})$ is
always an equilibrium of the system given Assumptions 1 and 2. ∎

We now consider the stability of the healthy state of the system and the
conditions needed to reach this equilibrium.

First, since $s_{i}+e_{i}+x_{i}+r_{i}=1$ for all $i\in[n]$, we can rewrite the
dynamics of (6) with respect to the exposed, infected, and recovered states as

$\displaystyle\begin{bmatrix}\dot{e}\\\ \dot{x}\\\
\dot{r}\end{bmatrix}=\underbrace{\begin{bmatrix}-\Sigma-\Gamma+\Phi&B(I-E-
X-R)&0\\\ \Sigma&-D-\Gamma+\Phi&0\\\
0&D&-A-\Gamma+\Phi\end{bmatrix}}_{\bar{H}}\begin{bmatrix}e\\\ x\\\
r\end{bmatrix},$ (13)

where $E=diag(e_{i})$, $X=diag(x_{i})$, and $R=diag(r_{i})$. The Jacobian
matrix of (13) evaluated at an equilibrium $(x^{*},e^{*},r^{*})$, denoted as
$J_{(e^{*},x^{*},r^{*})}$, is

$\displaystyle\begin{bmatrix}-\Sigma-\Gamma+\Phi-
BX^{*}&B(I-E^{*}-2X^{*}-R^{*})&-BX^{*}\\\ \Sigma&-D-\Gamma+\Phi&0\\\
0&D&-A-\Gamma+\Phi\end{bmatrix}.$ (14)

Further, we define the following matrix, the top 2x2 diagonal block of $J$
evaluated at $(\mathbf{0},\mathbf{0},\mathbf{0})$,

$U=\begin{bmatrix}-\Sigma-\Gamma+\Phi&B\\\
\Sigma&-D-\Gamma+\Phi\end{bmatrix},$ (15)

which yields the following conditions for the local stability and instability
of the healthy state.

######  Proposition 3.

If $s(U)<0$, then the healthy state of (6) is locally exponentially stable.

###### Proof.

Evaluating (14) at the healthy state equilibrium yields

$J_{(\mathbf{1},\mathbf{0},\mathbf{0},\mathbf{0})}=\begin{bmatrix}-\Sigma-\Gamma+\Phi&B&0\\\
\Sigma&-D-\Gamma+\Phi&0\\\ 0&D&-A-\Gamma+\Phi\end{bmatrix}.$ (16)

Since (16) is block lower triangular, we can examine the eigenvalues of each
block separately. By Assumption 1 we have that $\gamma_{i}=\sum_{j\neq
i}\phi_{ij}$, thus by the Gershgorin Circle Theorem we have that
$s(-A-\Gamma-\Phi)<0$. The upper block is equal to $U$, which is negative
definite by assumption. Thus, we have that
$s(J_{(\mathbf{1},\mathbf{0},\mathbf{0},\mathbf{0})})<0$, which satisfies
Lyapunov’s indirect method for determining the local stability of the healthy
state. ∎

######  Corollary 1.

If $s(U)>0$, then the healthy state of (6) is unstable.

###### Proof.

This result follows by the same method as the proof of Proposition 3, which
leverages (14) evaluated at the healthy state, shown in (16), which is block
lower triangular. Since the upper block is equal to $U$, if $s(U)>0$ then
$s(J_{(\mathbf{1},\mathbf{0},\mathbf{0},\mathbf{0})})>0$. Thus, by Lyapunov’s
indirect method, the healthy state is unstable. ∎

We now discuss the local stability of the healthy state for (13) with respect
to changes in the rate of flow between nodes.

######  Theorem 1.

If $J_{(\mathbf{1},\mathbf{0},\mathbf{0},\mathbf{0})}$, as defined in (16),
has distinct eigenvalues then, given Assumptions 1 and 2, there exists no
permissible perturbation to the population flows $\gamma_{i}$, for $i\in[n]$,
that will change the local stability of the healthy state of (6).

###### Proof.

Let $\gamma_{i}^{\prime}$ be some perturbation to $\gamma_{i},i\in[n]$ such
that Assumption 1 is satisfied and
$\gamma^{\prime}_{i}=\gamma_{i}+\theta_{i}$. We can then define the perturbed
system dynamics of (13) with the matrix

$\bar{H}^{\prime}=\begin{bmatrix}-\Sigma-\Gamma^{\prime}+\Phi^{\prime}&B(I-E-
X-R)&0\\\ \Sigma&-D-\Gamma^{\prime}+\Phi^{\prime}&0\\\
0&D&-A-\Gamma^{\prime}+\Phi^{\prime}\end{bmatrix},$ (17)

where $\Gamma^{\prime}=\Gamma+\Theta$, $\Theta=diag(\theta_{i})$, and
$\Phi^{\prime}=N^{-1}W(\Gamma+\Theta)N$. We can then separate the perturbed
system into

$\bar{H}^{\prime}=\bar{H}+\bar{E},$ (18)

where $\bar{H}$ is the original system dynamics defined in (13) and

$\bar{E}=\begin{bmatrix}N^{-1}W\Theta N-\Theta&0&0\\\ 0&N^{-1}W\Theta
N-\Theta&0\\\ 0&0&N^{-1}W\Theta N-\Theta\end{bmatrix}.$ (19)

Taking the Jacobian of (18) yields

$J^{\prime}_{(e^{*},x^{*},r^{*})}=J_{(e^{*},x^{*},r^{*})}+\bar{E}.$

We consider the healthy state equilibrium, which exists for all allowed
parameters of our system by Proposition 2. Since the eigenvalues of
$J_{{(\mathbf{1},\mathbf{0},\mathbf{0},\mathbf{0})}}$ are distinct, we have
that $J_{(\mathbf{1},\mathbf{0},\mathbf{0},\mathbf{0})}$ is diagonalizable by
some matrix $\bar{S}$, where
$J_{(\mathbf{1},\mathbf{0},\mathbf{0},\mathbf{0})}=\bar{S}\Lambda\bar{S}^{-1}$
and $\Lambda=diag(\lambda_{i})$. Thus, by Lemma 4 we have that if
$\lambda^{\prime}$ is an eigenvalue of $J^{\prime}$, then there exists some
eigenvalue $\lambda_{i}$ of $J$ such that

$|\lambda^{\prime}-\lambda_{i}|\leq\kappa(\bar{S})_{\infty}{||\bar{E}||}_{\infty}.$
(20)

However, by Assumption 1 we have that $\theta_{i}=\sum_{j\neq
i}\frac{N_{j}}{N_{i}}w_{ij}\theta_{j}$. Thus, ${||\bar{E}||}_{\infty}=0$ and
for any eigenvalue $\lambda^{\prime}$ of $J^{\prime}$ there exists some
eigenvalue $\lambda_{i},i\in[n]$ such that $\lambda^{\prime}=\lambda_{i}$.
Therefore, no perturbation of the travel flows that maintains Assumption 1,
and thus keeps the model well-defined, can change the local stability of the
healthy state of (6). ∎

###  4.2 Endemic State Analysis

We now present conditions for when an endemic state of the system, where
$z^{*}\gg\mathbf{0}$, will exist. For clarity, we define the following
matrices explicitly as a function of the full system state

$\displaystyle Q^{*}(z^{*})$
$\displaystyle=\begin{bmatrix}BX^{*}(z^{*})+\Gamma&0&0&0\\\
0&\Sigma+\Gamma&0&0\\\ 0&0&D+\Gamma&0\\\ 0&0&0&A+\Gamma\end{bmatrix}$ (21)
$\displaystyle M^{*}(z^{*})$ $\displaystyle=\begin{bmatrix}\Phi&0&0&A\\\
BX^{*}(z^{*})&\Phi&0&0\\\ 0&\Sigma&\Phi&0\\\ 0&0&D&\Phi\end{bmatrix},$ (22)

where
$X^{*}(z^{*})=diag\left(\begin{bmatrix}0&0&I_{n}&0\end{bmatrix}z^{*}\right).$

######  Theorem 2.

If $\inf_{t\geq t_{0}}s(-Q(t)+M(t))>0$, then there exists an endemic
equilibrium such that $z^{*}\gg{\color[rgb]{0,0,0}\mathbf{0}}$.

###### Proof.

Define the continuous map $f:(0,1]^{4n}\rightarrow[0,b]^{4n}$, where
$b\in\mathbb{R}_{\geq 1}$, given by

$f(z)=Q^{*}(z)^{-1}M^{*}(z)z.$ (23)

Since the domain of $f$ is $(0,1]^{4n}$, $z$ as an argument of $f$ satisfies
$z\gg 0$. Further, by Lemma 6, since $Q^{*}(z)^{-1}M^{*}(z)$ is an
irreducible, non-negative matrix for all $z\gg{\color[rgb]{0,0,0}\mathbf{0}}$,
by Lemma 2, there exists a $v\gg{\color[rgb]{0,0,0}\mathbf{0}}$ such that

$Q^{*}(z)^{-1}M^{*}(z)v=cv,$ (24)

where $c=\rho\big{(}Q^{*}(z)^{-1}M^{*}(z)\big{)}$. Since $\inf_{t\geq
t_{0}}s(-Q(t)+M(t))>0$ we have, by Proposition 1, that $c>1$. We can then find
an $\varepsilon>0$ such that for each $i\in[4n]$

$\displaystyle\varepsilon v_{i}\leq\frac{c-1}{c}$ (25)

and

$\varepsilon v_{i}\leq z_{i}.$ (26)

From (25), it follows that $1\leq\frac{c}{1+\varepsilon cv_{i}}$ and thus
$\varepsilon v_{i}\leq\frac{\varepsilon cv_{i}}{1+\varepsilon cv_{i}}$,
yielding

$\displaystyle\varepsilon v_{i}$
$\displaystyle\leq\frac{\big{(}Q^{*}(z)^{-1}M^{*}(z)\varepsilon
v\big{)}_{i}}{1+\big{(}Q^{*}(z)^{-1}M^{*}(z)\varepsilon v\big{)}_{i}}.$ (27)

We can show that the right hand side of (27) is less than or equal to
$f(\varepsilon v)_{i}=\big{(}Q^{*}(\varepsilon v)^{-1}M^{*}(\varepsilon
v)\varepsilon v\big{)}_{i}$ for each $i\in[4n]$ as follows. We can expand (23)
using (11) as

$\displaystyle f(z)=$ $\displaystyle(diag(Bx^{z})+\Gamma)^{-1}(Ar^{z}+\Phi
s^{z})$ (28a) $\displaystyle f(z)=$
$\displaystyle(\Sigma+\Gamma)^{-1}(diag(Bx^{z})s^{z}+\Phi e^{z})$ (28b)
$\displaystyle f(z)=$ $\displaystyle(D+\Gamma)^{-1}(\Sigma e^{z}+\Phi x^{z})$
(28c) $\displaystyle f(z)=$ $\displaystyle(A+\Gamma)^{-1}(Dx^{z}+\Phi r^{z})$
(28d)

where we can write
$z=\begin{bmatrix}{s^{z}}^{\top}&{e^{z}}^{\top}&{x^{z}}^{\top}&{r^{z}}^{\top}\end{bmatrix}^{\top}$.
Note that for (28b)-(28d), if we increase any value $z_{i}$ the output also
increases, i.e., (28b)-(28d) are monotonic. Thus, since $\varepsilon v_{i}\leq
z_{i}$, for any $f(\varepsilon v)_{i}$ in (28b)-(28d), we have $\varepsilon
v_{i}\leq f(\varepsilon v)_{i}$.

We now consider (28a), which is not monotonic, and show that $\varepsilon
v_{i}\leq f(\varepsilon v)_{i}$, for the last case (28a). To show this is
true, we write the vector $v$ as

$v=\begin{bmatrix}{v^{s}}^{\top}&{v^{e}}^{\top}&{v^{x}}^{\top}&{v^{r}}^{\top}\end{bmatrix}^{\top}$

and note that by (26) we have

$\displaystyle\beta_{i}\varepsilon v_{i}^{x}+\gamma_{i}$
$\displaystyle\leq\beta_{i}x_{i}+\gamma_{i},$ (29)

for every $i\in[n]$. To simplify notation, let

$a_{i}=(\alpha_{i}v^{r}_{i}+(\Phi v^{s})_{i}).$ (30)

We then have that

$\displaystyle\beta_{i}\varepsilon v_{i}^{x}+\gamma_{i}$
$\displaystyle\leq(\beta_{i}x_{i}+\gamma_{i})(1+(\beta_{i}x_{i}+\gamma_{i})^{-1}\varepsilon
a_{i})$ (31)

since $(\beta_{i}x_{i}+\gamma_{i})^{-1}\varepsilon a_{i}\geq 0$, which then
yields

$\displaystyle\frac{(\beta_{i}x_{i}+\gamma_{i})^{-1}}{1+(\beta_{i}x_{i}+\gamma_{i})^{-1}\varepsilon
a_{i}}$ $\displaystyle\leq(\beta_{i}\varepsilon v^{x}_{i}+\gamma_{i})^{-1}.$
(32)

Finally, multiplying both sides of (32) by $\varepsilon a_{i}$ we have

$\displaystyle\frac{(\beta_{i}x_{i}+\gamma_{i})^{-1}\varepsilon
a_{i}}{1+(\beta_{i}x_{i}+\gamma_{i})^{-1}\varepsilon a_{i}}$
$\displaystyle\leq(\beta_{i}\varepsilon v^{x}_{i}+\gamma_{i})^{-1}\varepsilon
a_{i}.$ (33)

Therefore, since by (27) we have

$\displaystyle\varepsilon
v_{i}^{s}\leq\frac{(\beta_{i}x_{i}+\gamma_{i})^{-1}\varepsilon(\alpha_{i}v^{r}_{i}+(\Phi
v^{s})_{i})}{1+(\beta_{i}x_{i}+\gamma_{i})^{-1}\varepsilon(\alpha_{i}v^{r}_{i}+(\Phi
v^{s})_{i})}$

and by (28a)

$\displaystyle f(\varepsilon v^{s})_{i}=(\beta_{i}\varepsilon
v^{x}_{i}+\gamma_{i})^{-1}\varepsilon(\alpha_{i}v^{r}_{i}+(\Phi v^{s})_{i}),$

we have that $\varepsilon v_{i}^{s}\leq f(\varepsilon v^{s})_{i}$, for all
$i\in[n]$, by (33). Thus, we have shown that $\varepsilon v_{i}\leq
f(\varepsilon v)_{i}$, for all $i\in[4n]$.

We now show that $f(z)$ is bounded by a $b\in{\rm I\\!R}$. Since
${\color[rgb]{0,0,0}\mathbf{0}}\ll z\leq{\color[rgb]{0,0,0}\mathbf{1}}$ we
have by Assumption 2 and (28d) that $f(z)$ must also be finite and we can
choose $b=\sup_{i\in[4n]}f(z)_{i}<\infty$. Thus, $f$ maps the convex compact
set $\mathcal{C}=\\{z\,\,|\,\,\varepsilon v\leq z\leq\mathbf{1}b\\}$ to
itself. By Brouwer’s fixed-point theorem, $f$ has a fixed point in
$\mathcal{C}$, which must be strictly positive. Let $z^{*}$ be this fixed
point, then $f(z^{*})=z^{*}$, i.e.

$z^{*}=Q^{*}(z^{*})^{-1}M^{*}(z^{*})z^{*}.$ (34)

Therefore, we have

$\displaystyle Q^{*}(z^{*})z^{*}$ $\displaystyle=M^{*}(z^{*})z^{*}$
$\displaystyle{\color[rgb]{0,0,0}\mathbf{0}}$
$\displaystyle=(M^{*}(z^{*})-Q^{*}(z^{*}))z^{*}$

and thus $z^{*}\gg{\color[rgb]{0,0,0}\mathbf{0}}$ is an equilibrium point of
(4). Further, by Lemma 5 we know that $z^{*}\in(0,1]^{4n}$. ∎

Note that since our method of proving existence of an endemic equilibrium for
(4) requires irreducibility we must include the susceptible state in our
analysis, causing the sufficient condition needed to prove existence to be
stronger than one that is time invariant. We now give conditions under which
the endemic equilibrium is unique.

######  Theorem 3.

Let there exist an endemic equilibrium for the system in (4). If
$\beta_{i}\geq\gamma_{i}$, for all $i\in[n]$, then the endemic equilibrium is
unique.

###### Proof.

Let $z$ and $y$ be two nonzero equilibria of (4) which can be divided into
their respective epidemic states as
$z=\begin{bmatrix}{s^{z}}^{\top}&{e^{z}}^{\top}&{x^{z}}^{\top}&{r^{z}}^{\top}\end{bmatrix}^{\top}$
and
$y=\begin{bmatrix}{s^{y}}^{\top}&{e^{y}}^{\top}&{x^{y}}^{\top}&{r^{y}}^{\top}\end{bmatrix}^{\top}$.
We can express $z$ in terms of $y$ for $i\in[n]$ as

$\displaystyle(s_{i}^{y}+\varepsilon_{s_{i}})+(e_{i}^{y}+\varepsilon_{e_{i}})+(x_{i}^{y}+\varepsilon_{x_{i}})+(r_{i}^{y}+\varepsilon_{r_{i}})$
(35) $\displaystyle=s_{i}^{z}+e_{i}^{z}+x_{i}^{z}+r_{i}^{z}=1,$

where

$\displaystyle\varepsilon_{s_{i}}+\varepsilon_{e_{i}}+\varepsilon_{x_{i}}+\varepsilon_{r_{i}}=0,$
(36)

for each $i\in[n]$. Since both $z\gg 0$ and $y\gg 0$ are fixed points of (4),
we can sum (4b), (4c), and (4d) to obtain

$\displaystyle\dot{e}_{i}^{z}+\dot{x}_{i}^{z}+\dot{r}_{i}^{z}$
$\displaystyle=\beta_{i}x_{i}^{z}{\color[rgb]{0,0,0}s_{i}^{z}}-\gamma_{i}(e_{i}^{z}+x_{i}^{z}+r_{i}^{z})-\alpha_{i}r_{i}^{z}+\sum_{j\neq
i}\phi_{ij}(e_{j}^{z}+x_{j}^{z}+r_{j}^{z})$
$\displaystyle=\beta_{i}(x_{i}^{y}+\varepsilon_{x_{i}})(s_{i}^{y}+\varepsilon_{s_{i}})-\alpha_{i}(r_{i}^{y}+\varepsilon_{r_{i}})-\gamma_{i}(e_{i}^{y}+x_{i}^{y}+r_{i}^{y}+\varepsilon_{e_{i}}+\varepsilon_{x_{i}}+\varepsilon_{r_{i}})$
$\displaystyle\,\,\,\,\,\,+\sum_{j\neq
i}\phi_{ij}(e_{j}^{y}+x_{j}^{y}+r_{j}^{y}+\varepsilon_{e_{j}}+\varepsilon_{x_{j}}+\varepsilon_{r_{j}})$
$\displaystyle=\beta_{i}x_{i}^{y}s_{i}^{y}-\gamma_{i}(e_{i}^{y}+x_{i}^{y}+r_{i}^{y})-\alpha_{i}r_{i}^{y}+\sum_{j\neq
i}\phi_{ij}(e_{j}^{y}+x_{j}^{y}+r_{j}^{y})$
$\displaystyle=\dot{e}_{i}^{y}+\dot{x}_{i}^{y}+\dot{r}_{i}^{y}=0,$

where $\phi_{ij}=\frac{N_{j}}{N_{i}}w_{ij}\gamma_{j}$. Thus, by taking the
difference of
$(\dot{e}_{i}^{z}+\dot{x}_{i}^{z}+\dot{r}_{i}^{z})-(\dot{e}_{i}^{y}+\dot{x}_{i}^{y}+\dot{r}_{i}^{y})$
we get

$\displaystyle\beta_{i}\big{(}x_{i}^{y}\varepsilon_{s_{i}}+s_{i}^{y}\varepsilon_{x_{i}}+\varepsilon_{x_{i}}\varepsilon_{s_{i}}\big{)}-\gamma_{i}(\varepsilon_{e_{i}}+\varepsilon_{x_{i}}+\varepsilon_{r_{i}})-\alpha_{i}\varepsilon_{r_{i}}+\sum_{j\neq
i}\phi_{ij}(\varepsilon_{e_{j}}+\varepsilon_{x_{j}}+\varepsilon_{r_{j}})$
$\displaystyle=0.$ (37)

We can express (37) for all $i\in[n]$ using matrix notation as follows. Define
$\varepsilon_{q}=[\varepsilon_{q_{1}},\dots,\varepsilon_{q_{n}}]^{\top}$,
where $q\in\\{s,e,x,r\\}$. Then, using (36) and (37) we have

$B(X^{y}\varepsilon_{s}+S^{y}\varepsilon_{x}+E_{x}\varepsilon_{s})+\Gamma\varepsilon_{s}-A\varepsilon_{r}-\Phi\varepsilon_{s}=\mathbf{0},$
(38)

where $X^{y}=diag(x_{i}^{y})$, $S^{y}=diag(s_{i}^{y})$, and
$E_{x}=diag(\varepsilon_{x_{i}})$. Further, we can put (38) into a block
matrix form as

$\displaystyle\underbrace{\begin{bmatrix}BX^{y}+E_{x}+\Gamma-\Phi&0&0\\\
0&BS^{y}&0\\\ 0&0&-A\end{bmatrix}}_{K}\begin{bmatrix}\varepsilon_{s}\\\
\varepsilon_{x}\\\ \varepsilon_{r}\end{bmatrix}=\mathbf{0}.$ (39)

Thus, the solution to the vector of perturbations to $s^{y},x^{y},r^{y}$ is
the nullspace of $K$. Note that if $K$ is full rank, then the only solution
for the perturbations must be
$\varepsilon_{s}=\varepsilon_{x}=\varepsilon_{r}=\mathbf{0}$. Considering the
rank of $K$, it is clear that the second and third block diagonals are full
rank since $y\gg 0$ and by Assumption 2, respectively. Further, by Lemma 3, we
have that if $BX^{y}+E_{x}+\Gamma-\Phi$ is diagonally dominant then it is also
full rank, or in other words, if

$|\beta_{i}x_{i}^{y}+\varepsilon_{x_{i}}+\gamma_{i}|>\sum_{j\neq
i}|\phi_{ij}|,\,\,\,\forall i\in[n],$ (40)

then $\det(BX^{y}+E_{x}+\Gamma-\Phi)\neq 0$.

We will now show that (40) holds for any choice of $\varepsilon_{x_{i}}$ if
$\beta_{i}\geq 1$ for all $i\in[n]$. Note that when $\varepsilon_{x_{i}}\geq
0$ we have that (40) is satisfied, since $\gamma_{i}=\sum_{j\neq i}\phi_{ij}$
for all $i\in[n]$ by Assumption 1, and $\beta_{i}x^{y}_{i}>0$ by Assumption 2.
We now consider when $\varepsilon_{x_{i}}<0$. Since
$x_{i}^{y}+\varepsilon_{x_{i}}=x_{i}^{z}>0$, we have that
$x_{i}^{y}>-\varepsilon_{x_{i}}$. Therefore, if $\varepsilon_{x_{i}}<0$, then
we have that $x_{i}^{y}>|\varepsilon_{x_{i}}|$. Thus, if $\beta_{i}\geq 1$ for
all $i\in[n]$, then (40) must hold for any feasible choice of
$\varepsilon_{x_{i}}$. Therefore, by (39) we have
$\varepsilon_{s}=\varepsilon_{x}=\varepsilon_{r}=\mathbf{0}$, which implies,
by (36), that $\varepsilon_{e}=\mathbf{0}$. Thus, by (35), $z=y$.

Now consider a system with an endemic equilibrium where $\exists i\in[n]$ such
that $\beta_{i}<1$. If $\beta_{i}\geq\gamma_{i}$ for all $i\in[n]$, there
exists some constant scalar $\eta\in\mathbb{R}$ such that by (6),

$\displaystyle\eta\left(Ar^{*}-(BX^{*}+\Gamma)s^{*}+N^{-1}W\Gamma
Ns^{*}\right)={\color[rgb]{0,0,0}\mathbf{0}}$
$\displaystyle\eta\left(BX^{*}s^{*}-(\Sigma+\Gamma)e^{*}+N^{-1}W\Gamma
Ne^{*}\right)={\color[rgb]{0,0,0}\mathbf{0}},$

where $\eta\beta_{i}\geq 1$ and $\eta\gamma_{i}\leq 1$ for all $i\in[n]$.
Thus, any endemic system with parameters where $\exists i\in[n]$ such that
$\beta_{i}<1$ can be mapped to one with the same endemic equilibrium where
$\beta_{i}\geq 1$ for all $i\in[n]$. ∎

######  Remark 1.

It should be noted that the condition $\beta_{i}\geq\gamma_{i},\forall
i\in[n],$ for uniqueness of the endemic equilibrium is only sufficient and
simulations show that an endemic equilibrium can still exist even when
$\exists i\in[n]$ such that $\beta_{i}<\gamma_{i}$.

The analysis throughout this section establishes several key points. First, we
find that short of a complete lockdown of infected nodes before spreading
begins, manipulating travel flows alone cannot change the long-term stability
of the healthy state of this system. Further, we have established conditions
under which the healthy state is exponentially stable, the healthy state is
unstable, an endemic equilibrium exists, and the endemic equilibrium is
unique.

##  5 Applications of Travel Flows Model

In this section, we shift our focus to constructing novel tools that will
enable us to leverage real-world data in Section 6. The proposed tools are 1)
learning the infection spread parameters given infection states and travel
flows data and 2) using travel flows data to predict the arrival time of a
disease to every node in the network given a disease origin. In this section,
we illustrate these tools using simulated data. We share the simulation code
used for this section on GitHub, which can be found in Butler (2023).

###  5.1 Parameter Identification

In this section we present a method for estimating the infection parameters
$\beta_{i},\sigma_{i},\delta_{i},\alpha_{i}$ for each node $i$, respectively,
given measurement data. Since infection data is generally collected at
discrete time intervals, we first discretize the system in (4) using Euler’s
method which yields

$\displaystyle s_{i}^{k+1}$
$\displaystyle=s_{i}^{k}+h\left(\alpha_{i}r_{i}^{k}-(\beta_{i}x_{i}^{k}+\gamma_{i}^{k})s_{i}^{k}+\sum_{j\neq
i}\frac{N_{j}}{N_{i}}w^{k}_{ij}\gamma^{k}_{j}s^{k}_{j}\right)$ (41a)
$\displaystyle e_{i}^{k+1}$
$\displaystyle=e_{i}^{k}+h\left(\beta_{i}x_{i}^{k}s_{i}^{k}-(\sigma_{i}+\gamma_{i}^{k})e_{i}^{k}+\sum_{j\neq
i}\frac{N_{j}}{N_{i}}w^{k}_{ij}\gamma^{k}_{j}e^{k}_{j}\right)$ (41b)
$\displaystyle x_{i}^{k+1}$
$\displaystyle=x_{i}^{k}+h\left(\sigma_{i}e_{i}^{k}-(\delta_{i}+\gamma_{i}^{k})x_{i}^{k}+\sum_{j\neq
i}\frac{N_{j}}{N_{i}}w^{k}_{ij}\gamma^{k}_{j}x^{k}_{j}\right)$ (41c)
$\displaystyle r_{i}^{k+1}$
$\displaystyle=r_{i}^{k}+h\left(\delta_{i}x_{i}^{k}-(\alpha_{i}+\gamma_{i}^{k})r_{i}^{k}+\sum_{j\neq
i}\frac{N_{j}}{N_{i}}w^{k}_{ij}\gamma^{k}_{j}r^{k}_{j}\right),$ (41d)

where $h>0$ is a sampling parameter and $k\in\mathbb{Z}_{\geq 0}$ is a given
time index. For this application, we assume a sampling parameter $h$ that is
small enough such that the system remains well-defined given the relative
values of the model parameters. For a more detailed discussion and analysis of
exactly when this condition holds, see our previous work in Butler et al.
(2021). Given $T$ samples of infection state data, we can estimate the
parameters for node $i$ as follows. Let

$\displaystyle\Delta
q_{i}=\begin{bmatrix}q^{1}_{i}-q^{0}_{i}+h(\gamma^{0}_{i}q^{0}_{i}-\sum_{j\neq
i}\frac{N_{j}}{N_{i}}w^{0}_{ij}\gamma^{0}_{j}q^{0}_{j})\\\ \vdots\\\
q^{T}_{i}-q^{T-1}_{i}+h(\gamma^{T-1}_{i}q^{T-1}_{i}-\sum_{j\neq
i}\frac{N_{j}}{N_{i}}w^{T-1}_{ij}\gamma^{T-1}_{j}q^{T-1}_{j})\end{bmatrix}$
(42)

where $q_{i}\in\\{s_{i},e_{i},x_{i},r_{i}\\}$, and define

$\displaystyle\Psi_{i}=h\begin{bmatrix}-s_{i}^{0}x_{i}^{0}&0&0&r_{i}^{0}\\\
\vdots&\vdots&\vdots&\vdots\\\ -s_{i}^{T-1}x_{i}^{T-1}&0&0&r_{i}^{T-1}\\\
s_{i}^{0}x_{i}^{0}&-e_{i}^{0}&0&0\\\ \vdots&\vdots&\vdots&\vdots\\\
s_{i}^{T-1}x_{i}^{T-1}&-e_{i}^{T-1}&0&0\\\ 0&e_{i}^{0}&-x_{i}^{0}&0\\\
\vdots&\vdots&\vdots&\vdots\\\ 0&e_{i}^{T-1}&-x_{i}^{T-1}&0\\\
0&0&x_{i}^{0}&-r_{i}^{0}\\\ \vdots&\vdots&\vdots&\vdots\\\
0&0&x_{i}^{T-1}&-r_{i}^{T-1}\\\ \end{bmatrix},$ (43)

then, we rewrite the system in (41) for a node $i$ given $T$ data samples of
data in terms of the spread parameters as

$\displaystyle\begin{bmatrix}\Delta s_{i}\\\ \Delta e_{i}\\\ \Delta x_{i}\\\
\Delta r_{i}\end{bmatrix}=\Psi_{i}\begin{bmatrix}\beta_{i}\\\ \sigma_{i}\\\
\delta_{i}\\\ \alpha_{i}\end{bmatrix}.$ (44)

Thus, we can solve for the parameters
$\beta_{i},\sigma_{i},\delta_{i},\alpha_{i}$ as

$\displaystyle\begin{bmatrix}\beta_{i}\\\ \sigma_{i}\\\ \delta_{i}\\\
\alpha_{i}\end{bmatrix}=\Psi_{i}^{\dagger}\begin{bmatrix}\Delta s_{i}\\\
\Delta e_{i}\\\ \Delta x_{i}\\\ \Delta r_{i}\end{bmatrix},$ (45)

where $\Psi_{i}^{\dagger}$ is the pseudo-inverse of $\Psi_{i}$. Thus, we have
that a solution to the spread parameters using (45) is unique if and only if
(43) is full column rank.

Using this method of parameter estimation, we can utilize real infection data
along with real travel flows data to determine the spreading rates for a given
disease outbreak, assuming that the model defined in Section 3 approximates
the behavior of the outbreak. We show an example of parameter identification
in Figure 2, where data from a 5-node simulated system with initial conditions
$s^{0}=[0.549,0.715,0.603,0.545,0.424]^{\top}$, $e^{0}=\mathbf{1}-s^{0}$, and
$x^{0}=r^{0}=\mathbf{0}$ is perturbed with Gaussian noise with zero mean and
standard deviation 0.01 and then used to learn spread the infection spread
parameters as outlined in this section, where the travel weight matrix used is
given in Table 1 and the spread parameters and RMSE of the learned parameters
are given in Table 2.

$w_{ij}$ | 1 | 2 | 3 | 4 | 5
---|---|---|---|---|---
1 | 0 | 0.212 | 0.275 | 0.25 | 0.212
2 | 0.249 | 0 | 0.26 | 0.299 | 0.338
3 | 0.246 | 0.198 | 0 | 0.204 | 0.178
4 | 0.285 | 0.29 | 0.259 | 0 | 0.272
5 | 0.22 | 0.299 | 0.206 | 0.247 | 0
Table 1: Travel connection matrix used in the simulated system shown in Figure 2. Node | 1 | 2 | 3 | 4 | 5 | RMSE
---|---|---|---|---|---|---
$\beta_{i}$ | 0.065 | 0.044 | 0.089 | 0.096 | 0.038 | 0.03
$\sigma_{i}$ | 0.079 | 0.053 | 0.057 | 0.093 | 0.007 | 0.01
$\delta_{i}$ | 0.001 | 0.001 | 0.008 | 0.008 | 0.009 | 0.002
$\alpha_{i}$ | 0.01 | 0.008 | 0.005 | 0.008 | 0.001 | 0.012
$\gamma_{i}$ | 0.002 | 0.002 | 0.002 | 0.002 | 0.005 | -
Table 2: Parameters used in the simulated system shown in Figure 2 along with
the RMSE of learned parameters. Figure 2: Simulated 5-node system from (41)
with parameters and travel matrix defined in Tables 2 and 1 where $n=5$,
$h=1$, and added mean zero Gaussian noise with standard deviation $0.01$ and
spread parameters chosen to induce an endemic state. Each solid line
represents the perturb simulated data for each node. Parameters from the noisy
system are identified and then re-simulated from the same initial condition,
where the predicted state trajectories for each node are shown by dashed lines
of the same color.

###  5.2 Effective Distance and Travel Flows

Given our notion of population flow and the definition of $w_{ij}$ being the
proportion of traveling members of the population at node $j$ flowing to node
$i$, we can use this proportional flow information to compute the effective
distance between any two nodes based on the most probable path an individual
would take between nodes. Further, we can then use the effective distance
between nodes to make predictions of the arrival time of an outbreak to a
given node based on the distance of that node from other infected nodes. A
similar notion of effective distance based on travel and using it to predict
disease arrival time is introduced and used in Brockmann and Helbing (2013),
where the effective distance between countries is computed using flight
traffic data. In our case, we apply this measure to inter-county travel data
in Minnesota during a period of the COVID-19 pandemic, which will be shown in
Section 6.

We compute the effective distance between nodes as follows. We interpret the
mobility matrix defined by the $w_{ij}$ as the transition matrix of a discrete
time homogeneous stochastic jump process, where at each time step a randomly
moving particle moves from node $j$ to node $i$ with probability $w_{ij}$. Let
$n_{k}$ be the location of this particle at time $k$. Now, consider the path
of a randomly moving particle after $L$ steps, forming the path
$\Xi=\\{n_{0},n_{1},\dots,n_{L}\\}$. If we fix any two endpoints in the graph
where $\Xi_{ij}$ is any path with $n_{0}=j$ and $n_{L}=i$, we can associate
any path taken between these points with a probability

$\displaystyle F(\Xi_{ij})={\color[rgb]{0,0,0}w_{in_{L-1}}\times\dots\times
w_{n_{1}j}}=\prod_{k=1}^{L}w_{n_{k}n_{k-1}}.$ (46)

Since the proposed notion of distance is additive, we establish a connection
between the probability of a given path and the effective distance of a given
path by employing the multiplicative property of logarithms

$D_{F(\Xi_{ij})}=-\log F(\Xi_{ij})=-\sum_{k=1}^{L}\log w_{n_{k}n_{k-1}},$ (47)

where $D_{F(\Xi_{ij})}$ is the log sum of the probabilities of each step taken
on the given path. Thus, since $w_{ij}\in[0,1]$ a low probability path will
equate to a high effective distance, with $D_{F(\Xi)}\rightarrow\infty$ as
$F(\Xi)\rightarrow 0$, and high probability paths result in lower effective
distance. We construct a new graph with edge weights between nodes being the
one-step negative log probabilities of traveling between nodes as

$d_{ij}=\begin{cases}0&i=j\\\ -\log w_{ij}&w_{ij}>0\\\
\infty&w_{ij}=0.\end{cases}$ (48)

Using this graph with edge weights between nodes $i$ and $j$ defined by (48),
we compute the effective distance as

$D_{ij}=\min_{\Xi_{ij}}D_{F(\Xi_{ij})}\geq 0,$ (49)

where $D_{ij}$ is the sum of the path lengths between nodes $j$ and $i$ that
yield the smallest distance, which can be found for any node $i\in[n]$ by
computing the minimum spanning tree for the distance graph from node $j$,
since $D_{F(\Xi_{ij})}=\sum_{k=1}^{L}d_{n_{k}n_{k-1}}$ describes the cost of
traveling any path between $i$ and $j$ in the graph. It is important to note
that since it is possible $d_{ij}\neq d_{ji}$, we also have that $D_{ij}\neq
D_{ji}$. Further, the shortest path to and from any two endpoints may not
necessarily contain the same intermediate points.

We can extend the notion of the most probable path and effective distance
further by considering the most probable path from a subgroup of nodes to any
given node. Since we are interested in the effective distance for an infection
to travel to uninfected nodes with respect to an infected group of nodes, we
can construct a node subgroup as follows

$\mathcal{X}=\big{\\{}i\in[n];x_{i}>p\big{\\}},$ (50)

where $x_{i}$ is the infection level at node $i$ and $p\in(0,1]$ is some
threshold at which we consider the node to be infected. We can then compute
the probability that an individual from node group $\mathcal{X}$ travels to
node $i$ as

$\tilde{w}_{i\mathcal{X}}=\frac{1}{\sum_{l\in\mathcal{X}}N_{l}}\sum_{j\in\mathcal{X}}N_{j}w_{ij}.$
(51)

Thus, we can construct another effective distance graph, with respect to a
group of nodes $\mathcal{X}$, using edge weights

$\tilde{d}_{ij}^{\mathcal{X}}=\begin{cases}0&i\in\mathcal{X}\text{ or }i=j\\\
-\log\tilde{w}_{i\mathcal{X}}&i\notin\mathcal{X},j\in\mathcal{X},\tilde{w}_{i\mathcal{X}}>0\\\
-\log w_{ij}&i,j\notin\mathcal{X},w_{ij}>0\\\
\infty&\text{otherwise},\end{cases}$ (52)

making the effective distance in the case from any node $j$ to every other
node $i\in[n]$, with respect to the node group ${\mathcal{X}}$,

$\tilde{D}_{ij}^{\mathcal{X}}=\min_{\Xi_{ij}}\big{(}\tilde{d}_{jn_{L-1}}^{\mathcal{X}}+\dots+\tilde{d}_{n_{1}i}^{\mathcal{X}}\big{)}.$
(53)

When a disease outbreak occurs, one important question for policymakers is how
long it will take for the disease to infect new nodes given past infection
data. Given travel data and knowledge of which nodes are currently infected,
we can use (53) to predict the arrival time of the disease to uninfected nodes
as follows. Let $T_{k}$ be the time of the $k$th arrival of the disease to an
uninfected node, and suppose $k>\tau$ disease arrivals have occurred, where
$\tau$ is the size of your training window:

  1. 1.

At time $T_{k}$, select training data that includes the disease arrivals from
time $T_{k-\tau}$ to $T_{k}$, where
$\tau=|\\{T_{k-\tau},T_{k-\tau+1},\dots,T_{k}\\}|$

  2. 2.

Compute the effective distance of each node in the training data to the
infected group $\mathcal{X}$ at time $T_{k-\tau}$

  3. 3.

Given the arrival times of the training data and the computed effective
distances to the infected group, fit a line to these training data

  4. 4.

Using the fitted line, predict the arrival time of the remaining uninfected
nodes given their effective distance to the current infected group
$\mathcal{X}$ at time $T_{k}$.

To ensure that this method always predicts sensible arrival times, as the
fitted line may predict arrival times that have already passed, a constant
shift is added to the predictions such that the next predicted arrival time is
greater than zero. This method differs from Brockmann and Helbing (2013) most
significantly in its ability to make updated predictions on the arrival time
of an outbreak, which depends on the current flow rate. Further, it
facilitates an updated notion on the distance of an outbreak from any given
node based on the likelihood of transmission from multiple source nodes as the
epidemic progresses, described by (52), rather than relying on the distance
from the original source. These factors contribute to a more accurate
prediction of arrival time in a more realistic evolving scenario.

We illustrate this algorithm by simulating the spread of a disease over a
network using our model as defined in Section 3. We construct our network to
emulate counties in Minnesota with travel between counties given by inter-
county travel data collected during part of the COVID-19 pandemic. A full
description of this data is given in Section 6. The system is then simulated
using infection parameters satisfying Assumption 2 and sampling parameter
$h=1$ (to emulate a daily resolution of data collection), where a node is
selected as the origin of the disease and predictions on the next nodes to be
infected are made. Using (49), we compute the effective distance of each node
from the origin of infection and plot the effective distance versus the
arrival time of the disease to each node, shown in Figure 3. Note the linear
relationship between effective distance and arrival time.

Figure 3: The arrival time of the disease to each node versus the effective
distance of each node in our simulated system from the disease origin,
computed by (49). The dashed line shows the linear fit of all data points.

In Figure 4 we show an example of using the shifting window algorithm to
predict the remaining infection arrival times with an arrival window of
$\tau=20$ disease arrival events. To evaluate the prediction error of our
shifting window method, we compute the root-mean-square (RMS) error of the
linear fit for all arrival times versus effective distance from the origin at
the start of the epidemic (with the linear fit shown in Figure 3) and compare
it with the average RMS error of the shifting window method that predicts the
next 10 arrival times whenever a new arrival occurs. This comparison yields an
RMS error of 4.59 timesteps (analogous to days since $h=1$) for the linear fit
of all arrival times versus effective distance from the disease origin (which
can only be evaluated after all data has been collected) versus an average RMS
error of 1.91 timesteps for the shifting window predicting the next 10 arrival
times whenever a new arrival occurs, which is a 58% reduction in prediction
error in addition to being usable in real-time to predict the next immediate
disease arrivals versus analysis after the fact.

(a) $T=T_{50}$

(b) $T=T_{60}$

(c) $T=T_{70}$

(d) $T=T_{80}$

Figure 4: Example of shifting prediction window of arrival times in the
simulated system at four different arrival times during the spreading process
($T=T_{50},T_{80}$) where training window $\tau=20$, with the data used to
make predictions (black circles), true arrival times for each node versus
effective distance (blue squares), and shifted best-fit predictions (red
triangles) shown at several arrival intervals, where the effective distances
are computed at each disease arrival using (53).

##  6 Case Study: Minnesota and COVID-19

We apply our model developed in Section 3 and methods described in Section 5
to data collected from the beginning of the 2020 COVID-19 pandemic in
Minnesota. Daily infection data at the county level is taken from usafacts.org
which contains the number of confirmed COVID-19 cases over the period of March
2020 to December 2020 MNd . To infer the infection states
$s_{i},e_{i},x_{i},r_{i}$ for each county $i$ from the confirmed cases data we
apply a constant shift that assumes confirmed individuals became exposed one
week before testing positive, at which point they are considered infected.
Then, we assume that infected individuals remain infected for one week after
testing positive, after which they are considered recovered. Finally, to
account for the rate at which individuals can lose immunity, we transfer
recovered individuals back into being susceptible after six weeks of being
recovered. For travel data, we worked with the Minnesota Department of
Transportation to obtain data collected by StreetLight str  that estimates the
number of trips between counties via anonymous geolocalization using
smartphones. This provides an estimate of the total number of trips made by
individuals between counties over a specified period of time. We choose a
weekly time scale in an effort to average out periodic behaviors such as
commuting and use this average to estimate the daily flow of individuals
between counties.

###  6.1 Parameter Identification

Using the inferred infection states for each county in MN along with the
travel flows computed from the collected traffic data, we can use (45) to
solve for the best-fit infection parameters using a convex solver and
constraining solutions to be strictly non-negative using our model defined in
Section 3. We then re-simulate the model using the same initial conditions and
learned parameters, with the results of the simulation compared with the
inferred infection states shown in Figure 5, where we show a subset of real
infection levels at five counties (solid lines) and the predicted infection
levels from the learned parameters (dashed lines). Note that although the
predicted infection peak of the model with learned parameters is much higher
than the actual data, we see that the peak infection time is more closely
approximated by the travel flows model, as shown in Table 3.

Regarding the error in magnitude prediction, there are several likely
contributing factors. First, in our estimation algorithm, we assume spread
parameters ($\beta_{i},\sigma_{i},\delta_{i},\alpha_{i}$) that are not time-
varying, which is most likely not the case in reality since lock-downs, social
distancing policies, medical interventions, and virus mutations can all change
the spread parameters over time, which violates this assumption. Additionally,
we are assuming that our rather naive method for state estimation (by choosing
a priori the recovering rate and the rate at which exposed individuals become
infected) yields the correct states for all counties and that all counties are
reporting infection cases in the same manner, which may not be true in all
cases. Thus, these combinations of structural assumptions and naivety in state
estimation are both likely large contributors to the magnitude error in state
prediction.

For comparison, we apply this parameter fit with another similar model,
modified from our previous work in Vrabac et al. (2022) to allow individuals
to become reinfected. The major difference between our model in Section 3 and
this modified model from Vrabac et al. (2022) is that travel data is used to
infer direct interaction between sub-populations rather than explicitly
modeling the flow of individuals between subpopulations. We then take the
best-fit parameters for this model given the infection and travel data and
simulate the model with the learned parameters, shown in Figure 6, where we
show a subset of real infection levels at the same five counties as Figure 5
(solid lines) and the predicted infection levels from the learned parameters
(dashed lines). Note that in this case, although the error is smaller in
magnitude than the prediction from the parameter fitting of the travel flows
model, the direct interaction model from Vrabac et al. (2022) fails to
reproduce the peaking behavior shown clearly in the inferred infection states
from the COVID-19 data and in Table 3.

| Infection RMSE | Peak Time RMSE
---|---|---
SEIRS (Travel Flows) | 0.023 | 4.15 Weeks
SEIRS (Vrabac et al. (2022)) | 0.007 | 21.21 Weeks
Table 3: The infection RMSE, defined as
$\sqrt{\frac{1}{nT}\sum_{i\in[n],k\in[T]}\left(x_{i}^{k}-\hat{x}_{i}^{k}\right)^{2}}$,
where $n$ is the number of nodes, $T$ is the number of samples, and
$\hat{x}_{i}^{k}$ is the predicted infection level at node $i$ for time step
$k$, and the peak time RMSE, defined as
$\sqrt{\frac{1}{n}\sum_{i\in[n]}\left(\operatorname*{arg\,max}_{k}(x_{i}^{k})-\operatorname*{arg\,max}_{k}(\hat{x}_{i}^{k})\right)^{2}}$,
for the parameter fitting predictions of the SEIRS travel flows model defined
in Section 3 (shown in Figure 5) and modified SEIRS model from Vrabac et al.
(2022) (shown in Figure 6), respectively, on the COVID-19 infection data.
Figure 5: Simulating (4) for select MN counties using traffic data to estimate
inter-county travel flows (dashed lines) with best-fit parameters learned from
the estimated epidemic states using COVID-19 infection data (solid lines).
Figure 6: Simulating the modified direct connection model from Vrabac et al.
(2022) for select MN counties using traffic data to estimate inter-county
travel flows (dashed lines) with best-fit parameters learned from the
estimated epidemic states using COVID-19 infection data (solid lines).

###  6.2 Predicting Disease Arrival Time

We use the methods described in Section 5.2 and the collected travel flows
data to make predictions on the arrival time of COVID-19 to counties in MN. In
Figure 7 we compare the arrival time of the first reported case with both the
effective distance and geographic distance of each county to the disease
origin (i.e., the county with the first reported case in MN), with the best-
fit line shown for each. We compute the R-value for each best-fit line to be
$0.53$ and $0.48$ for the effective distance and geographic distance,
respectively, suggesting that the effective distance is a better predictor of
arrival time, albeit by a small margin, than geographic distance when
considering all arrival times with respect to the disease origin. Finally, we
demonstrate the shifting window prediction method on the MN COVID-19 data in
Figure 8. We find an RMS error of close to 19 days for the linear fit of all
arrival times versus effective distance from the disease origin versus an
average RMS error of 4 days for the shifting window predicting the next 10
arrival times whenever a new arrival occurs, yielding a significant $79\%$
reduction in prediction error in addition to enabling real-time arrival
predictions.

Figure 7: The arrival time of the disease to each node versus the effective
distance (top) of each county in MN from the county with the first detected
COVID-19 cases during the 2020 pandemic, computed by (49), compared with the
geographic distance of each county (bottom) from the disease origin. The
dashed line shows the linear fit of all data points for each case
respectively.

(a) $T=T_{45}$

(b) $T=T_{55}$

(c) $T=T_{65}$

(d) $T=T_{75}$

Figure 8: Examples of the shifting prediction window of arrival times of
COVID-19 for counties in MN at four different arrival times during the
pandemic ($T=T_{45},T_{75}$) where training window $\tau=20$, with the data
used to make predictions (black circles), true arrival times for each county
versus effective distance (blue squares), and shifted best-fit predictions
(red triangles) shown at several arrival intervals, where the effective
distances is computed at each disease arrival using (53).

##  7 Conclusion

In this work, we have developed a model for simulating the spread of an SEIRS
epidemic process using networked population flows as the propagation
mechanism. We provide sufficient conditions under which the healthy state of
the system will be locally stable or unstable and show via model analysis that
there exists no valid perturbation to the population flows that will change
the local stability of any healthy state. These results suggest that
manipulating population flow alone is insufficient to achieve the long-term
eradication of the disease in the defined system where a loss of immunity is
present. Conversely, the result also suggests that perturbing the population
flows will not prevent a system from achieving disease eradication (i.e.,
causing the healthy state to become unstable if it is already locally stable).
Further, we provide sufficient conditions under which the system will enter a
unique endemic state. Although the analytical results show that travel control
is insufficient for affecting the long-term behavior of the epidemic, we show
that travel flows are a good predictor of disease arrival time, which can
provide valuable model-agnostic estimates of the arrival times for a given
disease that are useful in allocating a fixed amount of medical resources. To
this end, we then develop tools for applying our proposed model to data, such
as spreading parameter identification and disease arrival time prediction and
illustrate these tools using a case study of both travel and infection data
from the counties in MN during part of the COVID-19 pandemic.

There are several directions for continued research on the work presented in
this paper. Future directions for analysis included a study on the global
stability of the healthy state, as only local stability is explored in this
work, as well as a study of the stability of the endemic equilibrium for this
system. Additionally, in this work, we have assumed that no infections occur
during the transit of individuals between populations, and adding this element
to our travel flows model may provide additional meaningful insights into the
role of travel in spreading infectious diseases. Further, since the use of
effective distance to predict disease arrival time is independent of the
epidemic model being used, a more thorough analysis can be performed comparing
the usefulness of effective distance in predicting arrival time for other
epidemic models with travel flows (such as networked SIR, SIS, SAIR, etc.).

## References

  * Ruan et al. [2006] Shigui Ruan, Wendi Wang, and Simon A Levin.  The effect of global travel on the spread of SARS.  _Mathematical Biosciences & Engineering_, 3(1):205, 2006.
  * Tatem et al. [2006] Andrew J Tatem, David J Rogers, and Simon I Hay.  Global transport networks and infectious disease spread.  _Advances in Parasitology_ , 62:293–343, 2006.
  * Lowe et al. [2014] James Lowe, Phillip Gauger, Karen Harmon, Jianqiang Zhang, Joseph Connor, Paul Yeske, Timothy Loula, Ian Levis, Luc Dufresne, and Rodger Main.  Role of transportation in spread of porcine epidemic diarrhea virus infection, united states.  _Emerging Infectious Diseases_ , 20(5):872, 2014.
  * Findlater and Bogoch [2018] Aidan Findlater and Isaac I Bogoch.  Human mobility and the global spread of infectious diseases: a focus on air travel.  _Trends in Parasitology_ , 34(9):772–783, 2018.
  * Hollingsworth et al. [2006] T Déirdre Hollingsworth, Neil M Ferguson, and Roy M Anderson.  Will travel restrictions control the international spread of pandemic influenza?  _Nature Medicine_ , 12(5):497–499, 2006.
  * Tang et al. [2009] Ming Tang, Zonghua Liu, and Baowen Li.  Epidemic spreading by objective traveling.  _EPL (Europhysics Letters)_ , 87(1):18005, 2009.
  * Maital and Barzani [2020] Shlomo Maital and Ella Barzani.  The global economic impact of COVID-19: A summary of research.  _Samuel Neaman Institute for National Policy Research_ , 2020:1–12, 2020.
  * Kaye et al. [2021] Alan D Kaye, Chikezie N Okeagu, Alex D Pham, Rayce A Silva, Joshua J Hurley, Brett L Arron, Noeen Sarfraz, Hong N Lee, Ghali E Ghali, Jack W Gamble, et al.  Economic impact of COVID-19 pandemic on healthcare facilities and systems: International perspectives.  _Best Practice & Research Clinical Anaesthesiology_, 35(3):293–306, 2021.
  * Ahmad et al. [2020] Tauseef Ahmad, Mukhtiar Baig, and Jin Hui.  Coronavirus disease 2019 (COVID-19) pandemic and economic impact.  _Pakistan journal of medical sciences_ , 36(COVID19-S4):S73, 2020.
  * Kumar et al. [2020] Sandeep Kumar, V Maheshwari, J Prabhu, M Prasanna, P Jayalakshmi, P Suganya, Benjula Anbu Malar, and R Jothikumar.  Social economic impact of COVID-19 outbreak in India.  _International Journal of Pervasive Computing and Communications_ , 2020.
  * Kermack and McKendrick [1927] William Ogilvy Kermack and Anderson G McKendrick.  A contribution to the mathematical theory of epidemics.  _Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character_ , 115(772):700–721, 1927.
  * Brauer [2008] Fred Brauer.  _Compartmental Models in Epidemiology_ , pages 19–79.  Springer Berlin Heidelberg, Berlin, Heidelberg, 2008.
  * Prasse and Van Mieghem [2020] Bastian Prasse and Piet Van Mieghem.  Network reconstruction and prediction of epidemic outbreaks for general group-based compartmental epidemic models.  _IEEE Transactions on Network Science and Engineering_ , 7(4):2755–2764, 2020.
  * Fagnani and Zino [2017] Fabio Fagnani and Lorenzo Zino.  Time to extinction for the SIS epidemic model: New bounds on the tail probabilities.  _IEEE Transactions on Network Science and Engineering_ , 6(1):74–81, 2017.
  * Chen et al. [2019] Shanshan Chen, Michael Small, and Xinchu Fu.  Global stability of epidemic models with imperfect vaccination and quarantine on scale-free networks.  _IEEE Transactions on Network Science and Engineering_ , 7(3):1583–1596, 2019.
  * Xiang et al. [2021] Yue Xiang, Yonghong Jia, Linlin Chen, Lei Guo, Bizhen Shu, and Enshen Long.  COVID-19 epidemic prediction and the impact of public health interventions: A review of COVID-19 epidemic models.  _Infectious Disease Modelling_ , 6:324–342, 2021.
  * Rahmandad et al. [2021] Hazhir Rahmandad, Tse Yang Lim, and John Sterman.  Behavioral dynamics of COVID-19: Estimating underreporting, multiple waves, and adherence fatigue across 92 nations.  _System Dynamics Review_ , 37(1):5–31, 2021.
  * Wang [2020] Yuanhao Wang.  Evaluating the effectiveness of COVID-19 prevention and control measures based on SEIR model.  In _Proceedings of the 2020 International Conference on Communications, Information System and Computer Engineering (CISCE)_ , pages 143–146. IEEE, 2020.
  * Chen et al. [2020] Yi-Cheng Chen, Ping-En Lu, Cheng-Shang Chang, and Tzu-Hsuan Liu.  A time-dependent SIR model for COVID-19 with undetectable infected persons.  _IEEE Transactions on Network Science and Engineering_ , 7(4):3279–3294, 2020.
  * Nowzari et al. [2016] Cameron Nowzari, Victor M Preciado, and George J Pappas.  Analysis and control of epidemics: A survey of spreading processes on complex networks.  _IEEE Control Systems Magazine_ , 36(1):26–46, 2016.
  * Qian and Ukkusuri [2021] Xinwu Qian and Satish V Ukkusuri.  Connecting urban transportation systems with the spread of infectious diseases: A Trans-SEIR modeling approach.  _Transportation Research Part B: Methodological_ , 145:185–211, 2021.
  * Salehi et al. [2015] Mostafa Salehi, Rajesh Sharma, Moreno Marzolla, Matteo Magnani, Payam Siyari, and Danilo Montesi.  Spreading processes in multilayer networks.  _IEEE Transactions on Network Science and Engineering_ , 2(2):65–83, 2015.
  * Qu and Wang [2017] Bo Qu and Huijuan Wang.  SIS epidemic spreading with heterogeneous infection rates.  _IEEE Transactions on Network Science and Engineering_ , 4(3):177–186, 2017.
  * Kang et al. [2018] Huiyan Kang, Mengfeng Sun, Yajuan Yu, Xinchu Fu, and Bocheng Bao.  Spreading dynamics of an SEIR model with delay on scale-free networks.  _IEEE Transactions on Network Science and Engineering_ , 7(1):489–496, 2018.
  * Nadini et al. [2018] Matthieu Nadini, Alessandro Rizzo, and Maurizio Porfiri.  Epidemic spreading in temporal and adaptive networks with static backbone.  _IEEE Transactions on Network Science and Engineering_ , 7(1):549–561, 2018.
  * Li et al. [2021] Jian Li, Tao Xiang, and Linghui He.  Modeling epidemic spread in transportation networks: A review.  _Journal of Traffic and Transportation Engineering (English Edition)_ , 8(2):139–152, 2021.
  * Yang et al. [2012] Xu-Hua Yang, Bo Wang, Sheng-Yong Chen, and Wan-Liang Wang.  Epidemic dynamics behavior in some bus transport networks.  _Physica A: Statistical Mechanics and its Applications_ , 391(3):917–924, 2012.
  * Hackl and Dubernet [2019] Jürgen Hackl and Thibaut Dubernet.  Epidemic spreading in urban areas using agent-based transportation models.  _Future Internet_ , 11(4):92, 2019.
  * Paré et al. [2020] Philip E Paré, Carolyn L Beck, and Tamer Başar.  Modeling, estimation, and analysis of epidemics over networks: An overview.  _Annual Reviews in Control_ , 50:345–360, 2020.
  * Ye et al. [2020] Mengbin Ye, Ji Liu, Carlo Cenedese, Zhiyong Sun, and Ming Cao.  A network SIS meta-population model with transportation flow.  In _Proceedings of the IFAC World Congress_ , pages 2562–2567, 2020.
  * Yin et al. [2020] Qian Yin, Zhishuang Wang, Chengyi Xia, Matthias Dehmer, Frank Emmert-Streib, and Zhen Jin.  A novel epidemic model considering demographics and intercity commuting on complex dynamical networks.  _Applied Mathematics and Computation_ , 386:125517, 2020.
  * Brockmann and Helbing [2013] Dirk Brockmann and Dirk Helbing.  The hidden geometry of complex, network-driven contagion phenomena.  _Science_ , 342(6164):1337–1342, 2013.
  * Di Giamberardino et al. [2021] Paolo Di Giamberardino, Daniela Iacoviello, Federico Papa, and Carmela Sinisgalli.  A data-driven model of the COVID-19 spread among interconnected populations: Epidemiological and mobility aspects following the lockdown in Italy.  _Nonlinear Dynamics_ , 106:1239–1266, 2021.
  * Arino and van den Driessche [2003] Julien Arino and Pauline van den Driessche.  A multi-city epidemic model.  _Mathematical Population Studies_ , 10(3):175–193, 2003.
  * Liu and Stechlinski [2013] Xinzhi Liu and Peter Stechlinski.  Transmission dynamics of a switched multi-city model with transport-related infections.  _Nonlinear Analysis: Real World Applications_ , 14(1):264–279, 2013.
  * Knipl [2016] Diána Knipl.  Stability criteria for a multi-city epidemic model with travel delays and infection during travel.  _Electronic Journal of Qualitative Theory of Differential Equations_ , 2016(74):1–22, 2016.
  * Pujari and Shekatkar [2020] Bhalchandra S Pujari and Snehal Shekatkar.  Multi-city modeling of epidemics using spatial networks: Application to 2019-nCov (COVID-19) coronavirus in India.  _Medrxiv_ , pages 2020–03, 2020.
  * Arino and Van den Driessche [2006] Julien Arino and P Van den Driessche.  Disease spread in metapopulations.  _Fields Institute Communications_ , 48(1):1–13, 2006.
  * Wang and Li [2014] Lin Wang and Xiang Li.  Spatial epidemiology of networked metapopulation: An overview.  _Chinese Science Bulletin_ , 59(28):3511–3522, 2014.
  * Shao and Han [2022] Qi Shao and Dun Han.  Epidemic spreading in metapopulation networks with heterogeneous mobility rates.  _Applied Mathematics and Computation_ , 412:126559, 2022.
  * Zhu et al. [2021] Xuzhen Zhu, Yuxin Liu, Shengfeng Wang, Ruijie Wang, Xiaolong Chen, and Wei Wang.  Allocating resources for epidemic spreading on metapopulation networks.  _Applied Mathematics and Computation_ , 411:126531, 2021.
  * Vrabac et al. [2022] Damir Vrabac, Mingfeng Shang, Brooks Butler, Joseph Pham, Raphael Stern, and Philip E Paré.  Capturing the effects of transportation on the spread of COVID-19 with a multi-networked SEIR model.  _IEEE Control Systems Letters_ , 6:103–108, 2022.
  * Butler et al. [2021] Brooks Butler, Ciyuan Zhang, Ian Walter, Nishant Nair, Raphael Stern, and Philip E Paré.  The effect of population flow on epidemic spread: Analysis and control.  In _Proceedings of the 2021 60th IEEE Conference on Decision and Control (CDC)_ , pages 4260–4265, 2021.
  * Varga [1962] Richard S Varga.  _Iterative Analysis_.  Springer, 1962.
  * Horn and Johnson [2012] Roger A Horn and Charles R Johnson.  _Matrix Analysis_.  Cambridge University Press, 2012.
  * Shu et al. [2012] Hongying Shu, Dejun Fan, and Junjie Wei.  Global stability of multi-group SEIR epidemic models with distributed delays and nonlinear transmission.  _Nonlinear Analysis: Real World Applications_ , 13(4):1581–1592, 2012.
  * Hu et al. [2021] Maogui Hu et al.  Risk of Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2) Transmission Among Air Passengers in China.  _Clinical Infectious Diseases_ , 75(1):e234–e240, 09 2021.
  * Butler [2023] Brooks A. Butler.  Networked SEIRS with travel flows.  https://github.com/brooksbutler/SEIRS_PopulationFlows, 2023.
  * [49] Minnesota coronavirus cases and deaths.  https://usafacts.org/visualizations/coronavirus-{COVID-19}-spread-map/state/minnesota.  Accessed September 16, 2022.
  * [50] Street Light Data.  https://www.streetlightdata.com/.  Accessed November 2, 2022.