# Using the xr package on Overleaf

###### Abstract

We calculate epidemic thresholds and investigate the dynamics of a disease in
a networked metapopulation model. To study the specific role of mobility
levels and network geometry, we utilize the SIR-Network model and consider a
range of geometric structures. For _star-shaped_ networks where all nodes only
connect to a center, we obtain the same epidemic threshold formula as
previously found for fully connected networks in the case where all nodes have
the same infection rate except one. Next, we analyze _cycle-shaped_ networks
that yield different epidemic thresholds than star-shaped ones. We then
analyze more general classes of networks by combining the star, cycle, and
other structures, obtaining classes of networks with the same epidemic
threshold formulas. We present some conjectures on even more flexible networks
and complete our analysis by presenting simulations to explore the epidemic
dynamics for the different geometries.

Epidemic thresholds and disease dynamics in metapopulations: the role of
network geometry and human mobility.

Haridas K. Das111Department of Mathematics, Oklahoma State University,
Stillwater, OK 74078, e-mail: haridas.das$@$okstate.edu,222Department of
Mathematics, Dhaka University, Dhaka 1000, Bangladesh, e-mail:
<EMAIL_ADDRESS>and Lucas M. Stolerman333Department of Mathematics,
Oklahoma State University, Stillwater, OK 74078,
<EMAIL_ADDRESS>
##  1 Introduction

Epidemic transmission has become a significant public health concern in the
twenty-first century, with the recent COVID-19 pandemic servings as a
noteworthy example (?). It is widely acknowledged that globalization and
environmental factors have led to numerous infectious disease outbreaks (?, ?,
?). A recent study from Piret and Boivin (?) provided a comprehensive overview
of pandemics throughout history, noting that many of them occurred as a result
of increased global activities. Human mobility, in particular, is recognized
to play a significant role in shaping epidemic dynamics of respiratory (?, ?),
vector-borne (?, ?) and zoonotic diseases (?). The complex interplay between
disease transmission and contact patterns at multiple levels, from day-to-day
local interactions to international travels, has been a topic of intense study
in the past decades (?, ?, ?).

Historically, mathematical models have been used to explore the impacts of
human mobility in the dynamics of infectious diseases (?, ?, ?). The special
class of _metapopulation_ (or patch) models can include different
subpopulation sizes, movement rates, and location-specific parameters (?, ?,
?, ?, ?). In those models, the disease is assumed to spread in each
subpopulation, which can be divided into disease-related states (susceptible,
exposed, infected, among others). Sattenspiel and Dietz (?) introduced a
metapopulation demographic model with two distinct mobility patterns and
studied generalizations to describe a measles epidemic on the Caribbean island
of Dominica. Keeling et al. (?) investigated the role of individual identity
and its impact on epidemic dynamics. Other important contributions in the
field were made by Arino & Van Den Driessche (?) and Fulford et al. (?). More
recently, the SIR-network model was introduced in the context of dengue fever
epidemics (?). The simple structure of the SIR-Network model equations has
been used in a range of generalized systems, including models for COVID-19
epidemics with hospitalization and lockdown strategies (?), vaccination
behavior (?, ?), and spatially structured social dynamics (?).

Epidemic thresholds for infectious diseases are typically calculated with the
basic reproduction number $R_{0}$, the averaged cases generated by one
infected individual in a susceptible population (?, ?). In the past, several
studies have attempted to estimate $R_{0}$ for metapopulation models through
the next-generation matrix approach (?, ?, ?). Fulford, Roberts, and
Heesterbeek provided a methodology to compute $R_{0}$ for Tuberculosis in
Possums (?). Arino and Van Den Driessche (?) proposed a general multi-city SIS
model and computed bounds for $R_{0}$ with explicit formulas in particular
cases. Recently, more elaborated estimates of $R_{0}$ were proposed for
vector-borne diseases (?, ?). Despite the efforts to characterize epidemic
thresholds for disease-specific models with traditional methods, little
attention has been paid to the impact of the metapopulation network geometry
on these thresholds. For example, urban environments are intricately connected
through heterogeneous movement patterns, which can be described as complex
networks of transportation flow. These networks are far from equally
connected: they have big hubs, poorly connected nodes, and other nuanced
connections that can significantly impact the spread of infectious diseases.
Understanding the relationship between network geometry and epidemic dynamics
is thus critical for effective disease control and prevention efforts.

In this work, we calculated epidemic thresholds in a metapopulation model to
identify how the structure of connections between nodes and the intensity of
human movement may control or promote outbreaks. To this end, we built upon
previous work and established epidemic thresholds on various network
geometries using the SIR-network model (?). The model describes the dynamics
of a disease in a city, where people travel between neighborhoods. Each
neighborhood is represented by a node of a network and contains a population
divided into susceptible, infected, and recovered classes. To extend previous
results on fully connected networks, we assume that all nodes have the same
infection rate except one, which we call a _heterogeneous_ node. Our main
finding is the discovery of a class of networks with the same epidemic
thresholds as the fully connected network, in which the movement of people is
critical for epidemic control. Such class is led by _star-shaped_ networks
with the heterogeneous node in the center. In contrast, _cycle_ networks
exhibit different epidemic thresholds and thus are not part of the class. We
also performed numerical simulations to complement our theoretical analysis
and found significant differences between the disease dynamics depending on
the network structure. In particular, we observed how poorly connected
networks promoted delayed and damped epidemic peaks on nodes that are distant
from the heterogeneous nodes, while fully connected and star-shapes exhibited
more synchronized outbreaks.

The paper is organized as follows. First, we briefly review the SIR-network
model assumptions (Section 2) and the previously established epidemic
thresholds for fully connected networks (?) (Section 3). In section 4, we
introduce the new theoretical estimates for _star-shaped_ networks, _cycle_
networks, and the extended combinations of star-shaped and cycle (_cycle-
support_ and _star-cycle_ networks). We further establish more thresholds for
other networks such as _star-triangle_ and _star-background_ networks. We
complete our analysis with two conjectures about thresholds for the most
flexible network structures. In Section 5, we perform numerical simulations to
complement our theoretical results and investigate the epidemic dynamics
depending on the number of nodes (network size) and geometry. Finally, in
Section 6 we discuss our theoretical and numerical findings and their
practical importance.

##  2 The SIR-network model

In the SIR-network model (?), it is assumed that a city is divided into $M$
neighborhoods, which are the nodes of a network represented by the set
$V=\\{1,2,\ldots,M\\}$. The weighted edges of the network are the fractions
$\phi_{ij}\in[0,1]$ of residents moving daily between nodes $i$ and $j$. The
geometry of the network model can thus be set by the flux matrix
$\phi_{M\times M}=[\phi_{ij}]$, where $(i,j)\in V\times V$, or explicitly

$\phi_{M\times
M}=\begin{pmatrix}\phi_{11}&\phi_{12}&\phi_{13}&\ldots&\phi_{1M}\\\
\phi_{21}&\phi_{22}&\phi_{23}&\ldots&\phi_{2n}\\\
\vdots&\vdots&\vdots&\ddots&\vdots\\\
\phi_{M1}&\phi_{M2}&\phi_{M3}&\ldots&\phi_{MM}\end{pmatrix}.$ (1)

The flux matrix is stochastic by assumption since the population residents at
each node $i\in V$ is conserved; hence the fractions $\phi_{ij}$ (also called
fluxes) satisfy the following criterion:

$\sum_{j=1}^{M}\phi_{ij}=1.$

The populations of susceptible, infected, and removed individuals at node
$i\in V$ and time $t$ are given by $S_{i}(t),I_{i}(t),$ and $R_{i}(t)$
respectively. Each node $i$ has a total population $N_{i}$ and its
transmission rate $\beta_{i}$, which can be related to local transmission
probabilities and/or contact patterns driven by environmental factors, and
human behavior, among others. The model equations are given by

$\displaystyle\dot{S_{i}}(t)=-\sum_{j=1}^{M}\sum_{k=1}^{M}\beta_{j}\phi_{ij}S_{i}\frac{\phi_{kj}I_{k}}{N_{j}^{p}}$
(2)
$\displaystyle\dot{I_{i}}(t)=\sum_{j=1}^{M}\sum_{k=1}^{M}\beta_{j}\phi_{ij}S_{i}\frac{\phi_{kj}I_{k}}{N_{j}^{p}}-\gamma
I_{i}$ $\displaystyle\dot{R_{i}}(t)=\gamma I_{i},$

where $N_{j}^{p}=\sum_{k=1}^{M}\phi_{kj}N_{k}$ is the _present_ population at
node $j$ and $\gamma>0$ is the recovery rate for the disease. The basic
reproduction number $R_{0}$ is the largest eigenvalue of the next generation
matrix (NGM) (?, ?), which in this case can be calculated as

$\kappa=\frac{1}{\gamma}\begin{pmatrix}\sum_{j=1}^{M}\beta_{j}\phi_{1j}\phi_{1j}\frac{N_{1}}{N_{j}^{p}}&\sum_{j=1}^{M}\beta_{j}\phi_{1j}\phi_{2j}\frac{N_{1}}{N_{j}^{p}}&...&\sum_{j=1}^{M}\beta_{j}\phi_{1j}\phi_{Mj}\frac{N_{1}}{N_{j}^{p}}\\\
\sum_{j=1}^{M}\beta_{j}\phi_{2j}\phi_{1j}\frac{N_{2}}{N_{j}^{p}}&\sum_{j=1}^{M}\beta_{j}\phi_{2j}\phi_{2j}\frac{N_{2}}{N_{j}^{p}}&...&\sum_{j=1}^{M}\beta_{j}\phi_{2j}\phi_{Mj}\frac{N_{2}}{N_{j}^{p}}\\\
\vdots&\vdots&\ddots&\vdots\\\
\sum_{j=1}^{M}\beta_{j}\phi_{Mj}\phi_{1j}\frac{N_{M}}{N_{j}^{p}}&\sum_{j=1}^{M}\beta_{j}\phi_{Mj}\phi_{2j}\frac{N_{M}}{N_{j}^{p}}&...&\sum_{j=1}^{M}\beta_{j}\phi_{Mj}\phi_{Mj}\frac{N_{M}}{N_{j}^{p}}\\\
\end{pmatrix}.$ (3)

In the special case where the flux matrix is symmetric ($\phi_{ij}=\phi_{ji}$
for all $i$ and $j$ in $V$) and the populations at each node are equal
($N_{i}=N>0$ for all $i\in V$), we obtain the simplified present population

$N_{j}^{p}=\sum_{k=1}^{M}\phi_{kj}N_{k}=\left(\sum_{k=1}^{M}\phi_{kj}\right)N=N$

and the NGM becomes

$\kappa=\frac{1}{\gamma}\begin{pmatrix}\sum_{j=1}^{M}\beta_{j}\phi_{1j}\phi_{1j}&\sum_{j=1}^{M}\beta_{j}\phi_{1j}\phi_{2j}&...&\sum_{j=1}^{M}\beta_{j}\phi_{1j}\phi_{Mj}\\\
\sum_{j=1}^{M}\beta_{j}\phi_{2j}\phi_{1j}&\sum_{j=1}^{M}\beta_{j}\phi_{2j}\phi_{2j}&...&\sum_{j=1}^{M}\beta_{j}\phi_{2j}\phi_{Mj}\\\
\vdots&\vdots&\ddots&\vdots\\\
\sum_{j=1}^{M}\beta_{j}\phi_{Mj}\phi_{1j}&\sum_{j=1}^{M}\beta_{j}\phi_{Mj}\phi_{2j}&...&\sum_{j=1}^{M}\beta_{j}\phi_{Mj}\phi_{Mj}\\\
\end{pmatrix}.$ (4)

######  Remark 1.

Finding $R_{0}$ from the NGM approach is equivalent to performing a local
stability analysis of the disease-free equilibrium (DFE)
$(S^{*}_{i},I^{*}_{i},R^{*}_{i})=(N_{i},0,0)$ for $i\in V$. Substituting the
perturbation form $S_{i}(t)=S_{i}^{*}+\Delta S_{i}(t)$ and
$I_{i}(t)=I_{i}^{*}+\Delta I_{i}(t)$, where $i=1,\ldots,M$ into (Eq. 2) and
linearizing, we are able to obtain the following relation between the Jacobian
matrix $J$, and the next generation matrix $\kappa$

$J=\gamma(\kappa-I).$

To analyze the Jacobian matrix $J$ stability at the DFE, we find locally
asymptotic stable (LAS) DFE conditions when all eigenvalues have negative real
part or unstable DFE if at least one eigenvalue has a positive real part (?,
?).

##  3 Epidemic thresholds: previous results

This section briefly presents two theorems from Stolerman et al. (?). The
first theorem established the outbreak criteria for _homogeneous networks_ ,
where all infection rates are equal in each node.

######  Theorem (Homogeneous networks in (?)).

Let $\left(V,\phi_{M\times M}\right)$ be a network associated with the SIR-
network system (Eq. 2) in the case where all nodes have the same transmission
rate, $\beta_{j}=\beta_{0}\leavevmode\nobreak\ \leavevmode\nobreak\ \forall
j\in V$. Then basic reproduction number of the homogeneous network models is
$R_{0}^{\text{hom}}=\frac{\beta_{0}}{\gamma}$, and an epidemic outbreak occurs
if and only if $\beta_{0}>\gamma$.

The second theorem is devoted to fully connected networks in the case where a
single node has a different infection rate (heterogeneous node). The
particular flux matrix $\phi_{M\times M}$ in (Eq. 1) is given by

$\phi_{ij}=\begin{cases}\phi_{0}&\text{if $i\neq j$}\\\
1-(M-1)\phi_{0}&\text{if $i=j$},\end{cases}$ (5)

where $0<\phi_{0}\leq\frac{1}{M-1}$, since $\phi_{ij}>0$.

######  Theorem (Fully connected networks with one heterogeneous node (?)).

Let $\left(V,\phi_{M\times M}\right)$ be a fully connected network, where
$\phi_{ij}=\phi_{0}$ if $i\neq j$ and $\phi_{ii}=1-(M-1)\phi_{0}$. We suppose
that $\beta_{i}=\zeta\beta_{0}$ for $i\in V$ and $\beta_{k}=\beta_{0}$, for
all $k\in V$, where $k\neq i$. Under these conditions, if the background nodes
are stable, i.e., $R_{0}^{\text{hom}}<1$, there is a minimum value of $\zeta$
given by

$\zeta^{\text{crit}}=\left(\frac{M-(M-1)R_{0}^{\text{hom}}}{R_{0}^{\text{hom}}}\right)=1+\left(\frac{1-R_{0}^{\text{hom}}}{R_{0}^{\text{hom}}}\right)M.$

In the case $\frac{1}{R_{0}^{\text{hom}}}<\zeta<\zeta^{\text{crit}}$, there is
an interval

$\mathcal{I}=\left(\max\left\\{0,\frac{1}{M}-\tilde{\phi}\right\\},\min\left\\{\frac{1}{M-1},\frac{1}{M}+\tilde{\phi}\right\\}\right),$
(6)

where

$\displaystyle\tilde{\phi}=\frac{1}{M}\sqrt{\frac{M(1-R_{0}^{\text{hom}})-R_{0}^{\text{hom}}(\zeta-1)}{R_{0}^{\text{hom}}\left(M\zeta(1-R_{0}^{\text{hom}})-(\zeta-1)\right)}},$
(7)

such that if $\phi_{0}\in\mathcal{I}$, then the DFE is stable, and otherwise
it is unstable. Finally, when $\zeta<\frac{1}{R_{0}^{\text{hom}}}$, the DFE is
always stable.

When the infection rate $\beta_{0}$ is lower than the recovery rate $\gamma$,
theorem Theorem assures that the outbreak is controlled, which is the same
outbreak criteria of the SIR model. In theorem Theorem, Stolerman et al. (?)
_perturbed_ the homogeneous networks by adding a single node with a different
infection rate. In this case, the flux of people between nodes plays a vital
role in the epidemic criteria; there exists a critical value of the ratio
$\zeta=\frac{\beta_{M}}{\beta_{0}}$ such that there exists an interval of flux
values that control the epidemics. In section 4, we determine similar epidemic
thresholds for the SIR-network model in the different network geometries.

##  4 New epidemic thresholds

This section aims to introduce new classes of networks where the flux may
still play an important role in controlling epidemic outbreaks. Previously,
Stolerman et al. (?) considered the flux matrix (Eq. 5) to describe the _fully
connected_ networks where the fractions of movement between any two nodes are
given by a constant $\phi_{0}$. Therefore, the following question naturally
arises:

  * _What happens when the network is not fully connected?_

To answer that question and extend the results of Theorem 2 in (?), we explore
other classes of networks, still assuming that a single node has a different
infection rate. More precisely, given a network with $M$ nodes, we follow the
approach of (?) and define the positive quantities $\beta_{0}$ and $\zeta$
such that

$\beta_{i}=\begin{cases}\beta_{0}&\text{if $i=1,2,...,(M-1)$}\\\
\zeta\beta_{0}&\text{if $i=M$}\end{cases}$ (8)

Here $\zeta$ is a dimensionless multiplicative factor that modifies the
infection rate at the heterogeneous node $M$. In this study, we are
particularly interested in the case where all nodes are “stable” except one,
i.e, $\beta_{0}<\gamma$ and $\zeta\beta_{0}>\gamma$, so we may explore the
case where the single node $M$ can be unstable and thus promote an epidemic
outbreak in the whole network.

Fig 1 illustrates some of the networks we consider in this study. For _star-
shaped_ networks (subsection 4.1), we consider that the heterogeneous node $M$
(the one with a different infection rate) connects to all other nodes by the
same flux $\phi_{0}$, and all the outer nodes are not connected. These
networks mimic a city where a central node receives a daily flux of people. On
the other hand, _cycle-shaped_ networks (subsection 4.2), which can be
unidirectional or bidirectional, do not have such a central node. In this
case, all nodes are only connected to their immediate two neighbors,
potentially representing small towns where people usually commute short
distances. We also consider extensions of the cycle networks defining _cycle-
support_ and _star-cycle_ networks (subsection 4.3). Cycle support networks
are built by connecting the heterogeneous node $M$ in the cycle to the rest of
the nodes of the networks. On the other hand, for _star-cycle_ networks we
connect the immediate neighbors of the outer nodes in the star-shaped network.

Throughout this paper, we consider the SIR-network model with nonnegative and
symmetric flux matrix $\phi_{M\times M}$ (Eq. 1) and consider that all nodes
have the same population size, i.e, $N_{i}=N>0,\text{for}\
i=\leavevmode\nobreak\ 1,2,\ldots,M$. These assumptions allow us to solely
investigate the impact of the network geometry on the epidemic thresholds and
establish a comparison with previous results obtained in (?).

Figure 1: Different network geometries ($M=9$): All blue nodes have the same
infection rates $\beta_{0}>0$ except the red node, where
$\beta_{9}=\zeta\beta_{0}$. Here $\zeta>0$ is a dimensionless multiplicative
factor.

###  4.1 Star-shaped networks

We begin our analysis on star-shaped networks, where all nodes $i\in V$ are
uniquely connected to a central node $M$. A star-shaped network mimics a city
where a well-connected node (for instance, the downtown of a big city) has a
different infection rate and thus may influence the disease spread in the
network. The flux matrix $\phi_{M\times M}$ of a star-shaped network is
defined by

$\phi_{ij}=\begin{cases}\phi_{0}&\text{if $i=M,\leavevmode\nobreak\
\text{and}\leavevmode\nobreak\ j\in V\setminus M$}\leavevmode\nobreak\ \text{
or \leavevmode\nobreak\ if $j=M,\leavevmode\nobreak\
\text{and}\leavevmode\nobreak\ i\in V\setminus M$}\\\ 1-\phi_{0}&\text{if
$i=j,\leavevmode\nobreak\ \text{and}\leavevmode\nobreak\ \leavevmode\nobreak\
i,j\in V\setminus M$}\\\ 1-(M-1)\phi_{0}&\text{if $i=j=M$}\\\
0&\text{otherwise},\end{cases}$ (9)

where $0<\phi_{0}\leq\frac{1}{M-1}$ to ensure that all matrix elements are
nonnegative. We may also write the fluxes for Eq. 9 in the matrix format

$\phi_{M\times M}=\begin{pmatrix}1-\phi_{0}&0&0&\ldots&\phi_{0}\\\
0&1-\phi_{0}&0&\ldots&\phi_{0}\\\ 0&0&1-\phi_{0}&\ldots&\phi_{0}\\\
\vdots&\vdots&\vdots&\ddots&\vdots\\\
\phi_{0}&\phi_{0}&\phi_{0}&\ldots&1-(M-1)\phi_{0}\end{pmatrix}.$

For _star-shaped_ networks, we find the same minimum $\zeta$ (Eq. Theorem),
and interval $\mathcal{I}$ (Eq. 6) as the ones in the fully connected case
(Theorem 2 in (?)), where the flux between regions can control an epidemic
outbreak. We start with a proposition for $M=3$ to show a step-by-step
calculation using local stability analysis at the DFE.

######  Proposition 1.

Let $V=\left\\{1,2,3\right\\}$ be the nodes of the star-shaped network with
flux matrix $\phi_{3\times 3}$ be defined by (Eq. 9) for $M=3$, where
$\phi_{0}\in(0,\frac{1}{2}]$. For the infection rates $\beta_{i}$ (Eq. 8), we
assume that all nodes are stable except the center, i.e,
$R_{0}^{\text{hom}}=\frac{\beta_{0}}{\gamma}<1$ and
$\frac{\beta_{3}}{\gamma}=\frac{\zeta\beta_{0}}{\gamma}>1$. Then there exists
a critical value $\zeta^{\text{crit}}$ given by

$\zeta^{\text{crit}}=1+\left(\frac{1-R_{0}^{\text{hom}}}{R_{0}^{\text{hom}}}\right)3.$
(10)

such that if $\zeta<\zeta^{\text{crit}}$, there is an interval

$\mathcal{I}=\left(\max\left\\{0,\frac{1}{3}-\tilde{\phi}\right\\},\min\left\\{\frac{1}{2},\frac{1}{3}+\tilde{\phi}\right\\}\right),$
(11)

where

$\displaystyle\tilde{\phi}=\frac{1}{3}\sqrt{\frac{\gamma(3\gamma-2\beta_{0}-\zeta\beta_{0})}{\beta_{0}(\gamma+2\zeta\gamma-3\zeta\beta_{0})}}=\frac{1}{3}\sqrt{\frac{3(1-R_{0}^{\text{hom}})-R_{0}^{\text{hom}}(\zeta-1))}{R_{0}^{\text{hom}}\left(3\zeta(1-R_{0}^{\text{hom}})-(\zeta-1)\right)}},$
(12)

such that if $\phi_{0}\in\mathcal{I}$, then the DFE is stable, and otherwise
it is unstable.

###### Proof.

To establish the stability criteria for the DFE, we may perform a local
stability analysis of the DFE. The Jacobian matrix at the DFE corresponding to
the dynamical system of the $3\times 3$ star-shaped flux matrix is given by
the following expression:

$\displaystyle J_{3\times
3}=\begin{pmatrix}\beta_{1}(-1+\phi_{0})^{2}+\beta_{3}\phi_{0}^{2}-\gamma&\beta_{3}\phi_{0}^{2}&\phi_{0}\left(\beta_{1}+\beta_{3}-\beta_{1}\phi_{0}-2\beta_{3}\phi_{0}\right)\\\
\beta_{3}\phi_{0}^{2}&\beta_{2}(-1+\phi_{0})^{2}+\beta_{3}\phi_{0}^{2}-\gamma&\phi_{0}\left(\beta_{2}+\beta_{3}-\beta_{2}\phi_{0}-2\beta_{3}\phi_{0}\right)\\\
\phi_{0}\left(\beta_{1}+\beta_{3}-\beta_{1}\phi_{0}-2\beta_{3}\phi_{0}\right)&\phi_{0}\left(\beta_{2}+\beta_{3}-\beta_{2}\phi_{0}-2\beta_{3}\phi_{0}\right)&\beta_{3}(1-2\phi_{0})^{2}+(\beta_{1}+\beta_{2})\phi_{0}^{2}-\gamma\end{pmatrix}.$
(13)

Since symmetric matrices have only real eigenvalues, the eigenvalues of the
Jacobian matrix (Eq. 13) are real. We calculate the eigenvalues using the
symbolic packages of the software Mathematica, which are given by

$\displaystyle\lambda_{1}=\beta_{0}(-1+\phi_{0})^{2}-\gamma$
$\displaystyle\lambda_{2}=\frac{1}{2}\beta_{0}\left(1+\zeta-2\phi_{0}-4\zeta\phi_{0}+3\phi_{0}^{2}+6\zeta\phi_{0}^{2}+B-\gamma\right)$
$\displaystyle\lambda_{3}=\frac{1}{2}\beta_{0}\left(1+\zeta-2\phi_{0}-4\zeta\phi_{0}+3\phi_{0}^{2}+6\zeta\phi_{0}^{2}-B-\gamma\right)$

where
$B:=\sqrt{(1-2\phi_{0}+3\phi_{0}^{2})^{2}+\zeta^{2}\left(1-4\phi_{0}+6\phi_{0}^{2}\right)^{2}+2\zeta\left(-1+6\phi_{0}-\phi_{0}^{2}-24\phi_{0}^{3}+18\phi_{0}^{4}\right)}.$
Since $\lambda_{2}$ and $\lambda_{3}$ are real, the square root inside $B$
must be nonnegative for $\phi_{0}\in(0,\frac{1}{2}]$ and for any $\zeta$, $B$
is well defined and positive real number.

Next, we find the conditions for which all eigenvalues are negative (LAS DFE),
or at least one eigenvalue is positive (unstable DFE). First, since
$\phi_{0}\leq 1/2<1$, we have
$\lambda_{1}=\beta_{0}(-1+\phi_{0})^{2}-\gamma<\beta_{0}-\gamma<0$, since by
hypothesis $\frac{\beta_{0}}{\gamma}<1$. It is thus sufficient to check the
critical points of $\lambda_{2}$ and $\lambda_{3}$. However, since
$\lambda_{2}>\lambda_{3}$, it is thus sufficient to check the sign of
$\lambda_{2}$. It can be verified that $\phi_{0}=\frac{1}{3}$ is a critical
point of both eigenvalues $\lambda_{2}$, and $\lambda_{3}$, since
$\frac{\partial\lambda_{2}}{\partial\phi_{0}}(\phi_{0}=\frac{1}{3},\beta_{0},\zeta,\gamma)=0$,
and
$\frac{\partial\lambda_{3}}{\partial\phi_{0}}(\phi_{0}=\frac{1}{3},\beta_{0},\zeta,\gamma)=0$.
Moreover, we can compute the second derivatives at the critical point and
obtain

$\displaystyle\frac{\partial^{2}\lambda_{2}}{\partial\phi_{0}^{2}}(\phi_{0}=\frac{1}{3},\beta_{0},\zeta,\gamma)=\frac{12\beta_{0}(\zeta-1)^{2}}{2+\zeta}>0$
$\displaystyle\frac{\partial^{2}\lambda_{3}}{\partial\phi_{0}^{2}}(\phi_{0}=\frac{1}{3},\beta_{0},\zeta,\gamma)=\frac{54\zeta\beta_{0}}{2+\zeta}>0$

and so $\lambda_{2}$, and $\lambda_{3}$ have reached to the minimum value when
$\phi_{0}=\frac{1}{3}$. Finally,

$\displaystyle\lambda_{2}(\phi_{0}=\frac{1}{3},\beta_{0},\zeta,\gamma)=\frac{1}{6}\left(2\beta_{0}+\zeta\beta_{0}+\sqrt{(2+\zeta)^{2}}\beta_{0}-6\gamma\right)=\frac{1}{3}\left(\beta_{0}(2+\zeta)-3\gamma\right)$
$\displaystyle\lambda_{3}(\phi_{0}=\frac{1}{3},\beta_{0},\zeta,\gamma)=\frac{1}{6}\left(2\beta_{0}+\zeta\beta_{0}-\sqrt{(2+\zeta)^{2}}\beta_{0}-6\gamma\right)=-\gamma<0.$

In this way, if
$\lambda_{2}(\phi_{0}=\frac{1}{3},\beta_{0},\zeta,\gamma)=\frac{1}{3}\left(\beta_{0}(2+\zeta)-3\gamma\right)<0$,
then there exists an interval of $\phi_{0}$ values containing
$\phi_{0}=\frac{1}{3}$ such that all eigenvalues are negative and thus the DFE
is stable. The above condition yields the critical value $\zeta^{crit}$, i.e,

$\displaystyle\left(\beta_{0}(2+\zeta)-3\gamma\right)<0\Leftrightarrow\zeta<\frac{3\gamma-2\beta_{0}}{\beta_{0}}=1+\left(\frac{1-R_{0}^{\text{hom}}}{R_{0}^{\text{hom}}}\right)3:=\zeta^{\text{crit}}.$

In addition, the function
$\phi_{0}\mapsto\lambda_{2}(\phi_{0},\zeta,\beta_{0},\gamma)$ has two roots

$\displaystyle\phi_{c1}=\frac{1}{3}+\frac{1}{3}\sqrt{\frac{\gamma(3\gamma-2\beta_{0}-\zeta\beta_{0})}{\beta_{0}(\gamma+2\zeta\gamma-3\zeta\beta_{0})}}=\frac{1}{3}+\frac{1}{3}\sqrt{\frac{3(1-R_{0}^{\text{hom}})-R_{0}^{\text{hom}}(\zeta-1))}{R_{0}^{\text{hom}}\left(3\zeta(1-R_{0}^{\text{hom}})-(\zeta-1)\right)}},$
$\displaystyle\phi_{c2}=\frac{1}{3}-\frac{1}{3}\sqrt{\frac{\gamma(3\gamma-2\beta_{0}-\zeta\beta_{0})}{\beta_{0}(\gamma+2\zeta\gamma-3\zeta\beta_{0})}}=\frac{1}{3}-\frac{1}{3}\sqrt{\frac{3(1-R_{0}^{\text{hom}})-R_{0}^{\text{hom}}(\zeta-1))}{R_{0}^{\text{hom}}\left(3\zeta(1-R_{0}^{\text{hom}})-(\zeta-1)\right)}}.$

Therefore, taking into account the restriction
$\phi_{0}\in\left(0,\frac{1}{2}\right]$, the interval where there will be no
epidemic is given by
$\mathcal{I}(\zeta,R_{0}^{\text{hom}})=\left(\max\left\\{0,\frac{1}{3}-\tilde{\phi}\right\\},\min\left\\{\frac{1}{2},\frac{1}{3}+\tilde{\phi}\right\\}\right)$,
where

$\displaystyle\tilde{\phi}=\frac{1}{3}\sqrt{\frac{3(1-R_{0}^{\text{hom}})-R_{0}^{\text{hom}}(\zeta-1)}{R_{0}^{\text{hom}}\left(3\zeta(1-R_{0}^{\text{hom}})-(\zeta-1)\right)}}.$

It is worth noting that
$\frac{1}{R_{0}^{\text{hom}}}<\zeta<\zeta^{\text{crit}}$ guarantees that the
term inside of the above square root is nonnegative. Therefore, $\tilde{\phi}$
is a real number, and the interval $\mathcal{I}$ is well defined. This last
observation concludes the proof that the DFE is stable if
$\phi_{0}\in\mathcal{I}$, and otherwise unstable. ∎

The existence of a critical value $\zeta^{crit}$ suggests that the infection
rate at the central node plays a crucial role in determining the spread of the
epidemic. Specifically, when the value of $\zeta$ is sufficiently low, we
observe that the disease-free equilibrium (DFE) remains stable over a range of
flux values. On the other hand, if $\zeta>\zeta^{crit}$, then the DFE is
unstable for all flux values, a scenario where the central node promotes an
outbreak given the sufficiently high infection rate. Finally, if
$\zeta<\frac{1}{R_{0}^{\text{hom}}}$, an algebraic manipulation can be
performed to show that $\tilde{\phi}>\frac{1}{3}$ and thus
$\mathcal{I}(\zeta,R_{0}^{\text{hom}})=(0,\frac{1}{2})$, which means that the
epidemic is controlled for all flux values. This is consistent with our
intuition since in this case all nodes would be stable.

We also note that the threshold condition for the local stability of the DFE
can be written in terms of a critical infection rate
$\beta_{0}^{\text{crit}}$, which is given by

$\beta_{0}^{\text{crit}}=\gamma\left(\frac{3}{2+\zeta}\right).$ (14)

In section 5, we perform numerical simulations changing the $\beta_{0}$ and
observing the emergence of the $\phi_{0}$ intervals where the epidemic is
controlled.

Now we state and prove a general theorem for star-shaped networks of arbitrary
size. The particular case for $M=3$ is pedagogical since it highlights the
mathematical nature of such _flux-driven_ epidemic control. The main tool is
the classic local stability analysis, which is ultimately equivalent to the
next-generation matrix approach (?). In fact, the central argument of the
proof of the next theorem is the particular decomposition of the next-
generation matrix in the same format as found by Stolerman et al. (?).

######  Theorem 1.

Let $V=\left\\{1,2,3,\ldots,M\right\\}$ be the nodes of a _star-shaped_
network. We denote the infection rate $\beta_{i}=\beta_{0}$ at node $i$ for
all $i=1,2,3,\cdots,M-1$ and $\beta_{M}=\zeta\beta_{0}$. We also assume that
all the nodes are stable except the center $M$, i.e,
$R_{0}^{\text{hom}}=\frac{\beta_{0}}{\gamma}<1$ and
$\frac{\beta_{M}}{\gamma}=\frac{\zeta\beta_{0}}{\gamma}>1$. Then we have a
minimum $\zeta^{\text{crit}}$ given by

$\zeta^{\text{crit}}=1+\left(\frac{1-R_{0}^{\text{hom}}}{R_{0}^{\text{hom}}}\right)M.$
(15)

such that if $\frac{1}{R_{0}^{\text{hom}}}<\zeta<\zeta^{\text{crit}}$, there
exists an interval
$\mathcal{I}=\left(\max\left\\{0,\frac{1}{M}-\tilde{\phi}\right\\},\min\left\\{\frac{1}{M-1},\frac{1}{M}+\tilde{\phi}\right\\}\right)$,
where $\tilde{\phi}$ is given by (Eq. 7), such that if
$\phi_{0}\in\mathcal{I}$, then the DFE is stable, and otherwise it is
unstable. Finally, if $\zeta<\frac{1}{R_{0}^{\text{hom}}}$ i.e.
$\beta_{M}<\gamma$, then the DFE is always stable.

###### Proof.

We first find the next generation matrix (NGM) $\kappa=(\kappa_{ij})_{M\times
M}$, where
$\kappa_{ij}=\frac{1}{\gamma}\sum_{j=1}^{M}\beta_{j}\phi_{ij}\phi_{kj}\frac{N_{i}}{N_{j}^{p}}$
(Eq. 3). Similar to Section 2, we start from the premise that the populations
of each node are equal. i.e. $N_{i}=N>0$ for all $i\in V$. We also consider a
different infection rate in the central node (Eq. 8).

Now since we are assuming equal populations and symmetric flux matrix, the NGM
can be written as in Eq. 4). For a general star-shaped network, the NGM
entries can be explicitly calculated.

In fact, for any $u\in{1,2,3,4,5,...,M-1}$, and $v\in{1,2,3,4,5,...,M-1}$, and
$u=v$, then we have the following diagonal elements of the $\kappa$ are

$\displaystyle\kappa_{uu}=\sum_{j=1}^{M}\beta_{j}\phi_{uj}\phi_{uj}=\sum_{j=1}^{M}\beta_{j}\phi_{uj}^{2}$
$\displaystyle=\beta_{1}\phi_{u1}^{2}+\beta_{2}\phi_{u2}^{2}+\beta_{3}\phi_{u3}^{2}+\ldots+\beta_{u}\phi_{uu}^{2}+\ldots+\beta_{M}\phi_{uM}^{2}$
$\displaystyle=\beta_{0}.0^{2}+\beta_{0}.0^{2}+\beta_{0}.0^{2}+\ldots+\beta_{0}(1-\phi_{0})^{2}+\ldots+\zeta\beta_{0}\phi_{0}^{2}$
$\displaystyle=\beta_{0}\left(\zeta\phi_{0}^{2}+(1-\phi_{0})^{2}\right)$
$\displaystyle=\beta_{0}\left(p+c_{0}\right),$

and, if $u=v=M$,

$\displaystyle\kappa_{MM}=\sum_{j=1}^{M}\beta_{j}\phi_{Mj}^{2}=\beta_{1}\phi_{M1}^{2}+\beta_{2}\phi_{M2}^{2}+\beta_{3}\phi_{M3}^{2}+...+\beta_{M}\phi_{MM}^{2}$
$\displaystyle=\beta_{0}\phi_{0}^{2}+\beta_{0}\phi_{0}^{2}+\beta_{0}\phi_{0}^{2}+\ldots+\zeta\beta_{0}(1-(M-1)\phi_{0})^{2}$
$\displaystyle=\zeta\beta_{0}(1-(M-1)\phi_{0})^{2}+(M-1)\beta_{0}\phi_{0}^{2}$
$\displaystyle=(M-1)\beta_{0}\phi_{0}^{2}+\zeta\beta_{0}(1-(M-1)\phi_{0})^{2}$
$\displaystyle=\beta_{0}\left(M\phi_{0}^{2}-\phi_{0}^{2}+\zeta(1-2(M-1)\phi_{0}+(M-1)^{2}\phi_{0}^{2}\right)$
$\displaystyle=\beta_{0}\left(M\phi_{0}^{2}-\phi_{0}^{2}+\zeta-2\zeta
M\phi_{0}+2\zeta\phi_{0}+\zeta M^{2}\phi_{0}^{2}-2\zeta
M\phi_{0}^{2}+\zeta\phi_{0}^{2}\right)$
$\displaystyle=\beta_{0}\left(\zeta\phi_{0}^{2}+M\phi_{0}^{2}-\phi_{0}^{2}+\zeta-2\zeta
M\phi_{0}+2\zeta\phi_{0}+\zeta M^{2}\phi_{0}^{2}-2\zeta M\phi_{0}^{2}\right)$
$\displaystyle=\beta_{0}\left(\zeta\phi_{0}^{2}+\zeta-2\zeta(M-1)\phi_{0}+\phi_{0}^{2}\left(M-1-2M\zeta+M^{2}\zeta\right)\right)$
$\displaystyle=\beta_{0}\left(p+c_{1}\right).$

Finally, for any non-diagonal elements of the NGM $\kappa$, when
$u\in{1,2,3,4,5,...,M-1}$, and $v\in{1,2,3,4,5,...,M-1}$, and $u\not=v$ we
have

$\displaystyle\kappa_{uv}=\sum_{j=1}^{M}\beta_{j}\phi_{uj}\phi_{vj}$
$\displaystyle=\beta_{1}\phi_{u1}\phi_{v1}+\beta_{2}\phi_{u2}\phi_{v2}+\ldots+\beta_{u}\phi_{uu}\phi_{vu}+...+\beta_{v}\phi_{uv}\phi_{vu}...+\beta_{M}\phi_{uM}\phi_{vM}$
$\displaystyle=\beta_{0}0(0)+\beta_{0}.0+\ldots+\beta_{0}(1-\phi_{0}).0+\ldots+0+\zeta\beta_{0}(\phi_{0})\phi_{0}$
$\displaystyle=\beta_{0}(\zeta\phi_{0}^{2})$ $\displaystyle=\beta_{0}p,$

and for $u=1,2,3,4,\ldots,M-1$,

$\displaystyle\kappa_{uM}=\sum_{j=1}^{M}\beta_{j}\phi_{uj}\phi_{Mj}$
$\displaystyle=\beta_{1}\phi_{u1}\phi_{M1}+\beta_{2}\phi_{u2}\phi_{M2}+\beta_{3}\phi_{u3}\phi_{M3}+\ldots+\beta_{i}\phi_{uu}\phi_{Mu}+\ldots+\beta_{M}\phi_{uM}\phi_{MM}$
$\displaystyle=\beta_{0}0.\phi_{0}+\beta_{0}0.\phi_{0}+\beta_{0}0.\phi_{0}+\ldots+\beta_{0}(1-\phi_{0})\phi_{0}+...+\zeta\beta_{0}\phi_{0}(1-(M-1)\phi_{0})$
$\displaystyle=\beta_{0}[\phi_{0}-\phi_{0}^{2}+\zeta\left(\phi_{0}-(M-1)\phi_{0}^{2}\right)]$
$\displaystyle=\beta_{0}q,$

and also for $v=1,2,3,4,\ldots,M-1$,

$\begin{split}&\kappa_{Mv}=\sum_{j=1}^{M}\beta_{j}\phi_{Mj}\phi_{vj}\\\
&=\beta_{1}\phi_{M1}\phi_{v1}+\beta_{2}\phi_{M2}\phi_{v2}+\beta_{3}\phi_{M3}\phi_{v3}+\ldots+\beta_{i}\phi_{Mv}\phi_{vv}+\ldots+\beta_{M}\phi_{MM}\phi_{vM}\\\
&=\beta_{0}\phi_{0}.0+\beta_{0}\phi_{0}.0+\beta_{0}\phi_{0}.0+\ldots+\beta_{0}\phi_{0}(1-\phi_{0})+...+\zeta\beta_{0}\phi_{0}(1-(M-1)\phi_{0})\\\
&=\beta_{0}[\phi_{0}-\phi_{0}^{2}+\zeta\left(\phi_{0}-(M-1)\phi_{0}^{2}\right)]\\\
&=\beta_{0}q.\end{split}$

The resulting NGM can be thus written as

$\displaystyle\kappa=\frac{\beta_{0}}{\gamma}\begin{pmatrix}p+c_{0}&p&p&\ldots&q\\\
p&p+c_{0}&p&\ldots&q\\\ p&p&p+c_{0}&\ldots&q\\\
\vdots&\vdots&\vdots&\ddots&\vdots\\\ q&q&q&\ldots&p+c_{1}\\\
\end{pmatrix}=\frac{\beta_{0}}{\gamma}\begin{pmatrix}p&p&p&\ldots&q\\\
p&p&p&\ldots&q\\\ p&p&p&\ldots&q\\\ \vdots&\vdots&\vdots&\ddots&\vdots\\\
q&q&q&\ldots&p\\\
\end{pmatrix}+\frac{\beta_{0}}{\gamma}\begin{pmatrix}c_{0}\\\ c_{0}\\\
c_{0}\\\ \vdots\\\ c_{1}\\\ \end{pmatrix}I$ (16)

where $p=\zeta\phi_{0}^{2}$,
$q=\phi_{0}-\phi_{0}^{2}+\zeta\left(\phi_{0}-(M-1)\phi_{0}^{2}\right)$,
$c_{0}=(1-\phi_{0})^{2}$, and

$c_{1}=\zeta-2\zeta(M-1)\phi_{0}+\phi_{0}^{2}\left(M-1-2M\zeta+M^{2}\zeta\right).$

An alternative representation of the NGM can be done by writing

$\kappa=\frac{\beta_{0}}{\gamma}(\mathcal{P}+\mathcal{D}),$ (17)

where $\mathcal{P}=\frac{\beta_{0}}{\gamma}\begin{pmatrix}p&p&p&\ldots&q\\\
p&p&p&\ldots&q\\\ p&p&p&\ldots&q\\\ \vdots&\vdots&\vdots&\ddots&\vdots\\\
q&q&q&\ldots&p\\\ \end{pmatrix}$ and the diagonal matrix
$\mathcal{D}=\begin{pmatrix}c_{0}\\\ c_{0}\\\ c_{0}\\\ \vdots\\\ c_{1}\\\
\end{pmatrix}I$.

Since the next generation matrix $\kappa$ in (Eq. LABEL:NGM2) is symmetric,
the eigenvalues of $\kappa$ are all real. By doing elementary row or column
operations, we see that the rank of the matrix $\mathcal{P}$ is $2$. So all
submatrices bigger than $2\times 2$ have the determinant zero (see, e. g.
Friedberg (?)). As $\kappa$ in (Eq. LABEL:NGM2) is the sum of a low-rank
matrix and the diagonal matrix $\mathcal{D}$; therefore, we write the
characteristic polynomial of $\mathcal{P}+\mathcal{D}$ in the form

$(c_{0}-\lambda)^{M-2}[\lambda^{2}-(c_{0}+c_{1}+pM)\lambda+\left(c_{1}c_{0}+p(M-1)c_{1}+pc_{0}+(M-1)(p^{2}-q^{2})\right)]=0$

using the expansion (see, Collings (?)) for the eigenvalues $\lambda$. For the
square matrix, $\mathcal{P}+\mathcal{D}$, an eigenvector $V$ and eigenvalue
$\lambda$ satisfy the equation

$\left(\mathcal{P}+\mathcal{D}\right)V=\lambda V,$

and so for any scalar $R_{0}^{\text{hom}}$ we have

$R_{0}^{\text{hom}}\kappa
V=R_{0}^{\text{hom}}\left(\mathcal{P}+\mathcal{D}\right)V=R_{0}^{\text{hom}}\lambda
V.$

Therefore, $R_{0}^{\text{hom}}\lambda$ is the eigenvalue of the next
generation matrix $\kappa$. From the characteristic polynomial, the first
eigenvalue of $\kappa$ is
$\frac{\beta_{0}}{\gamma}c_{0}=R_{0}^{\text{hom}}c_{0}=R_{0}^{\text{hom}}(1-\phi_{0})^{2}$
with multiplicity $M-2$, and the other two eigenvalues are determined from the
quadratic term

$P(\lambda)=\lambda^{2}-(c_{0}+c_{1}+pM)\lambda+\left(c_{1}c_{0}+p(M-1)c_{1}+pc_{0}+(M-1)(p^{2}-q^{2})\right).$
(18)

Furthermore, we find the parameter conditions so that all eigenvalues of
$\kappa$ are within the unit circle. The first set of eigenvalues
automatically satisfies the condition

$R_{0}^{\text{hom}}(1-\phi_{0})^{2}<1,$ (19)

since $R_{0}^{\text{hom}}<1$ by assumption, and $\phi_{0}\in(0,\frac{1}{2}]$.
Now, the quadratic (Eq. 18) gives the other two eigenvalues lying within the
unit circle which can be tested by applying the Jury conditions (Murray (?),
page $507$). For a quadratic equation
$P(\lambda)=\lambda^{2}+a_{1}\lambda+a_{0}=0$, the Jury conditions tell us
$P(1)=1+a_{1}+a_{0}>0$, $P(-1)=1-a_{1}+a_{0}>0$, and $P(0)=a_{0}<1$.

As $R_{0}^{\text{hom}}\lambda:=\lambda^{\prime}$ is the eigenvalue of
$\kappa$, and so substituting
$\lambda=\frac{\lambda^{\prime}}{R_{0}^{\text{hom}}}$ in (Eq. 18) we have

$P(\lambda^{\prime})=\frac{1}{(R_{0}^{\text{hom}})^{2}}[(\lambda^{\prime})^{2}-a_{1}R_{0}^{\text{hom}}\lambda^{\prime}+a_{0}(R_{0}^{\text{hom}})^{2},$
(20)

where $a_{0}:=\left(c_{1}c_{0}+p(M-1)c_{1}+pc_{0}+(M-1)(p^{2}-q^{2})\right)$
and $a_{1}:=(c_{0}+c_{1}+pM)$. By applying the Jury conditions to (Eq. 20) we
obtain the following inequalities:

$\displaystyle\left\\{R_{0}^{\text{hom}}[\zeta\left(R_{0}^{\text{hom}}(M\phi_{0}-1)^{2}-(M-1)\phi_{0}(M\phi_{0}-2)-1\right)+\phi_{0}(2-M\phi_{0})-1]+1\right\\}>0,$
(21)
$\displaystyle\left\\{R_{0}^{\text{hom}}[\zeta\left(R_{0}^{\text{hom}}(M\phi_{0}-1)^{2}+(M-1)\phi_{0}(M\phi_{0}-2)+1\right)-\phi_{0}(2-M\phi_{0})+1]+1\right\\}>0,$
(22)

$\displaystyle(R_{0}^{\text{hom}})^{2}\zeta(M\phi_{0}-1)^{2}<1.$ (23)

The Jury conditions imply that (Eqs. 21-23) are identical to (Eqs. A.2-A.4) in
Stolerman et al. (?). Therefore, the layout of the proof follows the same
steps. ∎

Theorem 1 is not only mathematically remarkable, but it opens a new question:
_Is there a class of networks where the same critical $\zeta$ value and
interval $\mathcal{I}$ always emerge?_ A common feature of fully connected and
star-shaped networks is the connection that all nodes have with the node with
a different infection rate (the central node in the star). It is thus natural
to explore networks without this feature and observe if the results obtained
in Theorem 1 (and Theorem 2 in (?)) still hold. We start the next section with
the so-called _cycle-shaped_ network, perhaps the most natural example of a
network with no such highly connected nodes.

###  4.2 Cycle networks

In the cycle networks, the nodes only connect to their immediate neighbors.
These networks are a simple model for small cities or towns where people
usually commute short distances. For _unidirectional_ cycle networks, people
only travel from node $i$ to node $i+1$. This network can be represented by
the following flux matrix:

$\phi_{M\times M}=\begin{pmatrix}1-\phi_{0}&\phi_{0}&0&0&\ldots&0\\\
0&1-\phi_{0}&\phi_{0}&0&\ldots&0\\\ 0&0&1-\phi_{0}&\phi_{0}&\ldots&0\\\
\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\\
\phi_{0}&0&0&0&\ldots&1-\phi_{0}\end{pmatrix}.$ (24)

where $0<\phi_{0}<1$. On the other hand, _bidirectional_ cycle networks can be
summarized as follows by the geometry of the flux matrix:

$\phi_{M\times M}=\begin{pmatrix}1-2\phi_{0}&\phi_{0}&0&\ldots&0&\phi_{0}\\\
\phi_{0}&1-2\phi_{0}&\phi_{0}&\ldots&0&0\\\
0&\phi_{0}&1-2\phi_{0}&\ldots&0&0\\\
\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\\
\phi_{0}&0&0&\ldots&\phi_{0}&1-2\phi_{0}\end{pmatrix}.$ (25)

where $0<\phi_{0}\leq\frac{1}{2}$.

We start by exploring a particular case of a unidirectional cycle with $M=3$,
which indicates that unidirectional cycles are not in the class of networks
where the critical value $\zeta$ (Eq. Theorem) and interval $\mathcal{I}$ (Eq.
6) emerge. The proof of the following propositions can be found in the
appendix.

######  Proposition 2.

Let $V=\left\\{1,2,3\right\\}$ be the node of a _unidirectional cycle_ network
and the flux matrix $\phi_{3\times 3}$ be defined by (Eq. 24). We consider the
infection rates $\beta_{i}$ from (Eq. 8), where all the nodes are stable
except the third one, i.e, $R_{0}^{\text{hom}}<1$ and
$\frac{\zeta\beta_{0}}{\gamma}>1$. Then we have a critical value
$\zeta^{\text{crit}}$ given by

$\displaystyle\zeta^{\text{crit}}=\frac{1}{R_{0}^{\text{hom}}}\left(\frac{3R_{0}^{\text{hom}}-4}{R_{0}^{\text{hom}}-2}\right)$
(26)

such that if $\frac{1}{R_{0}^{\text{hom}}}<\zeta<\zeta^{\text{crit}}$, then
there exists an interval $\mathcal{I}$ around the point $\frac{1}{2}$ given by

$\mathcal{I}=\left(\max\left\\{0,\frac{1}{2}-\tilde{\phi}\right\\},\min\left\\{1,\frac{1}{2}+\tilde{\phi}\right\\}\right),$
(27)

where

$\displaystyle\tilde{\phi}=\frac{1}{2}\sqrt{\frac{\left(\zeta
R_{0}^{\text{hom}}(R_{0}^{\text{hom}}-2)+(4-3R_{0}^{\text{hom}})\right)}{R_{0}^{\text{hom}}\left(1+2\zeta-3\zeta
R_{0}^{\text{hom}}\right)}}$ (28)

such that if $\phi_{0}\in\mathcal{I}$, the DFE is stable, and otherwise, it is
unstable.

The $\zeta^{\text{crit}}$ found for _unidirectional cycle_ networks (Eq. 26)
are different from the _star-shaped_ networks (Eq. 10). The interval where
there will be no epidemic is given by (Eq. 27) in the _unidirectional cycle_
networks, which is different from the interval obtained for _star-shaped_
networks (Eq. 11). We conjecture that this difference also holds for $M>3$. In
section 5, we run numerical simulations for different network sizes and
observe that our conjecture is likely true.

If we connect all nodes of a cycle network to the last node of index $M$, we
can build different networks and further analyze the corresponding epidemic
thresholds. In the following section, we combine the bidirectional cycle-
shaped networks with the star-shaped networks.

###  4.3 Extension of heterogeneous cycle networks

This section aims to modify bidirectional cycle networks and investigate where
we find the minimum $\zeta$ (Eq. Theorem) and the interval $\mathcal{I}$ (Eq.
6). Specifically, we modify the bidirectional cycle by joining the
heterogeneous node $M$ directly to the other nodes across the cycle. In this
case, the adjacent nodes to the heterogeneous node are disjoint. We call this
network _cycle support_. We then extend these simple heuristic networks to the
_star-cycle_ networks by joining the heterogeneous node’s immediate nodes $1$
and $M-1$.

####  4.3.1 Cycle support networks

We introduce new networks called _cycle support_ , formed by the bidirectional
cycle networks and joining the heterogeneous node $M$ directly to the nodes
across the cycle. Their general flux matrix is given by:

$\phi_{M\times M}=\begin{pmatrix}1-2\phi_{0}&\phi_{0}&0&0&\ldots&0&\phi_{0}\\\
\phi_{0}&1-3\phi_{0}&\phi_{0}&0&\ldots&0&\phi_{0}\\\
0&\phi_{0}&1-3\phi_{0}&\phi_{0}&\ldots&0&\phi_{0}\\\
\vdots&\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\\
0&0&0&\ldots&\phi_{0}&1-2\phi_{0}&\phi_{0}\\\
\phi_{0}&\phi_{0}&\phi_{0}&\ldots&\phi_{0}&\phi_{0}&1-(M-1)\phi_{0}\end{pmatrix},$
(29)

where $0<\phi_{0}\leq\frac{1}{M-1}$, since $\phi_{ij}\geq 0$. We then obtain
the following proposition for _cycle support_ networks with $5$ nodes (the
proof can be found in the Appendix).

######  Proposition 3.

Let $V=\left\\{1,2,3,4,5\right\\}$ be a node of _cycle-supports_ network and
let the geometry of its flux matrix be defined by (Eq. 29), where
$\phi_{0}\in(0,\frac{1}{4}]$. The infection rate $\beta_{i}$ is defined in
(Eq. 8), and we also assume that all the outer nodes are stable
($R_{0}^{\text{hom}}=\frac{\beta_{0}}{\gamma}<1$) except the center node
$M=5$. Then the critical value of $\zeta$ called $\zeta^{\text{crit}}$ is
given by

$\zeta^{\text{crit}}=1+\left(\frac{1-R_{0}^{\text{hom}}}{R_{0}^{\text{hom}}}\right)5.$

Alternatively, the infection threshold of $\beta_{0}$,
$\beta_{0}^{\text{crit}}$, is given by

$\beta_{0}^{\text{crit}}=\gamma\left(\frac{5}{4+\zeta}\right).$ (30)

If $\frac{1}{R_{0}^{\text{hom}}}<\zeta<\zeta^{\text{crit}}$, then there is an
interval around the point $\frac{1}{5}$ given by
$\mathcal{I}=\left(\max\left\\{0,\frac{1}{5}-\tilde{\phi}\right\\},\min\left\\{\frac{1}{4},\frac{1}{5}+\tilde{\phi}\right\\}\right)$,
where

$\displaystyle\tilde{\phi}=\frac{1}{5}\sqrt{\frac{\gamma(5\gamma-4\beta_{0}-\zeta\beta_{0})}{\beta_{0}(\gamma+4\zeta\gamma-5\zeta\beta_{0})}}=\frac{1}{5}\sqrt{\frac{5(1-R_{0}^{\text{hom}})-R_{0}^{\text{hom}}(\zeta-1))}{R_{0}^{\text{hom}}\left(5\zeta(1-R_{0}^{\text{hom}})-(\zeta-1)\right)}},$

such that if $\phi_{0}\in\mathcal{I}$, the DFE is stable; and otherwise it is
unstable.

######  Remark 2.

We observe that if $\zeta<\frac{1}{R_{0}^{\text{hom}}}$, the DFE is always
stable. Note that we have a similar proof to $\\{1,\ldots,5\\}$, 3 for the
nodes $M=4,5,6,7$. On the other hand, simulations can be performed for every
node $M$.

####  4.3.2 Star-cycle networks

For star-cycle networks, we start from the heterogeneous node $M$ as the
center node of the network. Then the center node $M$ is connected to all other
outer nodes to create a _star-shaped_ network. Finally, the rest of the nodes
are in the form of bidirectional cycle networks called _star-cycle_ networks.
The geometry of the flux matrix is given by:

$\phi_{M\times
M}=\begin{pmatrix}1-3\phi_{0}&\phi_{0}&0&0&\ldots&\phi_{0}&\phi_{0}\\\
\phi_{0}&1-3\phi_{0}&\phi_{0}&0&\ldots&0&\phi_{0}\\\
0&\phi_{0}&1-3\phi_{0}&\phi_{0}&\ldots&0&\phi_{0}\\\
\vdots&\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\\
\phi_{0}&0&0&\ldots&\phi_{0}&1-3\phi_{0}&\phi_{0}\\\
\phi_{0}&\phi_{0}&\phi_{0}&\ldots&\phi_{0}&\phi_{0}&1-(M-1)\phi_{0}\end{pmatrix},$
(31)

where $0<\phi_{0}\leq\frac{1}{M-1}$, since $\phi_{ij}\geq 0$. Note that the
_cycle-support_ networks described in Section 4.3.1 are _star-cycle_ networks,
where the adjacent node of $M$ does not join with $M$. We then obtain the
following proposition for _star-cycle_ networks with $5$ nodes (the proof can
be found in the Appendix).

######  Proposition 4.

Let $V=\left\\{1,2,3,4,5\right\\}$ be the nodes of a _star-cycle_ network, and
the geometry of the corresponding network is defined by the flux (Eq. 31),
where $M=5,\leavevmode\nobreak\ \text{and}\leavevmode\nobreak\
\phi_{0}\in(0,\frac{1}{4}]$. The infection rate $\beta_{i}$ is defined in (Eq.
8), and we also assume that all the nodes are stable
($R_{0}^{\text{hom}}=\frac{\beta_{0}}{\gamma}<1$) except the node $M=5$. Then
the critical value of $\zeta$ called $\zeta^{\text{crit}}$ is given by

$\zeta^{\text{crit}}=1+\left(\frac{1-R_{0}^{\text{hom}}}{R_{0}^{\text{hom}}}\right)5.$
(32)

Alternatively, the infection threshold of $\beta_{0}$,
$\beta_{0}^{\text{crit}}$ is given by

$\beta_{0}^{\text{crit}}=\gamma\left(\frac{5}{4+\zeta}\right).$ (33)

If $\frac{1}{R_{0}^{\text{hom}}}<\zeta<\zeta^{\text{crit}}$, there is an
interval around the point $\frac{1}{5}$ given by
$\mathcal{I}=\left(\max\left\\{0,\frac{1}{5}-\tilde{\phi}\right\\},\min\left\\{\frac{1}{4},\frac{1}{5}+\tilde{\phi}\right\\}\right)$,
where

$\displaystyle\tilde{\phi}=\frac{1}{5}\sqrt{\frac{\gamma(5\gamma-4\beta_{0}-\zeta\beta_{0})}{\beta_{0}(\gamma+4\zeta\gamma-5\zeta\beta_{0})}}=\frac{1}{5}\sqrt{\frac{5(1-R_{0}^{\text{hom}})-R_{0}^{\text{hom}}(\zeta-1))}{R_{0}^{\text{hom}}\left(5\zeta(1-R_{0}^{\text{hom}})-(\zeta-1)\right)}},$

such that if $\phi_{0}\in\mathcal{I}$, the DFE is stable, and otherwise it is
unstable.

We observe that, if $\zeta<\frac{1}{R_{0}^{\text{hom}}}$, DFE is always
stable. Using symbolic calculations with the software Mathematica, we can also
perform a similar analysis for $M=4,5,6,7,9,11,13$. In the following
conjecture, we expand the Proposition 4 to an arbitrary size of _star-cycle_
networks.

######  Conjecture 1.

Let $V=\left\\{1,2,3,\ldots,M\right\\}$ be the nodes of a _star-cycle_ network
defined by the symmetric flux matrix (Eq. 31), where
$0<\phi_{0}\leq\frac{1}{M-1}$. We assume that the infection rate
$\beta_{i}=\beta_{0}$ at node $i$ for all $i=1,2,3,\cdots,M-1$ and
$\beta_{M}=\zeta\beta_{0}$. We also suppose that all the nodes are stable
$(R_{0}^{\text{hom}}=\frac{\beta_{0}}{\gamma}<1)$ except the center $M$. Then
we have a minimum $\zeta^{\text{crit}}$ given by (Eq. Theorem or 15). If
$\frac{1}{R_{0}^{\text{hom}}}<\zeta<\zeta^{\text{crit}}$, there exists an
interval
$\mathcal{I}=\left(\max\left\\{0,\frac{1}{M}-\tilde{\phi}\right\\},\min\left\\{\frac{1}{M-1},\frac{1}{M}+\tilde{\phi}\right\\}\right)$,
where $\tilde{\phi}$ is given by (Eq. 7), such that if
$\phi_{0}\in\mathcal{I}$, then the DFE is stable. Otherwise, it is unstable.
Finally, if $\zeta<\frac{1}{R_{0}^{\text{hom}}}$ i.e. $\beta_{M}<\gamma$, then
DFE is always stable.

###  4.4 More networks based on the star-shaped geometry

In the previous sections, we studied networks with different structures, but
assuming an equal flux $\phi_{0}$ between any two connected nodes.
Realistically, the fluxes between nodes could be different. In this section,
we explore more general networks with some of the structure of star-shaped
networks, but assuming the existence of two different flux values $\phi_{0}$
and $\phi_{1}$.

####  4.4.1 Star-triangle networks

We consider a _star-triangle_ network with an odd number of nodes $M$. The
nonnegative entries of its flux matrix $\phi_{M\times M}$ are defined by

$\phi_{ij}=\begin{cases}\phi_{0}&\text{if $i=M,\leavevmode\nobreak\
\text{and}\leavevmode\nobreak\ j\in V\setminus M$}\leavevmode\nobreak\ \text{
or \leavevmode\nobreak\ if $j=M,\leavevmode\nobreak\
\text{and}\leavevmode\nobreak\ i\in V\setminus M$}\\\
1-\phi_{0}-\phi_{1}&\text{if $i=j,\leavevmode\nobreak\
\text{and}\leavevmode\nobreak\ \leavevmode\nobreak\ i,j\in V\setminus M$}\\\
1-(M-1)\phi_{0}&\text{if $i=j=M$}\\\ \phi_{1}&\text{for all $i\in
V\setminus\\{M,M-1\\}\leavevmode\nobreak\ \&\leavevmode\nobreak\ j=i+1$}\\\
\phi_{1}&\text{for all $j\in V\setminus\\{M,M-1\\}\leavevmode\nobreak\
\&\leavevmode\nobreak\ i=j+1$}\\\ 0&\text{otherwise},\end{cases}$ (34)

where $\phi_{0}+\phi_{1}<1,0<\phi_{0}\leq\frac{1}{M-1}$, and $\phi_{1}>0$.
Alternatively, the explicit form of the flux matrix (Eq. 9) is given by

$\phi_{M\times
M}=\begin{pmatrix}1-\phi_{0}-\phi_{1}&\phi_{1}&\ldots&0&0&\phi_{0}\\\
\phi_{1}&1-\phi_{0}-\phi_{1}&\ldots&0&0&\phi_{0}\\\
\vdots&\vdots&\ddots&\vdots&\vdots&\vdots\\\
0&0&\ldots&1-\phi_{0}-\phi_{1}&\phi_{1}&\phi_{0}\\\
0&0&\ldots&\phi_{1}&1-\phi_{0}-\phi_{1}&\phi_{0}\\\
\phi_{0}&\phi_{0}&\ldots&\phi_{0}&\phi_{0}&1-(M-1)\phi_{0}\end{pmatrix}.$

In the following theorem, we find the exact minimum $\zeta$ and interval
$\mathcal{I}$ for _star-triangle_ networks as we do for _fully connected_
networks, where the flux between nodes can prevent an epidemic outbreak if it
lies in a specific interval. The proof can be found in the appendix.

######  Theorem 2.

Let $V=\left\\{1,2,3,\ldots,M\right\\}$ be the nodes of a star-triangle
network defined by the flux matrix (34), where
$\phi_{0}+\phi_{1}<1,0<\phi_{0}\leq\frac{1}{M-1}$, and $\phi_{1}>0$. We assume
that the infection rate $\beta_{i}=\beta_{0}$ at node $i$ for all
$i=1,2,3,\cdots,M-1$ and $\beta_{M}=\zeta\beta_{0}$. We also suppose that all
the nodes are stable $(R_{0}^{\text{hom}}=\frac{\beta_{0}}{\gamma}<1)$ except
the center node $M$. Then we have a minimum $\zeta^{\text{crit}}$ given by
(Eq. Theorem or 15). If
$\frac{1}{R_{0}^{\text{hom}}}<\zeta<\zeta^{\text{crit}}$, there exists an
interval
$\mathcal{I}=\left(\max\left\\{0,\frac{1}{M}-\tilde{\phi}\right\\},\min\left\\{\frac{1}{M-1},\frac{1}{M}+\tilde{\phi}\right\\}\right)$,
where $\tilde{\phi}$ is given by (Eq. 7), such that if
$\phi_{0}\in\mathcal{I}$, then the DFE is stable, and otherwise it is
unstable.

In the previous theorem, we note that if $\zeta<\frac{1}{R_{0}^{\text{hom}}}$,
i.e., $\beta_{M}<\gamma$, then DFE is always stable. Let us now generalize the
flux matrix (Eq. 34) accounting for the other nodes with different fluxes in
the following way:

$\phi_{M\times M}=\begin{cases}\phi_{2i-1,2i}=\phi_{i}&\text{ $1\leq
i\leq(\frac{M-1}{2})$}\\\ \phi_{2i,2i-1}=\phi_{i}&\text{ $1\leq
i\leq(\frac{M-1}{2})$}\\\ \phi_{M,i}=\phi_{(\frac{M+1}{2})}&\text{ $1\leq
i\leq M-1$}\\\ \phi_{i,M}=\phi_{(\frac{M+1}{2})}&\text{ $1\leq i\leq M-1$}\\\
\phi_{i,i}=1-\sum_{j=1}^{M}\phi_{i,j}&\text{ $1\leq i\leq M-1$
\leavevmode\nobreak\ and \leavevmode\nobreak\ $i\neq j$}\\\
\phi_{M,M}=1-(M-1)\phi_{(\frac{M+1}{2})}\\\ 0&\text{otherwise},\end{cases}$
(35)

where for all $\phi_{i,j}>0,i,j\in V$, and $M$ is an odd number. The fraction
of movement from the center node $M$ to the other nodes is defined as
$\phi_{(\frac{M+1}{2})}=\phi_{0}$ to be consistent with the other networks. We
prove the following result for this more general _star-triangle_ network for
the network size $M=11$ (the proof can be found in the Appendix).

######  Proposition 5.

Let $V=\left\\{1,2,3,4,\ldots,11\right\\}$ be the nodes of a star-triangle
network defined by the symmetric flux matrix (Eq. 35) and $M=11$ be an odd
number. We assume that the infection rate $\beta_{i}=\beta_{0}$ at node $i$
for all $i=1,2,3,\cdots,10$ and $\beta_{11}=\zeta\beta_{0}$. We also suppose
that all the nodes are stable
$(R_{0}^{\text{hom}}=\frac{\beta_{0}}{\gamma}<1)$ except the center node $M$.
Then the critical value of $\zeta$ called $\zeta^{\text{crit}}$ is given by

$\zeta^{\text{crit}}=1+\left(\frac{1-R_{0}^{\text{hom}}}{R_{0}^{\text{hom}}}\right)11.$
(36)

Alternatively, the infection threshold of $\beta_{0}$,
$\beta_{0}^{\text{crit}}$ is given by

$\beta_{0}^{\text{crit}}=\gamma\left(\frac{11}{10+\zeta}\right).$ (37)

If $\frac{1}{R_{0}^{\text{hom}}}<\zeta<\zeta^{\text{crit}}$, there exists an
interval
$\mathcal{I}=\left(\max\left\\{0,\frac{1}{11}-\tilde{\phi}\right\\},\min\left\\{\frac{1}{10},\frac{1}{11}+\tilde{\phi}\right\\}\right)$,
where $\tilde{\phi}$ is given by (Eq. 7), such that if
$\phi_{0}\in\mathcal{I}$, then the DFE is stable, and otherwise it is
unstable.

Also if $\zeta<\frac{1}{R_{0}^{\text{hom}}}$, i.e., $\beta_{M}<\gamma$, in
Proposition 5, then the DFE is always stable. Proposition 5 also holds for odd
$M$ values $M=3,5,7,9$ and higher. Our results for particular network sizes
can be thus generalized in the following conjecture.

######  Conjecture 2.

Let $V=\left\\{1,2,3,4,\ldots,M\right\\}$ be the nodes of a _star-triangle_
network defined by the symmetric flux matrix (Eq. 35), and let $M$ be an odd
number. We assume that the infection rate $\beta_{i}=\beta_{0}$ at node $i$
for all $i=1,2,3,\cdots,M-1$ and $\beta_{M}=\zeta\beta_{0}$. We also suppose
that all the nodes are stable
($R_{0}^{\text{hom}}=\frac{\beta_{0}}{\gamma}<1$) except the center node $M$.
Then we have a minimum $\zeta^{\text{crit}}$ given by (Eq. Theorem or 15). If
$\frac{1}{R_{0}^{\text{hom}}}<\zeta<\zeta^{\text{crit}}$, there exists an
interval
$\mathcal{I}=\left(\max\left\\{0,\frac{1}{M}-\tilde{\phi}\right\\},\min\left\\{\frac{1}{M-1},\frac{1}{M}+\tilde{\phi}\right\\}\right)$,
where $\tilde{\phi}$ is given by (Eq. 7), such that if
$\phi_{0}\in\mathcal{I}$, then the DFE is stable, otherwise, it is unstable.
Finally, if $\zeta<\frac{1}{R_{0}^{\text{hom}}}$, i.e., $\beta_{M}<\gamma$,
then the DFE is always stable.

####  4.4.2 Star-background networks

In Section 4.1, we studied _star-shaped_ networks and introduced theoretical
estimates for stability conditions for networks of arbitrary size $M$. In
Section 3, we reviewed a previous result for _fully connected_ networks
(Theorem Theorem). In this section, we define _star-background_ networks, from
which both fully connected and star-shaped networks can be obtained.
Specifically, we consider a heterogeneous node $M$ (the one with a different
infection rate) connecting to all other nodes with the same flux $\phi_{0}$,
and all the outer nodes are connected to all other nodes by a different flux
($\phi_{1}$). The corresponding flux matrix is given by

$\phi_{ij}=\begin{cases}\phi_{0}&\text{if $i=M,\leavevmode\nobreak\
\text{and}\leavevmode\nobreak\ j\in V$}\leavevmode\nobreak\ \text{ or
\leavevmode\nobreak\ if $j=M,\leavevmode\nobreak\
\text{and}\leavevmode\nobreak\ i\in V$}\\\ 1-\phi_{0}-(M-2)\phi_{1}&\text{if
$i=j,\leavevmode\nobreak\ \text{and}\leavevmode\nobreak\ \leavevmode\nobreak\
i,j\in V\setminus M$}\\\ 1-(M-1)\phi_{0}&\text{if $i=j=M$}\\\
\phi_{1}&\text{otherwise},\end{cases}$ (38)

where $\phi_{0}+(M-2)\phi_{1}<1,0<\phi_{0}\leq\frac{1}{M-1}$, and
$\phi_{1}>0$. The explicit form of the flux matrix (Eq. 38) is the following:

$\phi_{M\times
M}=\begin{pmatrix}1-\phi_{0}-(M-2)\phi_{1}&\phi_{1}&\phi_{1}&\ldots&\phi_{0}\\\
\phi_{1}&1-\phi_{0}-(M-2)\phi_{1}&\phi_{1}&\ldots&\phi_{0}\\\
\phi_{1}&\phi_{1}&1-\phi_{0}-(M-2)\phi_{1}&\ldots&\phi_{0}\\\
\vdots&\vdots&\vdots&\ddots&\vdots\\\
\phi_{0}&\phi_{0}&\phi_{0}&\ldots&1-(M-1)\phi_{0}\end{pmatrix}.$ (39)

The above flux matrix somehow connects both _fully connected_ and _star-
shaped_ networks. In fact, when $\phi_{1}=0$ we obtain _star-shaped_ networks.
On the other hand, when $\phi_{1}=\phi_{0}$, we obtain the _fully connected_
networks. _Star-background_ networks are thus an interesting “bridge” between
fully connected and _star-shaped_ networks. The following theorem presents the
epidemic thresholds for this class of networks (see appendix for the proof).

######  Theorem 3.

Let $V=\left\\{1,2,3,\ldots,M\right\\}$ be the nodes of a _star-background_
networks defined by the flux matrix (Eq. 38), where
$\phi_{0}+(M-2)\phi_{1}<1,0<\phi_{0}\leq\frac{1}{M-1}$, and $\phi_{1}>0$. We
assume that the infection rate $\beta_{i}=\beta_{0}$ at node $i$ for all
$i=1,2,3,\cdots,M-1$ and $\beta_{M}=\zeta\beta_{0}$. We also suppose that all
the nodes are stable ($R_{0}^{\text{hom}}=\frac{\beta_{0}}{\gamma}<1$) except
the center $M$. Then we have a minimum $\zeta^{\text{crit}}$ given by (Eq.
Theorem or 15). If $\frac{1}{R_{0}^{\text{hom}}}<\zeta<\zeta^{\text{crit}}$,
there exists an interval
$\mathcal{I}=\left(\max\left\\{0,\frac{1}{M}-\tilde{\phi}\right\\},\min\left\\{\frac{1}{M-1},\frac{1}{M}+\tilde{\phi}\right\\}\right)$,
where $\tilde{\phi}$ is given by (Eq. 7), such that if
$\phi_{0}\in\mathcal{I}$, then the DFE is stable, and otherwise it is
unstable. Finally, if $\zeta<\frac{1}{R_{0}^{\text{hom}}}$, i.e.,
$\beta_{M}<\gamma$, then the DFE is always stable.

######  Remark 3.

In particular, if $\phi_{1}=\phi_{0}$ in (Eq. 38), then the flux creates the
_fully connected_ networks described in Theorem Theorem, and the results of
Theorem 3 are identical to Theorem Theorem. Furthermore, $\phi_{1}=0$ in (Eq.
38) the network becomes _star-shaped_ (Fig. 1) network revealed in Theorem 1.
Then the results of the Theorem 3 are identical with Theorem 1.

####  4.4.3 Conjectures on more flexible network structures

To extend the previous section, we have constructed a general network diagram
representing a _star-class_ network. The network construction begins by
connecting a center node $M$ to outer nodes with flux value $\phi_{0}$ (by the
solid red arrows in Fig. 2). Then all the outer or external neighboring nodes
are joined (solid green arrows) with the new flux parameter $\phi_{1}$. Next,
we connect any other outer nodes with a gap of one node between them by the
flux $\phi_{2}$ (solid blue line). Then we connect all the external nodes with
a spacing of two nodes between them with the new flux parameter $\phi_{3}$
(maroon dashed line) and so on. Finally, the general flux of two
interconnecting outer nodes is $\phi_{p+1}$, where p is the number of nodes
between them.

Figure 2: A more flexible _star-class_ network of $M$ nodes: a center node is
connected to the outer nodes with a flux $\phi_{0}$. The outer notes with a
gap of $i$ nodes are connected with a different flux $\phi_{i}$.

The flux matrix for a network with an odd number $M$ of nodes can be written
as

$\phi_{M\times
M}=\setcounter{MaxMatrixCols}{15}\begin{pmatrix}\phi_{d}&\phi_{1}&\phi_{2}&\phi_{3}&\ldots&\phi_{\frac{M-1}{2}-1}&\phi_{\frac{M-1}{2}}&\phi_{\frac{M-1}{2}-1}&\phi_{\frac{M-1}{2}-2}&\ldots&\phi_{4}&\phi_{3}&\phi_{2}&\phi_{1}&\phi_{0}\\\
\phi_{1}&\phi_{d}&\phi_{1}&\phi_{2}&\ldots&\phi_{\frac{M-1}{2}-2}&\phi_{\frac{M-1}{2}-1}&\phi_{\frac{M-1}{2}}&\phi_{\frac{M-1}{2}-1}&\ldots&\phi_{5}&\phi_{4}&\phi_{3}&\phi_{2}&\phi_{0}\\\
\phi_{2}&\phi_{1}&\phi_{d}&\phi_{1}&\ldots&\phi_{\frac{M-1}{2}-3}&\phi_{\frac{M-1}{2}-2}&\phi_{\frac{M-1}{2}-1}&\phi_{\frac{M-1}{2}}&\ldots&\phi_{6}&\phi_{5}&\phi_{4}&\phi_{3}&\phi_{0}\\\
\phi_{3}&\phi_{2}&\phi_{1}&\phi_{d}&\ldots&\phi_{\frac{M-1}{2}-4}&\phi_{\frac{M-1}{2}-3}&\phi_{\frac{M-1}{2}-2}&\phi_{\frac{M-1}{2}-1}&\ldots&\phi_{7}&\phi_{6}&\phi_{5}&\phi_{4}&\phi_{0}\\\
\vdots&\vdots&\vdots&\vdots&\ddots&\vdots&\vdots&\vdots&\vdots&\ddots&\vdots&\vdots&\vdots&\vdots&\vdots\\\
\phi_{\frac{M-1}{2}-1}&\phi_{\frac{M-1}{2}-2}&\phi_{\frac{M-1}{2}-3}&\phi_{\frac{M-1}{2}-4}&\ldots&\phi_{d}&\phi_{1}&\phi_{2}&\phi_{3}&\ldots&\phi_{\frac{M-1}{2}-3}&\phi_{\frac{M-1}{2}-2}&\phi_{\frac{M-1}{2}-1}&\phi_{\frac{M-1}{2}}&\phi_{0}\\\
\phi_{\frac{M-1}{2}}&\phi_{\frac{M-1}{2}-1}&\phi_{\frac{M-1}{2}-2}&\phi_{\frac{M-1}{2}-3}&\ldots&\phi_{1}&\phi_{d}&\phi_{1}&\phi_{2}&\ldots&\phi_{\frac{M-1}{2}-4}&\phi_{\frac{M-1}{2}-3}&\phi_{\frac{M-1}{2}-2}&\phi_{\frac{M-1}{2}-1}&\phi_{0}\\\
\phi_{\frac{M-1}{2}-1}&\phi_{\frac{M-1}{2}}&\phi_{\frac{M-1}{2}-1}&\phi_{\frac{M-1}{2}-2}&\ldots&\phi_{2}&\phi_{1}&\phi_{d}&\phi_{1}&\ldots&\phi_{\frac{M-1}{2}-5}&\phi_{\frac{M-1}{2}-4}&\phi_{\frac{M-1}{2}-3}&\phi_{\frac{M-1}{2}-2}&\phi_{0}\\\
\phi_{\frac{M-1}{2}-2}&\phi_{\frac{M-1}{2}-1}&\phi_{\frac{M-1}{2}}&\phi_{\frac{M-1}{2}-1}&\ldots&\phi_{3}&\phi_{2}&\phi_{1}&\phi_{d}&\ldots&\phi_{\frac{M-1}{2}-6}&\phi_{\frac{M-1}{2}-5}&\phi_{\frac{M-1}{2}-4}&\phi_{\frac{M-1}{2}-3}&\phi_{0}\\\
\vdots&\vdots&\vdots&\vdots&\ddots&\vdots&\vdots&\vdots&\vdots&\ddots&\vdots&\vdots&\vdots&\vdots&\vdots\\\
\phi_{4}&\phi_{5}&\phi_{6}&\phi_{7}&\ldots&\phi_{\frac{M-1}{2}-3}&\phi_{\frac{M-1}{2}-4}&\phi_{\frac{M-1}{2}-5}&\phi_{\frac{M-1}{2}-6}&\ldots&\phi_{d}&\phi_{1}&\phi_{2}&\phi_{3}&\phi_{0}\\\
\phi_{3}&\phi_{4}&\phi_{5}&\phi_{6}&\ldots&\phi_{\frac{M-1}{2}-2}&\phi_{\frac{M-1}{2}-3}&\phi_{\frac{M-1}{2}-4}&\phi_{\frac{M-1}{2}-5}&\ldots&\phi_{1}&\phi_{d}&\phi_{1}&\phi_{2}&\phi_{0}\\\
\phi_{2}&\phi_{3}&\phi_{4}&\phi_{5}&\ldots&\phi_{\frac{M-1}{2}-1}&\phi_{\frac{M-1}{2}-2}&\phi_{\frac{M-1}{2}-3}&\phi_{\frac{M-1}{2}-4}&\ldots&\phi_{2}&\phi_{1}&\phi_{d}&\phi_{1}&\phi_{0}\\\
\phi_{1}&\phi_{2}&\phi_{3}&\phi_{4}&\ldots&\phi_{\frac{M-1}{2}}&\phi_{\frac{M-1}{2}-1}&\phi_{\frac{M-1}{2}-2}&\phi_{\frac{M-1}{2}-3}&\ldots&\phi_{3}&\phi_{2}&\phi_{1}&\phi_{d}&\phi_{0}\\\
\phi_{0}&\phi_{0}&\phi_{0}&\phi_{0}&\ldots&\phi_{0}&\phi_{0}&\phi_{0}&\phi_{0}&\ldots&\phi_{0}&\phi_{0}&\phi_{0}&\phi_{0}&\phi_{MM}\end{pmatrix}$
(40)

where
$\phi_{d}:=1-2(\phi_{0}+\phi_{1}+\ldots+\phi_{\frac{M-1}{2}-1})-\phi_{\frac{M-1}{2}}-\phi_{0}$
and $\phi_{MM}=1-(M-1)\phi_{0}$. The interested reader can find a similar
matrix for a network of an even number of nodes $M$ in the appendix.

We conjecture that the _star-class_ network has the same epidemic threshold as
the fully connected networks. The following proposition can be proved for
$M=9$ (the proof can be found in the appendix).

######  Proposition 6.

Let $V=\left\\{1,2,\ldots,9\right\\}$ be the nodes of a _general networks_
defined by the symmetric flux matrix (Eq. 40). We assume that the infection
rate $\beta_{i}=\beta_{0}$ at node $i$ for all $i=1,2,3,\ldots,9$ and
$\beta_{9}=\zeta\beta_{0}$. We also suppose that all the nodes are stable i.e.
$R_{0}^{\text{hom}}=\frac{\beta_{0}}{\gamma}<1$ except the center node $M$.
Then the critical value of $\zeta$ called $\zeta^{\text{crit}}$ is given by

$\zeta^{\text{crit}}=1+\left(\frac{1-R_{0}^{\text{hom}}}{R_{0}^{\text{hom}}}\right)9.$

Alternatively, the infection threshold of $\beta_{0}$,
$\beta_{0}^{\text{crit}}$ is given by

$\beta_{0}^{\text{crit}}=\gamma\left(\frac{9}{8+\zeta}\right).$

If $\frac{1}{R_{0}^{\text{hom}}}<\zeta<\zeta^{\text{crit}}$, and taking into
account the restriction $\phi_{ij}$, where $i,j=1,2,3,\ldots,9$, then the
interval where there will be no epidemic is given by
$\mathcal{I}(\zeta,R_{0}^{\text{hom}})=\left(\max\left\\{0,\frac{1}{9}-\tilde{\phi}\right\\},\min\left\\{\frac{1}{8},\frac{1}{9}+\tilde{\phi}\right\\}\right)$,
where $\tilde{\phi}$ is given by

$\displaystyle\tilde{\phi}=\frac{1}{9}\sqrt{\frac{9(1-R_{0}^{\text{hom}})-R_{0}^{\text{hom}}(\zeta-1))}{R_{0}^{\text{hom}}\left(9\zeta(1-R_{0}^{\text{hom}})-(\zeta-1)\right)}}$

such that if $\phi_{6}\in\mathcal{I}$, then the DFE is stable, and, otherwise,
it is unstable.

We present a conjecture which contains the findings in Proposition 6 for a
general network with $M$ nodes, where $M$ is an odd number.

######  Conjecture 3.

Let $V=\left\\{1,2,3,\ldots,M\right\\}$ be the nodes of a network with
different fluxes defined by the symmetric flux matrix (Eq. 40), where
$\phi_{ij}\geq 0$. We assume that the infection rate $\beta_{i}=\beta_{0}$ at
node $i$ for all $i=1,2,3,\cdots,M-1$ and $\beta_{M}=\zeta\beta_{0}$. We also
suppose that all the nodes are stable
($R_{0}^{\text{hom}}=\frac{\beta_{0}}{\gamma}<1$) except the center node $M$.
Then we have a minimum $\zeta^{\text{crit}}$ given by (Eq. Theorem or 15). If
$\frac{1}{R_{0}^{\text{hom}}}<\zeta<\zeta^{\text{crit}}$, there exists an
interval
$\mathcal{I}=\left(\max\left\\{0,\frac{1}{M}-\tilde{\phi}\right\\},\min\left\\{\frac{1}{M-1},\frac{1}{M}+\tilde{\phi}\right\\}\right)$,
where $\tilde{\phi}$ is given by (Eq. 7)), such that if
$\phi_{0}\in\mathcal{I}$, then the DFE is stable, and otherwise it is
unstable. Finally, if $\zeta<\frac{1}{R_{0}^{\text{hom}}}$ i.e.
$\beta_{M}<\gamma$, DFE is always stable.

To analyze a more general symmetric network, we may consider a flux matrix of
the form

$\phi=\begin{cases}\phi_{i,j+1}=\phi_{j}&\text{ $1\leq i\leq
M-1,\leavevmode\nobreak\ \text{and}\leavevmode\nobreak\ i\leq j\leq M-1$}\\\
\phi_{j+1,i}=\phi_{j}&\text{ $1\leq i\leq M-1,\leavevmode\nobreak\
\text{and}\leavevmode\nobreak\ i\leq j\leq M-1$}\\\
\phi_{i,i}=1-\sum_{j=1}^{M-1}\phi_{i,j}&\leavevmode\nobreak\ 1\leq i\leq M\\\
0&\text{ otherwise}\\\ \end{cases}$ (41)

where $\phi_{ij}>0,\leavevmode\nobreak\ \forall\leavevmode\nobreak\ i,j\in V$.
Here we define the fraction of movement from the center node $M$ to the other
nodes by $\phi_{M}=\phi_{0}$ to be consistent with the other networks. We
obtain a particular result for $M=7$ (see appendix for the proof) that we
extend as a second conjecture for this class of most flexible networks.

######  Conjecture 4.

Let $V=\left\\{1,2,3,\ldots,M\right\\}$ be the nodes of networks with the
different wights which are defined by the symmetric flux matrix (Eq. 41),
where $\phi_{ij}\geq 0$. We assume that the infection rate
$\beta_{i}=\beta_{0}$ at node $i$ for all $i=1,2,3,\cdots,M-1$ and
$\beta_{M}=\zeta\beta_{0}$. We also suppose that all the nodes are stable i.e.
$R_{0}^{\text{hom}}=\frac{\beta_{0}}{\gamma}<1$ except the center $M$. Then we
have a minimum $\zeta^{\text{crit}}$ given by (Eq. Theorem or 15). If
$\frac{1}{R_{0}^{\text{hom}}}<\zeta<\zeta^{\text{crit}}$, there exists an
interval
$\mathcal{I}=\left(\max\left\\{0,\frac{1}{M}-\tilde{\phi}\right\\},\min\left\\{\frac{1}{M-1},\frac{1}{M}+\tilde{\phi}\right\\}\right)$,
where $\tilde{\phi}$ is given by (Eq. 7)), such that if
$\phi_{0}\in\mathcal{I}$, then the DFEis stable, and otherwise it is unstable.
Finally, if $\zeta<\frac{1}{R_{0}^{\text{hom}}}$ i.e. $\beta_{M}<\gamma$, then
the DFE is always stable.

##  5 Numerical simulations

We perform numerical simulations to complement our analytical estimates. Our
goal is twofold: (i) to numerically verify the epidemic thresholds obtained
from our theoretical estimates and (ii) to explore the impact of the network
structures on the temporal dynamics of epidemic outbreaks. All simulations are
performed in MATLAB R2020a using the ode45 (4th/5th order Runge-Kutta-Fehlberg
method) function, keeping in mind the order stability area when running
simulations by numerically integrating the dynamics of our model.

### Epidemic thresholds and network structure

In section 4, we discovered a class of networks with the same epidemic
thresholds as the fully connected network. Here we explore the threshold
phenomenon using numerical simulations to determine parametric epidemic or
non-epidemic regions. More specifically, we investigate when the fluxes and
infection rates yielded an epidemic. Fig. 3 illustrates these results for
networks with nine nodes ($M=9$). We numerically determine that an outbreak
occurred when $\max_{t\in\tau}I_{\text{tot}}(t)\geq 2$, where
$I_{\text{tot}}(t)=\sum_{i=1}^{M}I_{i}(t)$ is obtained from the simulations.
In Fig. 3 a., we depict star-shaped, star-cycle, and a more general star-
background network. For those networks, each outer node (blue colored) has the
same transmission rate $\beta_{0}$, and the center node (red colored) has a
different transmission rate, $\beta_{9}=\zeta\beta_{0}$. We then plot
colormaps in the $\phi_{0}\times\beta_{0}$ plane to observe epidemic (yellow-
colored) and non-epidemic (blue-colored) parametric regions. Here $\phi_{0}$
represents the flux corresponding to the nodes connecting with the center and
$\beta_{0}$ is the infection rates at the outer nodes. The similarity among
the colormaps (including the same critical infection rate $\beta_{c}$) is
consistent with our analytical estimates. Fig 3 b highlights the difference in
the colormaps of cycle networks (uni and bi-directional), indicating that
cycles do not have the same epidemic threshold as fully connected or star-
shaped networks. A larger epidemic region for cycle networks suggests that a
low number of connections may promote outbreaks more easily, likely due to the
concentration of infected individuals in the troublesome node (red-colored in
the figure).

Figure 3: Epidemic (yellow colored) _vs._ non-epidemic (blue colored)
parametric regions depending on the flux between nodes and infection rates:
Panel (a) shows _star-shaped, star-cycle, and star-cycle-background_ networks
with their colormaps in the range of the flux $\phi_{0}=[0.01,0.12]$ and
transmission rate $\beta_{0}=[0,1]$, with the step sizes $\Delta\phi_{0}=0.01$
and $\Delta\beta=0.01$ for the network of size $M=9$. Panel (b) represents the
_fully connected, unidirectional, and bidirectional-cycle_ networks along with
the colormaps using the same information of $\phi_{0}$, $\beta_{0}$, and
$\phi_{1}=[0.01,0.12]$. The area above the threshold (shown by the red line)
reveals where the outbreak (yellow parametric region) occurs as predicted by
our analytical estimates. For comparison, in all panels we assume
$\beta_{9}=4\beta_{0}$ and $\phi_{0}\in[0.01,\frac{1}{8}]$, so star-shaped
networks have only positive fluxes. For the more general _star-background_
network connections, we consider two additional fluxes, $\phi_{1}=0.08$ and
$\phi_{2}=0.1$ respectively with flux $\phi_{0}=[0.01,\frac{1}{M-1}]$. For
these simulations, we also consider $\gamma=1$, and the initial conditions
$I_{2}(0)=1,I_{i}(0)=0$ for all $i\in V\setminus 2$;
$S_{2}(0)=9,999,S_{i}(0)=10,000$ for all $i\in V\setminus 2$, where
$V=\\{1,\ldots,9\\}$. We run simulations over the time range $\tau=[0,300]$
with step size $\Delta\tau=0.01$.

### Epidemic thresholds and network size

Cycle and star-class networks differ on the number of node connections, so it
is worth asking if the epidemic thresholds change significantly as the network
size $M$ increases. He numerically obtained $\beta_{0}$ threshold values by
plotting colormaps as in Fig 3 (see Fig 1 in the appendix for details). Fig 4
shows the network geometries and corresponding plots of
$M\times\beta_{0}^{\text{crit}}$ where $M$ is the number of nodes and
$\beta_{0}^{\text{crit}}$ is the minimum infection rate such that epidemics
can be controlled for certain values of the flux. We plot the critical
infection rates for different multiplicative factors $\zeta$, where
$\beta_{M}=\zeta\beta_{0}$. When $M$ increases, unidirectional and
bidirectional cycles exhibit epidemic thresholds that seem constant with
respect to the network size. A possible explanation for this phenomenon might
be the little importance of having more nodes in the cycle, given that the
epidemic likely occurs in the heterogeneous node and spreads to the
neighboring nodes only. On the other hand, star-class networks exhibit a
different phenomenon, where the critical infection rate converges to a fixed
number. In fact, if our conjecture 3 is true, then $\beta_{0}^{\text{crit}}$
for the star-class network is given by the expression

$\beta_{0}^{\text{crit}}=\left(\frac{M}{\zeta+M-1}\right)\gamma$

and thus $\beta_{0}^{\text{crit}}\to\gamma$ as $M\to\infty$, which agrees with
panel c in Fig 4 (in the simulations we consider $\gamma=1$). It is worth
noting that such _asymptotic_ critical value does not depend on the
multiplicative factor $\zeta=\frac{\beta_{M}}{\beta_{0}}$. In contrast, for
cycle networks $\beta_{0}^{\text{crit}}$ converges to different values
depending on $\zeta$, a feature that is also likely connected with the network
geometry.

Figure 4: Critical infection rates as network size $M$ increases: We plot
$\beta_{0}^{\text{crit}}$ for $M=1,2,...,50$ and different $\zeta$ values for
(a) unidirectional cycles, (b) bidirectional cycles and (c) the general star-
class network. Critical values are numerically obtained with colormaps as in
Fig 3 (see Fig 7 in the Appendix for details). For each simulation, we assume
that every node has the same infection rate $\beta_{0}$, except the center
node $M$, where the infection rate is $\beta_{M}=\zeta\beta_{0}$.

### Temporal dynamics of epidemic outbreaks

To better understand the temporal evolution of epidemic outbreaks on different
networks, we run numerical simulations for star-shaped, fully connected and
cycle (uni and bi-directional) networks and for three scenarios where flux
values are low, medium, and high. In Fig. 5 illustrates the chosen networks of
nine nodes ($M=9$) with stable blue colored nodes ($\beta_{0}<\gamma$) and the
red colored heterogeneous node ($\beta_{9}>\gamma$). We then plot the time
series of the infected population for flux values $0.01$, $0.06$, and $0.12$
at a single node $I_{2}$, the heterogeneous node $I_{9}$, and for all stable
nodes summed (represented by $I_{1}+I_{2}+...+I_{8}$). The infection dynamics
in star-shaped and fully connected networks (Fig. 5 a and b) are similar: when
the flux increases, the disease also spreads significantly from the center to
the other nodes. As a result, the sum of the infected individuals in nodes $1$
to $8$ increases, decreasing the epidemic at node $9$. In contrast, for
unidirectional and bidirectional cycles (Fig. 5 c and d) the epidemic peak
remains high in the heterogeneous node for all three flux values. The outbreak
spreads to the other nodes when the flux is high, but the peak still occurs
earlier compared to the dynamics at star-shaped and fully connected networks.
The difference between star-shaped and cycle networks is notable for low flux
values, as cycle networks exhibit a much more localized outbreak, which is
likely due to the lack of connections with the heterogeneous node.

Figure 5: Epidemic dynamics for low, medium, and high flux values:
Trajectories of the infected populations for star-shaped, fully connected, and
cycle networks, considering three scenarios of low, medium, and high fluxes.
Simulated epidemic individuals are obtained for time courses $\tau=[0,60]$
with the step size $\Delta\tau=0.01$, recovery rate $\gamma=1$, and infection
rate $\beta_{i}=\beta_{0}=0.95\leavevmode\nobreak\ \forall i\in V\setminus M$,
and $\beta_{M}=\zeta\beta_{0}$, where $\zeta=4$. Moreover, the initial
conditions are $I_{2}(0)=1,I_{i}(0)=0$ for all $i\in V\setminus 2$, and
$S_{2}(0)=9,999,S_{i}(0)=10,000$ for all $i\in V\setminus 2$.

### Epidemic peak intensity and timing across networks

We further investigate outbreak properties by plotting the maximum peak
intensity $I^{max}_{i}$ and the corresponding peak time $t(I^{max}_{i})$ for
each node $i$, considering star-shaped, fully connected and cycle networks
with nine nodes (Fig 6). For low flux values ($\phi_{0}=\phi_{1}=0.01$) and
for all four network structures, we observe a higher peak at the heterogeneous
node ($i=9$) in comparison to the other nodes (green dots). For the peak time
$t(I^{max}_{i})$, the network structure seemed to play a major role. In fact,
for both star-shaped and fully connected networks, $t(I^{max}_{i})$ had a
similar value of around 8 weeks for all stable nodes ($i=1,2,...,8$) being
lower for node 9 (around 6 weeks). In contrast, cycle networks exhibited
_delayed_ peaks for nodes that are distant from the heterogeneous node 9. For
example, nodes 4 and 5 reached their maximum number of infected only 30/40
weeks after the initial time. For higher flux values ($\phi_{0}=\phi_{1}=0.06$
and $\phi_{0}=\phi_{1}=0.12$), the peak intensities and timing were comparable
between stable and heterogeneous nodes for both star-shaped and fully
connected networks, thus indicating the role of mobility promoting the spread
of the disease. For cycle networks, we again observe delayed and damped peaks
for nodes at a higher distance of the heterogeneous node 9 due to the low
number of connections.

Figure 6: The intensity and timing of epidemic peaks: the maximum number of
infected individuals $I^{max}_{i}$ and peak time $t(I^{max}_{i})$ for each
node of star-shaped, fully connected, unidirectional, and bidirectional cycle
networks. When the flux is low, the center node significantly has the highest
peak in each network, and the epidemic occurs more quickly than the other
nodes. However, the disease spreads more rapidly in Star-shaped and fully
connected networks compared to cycle networks and shows the same dynamical
behavior in Star-shaped and fully connected networks. For cycle networks, the
highest peaks still occur at the center node for any fluxes. Results were
obtained by using the same parameter values described in Fig. 5.

##  6 Conclusions

In this study, we modeled the spread of a disease in a metapopulation, i.e., a
network of subpopulations (nodes) connected by the flux of individuals
(edges). We investigated epidemic thresholds using the SIR-network model in
the case where all subpopulations are equally sized and all nodes have the
same infection rate except one (the heterogeneous node). To determine outbreak
conditions, we utilized the classical next-generation matrix approach to find
explicit formulas for the intervals in which the flux of people may control
epidemics on different network geometries. Interestingly, we identified a
class of networks with the same epidemic thresholds as fully connected
networks, previously established in (?). This class is characterized by a
star-shaped structure, where a central heterogeneous node is equally connected
to all other nodes. On the other hand, cycle networks exhibited different
thresholds, suggesting that epidemics on poorly connected networks manifest
differently while still being influenced by infection rates and mobility
factors.

Our numerical simulations revealed distinct epidemic thresholds depending on
the network geometry. Fig 3 illustrates how the flux ($\phi_{0}$) may control
or promote outbreaks more prominently in star-shaped networks compared to
cycle networks. For cycles, the blue/yellow colormap division was primarily
driven by changes in the infection rates $\beta_{0}$. The epidemic regions for
cycles were larger than star-shaped networks, indicating that a highly
connected troublesome node may promote and control the disease spread.
Regarding temporal dynamics, cycle networks produced more localized outbreaks
at low flux values than star-shaped and fully connected networks. At high flux
values, epidemics on cycles also peaked faster, mainly in the heterogeneous
node. At the same time, star-shaped and fully connected networks could
“spread” the outbreak to the whole network.

Metapopulation models have played a significant role in exploring epidemic
thresholds and dynamics of infectious diseases on networks. Tocto-Erazo
recently explored the contribution of periodic human mobility on the dynamics
of Dengue fever in Hermosillo, Mexico (?). Also in the context of vector-borne
diseases, Iggidr and colleagues (?) explored the effects of population
circulation on strongly connected networks, estimating $R_{0}$ and analyzing
the stability of both disease-free and endemic equilibrium. Anzo-Hernández et
al. (?) conducted a detailed study on the influence of human mobility on
$R_{0}$ for vector-borne diseases, providing explicit formulas and numerical
simulations for various network configurations. The impact of vaccination game
theoretical approaches (?) and awareness with information spreading (?) has
also been explored along with mobility factors and network geometries to
provide insights on disease outbreaks and progression. In this work, we
adopted a different approach by systematically examining how mobility and
geometry shape epidemic thresholds for networks of increased complexity.
Unlike previous studies that focused on specific diseases or analyzed reduced-
node network structures, we established theoretical estimates that hold for
networks with an arbitrary number of nodes, including star-shaped, star-
background, and star-triangle networks. Our formulas provide explicit
parametric regions in which the flux of people can either control or promote
outbreaks. Additionally, our numerical simulations revealed critical
differences in the temporal dynamics of diseases based on the network
geometry.

The insights gained from this pave the way for potential future research
directions. Subsequent studies could delve into the role of population
heterogeneity by extending our findings to networks with subpopulations of
different sizes. For instance, it is reasonable to expect that a central node,
such as a city center or transportation hub, can be also highly populated.
Investigating the influence of highly populated central nodes could thus
provide valuable insights. Furthermore, future studies can relax the
assumption of symmetric flux matrices and explore more heterogeneous movement
patterns. While all networks considered in this work are theoretical, an
intriguing alternative study could involve extracting network features from
real-world mobility data. Finally, it is important to acknowledge that the
SIR-network model, originally conceived as a simplified model for dengue
fever, can generally describe diseases with direct transmission mechanisms.
However, it lacks specific elements such as loss of immunity, exposure
periods, and demographic changes. Therefore, an interesting avenue for further
exploration would be incorporating those factors in the SIR-network model and
examining their impact on epidemic thresholds and disease dynamics.

## Acknowledgements

The authors would like to acknowledge Cort Vanzant (OSU) and Dr. Pedro Maia
(UTA) for their careful reading and feedback on this work. The authors also
thank Dr. Anthony Kable (OSU) for the insightful comments, especially for
proposing the study of cycle-support networks. LMS thanks Dr. Stefanella
Boatto (Federal University of Rio de Janeiro) for initial discussions about
the SIR-network model and supportive mentorship during his undergraduate and
master’s studies.

## References

  * 1. World Health Organization.  Novel coronavirus (2019-ncov).  Retrieved April 11, 2023, from https://www.who.int/emergencies/diseases/novel-coronavirus-2019, 2023.
  * 2. Kate E Jones, Nikkita G Patel, Marc A Levy, Adam Storeygard, Deborah Balk, John L Gittleman, and Peter Daszak.  Global trends in emerging infectious diseases.  Nature, 451(7181):990–993, 2008.
  * 3. Stephen S Morse.  Factors in the emergence of infectious diseases.  Springer, 2001.
  * 4. S.C. Quinn and S. Kumar.  Health inequalities and infectious disease epidemics: a challenge for global health security.  Biosecurity and Bioterrorism: Biodefense Strategy, Practice, and Science, 12(5):263–273, 2014.
  * 5. J Piret and G Boivin.  Pandemics throughout history.  Frontiers in Microbiology, 11:631736, 2021.
  * 6. Vivek Charu, Scott Zeger, Julia Gog, Ottar N Bjørnstad, Stephen Kissler, Lone Simonsen, Bryan T Grenfell, and Cécile Viboud.  Human mobility and the spatial transmission of influenza in the united states.  PLoS computational biology, 13(2):e1005382, 2017.
  * 7. Stefano Merler and Marco Ajelli.  The role of population heterogeneity and human mobility in the spread of pandemic influenza.  Proceedings of the Royal Society B: Biological Sciences, 277(1681):557–565, 2010.
  * 8. Ben Adams and Durrell D Kapan.  Man bites mosquito: understanding the contribution of human movement to vector-borne disease dynamics.  PloS one, 4(8):e6763, 2009.
  * 9. Amy Wesolowski, Taimur Qureshi, Maciej F Boni, Pål Roe Sundsøy, Michael A Johansson, Syed Basit Rasheed, Kenth Engø-Monsen, and Caroline 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.
  * 10. Andrew M Kramer, J Tomlin Pulliam, Laura W Alexander, Andrew W Park, Pejman Rohani, and John M Drake.  Spatial spread of the west africa ebola epidemic.  Royal Society open science, 3(8):160294, 2016.
  * 11. James O Lloyd-Smith, Sebastian J Schreiber, P Ekkehard Kopp, and Wayne M Getz.  Superspreading and the effect of individual variation on disease emergence.  Nature, 438(7066):355–359, 2005.
  * 12. David E Bloom and Daniel Cadarette.  Infectious disease threats in the twenty-first century: strengthening the global response.  Frontiers in immunology, 10:549, 2019.
  * 13. Andrew J Tatem, David J Rogers, and Simon I Hay.  Global transport networks and infectious disease spread.  Advances in parasitology, 62:293–343, 2006.
  * 14. Cécile Viboud, Ottar N Bjørnstad, David L Smith, Lone Simonsen, Mark A Miller, and Bryan T Grenfell.  Synchrony, waves, and spatial hierarchies in the spread of influenza.  Science, 312(5772):447–451, 2006.
  * 15. Neil M Ferguson, Derek AT Cummings, Simon Cauchemez, Christophe Fraser, Steven Riley, Aronrag Meeyai, Sopon Iamsirithaworn, and Donald S Burke.  Strategies for containing an emerging influenza pandemic in southeast asia.  Nature, 437(7056):209–214, 2005.
  * 16. Vittoria Colizza, Alain Barrat, Marc Barthélemy, Alessandro Vespignani, and SIRIUS Friends.  Modeling the worldwide spread of pandemic influenza: baseline case and containment interventions.  PloS one, 2(12):e1316, 2006.
  * 17. L. Sattenspiel and K. Dietz.  A structured epidemic model incorporating geographic mobility among regions.  Mathematical Biosciences, 128(1-2):71–91, 1995.
  * 18. M. J. Keeling and T. D. Eames.  Networks and epidemic models.  J. R. Soc. Interface, 2:295–317, 2005.
  * 19. J. Arino and P. Van Den Driessche.  A multi-city epidemic model.  Mathematical Population Studies, 10(3):175–193, 2003.
  * 20. Glenn R Fulford, Michael G Roberts, and James AP Heesterbeek.  The metapopulation dynamics of an infectious disease: tuberculosis in possums.  Theoretical population biology, 61(1):15–29, 2002.
  * 21. Daniela Calvetti, Alexander P Hoover, Johnie Rose, and Erkki Somersalo.  Metapopulation network models for understanding, predicting, and managing the coronavirus disease covid-19.  Frontiers in Physics, 8:261, 2020.
  * 22. L. M. Stolerman, D. Coombs, and S. Boatto.  Sir-network model and its application to dengue fever.  SIAM Journal Applied Math, 75(6):2581–2609, 2015.
  * 23. Fabio Della Rossa, Davide Salzano, Anna Di Meglio, Francesco De Lellis, Marco Coraggio, Carmela Calabrese, Agostino Guarino, Ricardo Cardona-Rivera, Pietro De Lellis, Davide Liuzza, et al.  A network model of italy shows that intermittent regional strategies can alleviate the covid-19 epidemic.  Nature communications, 11(1):1–9, 2020.
  * 24. Sheryl Le Chang, Mahendra Piraveenan, and Mikhail Prokopenko.  The effects of imitation dynamics on vaccination behaviours in sir-network model.  International journal of environmental research and public health, 16(14):2477, 2019.
  * 25. Sheryl L Chang, Mahendra Piraveenan, and Mikhail Prokopenko.  Impact of network assortativity on epidemic and vaccination behaviour.  Chaos, solitons & fractals, 140:110143, 2020.
  * 26. Nathan Harding, Richard E Spinney, and Mikhail Prokopenko.  Population mobility induced phase separation in sis epidemic and social dynamics.  Scientific reports, 10(1):7646, 2020.
  * 27. Klaus Dietz.  The estimation of the basic reproduction number for infectious diseases.  Statistical methods in medical research, 2(1):23–41, 1993.
  * 28. Paul L Delamater, Erica J Street, Timothy F Leslie, Y Tony Yang, and Kathryn H Jacobsen.  Complexity of the basic reproduction number (r0).  Emerging infectious diseases, 25(1):1, 2019.
  * 29. Odo Diekmann, JAP Heesterbeek, and Michael G Roberts.  The construction of next-generation matrices for compartmental epidemic models.  Journal of the royal society interface, 7(47):873–885, 2010.
  * 30. Odo Diekmann and Johan Andre Peter Heesterbeek.  Mathematical epidemiology of infectious diseases: model building, analysis and interpretation, volume 5.  John Wiley & Sons, 2000.
  * 31. Aberrahman Iggidr, Gauthier Sallet, and Max O Souza.  On the dynamics of a class of multi-group models for vector-borne diseases.  Journal of Mathematical Analysis and applications, 441(2):723–743, 2016.
  * 32. Mayra R Tocto-Erazo, Daniel Olmos-Liceaga, and José A Montoya.  Effect of daily periodic human movement on dengue dynamics: The case of the 2010 outbreak in hermosillo, mexico.  Applied Mathematical Modelling, 97:559–567, 2021.
  * 33. Steven H. Strogatz.  Nonlinear Dynamics and Chaos: with Applications to Physics, Biology, Chemistry, and Engineering.  Westview Press, 2015.
  * 34. Lawrence Perko.  Differential equations and dynamical systems, volume 7.  Springer Science & Business Media, 2013.
  * 35. P. Van Den Driessche and J. Watmough.  Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission.  Mathematical Biosciences, 180:29–48, 2002.
  * 36. Stephen Friedberg, Arnold Insel, and Lawrence Spence.  Linear Algebra.  Pearson, 4th edition, 2019.
  * 37. B. J. Collings.  Characteristic polynomials by diagonal expansion.  The American Statistician, 37(3):233–235, 1983.
  * 38. James Dickson Murray.  Mathematical biology: I. An introduction.  Springer, 2002.
  * 39. Andres Anzo-Hernández, Beatriz Bonilla-Capilla, Jorge Velázquez-Castro, Moisés Soto-Bajo, and A Fraguela-Collar.  The risk matrix of vector-borne diseases in metapopulation networks and its relation with local and global r0.  Communications in Nonlinear Science and Numerical Simulation, 68:1–14, 2019.
  * 40. KM Ariful Kabir and Jun Tanimoto.  Evolutionary vaccination game approach in metapopulation migration model with information spreading on different graphs.  Chaos, Solitons & Fractals, 120:41–55, 2019.
  * 41. KM Ariful Kabir and Jun Tanimotoc.  Impact of awareness in metapopulation epidemic model to suppress the infected individuals for different graphs.  The European Physical Journal B, 92:1–16, 2019.

## Appendix

##  Appendix A A few proofs of Propositions and Theorems

The purpose of this appendix is to present a few proofs of the Propositions
and Theorems that are outlined in the main manuscript.

Proposition 2 Let $V=\left\\{1,2,3\right\\}$ be the node of a _unidirectional
cycle_ network and the flux matrix $\phi_{3\times 3}$ be defined by (Eq. 24).
We consider the infection rates $\beta_{i}$ from (Eq. 8), where all the nodes
are stable except the third one, i.e, $R_{0}^{\text{hom}}<1$ and
$\frac{\zeta\beta_{0}}{\gamma}>1$. Then we have a critical value
$\zeta^{\text{crit}}$ given by

$\displaystyle\zeta^{\text{crit}}=\frac{1}{R_{0}^{\text{hom}}}\left(\frac{3R_{0}^{\text{hom}}-4}{R_{0}^{\text{hom}}-2}\right)$
(42)

such that if $\frac{1}{R_{0}^{\text{hom}}}<\zeta<\zeta^{\text{crit}}$, then
there exists an interval $\mathcal{I}$ around the point $\frac{1}{2}$ given by

$\mathcal{I}=\left(\max\left\\{0,\frac{1}{2}-\tilde{\phi}\right\\},\min\left\\{1,\frac{1}{2}+\tilde{\phi}\right\\}\right),$
(43)

where

$\displaystyle\tilde{\phi}=\frac{1}{2}\sqrt{\frac{\left(\zeta
R_{0}^{\text{hom}}(R_{0}^{\text{hom}}-2)+(4-3R_{0}^{\text{hom}})\right)}{R_{0}^{\text{hom}}\left(1+2\zeta-3\zeta
R_{0}^{\text{hom}}\right)}}$ (44)

such that if $\phi_{0}\in\mathcal{I}$, the DFE is stable, and otherwise, it is
unstable.

###### Proof.

To establish the stability criteria for the DFE, we may perform a local
stability analysis of the DFE. The Jacobian matrix at the DFE corresponding to
the dynamical system of a $3\times 3$ star-shaped flux matrix is given by the
following expression:

$\displaystyle J_{3\times
3}=\begin{pmatrix}\beta_{1}(-1+\phi_{0})^{2}+\beta_{2}\phi_{0}^{2}-\gamma&-\beta_{2}(-1+\phi_{0})\phi_{0}&-\beta_{1}\left(-1+\phi_{0}\right)\phi_{0}\\\
-\beta_{2}(-1+\phi_{0})\phi_{0}&\beta_{2}\left(-1+\phi_{0}\right)^{2}+\beta_{3}\phi_{0}^{2}-\gamma&-\beta_{3}(-1+\phi_{0})\phi_{0}\\\
-\beta_{1}(-1+\phi_{0})\phi_{0}&-\beta_{3}(-1+\phi_{0})\phi_{0}&\beta_{3}(-1-\phi_{0})^{2}+\beta_{1}\phi_{0}^{2}-\gamma\end{pmatrix}.$
(45)

As symmetric matrices have only real eigenvalues, the eigenvalues of the
Jacobian matrix (Eq. 45) are real. For simplicity, we calculate the
eigenvalues using the symbolic packages of the software Mathematica, which are
given by:

$\displaystyle\lambda_{1}=(\beta_{0}-\gamma)-3\beta_{0}\phi_{0}(1-\phi_{0})$
$\displaystyle\lambda_{2}=-\gamma+\frac{1}{2}\beta_{0}\left(1+\zeta-\phi_{0}-2\zeta\phi_{0}+\phi_{0}^{2}+2\zeta\phi_{0}^{2}+B\right)$
$\displaystyle\lambda_{3}=-\gamma+\frac{1}{2}\beta_{0}\left(1+\zeta-\phi_{0}-2\zeta\phi_{0}+\phi_{0}^{2}+2\zeta\phi_{0}^{2}-B\right),$

where
$B:=\sqrt{(1-\phi_{0}+\phi_{0}^{2})^{2}+\zeta^{2}\left(1-2\phi_{0}+2\phi_{0}^{2}\right)^{2}+\zeta\left(-2+6\phi_{0}-2\phi_{0}^{2}-8\phi_{0}^{3}+4\phi_{0}^{4}\right)}.$
As the jacobian matrix $J_{3\times 3}$ (Eq.45) is a real symmetric matrix, the
eigenvalues are real. Since $\lambda_{2}$ and $\lambda_{3}$ are real, B must
be nonnegative for $\phi_{0}\in(0,1)$ and for any $\zeta$, $B$ is well-defined
and positive real number. We want to find conditions for which all eigenvalues
are negative (LAS DFE), or at least one eigenvalue is positive (unstable DFE).
We verify that $\lambda_{1}$ is always negative for all $\phi_{0}$, since by
hypothesis $\frac{\beta_{0}}{\gamma}<1$. Then we may check $\lambda_{2}$ and
$\lambda_{3}$, which are more complicated. However, a sufficient condition for
stability is the following: since $\lambda_{2}>\lambda_{3}$, then we do have
stability if the minimum value of $\lambda_{2}$ is negative. It can be easily
verified that $\frac{1}{2}$ is the critical point of the eigenvalues
$\lambda_{2}$, and $\lambda_{3}$, since
$\frac{\partial\lambda_{2}}{\partial\phi_{0}}(\phi_{0}=\frac{1}{2},\beta_{0},\zeta,\gamma)=0$,
and
$\frac{\partial\lambda_{3}}{\partial\phi_{0}}(\phi_{0}=\frac{1}{2},\beta_{0},\zeta,\gamma)=0$.

We find that

$\displaystyle\frac{\partial^{2}\lambda_{2}}{\partial\phi_{0}^{2}}(\phi_{0}=\frac{1}{2},\beta_{0},\zeta,\gamma)=1+2\zeta+\frac{3-16\zeta+4\zeta^{2}}{\sqrt{9-4\zeta+4\zeta^{2}}}>0$
$\displaystyle\frac{\partial^{2}\lambda_{3}}{\partial\phi_{0}^{2}}(\phi_{0}=\frac{1}{2},\beta_{0},\zeta,\gamma)=1+2\zeta-\frac{3-16\zeta+4\zeta^{2}}{\sqrt{9-4\zeta+4\zeta^{2}}}>0$

and so $\lambda_{2}$, and $\lambda_{3}$ have reached to the minimum value when
$\phi_{0}=\frac{1}{2}$. Moreover,

$\displaystyle\lambda_{1}(\phi_{0},\beta_{0},\zeta,\gamma)<0\leavevmode\nobreak\
\leavevmode\nobreak\ \forall\leavevmode\nobreak\ \leavevmode\nobreak\
\phi_{0}\in(0,1)$
$\displaystyle\lambda_{2}(\phi_{0}=\frac{1}{2},\beta_{0},\zeta,\gamma)=\frac{1}{8}\left(3+2\zeta+\sqrt{9-4\zeta+4\zeta^{2}}\right)\beta_{0}-\gamma$
$\displaystyle\lambda_{3}(\phi_{0}=\frac{1}{2},\beta_{0},\zeta,\gamma)=\frac{1}{8}\left(3+2\zeta-\sqrt{9-4\zeta+4\zeta^{2}}\right)\beta_{0}-\gamma.$

Since $\lambda_{2}>\lambda_{3}$, so it is enough to find the condition for
which
$\lambda_{2}(\phi_{0}=\frac{1}{2},\beta_{0},\zeta,\gamma)=\frac{1}{8}\left(3+2\zeta+\sqrt{9-4\zeta+4\zeta^{2}}\right)\beta_{0}-\gamma<0$.
This implies that

$\displaystyle\zeta<\frac{\gamma(3\beta_{0}-4\gamma)}{\beta_{0}(\beta_{0}-2\gamma)}=\frac{1}{R_{0}^{\text{hom}}}\left(\frac{3R_{0}^{\text{hom}}-4}{R_{0}^{\text{hom}}-2}\right):=\zeta^{\text{crit}}.$

In addition, the function
$\phi_{0}\mapsto\lambda_{2}(\phi_{0},\zeta,\beta_{0},\gamma)$ has roots

$\displaystyle\phi_{c1}=\frac{1}{2}+\frac{1}{2}\frac{\sqrt{-\beta_{0}(3\zeta\beta_{0}-\gamma-2\zeta\gamma)\left(\zeta\beta_{0}(\beta_{0}-2\gamma)+\gamma(-\beta_{0}+4\gamma)\right)}}{\beta_{0}(-3\zeta\beta_{0}+\gamma+2\zeta\gamma)}$
$\displaystyle=\frac{1}{2}+\frac{1}{2}\sqrt{\frac{\left(\zeta
R_{0}^{\text{hom}}(R_{0}^{\text{hom}}-2)+(4-3R_{0}^{\text{hom}})\right)}{R_{0}^{\text{hom}}\left(1+2\zeta-3\zeta
R_{0}^{\text{hom}}\right)}}$
$\displaystyle\phi_{c2}=\frac{1}{2}-\frac{1}{2}\frac{\sqrt{-\beta_{0}(3\zeta\beta_{0}-\gamma-2\zeta\gamma)\left(\zeta\beta_{0}(\beta_{0}-2\gamma)+\gamma(-\beta_{0}+4\gamma)\right)}}{\beta_{0}(-3\zeta\beta_{0}+\gamma+2\zeta\gamma)}$
$\displaystyle=\frac{1}{2}-\frac{1}{2}\sqrt{\frac{\left(\zeta
R_{0}^{\text{hom}}(R_{0}^{\text{hom}}-2)+(4-3R_{0}^{\text{hom}})\right)}{R_{0}^{\text{hom}}\left(1+2\zeta-3\zeta
R_{0}^{\text{hom}}\right)}}$

Notice that $\phi_{c1}$, and $\phi_{c2}$ both are well defined since
$\left(1+2\zeta-3\zeta R_{0}^{\text{hom}}\right)$, and $\zeta
R_{0}^{\text{hom}}(R_{0}^{\text{hom}}-2)+(4-3R_{0}^{\text{hom}})$ both are
negative and $R_{0}^{\text{hom}}>1$. Therefore, according to the restriction
$\phi_{0}\in\left(0,1\right)$, the interval where there will be no epidemic is
given by
$\mathcal{I}(\zeta,R_{0}^{\text{hom}})=\left(\max\left\\{0,\frac{1}{2}-\tilde{\phi}\right\\},\min\left\\{1,\frac{1}{2}+\tilde{\phi}\right\\}\right)$,
where

$\displaystyle\tilde{\phi}=\frac{1}{2}\sqrt{\frac{\left(\zeta
R_{0}^{\text{hom}}(R_{0}^{\text{hom}}-2)+(4-3R_{0}^{\text{hom}})\right)}{R_{0}^{\text{hom}}\left(1+2\zeta-3\zeta
R_{0}^{\text{hom}}\right)}}.$

Therefore, if $\phi_{0}\in\mathcal{I}$, then the DFE is stable and, otherwise,
it is unstable. ∎

Proposition 3

Let $V=\left\\{1,2,3,4,5\right\\}$ be a node of _cycle-supports_ network and
let the geometry of its flux matrix be defined by (Eq. 29), where
$\phi_{0}\in(0,\frac{1}{4}]$. The infection rate $\beta_{i}$ is defined in
(Eq. 8), and we also assume that all the outer nodes are stable
($R_{0}^{\text{hom}}=\frac{\beta_{0}}{\gamma}<1$) except the center node
$M=5$. Then the critical value of $\zeta$ called $\zeta^{\text{crit}}$ is
given by

$\zeta^{\text{crit}}=1+\left(\frac{1-R_{0}^{\text{hom}}}{R_{0}^{\text{hom}}}\right)5.$

Alternatively, the infection threshold of $\beta_{0}$,
$\beta_{0}^{\text{crit}}$, is given by

$\beta_{0}^{\text{crit}}=\gamma\left(\frac{5}{4+\zeta}\right).$ (46)

If $\frac{1}{R_{0}^{\text{hom}}}<\zeta<\zeta^{\text{crit}}$, then there is an
interval around the point $\frac{1}{5}$ given by
$\mathcal{I}=\left(\max\left\\{0,\frac{1}{5}-\tilde{\phi}\right\\},\min\left\\{\frac{1}{4},\frac{1}{5}+\tilde{\phi}\right\\}\right)$,
where

$\displaystyle\tilde{\phi}=\frac{1}{5}\sqrt{\frac{\gamma(5\gamma-4\beta_{0}-\zeta\beta_{0})}{\beta_{0}(\gamma+4\zeta\gamma-5\zeta\beta_{0})}}=\frac{1}{5}\sqrt{\frac{5(1-R_{0}^{\text{hom}})-R_{0}^{\text{hom}}(\zeta-1))}{R_{0}^{\text{hom}}\left(5\zeta(1-R_{0}^{\text{hom}})-(\zeta-1)\right)}},$

such that if $\phi_{0}\in\mathcal{I}$, the DFE is stable; and otherwise it is
unstable.

###### Proof.

To establish the stability criteria for the DFE in this case, we may perform a
local stability analysis of the DFE. The Jacobian matrix at the DFE
corresponding to the dynamical system of the $5\times 5$ cycle support flux
matrix (Eq. 29) is symmetric. As symmetric matrices have only real
eigenvalues, the eigenvalues of the above Jacobian matrix are real. For
simplicity, we calculate the eigenvalues using the symbolic packages of the
software Mathematica, which are given by

$\displaystyle\lambda_{1}=\beta_{0}(1-3\phi_{0})^{2}-\gamma$
$\displaystyle\lambda_{2}=\beta_{0}-\gamma-6\beta_{0}\phi_{0}+11\beta_{0}\phi_{0}^{2}-2\sqrt{2}\sqrt{\beta_{0}^{2}(1-3\phi_{0})^{2}\phi_{0}^{2}}$
$\displaystyle\lambda_{3}=\beta_{0}-\gamma-6\beta_{0}\phi_{0}+11\beta_{0}\phi_{0}^{2}+2\sqrt{2}\sqrt{\beta_{0}^{2}(1-3\phi_{0})^{2}\phi_{0}^{2}}$
$\displaystyle\lambda_{4}=-\gamma+\frac{1}{2}\beta_{0}\left(1+\zeta-2\phi_{0}-8\zeta\phi_{0}+5\phi_{0}^{2}+20\zeta\phi_{0}^{2}+B\right)$
$\displaystyle\lambda_{5}=-\gamma+\frac{1}{2}\beta_{0}\left(1+\zeta-2\phi_{0}-8\zeta\phi_{0}+5\phi_{0}^{2}+20\zeta\phi_{0}^{2}-B\right)$

where
$B:=\sqrt{(1-2\phi_{0}+5\phi_{0}^{2})^{2}+\zeta^{2}\left(1-8\phi_{0}+20\phi_{0}^{2}\right)^{2}+2\zeta\left(-1+10\phi_{0}-9\phi_{0}^{2}-80\phi_{0}^{3}+100\phi_{0}^{4}\right)}.$
Since all eigenvalues are real, B must be nonnegative for
$\phi_{0}\in(0,\frac{1}{4}]$, and for any $\zeta$, $B$ is a positive real
number. We want to find conditions for which all eigenvalues are negative (LAS
DFE), or at least one eigenvalue is positive (unstable DFE). One can see that
$\lambda_{1}$, $\lambda_{2}$, and $\lambda_{3}$ are negative for all
$\phi_{0}\in\left(0,\frac{1}{4}\right]$, since by hypothesis
$\frac{\beta_{0}}{\gamma}<1$. It is thus sufficient to check the critical
points of $\lambda_{4}$ and $\lambda_{5}$. However, a sufficient condition for
stability is the following: since $\lambda_{4}>\lambda_{5}$, it is thus
sufficient to check the stability of $\lambda_{4}$. It can be easily verified
that $\frac{1}{5}$ is the critical point of both eigenvalues $\lambda_{4}$,
and $\lambda_{5}$, since
$\frac{\partial\lambda_{4}}{\partial\phi_{0}}(\phi_{0}=\frac{1}{5},\beta_{0},\zeta,\gamma)=0$,
and
$\frac{\partial\lambda_{5}}{\partial\phi_{0}}(\phi_{0}=\frac{1}{5},\beta_{0},\zeta,\gamma)=0$.
We find that

$\displaystyle\frac{\partial^{2}\lambda_{4}}{\partial\phi_{0}^{2}}(\phi_{0}=\frac{1}{5},\beta_{0},\zeta,\gamma)=\frac{40\beta_{0}(\zeta-1)^{2}}{4+\zeta}>0$
$\displaystyle\frac{\partial^{2}\lambda_{5}}{\partial\phi_{0}^{2}}(\phi_{0}=\frac{1}{5},\beta_{0},\zeta,\gamma)=\frac{250\zeta\beta_{0}}{4+\zeta}>0$

and so $\lambda_{4}$ and $\lambda_{5}$ have reached to the minimum value when
$\phi_{0}=\frac{1}{5}$. Moreover,

$\displaystyle\lambda_{i}(\phi_{0},\beta_{0},\zeta,\gamma)<0\leavevmode\nobreak\
\leavevmode\nobreak\ \forall\leavevmode\nobreak\
\text{for}\leavevmode\nobreak\ \phi_{0}\in(0,\frac{1}{4}],\leavevmode\nobreak\
\text{where},i=1,2,3$
$\displaystyle\lambda_{4}(\phi_{0}=\frac{1}{5},\beta_{0},\zeta,\gamma)=\frac{1}{10}\left(4\beta_{0}+\zeta\beta_{0}+\sqrt{(4+\zeta)^{2}}\beta_{0}-10\gamma\right)=\frac{1}{5}\beta_{0}(4+\zeta)-\gamma$
$\displaystyle\lambda_{5}(\phi_{0}=\frac{1}{5},\beta_{0},\zeta,$

However, if $\frac{1}{5}\beta_{0}(4+\zeta)-\gamma<0$ and that implies

$\displaystyle\beta_{0}<\gamma\left(\frac{5}{4+\zeta}\right):=\beta_{0}^{\text{crit}}\leavevmode\nobreak\
\leavevmode\nobreak\ \text{and\leavevmode\nobreak\ alternatively}$
$\displaystyle\zeta<\frac{5\gamma-4\beta_{0}}{\beta_{0}}=\frac{5-4\frac{\beta_{0}}{\gamma}}{\frac{\beta_{0}}{\gamma}}=\frac{5-4R_{0}^{\text{hom}}}{R_{0}^{\text{hom}}}=1+\left(\frac{1-R_{0}^{\text{hom}}}{R_{0}^{\text{hom}}}\right)5:=\zeta^{\text{crit}}.$

In addition, the function
$\phi_{0}\mapsto\lambda_{4}(\phi_{0},\zeta,\beta_{0},\gamma)$ has roots

$\displaystyle\phi_{c1}=\frac{1}{5}+\frac{1}{5}\sqrt{\frac{\gamma(5\gamma-4\beta_{0}-\zeta\beta_{0})}{\beta_{0}(\gamma+4\zeta\gamma-5\zeta\beta_{0})}}=\frac{1}{5}+\frac{1}{5}\sqrt{\frac{5(1-R_{0}^{\text{hom}})-R_{0}^{\text{hom}}(\zeta-1))}{R_{0}^{\text{hom}}\left(5\zeta(1-R_{0}^{\text{hom}})-(\zeta-1)\right)}},$
$\displaystyle\phi_{c2}=\frac{1}{5}-\frac{1}{5}\sqrt{\frac{\gamma(5\gamma-4\beta_{0}-\zeta\beta_{0})}{\beta_{0}(\gamma+4\zeta\gamma-5\zeta\beta_{0})}}=\frac{1}{5}-\frac{1}{5}\sqrt{\frac{5(1-R_{0}^{\text{hom}})-R_{0}^{\text{hom}}(\zeta-1))}{R_{0}^{\text{hom}}\left(5\zeta(1-R_{0}^{\text{hom}})-(\zeta-1)\right)}}.$

Therefore, taking into account the restriction
$\phi_{0}\in\left(0,\frac{1}{4}\right]$, the interval where there will be no
epidemic is given by
$\mathcal{I}(\zeta,R_{0}^{\text{hom}})=\left(\max\left\\{0,\frac{1}{5}-\tilde{\phi}\right\\},\min\left\\{\frac{1}{4},\frac{1}{5}+\tilde{\phi}\right\\}\right)$,
where

$\displaystyle\tilde{\phi}=\frac{1}{5}\sqrt{\frac{5(1-R_{0}^{\text{hom}})-R_{0}^{\text{hom}}(\zeta-1))}{R_{0}^{\text{hom}}\left(5\zeta(1-R_{0}^{\text{hom}})-(\zeta-1)\right)}}.$

Therefore, if $\phi_{0}\in\mathcal{I}$, then the DFE is stable ∎

Proposition 4

Let $V=\left\\{1,2,3,4,5\right\\}$ be the nodes of a _star-cycle_ network, and
the geometry of the corresponding network is defined by the flux (Eq. 31),
where $M=5,\leavevmode\nobreak\ \text{and}\leavevmode\nobreak\
\phi_{0}\in(0,\frac{1}{4}]$. The infection rate $\beta_{i}$ is defined in (Eq.
8), and we also assume that all the nodes are stable
($R_{0}^{\text{hom}}=\frac{\beta_{0}}{\gamma}<1$) except the node $M=5$. Then
the critical value of $\zeta$ called $\zeta^{\text{crit}}$ is given by

$\zeta^{\text{crit}}=1+\left(\frac{1-R_{0}^{\text{hom}}}{R_{0}^{\text{hom}}}\right)5.$
(47)

Alternatively, the infection threshold of $\beta_{0}$,
$\beta_{0}^{\text{crit}}$ is given by

$\beta_{0}^{\text{crit}}=\gamma\left(\frac{5}{4+\zeta}\right).$ (48)

If $\frac{1}{R_{0}^{\text{hom}}}<\zeta<\zeta^{\text{crit}}$, there is an
interval around the point $\frac{1}{5}$ given by
$\mathcal{I}=\left(\max\left\\{0,\frac{1}{5}-\tilde{\phi}\right\\},\min\left\\{\frac{1}{4},\frac{1}{5}+\tilde{\phi}\right\\}\right)$,
where

$\displaystyle\tilde{\phi}=\frac{1}{5}\sqrt{\frac{\gamma(5\gamma-4\beta_{0}-\zeta\beta_{0})}{\beta_{0}(\gamma+4\zeta\gamma-5\zeta\beta_{0})}}=\frac{1}{5}\sqrt{\frac{5(1-R_{0}^{\text{hom}})-R_{0}^{\text{hom}}(\zeta-1))}{R_{0}^{\text{hom}}\left(5\zeta(1-R_{0}^{\text{hom}})-(\zeta-1)\right)}},$

such that if $\phi_{0}\in\mathcal{I}$, the DFE is stable, and otherwise it is
unstable.

###### Proof.

We perform a local stability analysis to establish the stability criteria for
DFE. The Jacobian matrix at DFE corresponding to the dynamical system of the
$5\times 5$ star-cycle flux matrix (Eq. 31) is symmetric. Since symmetric
matrices have only real eigenvalues, the eigenvalues of the Jacobian matrix
are real. For simplicity, we calculate the eigenvalues using the symbolic
packages of the software Mathematica, which are given by:

$\displaystyle\lambda_{1}=\beta_{0}(1-3\phi_{0})^{2}-\gamma\text{\leavevmode\nobreak\
with\leavevmode\nobreak\ multiplicity\leavevmode\nobreak\ 2}$
$\displaystyle\lambda_{2}=\beta_{0}(1-5\phi_{0})^{2}-\gamma$
$\displaystyle\lambda_{3}=-\gamma+\frac{1}{2}\beta_{0}\left(1+\zeta-2\phi_{0}-8\zeta\phi_{0}+5\phi_{0}^{2}+20\zeta\phi_{0}^{2}+B\right)$
$\displaystyle\lambda_{4}=-\gamma+\frac{1}{2}\beta_{0}\left(1+\zeta-2\phi_{0}-8\zeta\phi_{0}+5\phi_{0}^{2}+20\zeta\phi_{0}^{2}-B\right)$

where,
$B:=\sqrt{(1-2\phi_{0}+5\phi_{0}^{2})^{2}+\zeta^{2}\left(1-8\phi_{0}+20\phi_{0}^{2}\right)^{2}+2\zeta\left(-1+10\phi_{0}-9\phi_{0}^{2}-80\phi_{0}^{3}+100\phi_{0}^{4}\right)}$.

Since all eigenvalues are real and hence B must be nonnegative for
$\phi_{0}\in(0,\frac{1}{4}]$ and any $\zeta$, $B$ is a well-defined and
positive real number. We want to find conditions for which all eigenvalues are
negative (LAS DFE), or at least one eigenvalue is positive (unstable DFE). One
can see that both $\lambda_{1}$, and $\lambda_{2}$ are negative for all
$\phi_{0}\in\left(0,\frac{1}{4}\right]$, since by hypothesis
$\frac{\beta_{0}}{\gamma}<1$. It is thus sufficient to check the critical
points of the remaining eigenvalues. Then we may check $\lambda_{3}$ and
$\lambda_{4}$, which are more complicated. However, a sufficient condition for
stability is the following: since $\lambda_{3}>\lambda_{4}$, we have stability
if the minimum value of $\lambda_{3}$ is negative. It can be easily verified
that $\frac{1}{5}$ is the critical point of the eigenvalues $\lambda_{3}$, and
$\lambda_{4}$, since
$\frac{\partial\lambda_{3}}{\partial\phi_{0}}(\phi_{0}=\frac{1}{5},\beta_{0},\zeta,\gamma)=0$,
and
$\frac{\partial\lambda_{4}}{\partial\phi_{0}}(\phi_{0}=\frac{1}{5},\beta_{0},\zeta,\gamma)=0$.
We find that

$\displaystyle\frac{\partial^{2}\lambda_{3}}{\partial\phi_{0}^{2}}(\phi_{0}=\frac{1}{5},\beta_{0},\zeta,\gamma)=\frac{40\beta_{0}(\zeta-1)^{2}}{4+\zeta}>0$
$\displaystyle\frac{\partial^{2}\lambda_{4}}{\partial\phi_{0}^{2}}(\phi_{0}=\frac{1}{5},\beta_{0},\zeta,\gamma)=\frac{250\zeta\beta_{0}}{4+\zeta}>0$

and so $\lambda_{3}$ and $\lambda_{4}$ have reached to the minimum value when
$\phi_{0}=\frac{1}{5}$. Moreover,

$\displaystyle\lambda_{1}(\phi_{0},\beta_{0},\zeta,\gamma)<0\leavevmode\nobreak\
\leavevmode\nobreak\ \forall\leavevmode\nobreak\ \leavevmode\nobreak\
\phi_{0}\in(0,\frac{1}{4}]$
$\displaystyle\lambda_{2}(\phi_{0},\beta_{0},\zeta,\gamma)<0\leavevmode\nobreak\
\leavevmode\nobreak\ \forall\leavevmode\nobreak\ \leavevmode\nobreak\
\phi_{0}\in(0,\frac{1}{4}]$
$\displaystyle\lambda_{3}(\phi_{0}=\frac{1}{5},\beta_{0},\zeta,\gamma)=\frac{1}{10}\left(4\beta_{0}+\zeta\beta_{0}+\sqrt{(4+\zeta)^{2}}\beta_{0}-10\gamma\right)=\frac{1}{5}\beta_{0}(4+\zeta)-\gamma$
$\displaystyle\lambda_{4}(\phi_{0}=\frac{1}{5},\beta_{0},\zeta,\gamma)=\frac{1}{10}\left(4\beta_{0}+\zeta\beta_{0}-\sqrt{(4+\zeta)^{2}}\beta_{0}-10\gamma\right)=-\gamma<0$

However, if $\frac{1}{5}\beta_{0}(4+\zeta)-\gamma<0$ and that implies

$\displaystyle\beta_{0}<\gamma\left(\frac{5}{4+\zeta}\right):=\beta_{0}^{\text{crit}}\leavevmode\nobreak\
\leavevmode\nobreak\ \text{and\leavevmode\nobreak\
alternatively\leavevmode\nobreak\ we \leavevmode\nobreak\ have,}$
$\displaystyle\zeta<\frac{5\gamma-4\beta_{0}}{\beta_{0}}=\frac{5-4\frac{\beta_{0}}{\gamma}}{\frac{\beta_{0}}{\gamma}}=\frac{5-4R_{0}^{\text{hom}}}{R_{0}^{\text{hom}}}=1+\left(\frac{1-R_{0}^{\text{hom}}}{R_{0}^{\text{hom}}}\right)5:=\zeta^{\text{crit}}.$

In addition, the function
$\phi_{0}\mapsto\lambda_{3}(\phi_{0},\zeta,\beta_{0},\gamma)$ has roots

$\displaystyle\phi_{c1}=\frac{1}{5}+\frac{1}{5}\sqrt{\frac{\gamma(5\gamma-4\beta_{0}-\zeta\beta_{0})}{\beta_{0}(\gamma+4\zeta\gamma-5\zeta\beta_{0})}}=\frac{1}{5}+\frac{1}{5}\sqrt{\frac{5(1-R_{0}^{\text{hom}})-R_{0}^{\text{hom}}(\zeta-1))}{R_{0}^{\text{hom}}\left(5\zeta(1-R_{0}^{\text{hom}})-(\zeta-1)\right)}},$
$\displaystyle\phi_{c2}=\frac{1}{5}-\frac{1}{5}\sqrt{\frac{\gamma(5\gamma-4\beta_{0}-\zeta\beta_{0})}{\beta_{0}(\gamma+4\zeta\gamma-5\zeta\beta_{0})}}=\frac{1}{5}-\frac{1}{5}\sqrt{\frac{5(1-R_{0}^{\text{hom}})-R_{0}^{\text{hom}}(\zeta-1))}{R_{0}^{\text{hom}}\left(5\zeta(1-R_{0}^{\text{hom}})-(\zeta-1)\right)}}.$

Therefore, taking into account the restriction
$\phi_{0}\in\left(0,\frac{1}{4}\right]$, the interval where there will be no
epidemic is given by
$\mathcal{I}(\zeta,R_{0}^{\text{hom}})=\left(\max\left\\{0,\frac{1}{5}-\tilde{\phi}\right\\},\min\left\\{\frac{1}{4},\frac{1}{5}+\tilde{\phi}\right\\}\right)$,
where

$\displaystyle\tilde{\phi}=\frac{1}{5}\sqrt{\frac{5(1-R_{0}^{\text{hom}})-R_{0}^{\text{hom}}(\zeta-1))}{R_{0}^{\text{hom}}\left(5\zeta(1-R_{0}^{\text{hom}})-(\zeta-1)\right)}}.$

Therefore, if $\phi_{0}\in\mathcal{I}$, then the DFE is stable, and,
otherwise, it is unstable. ∎

Theorem 2

Let $V=\left\\{1,2,3,\ldots,M\right\\}$ be the nodes of a star-triangle
network defined by the flux matrix (34), where
$\phi_{0}+\phi_{1}<1,0<\phi_{0}\leq\frac{1}{M-1}$, and $\phi_{1}>0$. We assume
that the infection rate $\beta_{i}=\beta_{0}$ at node $i$ for all
$i=1,2,3,\cdots,M-1$ and $\beta_{M}=\zeta\beta_{0}$. We also suppose that all
the nodes are stable $(R_{0}^{\text{hom}}=\frac{\beta_{0}}{\gamma}<1)$ except
the center node $M$. Then we have a minimum $\zeta^{\text{crit}}$ given by
(Eq. Theorem or 15). If
$\frac{1}{R_{0}^{\text{hom}}}<\zeta<\zeta^{\text{crit}}$, there exists an
interval
$\mathcal{I}=\left(\max\left\\{0,\frac{1}{M}-\tilde{\phi}\right\\},\min\left\\{\frac{1}{M-1},\frac{1}{M}+\tilde{\phi}\right\\}\right)$,
where $\tilde{\phi}$ is given by (Eq. 7), such that if
$\phi_{0}\in\mathcal{I}$, then the DFE is stable, and otherwise it is
unstable.

###### Proof.

We first find the next generation matrix (NGM) $\kappa$ to check the stability
by finding the conditions on the parameters at the DFE. As the next generation
matrix $\kappa=(\kappa_{ij})_{M\times M}$, where
$\kappa_{ij}=\frac{1}{\gamma}\sum_{j=1}^{M}\beta_{j}\phi_{ij}\phi_{kj}\frac{N_{i}}{N_{j}^{p}}$
is given by (Eq. 3). Similarly to what we did in Section 2, we start from the
premise that the populations of each node are equal. As part of heterogeneous
network models, we consider a different infection rate in the central node
with (Eq. 8). Now for any node $j$ we get
$N_{j}^{p}=\sum_{k=1}^{M}\phi_{kj}N_{k}=(\sum_{k=1}^{M}\phi_{kj})N=1.N=N$,
since $\sum_{k=1}^{M}\phi_{kj}=1$, and therefore
$\frac{N_{j}}{N_{j}^{p}}=\frac{N}{N}=1$. Thus the NGM (Eq. 3) becomes (Eq. 4),
and applying the flux matrix (Eq. 34) we have

$\displaystyle\kappa=\frac{\beta_{0}}{\gamma}\begin{pmatrix}p+c_{0}&p+c_{2}&\ldots&p&p&q\\\
p+c_{2}&p+c_{0}&\ldots&p&p&q\\\ \vdots&\vdots&\ddots&\vdots&\vdots\\\
p&p&p+c_{0}&p+c_{2}&p&q\\\ p&p&p+c_{2}&p+c_{0}&p&q\\\
\vdots&\vdots&\vdots&\vdots&\ddots&q\\\ q&q&q&q&q&p+c_{1}\\\ \end{pmatrix},$
(49)

where
$p=\zeta\phi_{0}^{2},q=\phi_{0}(1-\phi_{0}+\zeta\left(1-(M-1)\phi_{0}\right))$,
$c_{0}=\phi_{1}^{2}+(-1+\phi_{0}+\phi_{1})^{2}$, and
$c_{1}=(M-1)\phi_{0}^{2}-\zeta\phi_{0}^{2}+\zeta\left(1-(M-1)\phi_{0}\right)^{2}$,
and $c_{2}=-2\phi_{1}\left(-1+\phi_{0}+\phi_{1}\right)$.

Since the next generation matrix, $\kappa$ in (Eq. 49 ) is symmetric and
positive definite, the eigenvalues of $\kappa$ are all real and positive.
Using Collings (22) method of expansion, the characteristic polynomial of
$\kappa$ is

$(c_{0}+c_{2}-\lambda)^{\frac{M-3}{2}}(c_{0}-c_{2}-\lambda)^{\frac{M-1}{2}}\left(\lambda^{2}-(c_{0}+c_{1}+pM)\lambda+\left(c_{0}c_{1}+c_{1}c_{2}+c_{0}p+(M-1)c_{1}p+c_{2}p+(M-1)(p^{2}-q^{2})\right)\right)=0.$
(50)

For the square matrix $\kappa$, and for the nonzero vector $V$, we have

$\kappa V=R_{0}^{\text{hom}}\lambda V.$ (51)

Therefore, $R_{0}^{\text{hom}}\lambda$ is the eigenvalue of the next
generation matrix $\kappa$. Therefore, eigenvalues of $\kappa$ are
$\frac{\beta_{0}}{\gamma}(c_{0}+c_{2})=R_{0}^{\text{hom}}(1-\phi_{0})^{2}$
with multiplicity $\frac{M-3}{2}$, also the eigenvalue of $\kappa$ is
$\frac{\beta_{0}}{\gamma}(c_{0}-c_{2})=R_{0}^{\text{hom}}(-1+\phi_{0}+2\phi_{1})^{2}$
with multiplicity $\frac{M-1}{2}$. Finally, the other two eigenvalues are
determined from the quadratic term

$p(\lambda)=\lambda^{2}-(c_{0}+c_{1}+c_{2}+pM)\lambda+\left(c_{1}c_{0}+c_{1}c_{2}+c_{0}p+p(M-1)c_{1}+c_{2}p+(M-1)(p^{2}-q^{2})\right).$
(52)

To determine the stability of the NGM $\kappa$ in (Eq. 49), we need to find
the conditions on parameters such that all of its eigenvalues are within the
unit circle, so that $R_{0}<1$. The eigenvalues automatically meet the
condition

$R_{0}^{\text{hom}}(1-\phi_{0})^{2}<1\leavevmode\nobreak\
\text{and}\leavevmode\nobreak\
R_{0}^{\text{hom}}(-1+\phi_{0}+2\phi_{1})^{2}<1$ (53)

since $R_{0}^{\text{hom}}<1$, and $1-(M-1)\phi_{0}>0$, and
$1-\phi_{0}-(M-2)\phi_{1}>0$ for all $M$. Even the condition (Eq. 53) holds
when $R_{0}^{\text{hom}}>1$, if $\phi_{0}$ near enough to $1$.

The quadratic (Eq. 52) gives the other two eigenvalues lying within the unit
circle, which can be tested by applying the Jury conditions (see Murray (?),
page $507$). Using the Jury condition, a quadratic equation
$P(\lambda)=\lambda^{2}+a_{1}\lambda+a_{0}=0$, these conditions are
$P(1)=1+a_{1}+a_{0}>0$, $P(-1)=1-a_{1}+a_{0}>0$, and $P(0)=a_{0}<1$.

As $R_{0}^{\text{hom}}\lambda:=\lambda^{\prime}$ is the eigenvalue of
$\kappa$,substituting $\lambda=\frac{\lambda^{\prime}}{R_{0}^{\text{hom}}}$ in
(Eq. 52) we have

$P(\lambda^{\prime})=\frac{1}{(R_{0}^{\text{hom}})^{2}}[(\lambda^{\prime})^{2}-a_{1}R_{0}^{\text{hom}}\lambda^{\prime}+a_{0}(R_{0}^{\text{hom}})^{2},$
(54)

where $a_{0}:=\left(c_{1}c_{0}+p(M-1)c_{1}+pc_{0}+(M-1)(p^{2}-q^{2})\right)$,
and $a_{1}:=(c_{0}+c_{1}+pM)$. Applying the Jury conditions to (Eq. 54), we
obtain the following inequalities:

$\displaystyle\left\\{R_{0}^{\text{hom}}[\zeta\left(R_{0}^{\text{hom}}(M\phi_{0}-1)^{2}-(M-1)\phi_{0}(M\phi_{0}-2)-1\right)+\phi_{0}(2-M\phi_{0})-1]+1\right\\}>0,$
(55)

$\displaystyle\left\\{R_{0}^{\text{hom}}[\zeta\left(R_{0}^{\text{hom}}(M\phi_{0}-1)^{2}+(M-1)\phi_{0}(M\phi_{0}-2)+1\right)-\phi_{0}(2-M\phi_{0})+1]+1\right\\}>0,$
(56)

$\displaystyle(R_{0}^{\text{hom}})^{2}\zeta(M\phi_{0}-1)^{2}<1.$ (57)

The Jury conditions imply that (Eqs. 55-57) are identical to (Eqs. A.2-A.4) in
Stolerman et al. (?). Therefore the rest of the proof follows.

∎

Proposition 5

Let $V=\left\\{1,2,3,4,\ldots,11\right\\}$ be the nodes of a star-triangle
network defined by the symmetric flux matrix (Eq. 35) and $M=11$ be an odd
number. We assume that the infection rate $\beta_{i}=\beta_{0}$ at node $i$
for all $i=1,2,3,\cdots,10$ and $\beta_{11}=\zeta\beta_{0}$. We also suppose
that all the nodes are stable
$(R_{0}^{\text{hom}}=\frac{\beta_{0}}{\gamma}<1)$ except the center node $M$.
Then the critical value of $\zeta$ called $\zeta^{\text{crit}}$ is given by

$\zeta^{\text{crit}}=1+\left(\frac{1-R_{0}^{\text{hom}}}{R_{0}^{\text{hom}}}\right)11.$
(58)

Alternatively, the infection threshold of $\beta_{0}$,
$\beta_{0}^{\text{crit}}$ is given by

$\beta_{0}^{\text{crit}}=\gamma\left(\frac{11}{10+\zeta}\right).$ (59)

If $\frac{1}{R_{0}^{\text{hom}}}<\zeta<\zeta^{\text{crit}}$, there exists an
interval
$\mathcal{I}=\left(\max\left\\{0,\frac{1}{11}-\tilde{\phi}\right\\},\min\left\\{\frac{1}{10},\frac{1}{11}+\tilde{\phi}\right\\}\right)$,
where $\tilde{\phi}$ is given by (Eq. 7), such that if
$\phi_{0}\in\mathcal{I}$, then the DFE is stable, and otherwise it is
unstable.

###### Proof.

To establish the stability criteria for DFE in this case, we may perform a
local stability analysis of the DFE. The Jacobian matrix at the DFE
corresponding to the dynamical system of the $11\times 11$ _star-triangle_
network and the flux matrix (Eq. 35) is symmetric. Since symmetric matrices
have only real eigenvalues, the eigenvalues of the above Jacobian matrix are
real. For simplicity, we calculate the eigenvalues using the symbolic packages
of the software Mathematica which are given by

$\displaystyle\lambda_{1}=-\gamma+\beta_{0}(-1-\phi_{6})^{2}\text{\leavevmode\nobreak\
with\leavevmode\nobreak\ multiplicity\leavevmode\nobreak\ 4}$
$\displaystyle\lambda_{2}=-\gamma+\beta_{0}(-1+2\phi_{1}+\phi_{6})^{2}$
$\displaystyle\lambda_{3}=-\gamma+\beta_{0}(-1+2\phi_{2}+\phi_{6})^{2}$
$\displaystyle\lambda_{4}=-\gamma+\beta_{0}(-1+2\phi_{3}+\phi_{6})^{2}$
$\displaystyle\lambda_{5}=-\gamma+\beta_{0}(-1+2\phi_{4}+\phi_{6})^{2}$
$\displaystyle\lambda_{6}=-\gamma+\beta_{0}(-1+2\phi_{5}+\phi_{6})^{2}$
$\displaystyle\lambda_{7}=-\gamma+\frac{1}{2}\beta_{0}\left(1+\zeta-2\phi_{6}-20\zeta\phi_{6}+11\phi_{6}^{2}+110\zeta\phi_{6}^{2}+B\right)$
$\displaystyle\lambda_{8}=-\gamma+\frac{1}{2}\beta_{0}\left(1+\zeta-2\phi_{6}-20\zeta\phi_{6}+11\phi_{6}^{2}+110\zeta\phi_{6}^{2}-B\right)$

where
$B:=\sqrt{(1-2\phi_{6}+11\phi_{6}^{2})^{2}+\zeta^{2}\left(1-20\phi_{6}+110\phi_{6}^{2}\right)^{2}+2\zeta\left(-1+22\phi_{6}-81\phi_{6}^{2}-440\phi_{6}^{3}+1210\phi_{6}^{4}\right)}$.
Since all eigenvalues are real and B must be nonnegative for
$\phi_{0}\in(0,\frac{1}{10}]$ and any $\zeta$, $B$ is a positive real number.
We want to find conditions for which all eigenvalues are negative (LAS DFE),
or at least one eigenvalue is positive (unstable DFE). Because $\phi_{ij}>0$,
we can see that all $\lambda_{i}$’s are negative for all $i=1,2,\ldots,6$. It
is thus sufficient to check the critical points of $\lambda_{7}$ and
$\lambda_{8}$. However, a sufficient condition for their stability is the
following: since $\lambda_{7}>\lambda_{8}$, it is thus sufficient to check the
stability of $\lambda_{7}$. It can be easily verified that $\frac{1}{11}$ is
the critical point of both eigenvalues $\lambda_{7}$ and $\lambda_{8}$, since
$\frac{\partial\lambda_{7}}{\partial\phi_{6}}(\phi_{6}=\frac{1}{11},\beta_{6},\zeta,\gamma)=0$,
and
$\frac{\partial\lambda_{8}}{\partial\phi_{6}}(\phi_{6}=\frac{1}{11},\beta_{0},\zeta,\gamma)=0$.
We can find that

$\displaystyle\frac{\partial^{2}\lambda_{7}}{\partial\phi_{6}^{2}}(\phi_{6}=\frac{1}{11},\beta_{0},\zeta,\gamma)=\frac{220\beta_{0}(\zeta-1)^{2}}{10+\zeta}>0$
$\displaystyle\frac{\partial^{2}\lambda_{8}}{\partial\phi_{6}^{2}}(\phi_{6}=\frac{1}{11},\beta_{0},\zeta,\gamma)=\frac{2662\zeta\beta_{0}}{10+\zeta}>0$

and so $\lambda_{7}$ and $\lambda_{8}$ have reached to their minimum value
when $\phi_{6}=\frac{1}{11}$. Moreover,

$\displaystyle\lambda_{i}(\phi_{6},\beta_{0},\zeta,\gamma)<0\leavevmode\nobreak\
\leavevmode\nobreak\ \text{for \leavevmode\nobreak\ the \leavevmode\nobreak\
choice\leavevmode\nobreak\ of}\leavevmode\nobreak\ \leavevmode\nobreak\
\phi_{i},\text{where}\leavevmode\nobreak\ i=1,2,\ldots,6$
$\displaystyle\lambda_{7}(\phi_{6}=\frac{1}{11},\beta_{0},\zeta,\gamma)=\frac{1}{22}\left(10\beta_{0}+\zeta\beta_{0}+\sqrt{(10+\zeta)^{2}}\beta_{0}-22\gamma\right)=\frac{1}{11}\beta_{0}(10+\zeta)-\gamma$
$\displaystyle\lambda_{8}(\phi_{6}=\frac{1}{11},\beta_{0},\zeta,\gamma)=\frac{1}{10}\left(10\beta_{0}+\zeta\beta_{0}-\sqrt{(4+\zeta)^{2}}\beta_{0}-10\gamma\right)=-\gamma<0$

However, if $\frac{1}{11}\beta_{0}(10+\zeta)-\gamma<0$ that implies

$\displaystyle\beta_{0}<\gamma\left(\frac{11}{10+\zeta}\right):=\beta_{0}^{\text{crit}}\leavevmode\nobreak\
\leavevmode\nobreak\ \text{and}$
$\displaystyle\zeta<\frac{11\gamma-10\beta_{0}}{\beta_{0}}=\frac{11-10\frac{\beta_{0}}{\gamma}}{\frac{\beta_{0}}{\gamma}}=\frac{11-10R_{0}^{\text{hom}}}{R_{0}^{\text{hom}}}=1+\left(\frac{1-R_{0}^{\text{hom}}}{R_{0}^{\text{hom}}}\right)11:=\zeta^{\text{crit}}.$

In addition, the function
$\phi_{6}\mapsto\lambda_{7}(\phi_{6},\zeta,\beta_{0},\gamma)$ has roots

$\displaystyle\phi_{c1}=\frac{1}{11}+\frac{1}{11}\sqrt{\frac{\gamma(11\gamma-10\beta_{0}-\zeta\beta_{0})}{\beta_{0}(\gamma+10\zeta\gamma-11\zeta\beta_{0})}}=\frac{1}{11}+\frac{1}{11}\sqrt{\frac{11(1-R_{0}^{\text{hom}})-R_{0}^{\text{hom}}(\zeta-1))}{R_{0}^{\text{hom}}\left(11\zeta(1-R_{0}^{\text{hom}})-(\zeta-1)\right)}},$
$\displaystyle\phi_{c2}=\frac{1}{11}-\frac{1}{11}\sqrt{\frac{\gamma(11\gamma-10\beta_{0}-\zeta\beta_{0})}{\beta_{0}(\gamma+10\zeta\gamma-11\zeta\beta_{0})}}=\frac{1}{11}-\frac{1}{11}\sqrt{\frac{11(1-R_{0}^{\text{hom}})-R_{0}^{\text{hom}}(\zeta-1))}{R_{0}^{\text{hom}}\left(11\zeta(1-R_{0}^{\text{hom}})-(\zeta-1)\right)}}.$

Therefore, according to the restriction
$\phi_{6}\in\left(0,\frac{1}{10}\right]$, the interval where there will be no
epidemic is given by
$\mathcal{I}(\zeta,R_{0}^{\text{hom}})=\left(\max\left\\{0,\frac{1}{11}-\tilde{\phi}\right\\},\min\left\\{\frac{1}{10},\frac{1}{11}+\tilde{\phi}\right\\}\right)$,
where

$\displaystyle\tilde{\phi}=\frac{1}{11}\sqrt{\frac{11(1-R_{0}^{\text{hom}})-R_{0}^{\text{hom}}(\zeta-1))}{R_{0}^{\text{hom}}\left(5\zeta(1-R_{0}^{\text{hom}})-(\zeta-1)\right)}}.$

Therefore, if $\phi_{6}\in\mathcal{I}$, then the DFE is stable, and otherwise,
it is unstable. ∎

Theorem 3 Let $V=\left\\{1,2,3,\ldots,M\right\\}$ be the nodes of a _star-
background_ networks defined by the flux matrix (Eq. 38), where
$\phi_{0}+(M-2)\phi_{1}<1,0<\phi_{0}\leq\frac{1}{M-1}$, and $\phi_{1}>0$. We
assume that the infection rate $\beta_{i}=\beta_{0}$ at node $i$ for all
$i=1,2,3,\cdots,M-1$ and $\beta_{M}=\zeta\beta_{0}$. We also suppose that all
the nodes are stable ($R_{0}^{\text{hom}}=\frac{\beta_{0}}{\gamma}<1$) except
the center $M$. Then we have a minimum $\zeta^{\text{crit}}$ given by (Eq.
Theorem or 15). If $\frac{1}{R_{0}^{\text{hom}}}<\zeta<\zeta^{\text{crit}}$,
there exists an interval
$\mathcal{I}=\left(\max\left\\{0,\frac{1}{M}-\tilde{\phi}\right\\},\min\left\\{\frac{1}{M-1},\frac{1}{M}+\tilde{\phi}\right\\}\right)$,
where $\tilde{\phi}$ is given by (Eq. 7), such that if
$\phi_{0}\in\mathcal{I}$, then the DFE is stable, and otherwise it is
unstable. Finally, if $\zeta<\frac{1}{R_{0}^{\text{hom}}}$, i.e.,
$\beta_{M}<\gamma$, then the DFE is always stable.

###### Proof.

We first find the next generation matrix (NGM) $\kappa$ to check the stability
by finding the conditions on the parameters at the DFE. As the next generation
matrix $\kappa=(\kappa_{ij})_{M\times M}$, where
$\kappa_{ij}=\frac{1}{\gamma}\sum_{j=1}^{M}\beta_{j}\phi_{ij}\phi_{kj}\frac{N_{i}}{N_{j}^{p}}$
is given by (Eq. 3). Similar to Section 2, we start from the premise that the
populations of each node are equal. i.e. $N_{i}=N>0$ for all $i\in V$. As part
of heterogeneous network models, we consider a different infection rate in the
central node with the (Eq. 8).

Now for any node $j$ we get
$N_{j}^{p}=\sum_{k=1}^{M}\phi_{kj}N_{k}=(\sum_{k=1}^{M}\phi_{kj})N=1.N=N$,
since $\sum_{k=1}^{M}\phi_{kj}=1$, and therefore
$\frac{N_{j}}{N_{j}^{p}}=\frac{N}{N}=1$. Thus the NGM (Eq. 3) becomes (Eq. 4),
and applying the flux matrix (Eq. 38) we obtain

$\displaystyle\kappa=\frac{\beta_{0}}{\gamma}\begin{pmatrix}p+c_{0}&p&p&\ldots&q\\\
p&p+c_{0}&p&\ldots&q\\\ p&p&p+c_{0}&\ldots&q\\\
\vdots&\vdots&\vdots&\ddots&\vdots\\\ q&q&q&\ldots&p+c_{1}\\\
\end{pmatrix}=\frac{\beta_{0}}{\gamma}\begin{pmatrix}p&p&p&\ldots&q\\\
p&p&p&\ldots&q\\\ p&p&p&\ldots&q\\\ \vdots&\vdots&\vdots&\ddots&\vdots\\\
q&q&q&\ldots&p\\\
\end{pmatrix}+\frac{\beta_{0}}{\gamma}\begin{pmatrix}c_{0}\\\ c_{0}\\\
c_{0}\\\ \vdots\\\ c_{1}\\\ \end{pmatrix}$ (60)

which is equal to

$\kappa=\frac{\beta_{0}}{\gamma}(\mathcal{P}+\mathcal{D})$ (61)

where $p=\zeta\phi_{0}^{2}+\phi_{1}(2-2\phi_{0}+\phi_{1}-M\phi_{1})$,
$q=\phi_{0}(1-\phi_{0}+\zeta\left(1-(M-1)\phi_{0}\right))$,
$c_{0}=(1-\phi_{0}+(M-1)\phi_{1})^{2}$, and
$c_{1}=(m-1)\phi_{0}^{2}-\zeta\phi_{0}^{2}+\zeta\left(\phi_{0}(m-1)-1\right)^{2}+\phi_{1}\left(2\phi_{0}+(m-1)\phi_{1}-2\right)$.
We find the entries of the NGM in the following way: For
$u\in{1,2,3,4,5,...,M-1}$, and $v\in{1,2,3,4,5,...,M-1}$, and $u=v$, the
diagonal elements of the NGM $\kappa$ are

$\displaystyle\kappa_{uu}=\sum_{j=1}^{M}\beta_{j}\phi_{uj}\phi_{uj}=\sum_{j=1}^{M}\beta_{j}\phi_{uj}^{2}$
(62)
$\displaystyle=\beta_{1}\phi_{u1}^{2}+\beta_{2}\phi_{u2}^{2}+\beta_{3}\phi_{u3}^{2}+\ldots+\beta_{u}\phi_{uu}^{2}+\ldots+\beta_{M}\phi_{uM}^{2}$
$\displaystyle=\beta_{0}.\phi_{1}^{2}+\beta_{0}.\phi_{1}^{2}+\beta_{0}.\phi_{1}^{2}+\ldots+\beta_{0}(1-\phi_{0}-(M-2)\phi_{1})^{2}+\ldots+\zeta\beta_{0}\phi_{0}^{2}$
$\displaystyle=\beta_{0}\left(\zeta\phi_{0}^{2}+\phi_{1}^{2}(M-2)+(1-\phi_{0}-(M-2)\phi_{1})^{2}\right)$
$\displaystyle=\beta_{0}\left(p+c_{0}\right),$

and, if $u=v=M$,

$\displaystyle\kappa_{MM}=\sum_{j=1}^{M}\beta_{j}\phi_{Mj}^{2}\ldots$
$\displaystyle=\beta_{1}\phi_{M1}^{2}+\beta_{2}\phi_{M2}^{2}+\beta_{3}\phi_{M3}^{2}+...+\beta_{M}\phi_{MM}^{2}$
$\displaystyle=\beta_{0}\phi_{0}^{2}+\beta_{0}\phi_{0}^{2}+\beta_{0}\phi_{0}^{2}+\ldots+\zeta\beta_{0}(1-(M-1)\phi_{0})^{2}$
$\displaystyle=\zeta\beta_{0}(1-(M-1)\phi_{0})^{2}+(M-1)\beta_{0}\phi_{0}^{2}$
$\displaystyle=(M-1)\beta_{0}\phi_{0}^{2}+\zeta\beta_{0}(1-(M-1)\phi_{0})^{2}$
$\displaystyle=\beta_{0}\left(M\phi_{0}^{2}-\phi_{0}^{2}+\zeta(1-2(M-1)\phi_{0}+(M-1)^{2}\phi_{0}^{2}\right)$
$\displaystyle=\beta_{0}\left(M\phi_{0}^{2}-\phi_{0}^{2}+\zeta-2\zeta
M\phi_{0}+2\zeta\phi_{0}+\zeta M^{2}\phi_{0}^{2}-2\zeta
M\phi_{0}^{2}+\zeta\phi_{0}^{2}\right)$
$\displaystyle=\beta_{0}\left(\zeta\phi_{0}^{2}+M\phi_{0}^{2}-\phi_{0}^{2}+\zeta-2\zeta
M\phi_{0}+2\zeta\phi_{0}+\zeta M^{2}\phi_{0}^{2}-2\zeta M\phi_{0}^{2}\right)$
$\displaystyle=\beta_{0}\left(\zeta\phi_{0}^{2}+\zeta-2\zeta(M-1)\phi_{0}+\phi_{0}^{2}\left(M-1-2M\zeta+M^{2}\zeta\right)\right)$
$\displaystyle=\beta_{0}\left(p+c_{1}\right).$

Finally, for any non-diagonal elements of the NGM $\kappa$, when
$u\in{1,2,3,4,5,...,M-1}$, and $v\in{1,2,3,4,5,...,M-1}$, and $u\not=v$ we
have

$\displaystyle\kappa_{uv}=\sum_{j=1}^{M}\beta_{j}\phi_{uj}\phi_{vj}$
$\displaystyle=\beta_{1}\phi_{u1}\phi_{v1}+\beta_{2}\phi_{u2}\phi_{v2}+\ldots+\beta_{u}\phi_{uu}\phi_{vu}.+\beta_{v}\phi_{uv}\phi_{vv}...+\beta_{M}\phi_{uM}\phi_{vM}$