Adaptive Susceptibility and Heterogeneity in Contagion Models on Networks

TL;DR

This paper introduces a generalized network epidemic model, SIRI, accounting for adaptive susceptibility and heterogeneity, revealing complex dynamics like bistability and resurgent epidemics, with implications for control strategies.

Contribution

It develops a novel SIRI model incorporating adaptive susceptibility and heterogeneity, providing necessary and sufficient conditions for various epidemic regimes, including bistability.

Findings

Identification of four dynamic regimes: infection-free, epidemic, endemic, bistable.

Bistable regime can cause rapid epidemic resurgence after apparent control.

Model aids in designing effective containment strategies.

Abstract

Contagious processes, such as spread of infectious diseases, social behaviors, or computer viruses, affect biological, social, and technological systems. Epidemic models for large populations and finite populations on networks have been used to understand and control both transient and steady-state behaviors. Typically it is assumed that after recovery from an infection, every agent will either return to its original susceptible state or acquire full immunity to reinfection. We study the network SIRI (Susceptible-Infected-Recovered-Infected) model, an epidemic model for the spread of contagious processes on a network of heterogeneous agents that can adapt their susceptibility to reinfection. The model generalizes existing models to accommodate realistic conditions in which agents acquire partial or compromised immunity after first exposure to an infection. We prove necessary and sufficient conditions on model parameters and network structure that distinguish four dynamic regimes: infection-free, epidemic, endemic, and bistable. For the bistable regime, which is not accounted for in traditional models, we show how there can be a rapid resurgent epidemic after what looks like convergence to an infection-free population. We use the model and its predictive capability to show how control strategies can be designed to mitigate problematic contagious behaviors.

Tables

Table 1. TABLE I: Network SIRI model cases

Case	Parameter Value	Equivalent Model
1	$D=\overline{\mathrm{𝟎}}$	SI
2	$\hat{B}=\overline{\mathrm{𝟎}}$	SIR
3	$B=\hat{B}$	SIS
4	$B\succ \hat{B}\succ \overline{\mathrm{𝟎}}$	Partial Immunity
5	$\hat{B}\succ B\succ \overline{\mathrm{𝟎}}$	Compromised Immunity
6	Otherwise	Mixed Immunity

Table 2. TABLE II: Applicability of Theorem  1 statements                                                                                                                                    to the Network SIRI model cases

Case	Model	A	B	C	D	E	F
1	SI
2	SIR	✓			✓		✓
3	SIS	✓	✓	✓
4	Partial	✓	✓		✓	✓	✓
5	Comprom.	✓	✓		✓	✓	✓
6	Mixed	✓	✓		✓	✓	✓

Table 3. TABLE III: Behavioral Regimes of the Network SIRI Cases

Case	Model	Inf.-Free	Endemic	Epidemic	Bistable
1	SI		✓
2	SIR	✓		✓
3	SIS	✓	✓
4	Partial	✓	✓	✓
5	Comprom.	✓	✓		✓
6	Mixed	✓	✓	✓	✓

Equations

$${\lambda }_{max}(K_{min})<{\lambda }_{max}(K)<{\lambda }_{max}(K_{max}).$$

$${\dot{p}}_j^S=$$

$${\dot{p}}_j^I=$$

$${\dot{p}}_j^R=$$

$${\dot{p}}_j^S=$$

$${\dot{p}}_j^I=$$

$${\dot{p}}^S$$

$${\dot{p}}^I$$

$$B^{\ast }(p^S)=(\mathbb{I}-P^S)\hat{B}+P^{S}B$$

$$B$$

$$\hat{B}$$

$$D$$

$${\bar{B}}_{max}=[\max ({\beta }_{jk},{\hat{\beta }}_{jk})],{\bar{B}}_{min}=[\min ({\beta }_{jk},{\hat{\beta }}_{jk})].$$

$$p_j^{I\ast }=\frac{{\sum }_{k=1}^N{\hat{\beta }}_{jk}{p}_k^{I\ast }}{{\delta }_j+{\sum }_{k=1}^N{\hat{\beta }}_{jk}{p}_k^{I\ast }}.$$

$$0=(\hat{B}-D-P^{I\ast }\hat{B})p^{I\ast }=(\hat{B}-D-\text{diag}(\hat{B}{p}^{I\ast }))p^{I\ast }.$$

$${\dot{p}}^I=(B-D)p^I-P^{I}B{p}^I.$$

$${\dot{p}}^S$$

$${\dot{p}}^I$$

$$B^{\ast }(p)=(\mathbb{I}-P)\hat{B}+PB,$$

$$R(p)=\rho (B^{\ast }(p)D^{-1}).$$

$$R_0=\rho (B{D}^{-1}),R_1=\rho (\hat{B}{D}^{-1}),$$

$$R_{max}=\underset{0⪯p⪯1}{\max }\rho (B^{\ast }(p)D^{-1}),R_{min}=\underset{0⪯p⪯1}{\min }\rho (B^{\ast }(p)D^{-1}).$$

$${\bar{R}}_{min}\le R_{min}\le {\lambda }_{max}(B^{\ast }(p^S)D^{-1})\le R_{max}\le {\bar{R}}_{max}$$

$$J_{EE}=\left[\begin{array}{cc}-\text{diag}(B{p}^{I\ast })& \bar{0}\\ \text{diag}((B-\hat{B})p^{I\ast })& J_a\end{array}\right],$$

$$J_a{p}^{I\ast }=-\text{diag}(\hat{B}{p}^{I\ast })p^{I\ast }\ll 0,$$

$$J_{\mathcal{M}}(x)=\left[\begin{array}{cc}\bar{0}& -P^{S\ast }B\\ \bar{0}& J_T(p^{S\ast })\end{array}\right],$$

$$\{\begin{array}{ll}{p}_j^{\prime }>-\frac{{\sum }_{j=1}^N{\hat{\beta }}_{jk}}{{\sum }_{j=1}^N({\beta }_{jk}-{\hat{\beta }}_{jk})}& \mathrm{if}{\sum }_{\mathrm{j}=1}^{\mathrm{N}}({\beta }_{\mathrm{jk}}-{\hat{\beta }}_{\mathrm{jk}})\ge 0,\\ p_j^{\prime }<\frac{{\sum }_{j=1}^N{\hat{\beta }}_{jk}}{{\sum }_{j=1}^N({\hat{\beta }}_{jk}-{\beta }_{jk})}& \mathrm{if}{\sum }_{\mathrm{j}=1}^{\mathrm{N}}({\beta }_{\mathrm{jk}}-{\hat{\beta }}_{\mathrm{jk}})\le 0.\end{array}$$

$$\frac{\partial \Lambda }{\partial p_j}(p)=w_j\sum _{k=1}^N({\beta }_{jk}-{\hat{\beta }}_{jk})v_k$$

$$\nabla \Lambda (p)=\text{diag}(w)(B-\hat{B})v.$$

$$w^T\frac{\partial J_T}{\partial p_j}v+w^T{J}_T\frac{\partial v}{\partial p_j}=w^{T}v\frac{\partial \Lambda }{\partial p_j}+w^T\Lambda \frac{\partial v}{\partial p_j},$$

$$z_{kl}=\{\begin{array}{ll}{\beta }_{kl}-{\hat{\beta }}_{kl}& k=j\text{and}(k,l)\in \mathcal{E},\\ 0& \text{otherwise}.\end{array}$$

$$\frac{\partial \Lambda }{\partial p_j}(p)=w^T\frac{\partial J_T}{\partial p_j}(p)v=w_j\sum _{k=1}^N({\beta }_{jk}-{\hat{\beta }}_{jk})v_k.$$

$$J_T(p)v_0$$

$$=(\hat{B}-D)v_0+\text{diag}(p)(B-\hat{B})v_0$$

Full text

I Introduction

Contagious processes affect biological, social, and technological network systems. In biology, a concern is the spread of disease across a population [1]. In social networks, the spread of information, behaviors, or cultural norms has important implications for decisions about politics, the environment, health care, etc. In many contexts a spike in the “infected” population is detrimental, e.g., in the spread of misinformation [2]. In other cases, the spike can be crucial, e.g., in the spread of safety instructions in an emergency [3]. In technological networks, like computer networks and mobile sensor networks, the spread of a virus can lead to disruptions, while the spread of information on a changing environment can be critical to a successful mission. To understand and control the dynamics of contagion, models with sufficiently good predictive capability are warranted.

Epidemic models have been successfully used to study contagious processes in a large number of systems, ranging from the spread of infectious diseases on populations [4, 5] and memes on social networks [6] to the evolution of riots [7] and power grid failures [8]. The wide applicability of epidemic models has led to an increase in recent years in the number of studies in the physics and control communities focusing on theoretical epidemic models for the propagation of contagious processes on networks [9, 10, 11, 12, 13]. Studies are typically based on the SIS (Susceptible-Infected-Susceptible) model [14, 15, 16, 11, 12, 13], in which every recovered individual experiences no change to its susceptibility to the infection after recovery, or on the SIR (Susceptible-Infected-Recovered) model [17, 13], in which recovered individuals gain full immunity to the infection. The SIS model captures endemic behaviors in which there is a constant fraction of the population that is infected at steady-state, while the SIR model captures epidemic behaviors in which there is a rapid increase in the number of infected individuals that eventually subsides.

Although the SIR and SIS models are useful, they fall short in accounting for what can happen in many real-world systems when agents adapt their susceptibility in more general ways after their first exposure to the infection. For example, in the case of infectious diseases, the susceptibility of individuals to the infection can decrease after a first exposure, resulting in partial immunity as in the case of influenza [18]. Or susceptibility can increase, resulting in compromised immunity as in the case of tuberculosis in particular populations [19]. In the spread of social behaviors, the susceptibility of individuals to the “infection” might decrease (increase) as a result of a negative (positive) past experience that decreases (increases) the propensity of an individual to engage in the behavior.

For example, in psychology there is a well-established “inoculation theory” [20, 21, 2], which has shown that an effective way to combat the spread of misinformation is to pre-expose people to misinformation so that they develop at least partial cognitive immunity to subsequent contact with misinformation. Examples also abound in social animal behavior. For example, desert harvester ants regulate foraging for seeds by means of a contagious process in which successful ants returning to the nest motivate available ants in the nest to go out and forage [22]. Resilience of foraging rates to temperature and humidity changes outside the nest has been attributed to adaptation of susceptibility to the “infection” by available ants that have already been exposed to outside conditions [23].

In technological systems, such as systems comprised of interconnected mechanical parts, susceptibility to cascading failures can increase the second time around when parts become worn or compromised the first time. When there is the opportunity for learning and design, such as in a network of autonomous robots, protocols can be designed so that agents modulate their response to an “infected” agent based on what they have learned from previous interactions.

In the present paper, we analyze the role of adaptive susceptibility in the spread of a contagious process over a network of heterogeneous agents using the network SIRI (Susceptible-Infected-Recovered-Infected) epidemic model. The network SIRI model describes susceptible agents that become infected from contact with infected agents, infected agents that recover, and recovered agents that become reinfected from contact with infected agents. Susceptibility is adaptive when the rates of infection and reinfection differ. An agent acquires partial immunity (compromised immunity) to its neighbor, if after recovering from a first-time infection, it experiences a decrease (increase) in susceptibility to that neighbor. Every agent may have a different rate of recovery and different rates of susceptibility to infection and reinfection to each of its neighbors. SIS and SIR are special cases of SIRI.

To the best of our knowledge we provide the first rigorous and comprehensive analysis of the network SIRI model. Models that consider reinfection usually only consider the case of partial immunity [24, 25, 26]. In [9] the authors studied a discrete-time mean-field model with global recovery, infection, and reinfection rates, and showed through numerical simulations that when the reinfection rate is larger than the infection rate, there is an abrupt transition from an infection-free to an endemic steady state. We formalize this observation by proving new results on the existence of a bistable regime in which a critical manifold of initial conditions separates solutions for which the infection dies out from solutions for which the infection spreads.

Our analysis of the network SIRI model generalizes our novel results on the SIRI model of a well-mixed population [27]. Our contributions are as follows. First, we introduce the network SIRI model over strongly connected digraphs and present a rigorous stability analysis. We show there can exist only a set of non-isolated infection-free equilibria (IFE) and a stable isolated endemic equilibrium (EE), and we prove conditions on the graph structure and system parameters that determine the stability of the IFE. Second, we prove that the model exhibits the same four distinct behavioral regimes observed in the well-mixed SIRI model: infection-free, epidemic, endemic, and bistable. We show how the four behavioral regimes are characterized by four numbers that generalize, to the network setting, the two reproduction numbers of [27] and generalize previous results for the network SIS and SIR models in [14, 15, 28, 17, 13]. Third, we prove features of the geometry of solutions near the IFE in the epidemic and bistable regimes that dictate both transient and steady-state possibilities. For the bistable regime, which is not accounted for by the SIS and SIR models, we show how solutions can exhibit a resurgent epidemic in which an initial infection appears to die out for an arbitrarily long period of time before abruptly resurging to an endemic steady state. Finally, we show how our results can be used to design control strategies that guarantee the eradication or spread of the infection through a network of heterogeneous agents with adaptive susceptibility.

In Section II we present notation and well-known results used in the paper. In Section III we introduce the network SIRI model and its classification into six cases. In Section IV we study the equilibria and define new notions of reproduction numbers. In Section V we analyze stability, and in Section VI we prove our main result on conditions for the four behavioral regimes. In Section VII we examine the geometry and behavior of solutions in the epidemic and bistable regimes. We apply our theory to control laws that guarantee desired steady-state behavior in Section VIII. We make final remarks in Section IX.

II Mathematical Preliminaries

II-A Notation

We denote the $j$-th entry of $x\in\mathbb{R}^{N}$ as $x_{j}$, and the $(j,k)$-th entry of $M\in\mathbb{R}^{N\times N}$ as $m_{jk}$. We define $e^{j}\in\mathbb{R}^{N}$, $j=1,\dots,N$ as the standard basis vectors. We define $\mathbf{0}\in\mathbb{R}^{N}$ as the zero vector, $\mathbf{1}\in\mathbb{R}^{N}$ as the vector with every entry 1, $\bar{\mathbf{0}}\in\mathbb{R}^{N\times N}$ as the zero square matrix, and $\mathbb{I}\in\mathbb{R}^{N\times N}$ as the identity matrix. We let $\textrm{diag}(x)\in\mathbb{R}^{N\times N}$ be the diagonal matrix with entries given by the entries of $x\in\mathbb{R}^{N}$.

For any vectors $x,y\in\mathbb{R}^{N}$, we write $x\gg y$ if $x_{j}>y_{j}$ for all $j$, $x\succ y$ if $x_{j}\geq y_{j}$ for all $j$, but $x\neq y$, and $x\succeq y$ if $x_{j}\geq y_{j}$ for all $j$. Similarly, for any two matrices $M,Q\in\mathbb{R}^{N\times N}$ we write $M\gg Q$ if $m_{jk}>q_{jk}$, $M\succ Q$ if $m_{jk}\geq q_{jk}$, for any $j,k$, but $M\neq Q$, and $M\succeq Q$ if $m_{jk}\geq q_{jk}$ for any $j,k$.

A square matrix $M$ is Hurwitz (stable) if it has no eigenvalue with positive or zero real part, and it is unstable if at least one of its eigenvalues has positive real part. A real square matrix $M$ is Metzler if $m_{jk}\geq 0$ for $j\neq k$. We denote the spectrum of a square matrix $M$ as $\lambda(M)=\{\lambda_{1},\lambda_{2},\dots,\lambda_{N}\}$, its spectral radius as \rho(M)=\max{\{|\lambda_{j}|\,\big{|}\,\lambda_{j}\in\lambda(M)\}} and its leading eigenvalue as $\lambda_{max}(M)=\operatorname*{argmax}_{\lambda_{j}\in\lambda(M)}|\lambda_{j}|$.

A weighted digraph $\mathcal{G=(V,E)}$ consists of a set of nodes $\mathcal{V}$ and a set of edges $\mathcal{E}\subseteq\mathcal{V}\times\mathcal{V}$. Each edge $(j,k)\in\mathcal{E}$ from node $j\in\mathcal{V}$ to node $k\in\mathcal{V}$ has an associated weight $a_{jk}>0$ with $a_{jk}=0$ if $(j,k)\not\in\mathcal{E}$. The set of neighbors of node $j$ is $\mathcal{N}_{j}=\{k\in\mathcal{V}|(j,k)\in\mathcal{E}\}$. $\mathcal{G}$ is strongly connected if there exists a directed path from any node $j\in\mathcal{V}$ to any other node $k\in\mathcal{V}$. The adjacency matrix $A=\{a_{jk}\}$ of $\mathcal{G}$ is irreducible if $\mathcal{G}$ is strongly connected. The degree of node $j$ is $d_{j}=\sum_{k=1}^{N}a_{jk}$. $\mathcal{G}$ is $d$-regular if $d_{j}=d$ for all $j=1,\ldots,N$.

II-B Properties of Metzler Matrices

We use well-known properties of Metzler matrices, which we summarize in the following three propositions (see [29] Theorem 6.2.3, [30] Theorem 11 and 17, and [31] Ch. 8).

Proposition 1**.**

Let $K$ be a Metzler matrix. Then,

$\lambda_{max}(K)\in\mathbb{R}$. If $K$ is irreducible, $\lambda_{max}(K)$ has multiplicity one. 2. 2.

Let $w^{T}$ and $v$ be left and right eigenvectors corresponding to $\lambda_{max}(K)$. Then, $w,v\succeq\mathbf{0}$. If $K$ is irreducible, then $w,v\gg\mathbf{0}$, and every other eigenvector of $K$ has at least one negative entry. 3. 3.

*Let $K_{min},K_{max}$ be irreducible Metzler matrices where *

$K_{min}\prec K\prec K_{max}$, then

[TABLE]

Proposition 2**.**

Let $K$ be a Metzler matrix. Then, the following statements are equivalent:

K is Hurwitz. 2. 2.

There exists a vector $v\gg\mathbf{0}$ such that $Kv\ll\mathbf{0}$. 3. 3.

There exists a vector $w\gg\mathbf{0}$ such that $w^{T}K\ll\mathbf{0}$.

Definition II.1** (Regular Splitting).**

Let $K$ be a Metzler matrix. $K=T+U$ is a regular splitting of $K$ if $T\succeq\bar{\mathbf{0}}$ and $U$ is a Hurwitz Metzler matrix.

Proposition 3**.**

Let $K$ be a Metzler matrix and let $K=T+U$ be a regular splitting. Then,

$\lambda_{max}(K)<0$* if and only if $\rho(-TU^{-1})<1$.* 2. 2.

$\lambda_{max}(K)=0$* if and only if $\rho(-TU^{-1})=1$.* 3. 3.

$\lambda_{max}(K)>0$* if and only if $\rho(-TU^{-1})>1$.*

II-C Properties of Gradient Systems

A gradient system on an open set $\Omega\subseteq\mathbb{R}^{N}$ is a system of the form $\dot{\zeta}=-\nabla V(\zeta)$ where $\zeta(t)\in\Omega$, $V\in C^{2}(\Omega)$ is the potential function, and $\nabla V=[\partial V/\partial\zeta_{1},\dots,\partial V/\partial\zeta_{N}]$ is the gradient of $V$ with respect to $\zeta$. The level surfaces of $V$ are the subsets $V_{c}=\{V^{-1}(c)\in\Omega\,|\,c\in\mathbb{R}\}$. A point $\zeta_{0}\in\Omega$ is a regular point if $\nabla V(\zeta_{0})\neq\mathbf{0}$ and a critical point if $\nabla V(\zeta_{0})=\mathbf{0}$. If $\nabla V(\zeta)\neq\mathbf{0}$ for all $\zeta\in V_{c}$, then $c$ is a regular value for $V$.

Proposition 4** (Properties of Gradient Systems [32, 33]).**

Consider the gradient system $\dot{\zeta}=-\nabla V(\zeta)$ where $V\in C^{2}(\Omega)$, $\zeta(t)\in\Omega\subseteq\mathbb{R}^{N}$. Then,

$V(\zeta)$* is a Lyapunov function of the gradient system. Moreover, $\dot{V}(\zeta)=0$ if and only if $\zeta$ is an equilibrium.* 2. 2.

The critical points of $V$ are the system equilibria. 3. 3.

If $c$ is a regular value for $V$, then the surface set $V_{c}$ forms an $N-1$ dimensional surface in $\Omega$ and the vector field is perpendicular to $V_{c}$. 4. 4.

At every point $\zeta\in\Omega$, the directional derivative along $w\in\mathbb{R}^{N}$ is given by $D_{w}V(\zeta)=w^{T}\nabla V(\zeta)$. 5. 5.

Let $\zeta_{0}$ be an $\alpha$-limit point or an $\omega$-limit point of a solution of the gradient system. Then $\zeta_{0}$ is an equilibrium. 6. 6.

The linearized system at any equilibrium has only real eigenvalues. No periodic solutions are possible.

III Network SIRI Model Dynamics

In this section we present the network SIRI model dynamics, which represent a contagious process with reinfection in a population of $N$ agents. Consider a strongly connected digraph $\mathcal{G=(V,E)}$ with adjacency matrix $A$, where each node in $\mathcal{V}$ represents an agent. The state of each agent $j$ is given by the random variable $X_{j}(t)\in\{S,I,R\}$, where $S$ is “susceptible”, $I$ is “infected”, and $R$ is “recovered”. Let transitions between states for each agent be independent Poisson processes with rates defined as follows. Susceptible agent $j$ becomes infected through contact with infected neighbor $k$ at the rate $\beta_{jk}\geq 0$. We assume $\beta_{jk}=0$ if and only if $a_{jk}=0$. Infected agent $j$ recovers from the infection at the rate $\delta_{j}\geq 0$. Recovered agent $j$ becomes reinfected through contact with infected neighbor $k$ at the rate $\hat{\beta}_{jk}\geq 0$, where $\hat{\beta}_{jk}=0$ if $a_{jk}=0$. These transitions are summarized as

\schemestart

$S_{j}+I_{k}$\arrow-¿[$\beta_{jk}$]$I_{j}+I_{k}$\arrow¡-[$\hat{\beta}_{jk}$]$R_{j}+I_{k}$\schemestop

\schemestart

$I_{j}$\arrow-¿[$\delta_{j}$]$R_{j}$.\schemestop

The dynamics are described by a continuous-time Markov chain, where the probability that an agent transitions state at time $t$ can depend on the state of its neighbors at time $t$. Thus, the dimension of the state space can be as large as $3^{N}$.

To reduce the size of the state space, we use an individual-based mean-field approach (IBMF) [9, 15]. This approach assumes that the state of every node is statistically independent from the state of its neighbors. The approximation reduces the state of every agent $j$ to the probabilities $p_{j}^{S}(t)$, $p_{j}^{I}(t)$, and $p_{j}^{R}(t)$ of agent $j$ being in state $S$, $I$, and $R$, respectively, at time $t\geq 0$. Since at every time $t\geq 0$, these probabilities sum to 1, the state of every agent $j$ evolves on the 2-simplex $\Delta:=\{(p_{j}^{S},p_{j}^{I},p_{j}^{R})\in[0,1]^{3}|\ p_{j}^{S}+p_{j}^{I}+p_{j}^{R}=1\}$. The reduced state space corresponds to $N$ copies of $\Delta$, denoted $\Delta_{N}$, which has dimension $2N$.

The dynamics retain the full topological structure of the network encoded in the infection and reinfection rates $\beta_{jk}$ and $\hat{\beta}_{jk}$, which depend on the entries of the adjacency matrix $A$. We refer the reader to [28] for a detailed derivation of the individual mean-field approximation for the SIS model, and to [16, 34] for a discussion and numerical exploration of the accuracy of mean-field approximations in network dynamics.

Under the individual mean-field approximation, the dynamics of the network SIRI model on $\Delta_{N}$ are given by

[TABLE]

where $p_{j}^{S}(0)+p_{j}^{I}(0)+p_{j}^{R}(0)=1$ for all $j$.

We can reduce the number of equations from $3N$ to $2N$ by using the substitution $p_{j}^{R}=1-p_{j}^{S}-p_{j}^{I}$, for all $j$, in (III):

[TABLE]

where $p_{j}^{S}(0)+p_{j}^{I}(0)=1-p_{j}^{R}(0)$ for all $j$.

The dynamics can be written in matrix form where $p^{\Omega}=[p_{1}^{\Omega},\cdots,p_{N}^{\Omega}]^{T}$ and $P^{\Omega}=\textrm{diag}(p^{\Omega})$ for $\Omega\in\{S,I\}$:

[TABLE]

where

[TABLE]

and

[TABLE]

Further, we define

[TABLE]

The network SIRI model dynamics provide sufficient richness to describe a family of models which can be classified into six different cases (summarized in Table I):

• •

Case 1 (SI): When $D=\bar{\mathbf{0}}$ the network SIRI model specializes to the network SI model.

• •

Case 2 (SIR): When $\hat{B}=\bar{\mathbf{0}}$, the network SIRI model specializes to the network SIR model.

• •

Case 3 (SIS): When $B=\hat{B}$ the network SIRI model specializes to the network SIS model with $p^{S}+p^{R}\mapsto p^{S}$.

• •

Case 4 (Partial Immunity): When ${B}\succ\hat{B}\succ\bar{\mathbf{0}}$, every recovered agent acquires partial (or no) immunity to each of its infected neighbors.

• •

Case 5 (Compromised Immunity): When $\hat{B}\succ B\succ\bar{\mathbf{0}}$, every recovered agent acquires compromised (or no) immunity to each of its infected neighbors.

• •

Case 6 (Mixed Immunity): Models not in Cases 1-5. Notably, there is at least one pair of edges $(j,k)$ and $(l,m)$ such that $\beta_{jk}\geq\hat{\beta}_{jk}$ and $\beta_{lm}<\hat{\beta}_{lm}$. We classify mixed immunity into two sub-cases:

• –

Case 6a (Weak Mixed Immunity): For every agent $j$, $\beta_{jk}-\hat{\beta}_{jk}\geq 0$ for all $k\in\mathcal{N}_{j}$ or $\beta_{jk}-\hat{\beta}_{jk}\leq 0$ for all $k\in\mathcal{N}_{j}$.

• –

Case 6b (Strong Mixed Immunity): Mixed immunity that is not weak.

IV Equilibria and Reproduction Numbers

In this section we analyze the equilibria of the network SIRI model dynamics and define the notion of basic and extreme basic reproduction numbers. We denote the value of $p^{S}$ and $p^{I}$ at equilibrium as $p^{S*}$ and $p^{I*}$, respectively.

IV-A Equilibria

Proposition 5**.**

The only equilibria of the network SIRI model (3) are an invariant set of infection-free equilibria (IFE) and one or more isolated endemic equilibria (EE). The IFE set is defined as $\mathcal{M}=\{(p^{S*},\mathbf{0})\in\Delta_{N}|\ \mathbf{0}\preceq p^{S*}\preceq\mathbf{1}\}$, corresponding to all equilibria in which $p^{I*}=\mathbf{0}$, i.e., $p^{S*}+p^{R*}=\mathbf{1}$. The EE are defined as equilibria where $p^{I*}\succ\mathbf{0}$ satisfies

[TABLE]

If $\hat{B}$ is irreducible then, for every EE, $p^{S*}=\mathbf{0}$ and $p^{I*}\gg\mathbf{0}$.

Proof.

Setting $\dot{p}^{S}=\mathbf{0}$ in (III), we get $P^{S*}Bp^{I*}=\mathbf{0}$. Since $\mathcal{G}$ is strongly connected and $B$ preserves the connectivity of $A$, then for every agent $j$ we must have $p_{j}^{S*}=0$ or $\sum_{\mathcal{N}_{j}}p_{k}^{I*}=0$. Moreover, since $\hat{B}$ has a zero at every entry where $B$ has a zero, it follows that $P^{S*}\hat{B}p^{I*}=\mathbf{0}$.

Setting $\dot{p}^{I}=\mathbf{0}$ in (III) and using $P^{S*}Bp^{I*}=P^{S*}\hat{B}p^{I*}=\mathbf{0}$ we get

[TABLE]

One solution is the invariant set $\mathcal{M}=\{(p^{S*},\mathbf{0})\in\Delta_{N}|\ \mathbf{0}\preceq p^{S*}\preceq\mathbf{1}\}$. The only other solutions are isolated equilibria $p^{I*}\succ\mathbf{0}$ satisfying (5).

If $\hat{B}$ is irreducible, $\beta_{jk}>0$ for any $(j,k)\in\mathcal{E}$, and if $p_{k}^{I*}>0$ for any $k\in\mathcal{V}$, then by (5) $p_{j}^{I*}>0$ for any $j$ where $k\in\mathcal{N}_{j}$. So $p_{i}^{I*}>0$ for any $i$ where $j\in\mathcal{N}_{i}$. This argument can be recursively applied until all nodes in $\mathcal{G}$ are covered. Since $P^{S*}Bp^{I*}=\mathbf{0}$, $p^{I*}\gg\mathbf{0}$ implies $p^{S*}=\mathbf{0}$. ∎

Proposition 6**.**

The boundary of $\mathcal{M}$ is $\partial\mathcal{M}=\{x=(p^{S*},\mathbf{0})\in\mathcal{M}\,|\,\exists j,\,p_{j}^{S*}\in\{0,1\}\}$. The corner set of $\mathcal{M}$ is $\hat{\mathcal{M}}=\{x=(p^{S*},\mathbf{0})\in\partial\mathcal{M}\,|\,p_{j}^{S*}\in\{0,1\},\,\forall j\}$. The interior of $\mathcal{M}$ is $\textrm{int}(\mathcal{M})=\mathcal{M}\setminus\partial\mathcal{M}$.

Proof.

The proof follows from the definition of $\mathcal{M}$. ∎

Remark 1*.*

For $\hat{B}$ irreducible, the equilibria of the network SIRI model are equivalent to the equilibria of the network SIS model (Case 3), where $B=\hat{B}$. This follows since any equilibrium of (III) satisfies (6), and therefore is also an equilibrium of the network SIS dynamics [4, 28, 14]:

[TABLE]

For the network SI model (Case 1), the only equilibrium is a unique EE with $p^{I*}=\mathbf{1}$ and $p^{S*}=\mathbf{0}$. For the network SIR model (Case 2), the only equilibria are the IFE set $\mathcal{M}$.

Remark 2*.*

For initial conditions $p^{S}(0)=\mathbf{1}-p^{I}(0)$, the network SIRI dynamics (III) initially behave as the network SIS model (7) with infection matrix $B$. As agents become exposed to the infection for the first time, the dynamics transition to network SIS dynamics with infection matrix $\hat{B}$.

In the remainder of this paper, we assume $D$ is nonsingular and $\hat{B}$ is irreducible; thus every EE is strong since $p^{I*}\gg\mathbf{0}$. The generalization to reducible $\hat{B}$ is straightforward.111The graph $\mathcal{G}_{\hat{B}}$ with reducible adjacency matrix $\hat{B}$ is weakly connected or disconnected. If $\mathcal{G}_{\hat{B}}$ is weakly connected, the adjacency matrix of $\mathcal{G}_{\hat{B}}$ can be written as an upper block triangular matrix with $K$ diagonal irreducible blocks that describe the $K$ strongly connected subgraphs of $\mathcal{G}$ [35]. If $\mathcal{G}_{\hat{B}}$ is disconnected, it is sufficient to study each connected subgraph of $\mathcal{G}_{\hat{B}}$.

IV-B Basic Reproduction Numbers

The basic reproduction number is a key concept in epidemiology, defined as the expected number of new cases of infection caused by a typical infected individual in a population of susceptible individuals [1, 36, 37]. For many deterministic epidemiological models, an infection can invade and persist in a fully susceptible population if and only if $R_{0}>1$ [38]. For the SIS model in a well-mixed population in which all agents have the same infection rate $\beta$ and recovery rate $\delta$, the basic reproduction number, $R_{0}=\beta/\delta$, governs the steady-state behavior of solutions. If $R_{0}\leq 1$, the infection eventually dies out. If $R_{0}>1$, the infection spreads through the population, converging to an endemic steady-state solution with a fixed fraction of the population in the infected state. In the network SIS model (7), the basic reproduction number is ${R_{0}=\rho(B\Gamma^{-1})}$ [37, 14, 39]. In this case, if $R_{0}\leq 1$, solutions reach the IFE set $\mathcal{M}$ as $t\to\infty$ while if $R_{0}>1$, solutions reach the unique EE (5) as $t\to\infty$ [4, 14].

In previous work [27] we proved that, in well-mixed settings, the transient and steady-state behavior of solutions in the SIRI model depend on two numbers $R_{0}$ and $R_{1}$, corresponding to the basic reproduction number for a population of susceptible individuals and for a population of recovered individuals, respectively. Here we extend the definition of $R_{0}$ and $R_{1}$ in [27] to network topologies and introduce the notion of extreme basic reproduction numbers.

Definition IV.1** (Basic and Extreme Basic Reproduction Numbers).**

Consider the following system, which is a modification to dynamics (III):

[TABLE]

where

[TABLE]

$\mathbf{0}\preceq p\preceq\mathbf{1}$* is constant, and $P=\textrm{diag}(p)$. Let $R(p)$ be the basic reproduction number for (IV.1). We distinguish $R(p)$ for four different values of $p$ as follows:*

• •

Basic infection reproduction number* $R_{0}$ is the reproduction number for (IV.1) with $p=\mathbf{1}$.*

• •

Basic reinfection reproduction number* $R_{1}$ is the reproduction number for (IV.1) with $p=\mathbf{0}$.*

• •

Maximum basic reproduction number* $R_{max}$ is the reproduction number for (IV.1) with $p=\operatorname*{argmax}_{\mathbf{0}\preceq q\preceq\mathbf{1}}R(q)$.*

• •

Minimum basic reproduction number* $R_{min}$ is the reproduction number for (IV.1) with $p=\operatorname*{argmin}_{\mathbf{0}\preceq q\preceq\mathbf{1}}R(q)$.*

Proposition 7** (Spectral Radius Formulas for Reproduction Numbers).**

The basic reproduction number $R(p)$ of (IV.1) for $\mathbf{0}\preceq p\preceq\mathbf{1}$ can be computed as

[TABLE]

Therefore,

[TABLE]

[TABLE]

Proof.

The linear term of the dynamics of $p^{I}$ in (IV.1) is $K=B^{*}(p)-D$. Thus, $K$ is Metzler and by Definition II.1 $K=T+U$ is a regular splitting where $T=B^{*}(p)$ and $U=-D$. Following [39], $R(p)=\rho(-TU^{-1})=\rho(B^{*}(p)D^{-1})$. ∎

Proposition 8** (Reproduction number ordering).**

Let $\bar{R}_{max}=\rho(\bar{B}_{max}D^{-1})$ and $\bar{R}_{min}=\rho(\bar{B}_{min}D^{-1})$. Then,

[TABLE]

for any $\mathbf{0}\preceq p^{S}\preceq\mathbf{1}$. If $B\succeq\hat{B}$, then $\bar{R}_{max}=R_{max}=R_{0}$ and $\bar{R}_{min}=R_{min}=R_{1}$. If $\hat{B}\succeq B$, then $\bar{R}_{max}=R_{max}=R_{1}$, and $\bar{R}_{min}=R_{min}=R_{0}$.

Proof.

Any matrix with nonnegative entries is Metzler. Thus, $Y(p^{S})=B^{*}(p^{S})D^{-1}$ is an irreducible Metzler matrix since $B\succ\bar{\mathbf{0}}$ and $\hat{B}\succeq\bar{\mathbf{0}}$ are irreducible. By Proposition 1, $\lambda_{max}(Y(p^{S}))$ increases (decreases) as any entry in $Y(p^{S})$ increases (decreases). Since every non-zero entry of $Y(p^{S})$ is a scaled convex sum of $\beta_{jk}$ and $\hat{\beta}_{jk}$, it follows that $\bar{B}_{min}D^{-1}\preceq Y(p^{S})\preceq\bar{B}_{max}D^{-1}$ for any $\mathbf{0}\preceq p^{S}\preceq\mathbf{1}$. Consequently, (9) holds for any $\mathbf{0}\preceq p^{S}\preceq\mathbf{1}$. If $B\succeq\hat{B}$, then $\bar{B}_{max}=B$ and $\bar{B}_{min}=\hat{B}$. If $\hat{B}\succeq B$, then $\bar{B}_{max}=\hat{B}$ and $\bar{B}_{min}=B$. The stated results then follow from the definitions of the reproduction numbers. ∎

V Stability of Equilibria

In this section we prove conditions for the local stability of the EE and of points in the IFE set $\mathcal{M}$.

V-A Stability of the Endemic Equilibria

Proposition 9**.**

The network SIRI dynamics (III) have a unique EE, given by (5), which exists if and only if $R_{1}>1$. When it exists, the EE is locally stable.

Proof.

Since $B$ and $\hat{B}$ are irreducible, $p^{S*}=\mathbf{0}$. Following the proof of Proposition 5 and the argument of Remark 1, the equilibria of the network SIRI model (III) are equivalent to the equilibria of the network SIS model (7) where $B=\hat{B}$. Thus, the proof of existence and uniqueness of the EE of (III), if and only if $R_{1}>1$, follows the proof in Section 2.2 of [14] for existence and uniqueness of the EE of (7) where $B=\hat{B}$.

To prove local stability, we compute the Jacobian of (III) at the unique EE (5):

[TABLE]

where $J_{a}=\hat{B}-D-P^{I*}\hat{B}-\textrm{diag}(\hat{B}p^{I*})$. Since $-\textrm{diag}(Bp^{I*})$ is Hurwitz, showing that the EE is locally stable is equivalent to showing that the Metzler matrix $J_{a}$ is Hurwitz.

By (6), $(\hat{B}-D-P^{I*}\hat{B})p^{I*}=0$, and thus

[TABLE]

where the inequality follows from $p^{I*}\gg\mathbf{0}$. By Proposition 2 we conclude that $J_{a}$ is Hurwitz. ∎

V-B Stability of Infection-Free Equilibria

In this section we prove results on the stability of the IFE set $\mathcal{M}$. The equilibria in $\mathcal{M}$ are non-hyperbolic: the Jacobian of (III) at a point $x\in\mathcal{M}$ has $N$ zero eigenvalues corresponding to the $N$-dimensional space tangent to $\mathcal{M}$. The remaining $N$ eigenvalues are called transverse as they correspond to the $N$-dimensional space transverse to $\mathcal{M}$. The Shoshitaishvili Reduction Principle [40], which extends the Hartman-Grobman Theorem to non-hyperbolic equilibria, can be used to study the local stability of points in $\mathcal{M}$ in terms of the transverse eigenvalues of the Jacobian and the dynamics on the center manifold. We show how the irreducibility of $B$ and $\hat{B}$ imply that the behavior of solutions in $\Delta_{N}$ close to a point $x\in\mathcal{M}$ depends only on the sign of the leading transverse eigenvalue of the Jacobian at $x$.

Throughout the rest of this paper, we consider the topological space $\Delta_{N}$ as a subspace of $\mathbb{R}^{2N}$. This allows us to study points in $\partial\mathcal{M}$ and in $\textrm{int}(\mathcal{M})$ simultaneously. In $\mathbb{R}^{2N}$, the invariant set $\mathcal{M}$ of IFE points becomes a subset of the invariant manifold of equilibria $\mathcal{M}^{\prime}=\{(p,\mathbf{0})|p\in\mathbb{R}^{N}\}$.

Lemma 1** (Local Stability of Points in the IFE set $\mathcal{M}$).**

Let $x=(p^{S*},\mathbf{0})\in\mathcal{M}$. Let $J_{\mathcal{M}}(x)$ be the Jacobian of (III) at $x$ and $\lambda_{Tmax}(J_{\mathcal{M}}(x))$ the leading transverse eigenvalue of $J_{\mathcal{M}}(x)$. Then, $\lambda_{Tmax}(J_{\mathcal{M}}(x))\in\mathbb{R}$ and the following hold.

• •

Suppose $\lambda_{Tmax}(J_{\mathcal{M}}(x))<0$. Then, $x$ is locally stable. I.e., given a neighborhood $U$ of $x$ on $\mathcal{M}^{\prime}$ such that $\lambda_{Tmax}(J_{\mathcal{M}}(u))<0$ for all $u\in U$, there exists $V\subset\Delta_{N}$ and $x\in V$ such that any solution starting in $V$ converges exponentially to a point in $U\cap\Delta_{N}$.

• •

Suppose $\lambda_{Tmax}(J_{\mathcal{M}}(x))>0$. Then, $x$ is unstable. I.e., there exists $W\subset\Delta_{N}$ and $x\in W$, such that any solution starting in $W$ leaves $W$.

Proof.

For an arbitrary point $x=(p^{S*},\mathbf{0})\in\mathcal{M}$,

[TABLE]

where $J_{T}(p^{S*})=B^{*}(p^{S*})-D$. The $N$ transverse eigenvalues of $J_{\mathcal{M}}(x)$ are the eigenvalues of $J_{T}(p^{S*})$ and so $\lambda_{Tmax}(J_{\mathcal{M}}(x))=\lambda_{max}(J_{T}(p^{S*}))$. The matrix $J_{T}(p^{S*})$ is Metzler irreducible since $B$ and $\hat{B}$ are Metzler irreducible. By Proposition 1 $\lambda_{max}(J_{T}(p^{S*}))\in\mathbb{R}$.

Consider an arbitrary point $x^{\prime}=(p^{\prime},\mathbf{0})\in\mathcal{M}^{\prime}\setminus\mathcal{M}$. The Jacobian of (III) at $x^{\prime}$ takes on the same form as (12), and $J_{T}(p^{\prime})$ is Metzler irreducible if every entry of $p^{\prime}$ satisfies

[TABLE]

Since $B,\hat{B}$ are irreducible, $\mid\sum_{j=1}^{N}\hat{\beta}_{jk}/\sum_{j=1}^{N}(\beta_{jk}-\hat{\beta}_{jk})\mid>1$ for all $j$. So, for any $x\in\partial\mathcal{M}$, there exists a neighborhood $\bar{U}$ of $x$ on $\mathcal{M}^{\prime}$ such that $J_{T}(\bar{u}^{\prime})$ is Metzler irreducible for every $\bar{u}=(\bar{u}^{\prime},\mathbf{0})\in\bar{U}$. By Proposition 1 $\lambda_{max}(J_{T}(p^{\prime}))\in\mathbb{R}$.

Let $U$ be a neighborhood of $x$ on $\mathcal{M}^{\prime}$ such that $\lambda_{Tmax}(J_{\mathcal{M}}(u))$ has the same sign as $\lambda_{Tmax}(J_{\mathcal{M}}(x))$ for all $u\in U$. Then, $\lambda_{max}(J_{T}(u^{\prime}))$ has the same sign as $\lambda_{max}(J_{T}(p^{S*}))$ for all $u=(u^{\prime},\mathbf{0})\in U$. By Proposition 1 every left and right eigenvector of every eigenvalue of $J_{T}(u^{\prime})$, other than $\lambda_{max}(J_{T}(u^{\prime}))$, contains at least one negative entry. Thus, for any $\bar{u}\in U\cap\Delta_{N}$, the eigenvector corresponding to $\lambda_{Tmax}(J_{\mathcal{M}}(\bar{u}))$ lies in $\Delta_{N}$, and the eigenvectors corresponding to the other $N-1$ transverse eigenvalues lie outside $\Delta_{N}$.

If $\lambda_{max}(J_{T}(p^{S*}))<0$, then every transverse eigenvalue of $J_{\mathcal{M}}(x)$ has negative real part. By the Shoshitaishvili Reduction Principle [40], there exists a neighborhood $V^{\prime}\in\mathbb{R}^{2N}$ of $x$ that is positively invariantly foliated by a family of stable manifolds corresponding to the family of stationary solutions in $U$ (see [41, 42, 43]), each stable manifold spanned by the (generalized) eigenvectors associated with the $N$ negative transverse eigenvalues of $J_{\mathcal{M}}(u)$. Let $V=V^{\prime}\cap\Delta_{N}$. Then $V\subset\Delta_{N}$ is positively invariantly foliated by a family of stable manifolds. The invariance of $\Delta_{N}$ implies each of these stable manifolds corresponds to a point $\bar{u}\in U\cap\mathcal{M}$. Thus, any solution starting in $V$ converges exponentially along a stable manifold to the corresponding stationary solution in $U\cap\mathcal{M}$.

If $\lambda_{max}(J_{T}(p^{S*}))>0$, then there is at least one transverse eigenvalue of $J_{\mathcal{M}}(x)$ with positive real part. The trace of $J_{T}(p^{S*})$ is negative for any ${\mathbf{0}\preceq p^{S*}\preceq\mathbf{1}}$, so the sum of the eigenvalues of $J_{T}(p^{S*})$ is always negative and $J_{\mathcal{M}}(x)$ has at least one transverse eigenvalue with negative real part. By the Shoshitaishvili Reduction Principle [40] there exists a neighborhood $W^{\prime}\subset\mathbb{R}^{2N}$ of $x$ that is positively invariantly foliated by a family of stable, unstable, and possibly center manifolds corresponding to the family of stationary solutions in $U$. Let $W=W^{\prime}\cap\Delta_{N}$. Then, the stable and center manifolds of each stationary solution $\bar{u}\in W\cap\mathcal{M}$ lie outside $\Delta_{N}$. Thus, no solution starting in $W$ can remain in $W$ for all time, i.e., any solution starting in $W$ leaves $W$. ∎

Definition V.1** (Stable, unstable, and center IFE subsets).**

The stable IFE subset is $\mathcal{M}_{-}=\{x\in\mathcal{M}\,|\,\lambda_{Tmax}(J_{\mathcal{M}}(x))<0\}$. The unstable IFE subset is $\mathcal{M}_{+}=\{x\in\mathcal{M}\,|\,\lambda_{Tmax}(J_{\mathcal{M}}(x))>0\}$. The center IFE subset is $\mathcal{M}_{0}=\{x\in\mathcal{M}\,|\,\lambda_{Tmax}(J_{\mathcal{M}}(x))=0\}$.

Proposition 10**.**

${\mathcal{M}_{-}\cup\mathcal{M}_{+}\cup\mathcal{M}_{0}=\mathcal{M}}$. Every point in $\mathcal{M}_{-}$ is locally stable and every point in $\mathcal{M}_{+}$ is unstable.

Proof.

This follows from Definition V.1 and Lemma 1. ∎

We now state the first two theorems of the paper, which relate the extreme basic reproduction numbers $R_{max}$ and $R_{min}$ to the stable, unstable, and center subsets of the IFE set $\mathcal{M}$.

Theorem 1** (Stability of the IFE set $\mathcal{M}$).**

1. (A)

If $R_{max}<1$, then $\mathcal{M}_{-}=\mathcal{M}$. 2. (B)

If $R_{min}>1$, then $\mathcal{M}_{+}=\mathcal{M}$. 3. (C)

If $R_{max}=R_{min}=1$, then $\mathcal{M}_{0}=\mathcal{M}$. 4. (D)

If $R_{min}<R_{max}=1$, then $\mathcal{M}_{-}=\mathcal{M}\setminus\mathcal{M}_{0}$ and $\mathcal{M}_{0}\subset\partial\mathcal{M}$. 5. (E)

If $R_{max}>R_{min}=1$, then $\mathcal{M}_{+}=\mathcal{M}\setminus\mathcal{M}_{0}$ and $\mathcal{M}_{0}\subset\partial\mathcal{M}$. 6. (F)

If ** $R_{max}>1>R_{min}$**, then $\mathcal{M}_{-},\mathcal{M}_{+},\mathcal{M}_{0}\neq\emptyset$ and each subset consists of $n_{-},n_{+},n_{0}$ connected sets, respectively. Each of the center connected sets $\mathcal{M}_{0}^{j}$, $j=1,\dots,n_{0}$, is an $N-1$-dimensional smooth hypersurface with boundary $\partial\mathcal{M}_{0}^{j}\subset\partial\mathcal{M}$. Each $\mathcal{M}_{0}^{j}$ separates an $N$-dimensional stable connected hypervolume from an $N$-dimensional unstable connected hypervolume.

Remark 3*.*

Theorem 1 applies to the six different cases of the network SIRI model as follows (see Table II). A applies to Cases 2, 3, 4, 5, and 6. B applies to Cases 3, 4, 5, and 6. C applies to Case 3. D applies to Cases 2, 4, 5, and 6. E applies to Cases 4, 5, and 6. F applies to Cases 2, 4, 5, and 6. We specialize F in Theorem 2 to provide the key to characterizing global behavior in Cases 2, 4, 5, and 6a.

Theorem 2** (Uniqueness of stable, unstable, and center subsets).**

If $R_{max}>1>R_{min}$, then for Case 2 (SIR), Case 4 (partial immunity), Case 5 (compromised immunity), and Case 6a (weak mixed immunity), $\mathcal{M}_{0}$ consists of a unique $(N-1)$-dimensional hypersurface with boundary $\partial\mathcal{M}_{0}\subset\partial\mathcal{M}$ dividing $\mathcal{M}$ into $\mathcal{M}_{-}$ and $\mathcal{M}_{+}$.

Remark 4*.*

We conjecture that Theorem 2 can be extended to Case 6b. Extensive computations of $\mathcal{M}_{0}$, $\mathcal{M}_{-}$, and $\mathcal{M}_{+}$, for an $N=3$ agent network, with different network configurations and parameter values, consistently show a unique connected surface $\mathcal{M}_{0}$ dividing $\mathcal{M}$ into $\mathcal{M}_{-}$ and $\mathcal{M}_{+}$.

The proofs of Theorems 1 and 2 make use of the following definition and lemmas.

Lemma 2** (Neighborhood $E\subset\mathcal{M}^{\prime}$ of $\mathcal{M}$).**

Let $E$ be the union of $\mathcal{M}$ and the neighborhoods $\bar{U}\subset\mathcal{M}^{\prime}$ of every $\bar{x}\in\partial\mathcal{M}$ described in the proof of Lemma 1. Then, $E\subset\mathcal{M}^{\prime}$ is a neighborhood of $\mathcal{M}$.

Lemma 3** ($J_{T}$ and $\Lambda$).**

Let $J_{T}(p)=B^{*}(p)-D$ for $(p,\mathbf{0})\in E$ and $\Lambda:E\to\mathbb{R}$, $(p,\mathbf{0})\mapsto\lambda_{max}(J_{T}(p))$. For ease of notation we use $\Lambda(p)$ for $\Lambda(p,\mathbf{0})$. Then $\Lambda(p)=\lambda_{Tmax}(J_{\cal M}(p,\mathbf{0}))$.

Definition V.2** (Stubborn agents).**

An agent $j\in\mathcal{V}$ is stubborn if $(\beta_{jk}-\hat{\beta}_{jk})=0$ for all $k\in\mathcal{N}_{j}$.

Lemma 4** (IFE subsets as level surfaces of $\Lambda$).**

*Let $\Lambda_{c}=\{(p,\mathbf{0})\in E\,|\,\Lambda(p)=c\}$ be the level surface of $\Lambda$ on $E\subset\mathcal{M}^{\prime}$ corresponding to $c\in\mathbb{R}$. Then, $\mathcal{M}_{0}=\Lambda_{0}\cap\mathcal{M}$, $\mathcal{M}_{-}=\bigcup_{c<0}\Lambda_{c}\cap\mathcal{M}$ and $\mathcal{M}_{+}=\bigcup_{c>0}\Lambda_{c}\cap\mathcal{M}$. *

Proof.

This follows from Definition V.1 and Lemma 3. ∎

Lemma 5** (Gradient of $\Lambda$).**

For $\tilde{x}=(p,\mathbf{0})\in E$, let $w^{T},v\in\mathbb{R}^{N}$ be left and right eigenvectors of $J_{T}(p)$ for $\Lambda(p)$. Then, $\Lambda$ is smooth on $E$, i.e., $\Lambda(\cdot)\in C^{\infty}(E)$, with partial derivatives

[TABLE]

and gradient

[TABLE]

In addition, the following hold.

• •

If $B=\hat{B}$, all points $\tilde{x}\in E$ are critical points of $\Lambda$.

• •

If $j$ is a stubborn agent, $(\partial\Lambda/\partial p_{j})(p)=0$ for all $\tilde{x}\in E$.

• •

If $B\succ\hat{B}$, and there are no stubborn agents in $\mathcal{G}$, $\nabla\Lambda(p)\gg\mathbf{0}$ for all $\tilde{x}\in E$, and $\Lambda$ has no critical points.

• •

If $\hat{B}\succ B$, and there are no stubborn agents in $\mathcal{G}$, $\nabla\Lambda(p)\ll\mathbf{0}$ for all $\tilde{x}\in E$ and $\Lambda$ has no critical points.

• •

If $B\neq\hat{B}$ then either $\Lambda$ has no critical points or all points $\tilde{x}\in E$ are critical points of $\Lambda$.

Proof.

Let $\tilde{x}=(p,\mathbf{0})\in E$. By the proof of Lemma 1, $J_{T}(p)$ is Metzler irreducible. By Proposition 1, $\lambda_{max}(J_{T}(p))\in\mathbb{R}$ and has multiplicity one. Thus by [44, p. 66-67], $\Lambda(\cdot)\in C^{\infty}(E)$.

Differentiating the right-eigenvector equation $J_{T}(p)v=\Lambda(p)v$ with respect to $p_{j}$, and premultiplying by $w^{T}$ we get

[TABLE]

where all terms are evaluated at $p$, $\partial J_{T}/\partial p_{j}(p)=\{z_{kl}\}$ and

[TABLE]

Because $\Lambda(p)$ has multiplicity one, we can always pick $w$ such that $w^{T}v=1$. Using $w^{T}J_{T}(p)=w^{T}\Lambda(p)$, we obtain

[TABLE]

In vector form, (15) becomes $\nabla\Lambda(p)=\textnormal{diag}(w)(B-\hat{B})v.$

By Proposition 1 $w,v\gg\mathbf{0}$. If $B=\hat{B}$, then $\nabla\Lambda(p)=\mathbf{0}$ for all $\tilde{x}\in E$. If $j$ is a stubborn agent, then $\beta_{jk}-\hat{\beta}_{jk}=0$ for all $k\in\mathcal{N}_{j}$. Therefore, $(\partial\Lambda/\partial p_{j})(p)=0$ for all $\tilde{x}\in E$. If $B\succ\hat{B}$ and there are no stubborn agents in $\mathcal{G}$, then $\sum_{k=1}^{N}(\beta_{jk}-\hat{\beta}_{jk})v_{k}>0$ for any $j$. Therefore, $\nabla\Lambda(p)\gg\mathbf{0}$ for all $\tilde{x}\in E$. A similar argument is made if $\hat{B}\succ B$, and there are no stubborn agents in $\mathcal{G}$.

By Proposition 4 and (14), $\tilde{x}$ is a critical point of $\Lambda$ if and only if $(B-\hat{B})v=\mathbf{0}$. Let $\tilde{x}_{0}=(p_{0},\mathbf{0})\in E$ be a critical point, then $(B-\hat{B})v_{0}=\mathbf{0}$. Thus, for any $p$, we have

[TABLE]

In particular, $(\hat{B}-D)v_{0}=\Lambda(p_{0})v_{0}$. Hence, $(\Lambda(p_{0}),v_{0})$ is an eigenpair of $\hat{B}-D$ and by (V-B) also an eigenpair of $J_{T}(p)$ for any $p$. Since $\Lambda(p)$ has multiplicity one and $v$ is the only eigenvector satisfying $v\gg\mathbf{0}$, then $(\Lambda(p),v)=(\Lambda(p_{0}),v_{0})$ and $\tilde{x}\in E$ is a critical point of $\Lambda$. ∎

Lemma 6** ($\Lambda_{c}$ with stubborn agents).**

If agents $1,\dots,k$, $k<N$, are stubborn, then at any $(p,\mathbf{0})\in E$, the subspace spanned by $\{e^{1},\dots,e^{k}\}$ is tangent to the level surface $\Lambda_{c}$ with $p\in\Lambda_{c}$.

Proof.

By Lemma 5, for any stubborn agent $j$, $(e^{j})^{T}\nabla\Lambda(p)=0$ for any $(p,\mathbf{0})\in E$. By Proposition 4, $e^{j}$ is tangent to the level surface $\Lambda_{c}$ with $p\in\Lambda_{c}$. ∎

Lemma 7** (Maximum and minimum values of $\Lambda$ on $\mathcal{M}$).**

$\Lambda$* achieves its global maximum and minimum $c_{max},c_{min}$ on $\mathcal{M}$ at one or more points in $\partial\mathcal{M}$. In Cases 2, 4, 5, and 6a, there exist unique corner points $(p^{max},\mathbf{0})\in\hat{\mathcal{M}}$, $(p^{min},\mathbf{0})\in\hat{\mathcal{M}}$ such that $\Lambda(p^{max})=c_{max}$, $\Lambda(p^{min})=c_{min}$, $R_{max}=\rho(B^{*}(p^{max})D^{-1})$, and $R_{min}=\rho(B^{*}(p^{min})D^{-1})$. Moreover, $p^{max}$ and $p^{min}$ are the respective unique global maximum and minimum points of $\Lambda$ on $\mathcal{M}$ if and only if there are no stubborn agents in $\mathcal{G}$.*

Proof.

The case in which every point in $E$ is a critical point of $\Lambda$ is trivial. Assume $\Lambda$ has no critical points in $E$. Since $\mathcal{M}\subset\mathcal{M}^{\prime}$ is a compact set and $\Lambda$ is continuous, then by the Extreme Value Theorem, $\Lambda$ achieves $c_{max}$ and $c_{min}$ at one or more points in $\partial\mathcal{M}$. Let $(p^{max},\mathbf{0})\in\partial\mathcal{M}$ and $(p^{min},\mathbf{0})\in\partial\mathcal{M}$ be points such that $\Lambda(p^{max})=c_{max}$ and $\Lambda(p^{min})=c_{min}$.

Let $(p^{S},\mathbf{0})\in\mathcal{M}$. Assume there are no stubborn agents. The components $(J_{T}(p^{S}))_{jk}=(1-p_{j}^{S})\hat{\beta}_{jk}+p_{j}^{S}\beta_{jk}$, $j\neq k$, are maximized and minimized at $p_{j}^{S}=0$ or $p_{j}^{S}=1$ for all $j$. If $B\succ\hat{B}$ (Cases 2 and 4), the entries $(J_{T}(p^{S}))_{jk}$ are maximized at $p_{j}^{S}=1$ and minimized at $p_{j}^{S}=0$ for all $j$. It follows from Proposition 1 that $p^{max}=\mathbf{1}\in\hat{\mathcal{M}}$ and $p^{min}=\mathbf{0}\in\hat{\mathcal{M}}$. Using a similar argument if $\hat{B}\succ B$ (Case 5) $p^{max}=\mathbf{0}\in\hat{\mathcal{M}}$ and $p^{min}=\mathbf{1}\in\hat{\mathcal{M}}$. Further, in Case 6a with no stubborn agents, $\beta_{jk}-\hat{\beta}_{jk}>0$ for all $k\in\mathcal{N}_{j}$ or $\beta_{jk}-\hat{\beta}_{jk}<0$ for all $k\in\mathcal{N}_{j}$. Thus, similarly, $p_{j}^{max},p_{j}^{min}\in\{0,1\}$ for all $j$.

For any $p^{S}\neq p^{max},p^{min}$, it follows that $B^{*}(p^{max})D^{-1}\succ B^{*}(p^{S})D^{-1}\succ B^{*}(p^{min})D^{-1}$. By Proposition 1, $R_{max}=\rho(B^{*}(p^{max})D^{-1})$ and $R_{min}=\rho(B^{*}(p^{min})D^{-1})$.

Now suppose there are stubborn agents $1,\dots,k$, $k<N$. Let $v_{e}\in{\rm span}\{e^{1},\dots,e^{k}\}$. Then, by Lemma 6, if $p^{max}\in\Lambda_{c_{max}}$ then $p^{max}+v_{e}\in\Lambda_{c_{max}}$ as long as $(p^{max}+v_{e},\mathbf{0})\in E$. Similarly, if $p^{min}\in\Lambda_{c_{min}}$ then $p^{min}+v_{e}\in\Lambda_{c_{min}}$ as long as $(p^{min}+v_{e},\mathbf{0})\in E$. ∎

We now present the proof for Theorems 1 and 2.

Proof of Theorem 1.

To prove A, assume $R_{max}<1$. By Proposition 3, $\lambda_{Tmax}(J_{\mathcal{M}}(x))<0$ for all $x\in\mathcal{M}$. Therefore, $\mathcal{M}_{-}=\mathcal{M}$, $\mathcal{M}_{+}=\emptyset$, and $\mathcal{M}_{0}=\emptyset$. The proof of B and C follow similarly. To prove D, assume $R_{min}<R_{max}=1$. By Proposition 3, $\max_{x\in\mathcal{M}}\lambda_{Tmax}(J_{\mathcal{M}}(x))=0$. Therefore, $\mathcal{M}_{0}\neq\emptyset$ and $\mathcal{M}_{-}=\mathcal{M}\setminus\mathcal{M}_{0}$. By Lemma 7, $\mathcal{M}_{0}\subset\partial\mathcal{M}$. The proof of E follows similarly.

To prove F, assume $R_{max}>1>R_{min}$. By Proposition 3 and Lemma 7, there exist points $x_{max}=(p_{max}^{S},\mathbf{0})$ and $x_{min}=(p_{min}^{S},\mathbf{0})$ in $\partial\mathcal{M}$ such that $\lambda_{Tmax}(J_{\mathcal{M}}(x_{max}))>0$ and $\lambda_{Tmax}(J_{\mathcal{M}}(x_{min}))<0$. By the continuity of $\Lambda$ on $E$ and Lemma 4, it follows that $\mathcal{M}_{-},\mathcal{M}_{+},\mathcal{M}_{0}\neq\emptyset$, with each subset consisting of $n_{-},n_{+}$, and $n_{0}$ connected sets, respectively, where each of the $n_{0}$ center sets separates $n_{-}$ stable connected sets from $n_{+}$ unstable connected sets.

Since, by Lemma 5, $\Lambda$ has no critical points in $E$, it follows that $c$ is a regular value of $\Lambda$ for any $c\in\mathbb{R}$. Hence, by the Implicit Function Theorem, every center connected set $M_{0}$ is an $(N-1)$-dimensional smooth hypersurface, and every stable and unstable connected set is an $N$-dimensional hypervolume.

Since $\Lambda$ has no critial points in $E$, the gradient dynamics of $\Lambda$ have no equilibria in $\mathcal{M}\subset E$. Thus, no center connected set $M_{0}$ in $\mathcal{M}_{0}$ is compact in $\textrm{int}(\mathcal{M})$ in any direction, since by Proposition 4 there cannot be any $\alpha$-limit points or $\omega$-limit points in $\mathcal{M}$. Thus, $\partial M_{0}$ in $\mathcal{M}$ must be contained in $\partial\mathcal{M}$. ∎

Proof of Theorem 2.

Let $R_{max}>1>R_{min}$. In Cases 2, 4, 5 and 6a, by Lemma 7, there exists a unique corner point $x^{min}=(p^{min},\mathbf{0})\in\hat{\mathcal{M}}$ where $\Lambda(p^{min})\leq\Lambda(p^{S})$ for any $x=(p^{S},\mathbf{0})\in\mathcal{M}$, and $R_{min}=\rho(B^{*}(p^{min})D^{-1})$. By Definition IV.1 and Proposition 3, it follows that $\Lambda(p^{min})<0$.

Given $p^{S}$, there exists a point $(p^{end},\mathbf{0})\in\partial\mathcal{M}$ such that the line connecting $p^{min}$ to $p^{end}$ passes through $p^{S}$. Each point on the line is parameterized by $r\in[0,1]$ as follows: $s(r,p^{end})=\{(p(r),\mathbf{0})\in\mathcal{M}\,|\,p(r)=(1-r)p^{min}+rp^{end}\}$.

Let $\bar{p}=p^{end}-p^{min}$. Then, $\bar{p}$ is tangent to $s(r,p^{end})$, and $\bar{p}_{j}\geq 0$ if $p^{min}_{j}=0$ and $\bar{p}_{j}\leq 0$ if $p^{min}_{j}=1$, for all $j$. By Proposition 4, the directional derivative of $\Lambda$ at $s(r,p^{end})$ in the direction $\bar{p}$ is $D_{\bar{p}}\Lambda(s(r,p^{end}))=\bar{p}^{T}\nabla\Lambda(s(r,p^{end}))$.

By Lemma 7, if $B\succ\hat{B}$ (Cases 2 and 4), then $p^{min}=\mathbf{0}$ and $\bar{p}\succ\mathbf{0}$. By Lemma 5, $\nabla\Lambda(s(r,p^{end}))\succ 0$ and so $\bar{p}^{T}\nabla\Lambda(s(r,p^{end}))\geq 0$ for any $r\in[0,1]$. Similarly, for Cases 5 and 6a, $\bar{p}^{T}\nabla\Lambda(s(r,p^{end}))\geq 0$ for all $r\in[0,1]$.

If $\bar{p}^{T}\nabla\Lambda(s(r,p^{end}))=0$, $\Lambda$ is constant along the line $s(r,p^{end})$, $r\in[0,1]$, and the line describes a level surface $\Lambda_{c}$. Moreover, since $\Lambda(p^{min})<0$, by Lemma 4 $c\neq 0$, i.e., $\mathcal{M}_{0}$ does not intersect $\Lambda_{c}$. For all other lines, $\bar{p}^{T}\nabla\Lambda(s(r,p^{end}))>0$, which implies that $\Lambda$ is strictly increasing from a negative value at the corner $(p^{min},\mathbf{0})$ to the value at $(p^{end},\mathbf{0})\in\partial\mathcal{M}$. By Theorem 1, $\partial\mathcal{M}_{0}\subset\partial\mathcal{M}$. Thus, there is only the possibility of a single crossing of $\Lambda_{0}$ on each of these lines. By Lemma 4 there is a unique center connected hypersurface. ∎

VI Reproduction Numbers Predict Behavioral Regimes

In this section we prove our third theorem, which provides conditions that determine whether solutions of (III) converge to a point in the IFE set $\mathcal{M}$ or to the EE as $t\to\infty$. We show that the basic and extreme basic reproduction numbers $R_{0},R_{1},R_{min},R_{max}$ distinguish four behavioral regimes in the network SIRI model, each characterized by qualitatively different transient and steady-state behaviors.

VI-A The $\omega$-limit Set of Solutions

The components of $p^{S}$ decrease monotonically along solutions of (III). Here we show that this monotonicity implies that all solutions either converge to a point in $\mathcal{M}$ or to the EE as $t\to\infty$. Moreover, this means that when the EE is not an equilibrium of the dynamics, the infection cannot survive in the network and all solutions reach an IFE point in $\mathcal{M}$. The results in this section are valid even if $\hat{B}$ is not irreducible.

Lemma 8**.**

Let $y(t,y_{0})=(p^{S}(t),p^{I}(t))$ be the solution of (III) with initial condition $y_{0}\in\Delta_{N}$. Then the following hold:

• •

Every point in the $\omega$-limit set $\Omega(y_{0})$ of $y(t,y_{0})$ is an equilibrium of (III).

• •

$y(t,y_{0})$* converges to a point in $\mathcal{M}$ as $t\to\infty$ if $R_{1}\leq 1$.*

• •

$y(t,y_{0})$* converges to a point in $\mathcal{M}$ or to the EE as $t\to\infty$ if $R_{1}>1$.***

Proof.

By invariance of $\Delta_{N}$, any solution $y=y(t,y_{0})$ of (III) with initial condition $y_{0}\in\Delta_{N}$ is bounded and stays in $\Delta_{N}$ for $t\geq 0$. Therefore its $\omega$-limit set $\Omega(y_{0})$ is a nonempty, compact, invariant set, and $y$ approaches $\Omega(y_{0})$ as $t\to\infty$ (see Lemma 4.1 in [45]). Let $V=\mathbf{1}^{T}p^{S}$, then $\dot{V}=-(p^{S})^{T}Bp^{I}\leq 0$ in $\Delta_{N}$. By LaSalle’s Invariance Principle [45], $y$ approaches the largest invariant set $W$ in the set $L=\{(p^{S},p^{I})\in\Delta_{N}\,|\,\dot{V}=0\}$. It follows that, on $L$, $P^{S}Bp^{I}=\mathbf{0}$ which implies $\dot{p}^{S}=\mathbf{0}$. Moreover, since $\hat{B}$ has a zero at every entry where $B$ has a zero, it follows that $P^{S}\hat{B}p^{I}=\mathbf{0}$ on $L$. This in turn implies that the $p^{I}$ dynamics on $L$ are given by the network SIS dynamics (7). Since solutions of (7) either converge to the IFE point $p^{I*}=\mathbf{0}$ or to the EE (5) as $t\to\infty$ [14], it follows that every invariant set of $L$ consists only of equilibria of (III) (see Remark 1). By Proposition 5, $W$ is the union of the IFE set $\mathcal{M}$ and the EE (5). Furthermore, since all $\omega$-limit points are equilibria, $\Omega(y_{0})$ contains a single point corresponding to either a point in the IFE set $\mathcal{M}$ or the EE. If $R_{1}\leq 1$ it follows from Propositions 5 and 9 that the IFE set $\mathcal{M}$ comprises the only equilibria of (III). Therefore, $y$ converges to a point in $\mathcal{M}$ as $t\to\infty$. If $R_{1}>1$, $y$ converges to a point in $\mathcal{M}$ or the EE as $t\to\infty$. ∎

VI-B Behavioral Regimes

We now state the third theorem of the paper, which shows how the four reproduction numbers distinguish the four behavioral regimes: infection-free, endemic, epidemic, and bistable. We interpret and illustrate in Figure 1. The proof follows.

Theorem 3** (Behavioral Regimes).**

Let $p^{I}(0)\succ\mathbf{0}$, and $w^{T}_{m}$ be the leading left-eigenvector associated with $R_{max}$. Then the network SIRI model (III) exhibits four qualitatively distinct behavioral regimes:

Infection-Free Regime: If $R_{max}\leq 1$ the following hold:

1. (a)

All solutions converge to a point in $\mathcal{M}$ as $t\to\infty$. 2. (b)

If $B\succeq\hat{B}$ or $\hat{B}\succeq B$, the weighted infected average $p^{I}_{\rm avg}(t)=w^{T}_{m}D^{-1}p^{I}(t)$ decays monotonically to zero. 2. 2.

Endemic Regime: If $R_{min}>1$, all solutions converge to the EE as $t\to\infty$. 3. 3.

Epidemic Regime: If $R_{max}>1>R_{min}$, and $R_{1}\leq 1$, the following hold:

1. (a)

All solutions converge to a point in $\mathcal{M}$ as $t\to\infty$. 2. (b)

There exists $H\subset\Delta_{N}$ and $H\supset\mathcal{M}_{+}$ that is foliated by families of heteroclinic orbits, each orbit connecting two points in $\mathcal{M}$. 4. 4.

Bistable Regime: If $R_{max}>1>R_{min}$ and $R_{1}>1$, then, depending on the initial conditions, solutions converge to a point in $\mathcal{M}$ or to the EE as $t\to\infty$.

Table III summarizes which regimes of Theorem 3 exist for each of the six cases of the network SIRI model. Figure 1 illustrates the four regimes of Theorem 3 near the IFE set $\mathcal{M}$ when $B\succeq\hat{B}$ or $\hat{B}\succeq B$.

In Case 2 (SIR), since $\hat{B}=\bar{\mathbf{0}}$, $R_{max}=R_{0}$ and $R_{min}=R_{1}=0$. So only the infection-free and epidemic regimes are possible. This corresponds to the line $R_{1}=0$ in Figure 1.

In Case 3 (SIS), since $B=\hat{B}$, $R_{max}=R_{0}=R_{1}=R_{min}$. So only the infection-free and endemic regimes are possible. This corresponds to the line $R_{1}=R_{0}$ in Figure 1.

In Case 4 (partial immunity), since $B\succ\hat{B}$, by Proposition 8, $R_{max}=R_{0}$ and $R_{min}=R_{1}$. So only the infection-free, endemic, and epidemic regimes are possible. This corresponds to the region $R_{1}<R_{0}$ in Figure 1.

In Case 5 (compromised immunity), since $\hat{B}\succ B$, by Proposition 8, $R_{max}=R_{1}$ and $R_{min}=R_{0}$. So only the infection-free, endemic, and bistable regimes are possible. This corresponds to the region $R_{1}>R_{0}$ in Figure 1.

In Case 6 (mixed immunity), all four regimes are possible.

The $N$-dimensional set $\mathcal{M}$ is illustrated in Figure 1 as a plane ($N=2$) for ease of visualization. The blue region represents $\mathcal{M}_{-}$ (the set of stable points in $\mathcal{M}$) and the red region represents $\mathcal{M}_{+}$ (the set of unstable points in $\mathcal{M}$). Black arrows illustrate the flow of solutions near $\mathcal{M}$. Theorems 1 and 2 prove which regions of $\mathcal{M}$ exist in each of the regimes. In the infection-free regime, $\mathcal{M}=\mathcal{M}_{-}$ and all solutions converge to the IFE set $\mathcal{M}$. In the endemic regime, $\mathcal{M}=\mathcal{M}_{+}$ and all solutions converge to the EE. In Cases 2 and 4 in the epidemic regime and in Case 5 in the bistable regime, $\mathcal{M}_{-}$ and $\mathcal{M}_{+}$ both exist and there is a unique hypersurface $\mathcal{M}_{0}$ shown as a black dashed line separating the two.

In Section VII we study the geometry of solutions near $\mathcal{M}$ and the stable manifold (green) and unstable manifold (magenta) of $\mathcal{M}_{0}$ in the epidemic regime of Cases 2 and 4 and the bistable regime of Case 5. These manifolds are included in Figure 1 and help illustrate how solutions can flow.

Proof of Theorem 3.

To prove 1, let $R_{max}\leq 1$. Then by definition $R_{1}\leq 1$. By Lemma 8, all solutions converge to a point in $\mathcal{M}$ as $t\to\infty$. By (III), the dynamics of $p^{I}_{\rm avg}$ are

[TABLE]

The inequality follows from (4), and the last equality follows from $\rho(D^{-1}\bar{B}_{max})=\rho(\bar{B}_{max}D^{-1})$. Since $\hat{B}$ is irreducible and by Proposition 1 $w_{m}\gg\mathbf{0}$, the nonlinear term $w_{m}^{T}D^{-1}P^{I}\hat{B}p^{I}$ is nonnegative. And since $B\succeq\hat{B}$ or $\hat{B}\succeq B$, by Proposition 8, $\bar{R}_{max}=R_{max}$. Therefore, if $R_{max}<1$, then $\dot{p}^{I}_{\rm avg}<0$ and $p^{I}_{\rm avg}$ decays monotonically to zero as $t\to\infty$.

If $R_{max}=1$, $\dot{p}^{I}_{\rm avg}\leq 0$ with equality holding only at points in $\Sigma=\{\mathbf{0}\preceq p^{I}\preceq\mathbf{1}\,|\,p^{I}_{j}>0{\rm~{}implies~{}}p_{k}^{I}=0,k\in\mathcal{N}_{j}\}$.

At any point in $\Sigma$, the dynamics of node $j$ where $p^{I}_{j}>0$ reduce to $\dot{p}^{I}_{j}=-\delta_{j}p_{j}^{I}$. Thus, no solution can stay in $\Sigma$ except for the trivial solution $p^{I}=\mathbf{0}$. By LaSalle’s Invariance Principle, $p^{I}_{\rm avg}$ decays monotonically to zero as $t\to\infty$.

To prove 2, let $R_{min}>1$. Then by definition $R_{1}>1$. By Theorem 1, all points in $\mathcal{M}$ are unstable. Therefore, no non-trivial solution can converge to a point in $\mathcal{M}$. By Lemma 8, it follows that the $\omega$-limit set of all solutions with $p^{I}(0)\succ\mathbf{0}$ is the EE. Therefore, all solutions converge to the EE as $t\to\infty$.

To prove 3, let $R_{max}>1>R_{min}$, and $R_{1}\leq 1$. By Theorem 1, $\mathcal{M}_{-},\mathcal{M}_{+},\mathcal{M}_{0}\neq\emptyset$. By Lemma 8, all solutions converge to a point in $\mathcal{M}$ as $t\to\infty$. By the proof of Lemma 1, the unstable manifold of any unstable point $x\in\mathcal{M}_{+}$, lies partially or entirely in $\Delta_{N}$. Let $y=y(t,y_{0})$ be a solution of (III) with initial condition $y_{0}$ on the unstable manifold of $x$. Then, $y$ converges to $x$ as $t\to-\infty$. By Lemma 8, $y$ converges to a point $x^{\prime}\in\mathcal{M}$, $x^{\prime}\neq x$, as $t\to\infty$. Thus, $y$ forms a heteroclinic orbit. Let $H\subset\Delta_{N}$ be the union of the unstable manifolds of all points $x\in\mathcal{M}_{+}$. Then, every solution in $H$ forms a heteroclinic orbit connecting two points in $\mathcal{M}$.

To prove 4, let $R_{1}>1$. By Proposition 9, the EE exists and it is locally stable. Therefore any solution in the region of attraction of the EE converges to the EE as $t\to\infty$. By Lemma 1, for any locally stable point $x\in\mathcal{M}_{-}$, there exists $V\subset\Delta_{N}$ and $x\in V$ such that any solution starting in $V$ converges to a point in $V\cap\mathcal{M}_{-}$ at an exponential rate. ∎

VII Bistable and Epidemic Regimes

VII-A Geometry of Solutions near $\mathcal{M}$

In this section we examine the geometry of solutions near the IFE set $\mathcal{M}$ in the epidemic regime for Case 2 (SIR) and Case 4 (partial immunity), and the bistable regime for Case 5 (compromised immunity). The bistable regime for the network SIRI model, which doesn’t exist for the well-studied SIS and SIR models, generalizes that proved for the well-mixed SIRI model studied in [27].

Definition VII.1** (Transversal crossing of $\Lambda_{c}$).**

Let $y(t)=(p^{S}(t),p^{I}(t))\in\Delta_{N}$, $t\geq 0$, be a solution of (III). We say that $y$ crosses $\Lambda_{c}$ transversally if $p^{S}$, the projection of $y$ onto $\mathcal{M}$, crosses $\Lambda_{c}$ transversally. This holds if there exists a time $t^{\prime}>0$ and $m\in\Lambda_{c}$ such that $p^{S}(t^{\prime})=m$ and $\dot{p}^{S}(t^{\prime})^{T}\nabla\Lambda(m)\neq 0$.

Proposition 11** (Transversal crossing direction).**

Let $y(t)=(p^{S}(t),p^{I}(t))\in\Delta_{N}$, $t\geq 0$, be a solution of (III) that crosses $\Lambda_{c}$ transversally at the point $(m,\mathbf{0})\in\mathcal{M}$ and time $t=t^{\prime}$, where $c=\Lambda(m)$. If $\dot{p}^{S}(t^{\prime})^{T}\nabla\Lambda(m)<0$, then $\Lambda$ decreases as $p^{S}$ crosses $\Lambda_{c}$, and if $\dot{p}^{S}(t^{\prime})^{T}\nabla\Lambda(m)>0$, then $\Lambda$ increases as $p^{S}$ crosses $\Lambda_{c}$. Suppose $(m,\mathbf{0})\in\mathcal{M}_{0}$, then if $\dot{p}^{S}(t^{\prime})^{T}\nabla\Lambda(m)<0$, $p^{S}$ crosses $\mathcal{M}_{0}$ from $\mathcal{M}_{+}$ to $\mathcal{M}_{-}$ and if $\dot{p}^{S}(t^{\prime})^{T}\nabla\Lambda(m)>0$, then $p^{S}$ crosses $\mathcal{M}_{0}$ from $\mathcal{M}_{-}$ to $\mathcal{M}_{+}$.

Proof.

Since $\dot{p}^{S}=-P^{S}Bp^{I}$, then $\dot{p}^{S}(t^{\prime})^{T}\nabla\Lambda(m)$ is the derivative of $\Lambda$ at $m$ along solutions of (III) (see Proposition 4). If $\dot{p}^{S}(t^{\prime})^{T}\nabla\Lambda(m)<0$, $\Lambda$ decreases as $p^{S}(t)$ crosses $\Lambda_{c}$ and if $\dot{p}^{S}(t^{\prime})^{T}\nabla\Lambda(m)>0$, $\Lambda$ increases as $p^{S}(t)$ crosses $\Lambda_{c}$. Suppose $(m,\mathbf{0})\in\mathcal{M}_{0}$. By Lemma 4, if $\dot{p}^{S}(t^{\prime})^{T}\nabla\Lambda(m)<0$, $p^{S}$ crosses $\mathcal{M}_{0}$ from $\mathcal{M}_{+}$ to $\mathcal{M}_{-}$. Similarly, if $\dot{p}^{S}(t^{\prime})^{T}\nabla\Lambda(m)>0$, $p^{S}$ crosses $\mathcal{M}_{0}$ from $\mathcal{M}_{-}$ to $\mathcal{M}_{+}$. ∎

Theorem 4** (Transversality of solutions).**

Consider Cases 2 and 4 in the epidemic regime ($R_{0}>1$, $R_{1}<1$) and Case 5 in the bistable regime ($R_{0}<1$, $R_{1}>1$). Assume no stubborn agents. Let $y(t)=(p^{S}(t),p^{I}(t))\in\Delta_{N}$, $t\geq 0$, be a solution of (III) for which there exists a time $t^{\prime}>0$ and $(m,\mathbf{0})\in\textrm{int}(\mathcal{M})$, such that $p^{S}(t^{\prime})=m$ and $p^{I}(t^{\prime})\succ\mathbf{0}$. Let $c=\Lambda(m)$ so that $m\in\Lambda_{c}$. Then $y$ crosses $\Lambda_{c}$ transversally. Suppose $(m,\mathbf{0})\in\textrm{int}(\mathcal{M}_{0})$. In the epidemic regime of Cases 2 and 4, $p^{S}$ crosses $\mathcal{M}_{0}$ from $\mathcal{M}_{+}$ to $\mathcal{M}_{-}$, and the stable and unstable manifolds of $\mathcal{M}_{0}$ lie outside $\Delta_{N}$. In the bistable regime of Case 5, $p^{S}$ crosses $\mathcal{M}_{0}$ from $\mathcal{M}_{-}$ to $\mathcal{M}_{+}$, and the stable and unstable manifolds of $\mathcal{M}_{0}$ lie inside $\Delta_{N}$.

Proof.

Assume no stubborn agents. Let $t^{\prime}>0$ such that $p^{S}(t^{\prime})=m$, $p^{I}(t^{\prime})\succ\mathbf{0}$ and $(m,\mathbf{0})\in\textrm{int}(\mathcal{M}$). By Definition 6, $m\gg\mathbf{0}$. Let $c=\Lambda(m)$. So $\dot{p}^{S}(t^{\prime})=-\textrm{diag}(m)Bp^{I}(t^{\prime})\prec\mathbf{0}$. By Lemma 5, if $B\succ\hat{B}$ (Cases 2 and 4) then $\nabla\Lambda\gg\mathbf{0}$ and if $\hat{B}\succ B$ (Case 5) then $\nabla\Lambda\ll\mathbf{0}$. Thus $\dot{p}^{S}(t^{\prime})^{T}\nabla\Lambda(m)\neq 0$, and so, by Proposition 11, $y$ crosses $\Lambda_{c}$ transversally. Suppose $(m,\mathbf{0})\in\textrm{int}(\mathcal{M}_{0}$). Then, by Proposition 11, $p^{S}$ crosses from $\mathcal{M}_{+}$ to $\mathcal{M}_{-}$ in Cases 2 and 4, and from $\mathcal{M}_{-}$ to $\mathcal{M}_{+}$ in Case 5. By continuity of solutions with respect to initial conditions, it follows that the stable and unstable manifolds of $\mathcal{M}_{0}$ must lie outside $\Delta_{N}$ in the epidemic regime of Cases 2 and 4 and inside $\Delta_{N}$ in the bistable regime of Case 5, as illustrated in Figure 1. ∎

Corollary 1**.**

In the epidemic regime of Cases 2 and 4, every heteroclinic orbit in $\Delta_{N}$ connects a point in $\mathcal{M}_{+}$ to a point in $\mathcal{M}_{-}$.

Proof.

Since by Theorem 4 the stable manifold of any point in $\mathcal{M}_{0}$ lies outside $\Delta_{N}$, it follows by Theorem 3 that the orbit connects a point in $\mathcal{M}_{+}$ to a point in $\mathcal{M}_{-}$. ∎

Corollary 2**.**

Consider Case 5 in the bistable regime. Let $y(t)=(p^{S}(t),p^{I}(t))$ be a solution of (III). Then it holds that

• •

If $y$ crosses $\mathcal{M}_{0}$ transversally or $(p^{S}(0),\mathbf{0})\in\mathcal{M}_{+}$, then $y$ converges to the EE as $t\to\infty$. Moreover, the EE lies on the unstable manifold of $\mathcal{M}_{0}$.

• •

The stable manifold of $\mathcal{M}_{0}$ intersects the boundary of $\Delta_{N}$ where $p^{S}=\mathbf{1}-p^{I}$.

Proof.

If $y$ crosses $\mathcal{M}_{0}$ transversally at $t=t^{\prime}$, then, by Theorem 4, $p^{S}$ crosses every level surface $\Lambda_{c}$ in $\mathcal{M}$ transversally and crosses $\mathcal{M}_{0}$ from $\mathcal{M}_{-}$ to $\mathcal{M}_{+}$. By Proposition 11, $\Lambda$ strictly increases along $p^{S}$. It follows by Lemma 4 that $(p^{S},\mathbf{0})\in\mathcal{M}_{+}$ for all $t>t^{\prime}$. Since by Lemma 1 $y$ cannot converge to a point in the IFE subset $\mathcal{M}_{+}$, it follows by Lemma 8 that $y$ converges to the EE as $t\to\infty$. The same argument holds if $(p^{S}(0),\mathbf{0})\in\mathcal{M}_{+}$. It follows that the EE lies on the unstable manifold of $\mathcal{M}_{0}$.

Consider any point $y=(p^{S},p^{I})\in\Delta_{N}$ on the stable manifold of $\mathcal{M}_{0}$. Then, as $t\to-\infty$ and $y(t)\in\Delta_{N}$, $y$ remains on the stable manifold of $\mathcal{M}_{0}$ and $(p^{S},\mathbf{0})$ remains in $\mathcal{M}_{-}$. By Lemma 1, points in $\mathcal{M}_{-}$ have no unstable manifold in $\Delta_{N}$ so the stable manifold of $\mathcal{M}_{0}$ cannot intersect $\mathcal{M}_{-}$. Instead, because the components of $p^{S}$ increase monotonically as $t\to-\infty$, $y$ intersects $\partial\Delta_{N}$ where $p^{S}=\mathbf{1}-p^{I}$. ∎

The locations of the stable and unstable manifolds of $\mathcal{M}_{0}$, as proved in Theorem 4, are illustrated in Figure 1: outside $\Delta_{N}$ in the epidemic regime and inside $\Delta_{N}$ in the bistable regime. The figure shows the heteroclinic orbits proved in Corollary 1 for the epidemic regime. The solutions along the heteroclinic orbits cross $\mathcal{M}_{0}$ transversally, with their projection onto $\mathcal{M}$ crossing from $\mathcal{M}_{+}$ to $\mathcal{M}_{-}$, as proved in Theorem 4.

Figure 1 also shows the local flow in the bistable regime, as proved in Corollary 2. These solutions also cross $\mathcal{M}_{0}$ transversally, but with their projection onto $\mathcal{M}$ crossing from $\mathcal{M}_{-}$ to $\mathcal{M}_{+}$, as proved in Theorem 4.

VII-B Bistability and Resurgent Epidemic

In this section we examine the bistable regime of Theorem 3, which exists for Case 5 (compromised immunity) and for Case 6 (mixed immunity). We show how solutions in this regime can exhibit a resurgent epidemic in which an initial infection appears to die out for an arbitrarily long period of time, but then abruptly and surprisingly resurges to the EE.

Conditions on the initial state that predict a resurgent epidemic were proved for the well-mixed SIRI model in [27]. Here we compute a critical condition on the initial state in the special case that $\mathcal{G}$ is a $d$-regular digraph, i.e., every agent $j$ in $\mathcal{G}$ has degree $d_{j}=d$, and every node has the same initial state. We then illustrate numerically for a more general digraph with compromised immunity in Figure 2.

Consider a $d$-regular digraph with global recovery, infection, and reinfection rates: $D=\delta\mathbb{I}$, $B=\beta A$, and $\hat{B}=\hat{\beta}A$, where $A$ is the adjacency matrix. The network SIRI dynamics (2) are

[TABLE]

where we have used the identity $\sum_{k=1}^{N}a_{jk}p^{I}_{k}=p^{I}_{j}d+\sum_{k=1}^{N}a_{jk}(p^{I}_{k}-p^{I}_{j})$. Let $p^{I}(0)=p_{ic}\mathbf{1}$. Then (18) reduce to

[TABLE]

(VII-B) describes identical and uncoupled dynamics for every agent $j$, which are equivalent to the dynamics of the well-mixed SIRI model [27] with infection rate $\beta d$ and reinfection rate $\hat{\beta}d$. Following [27], we find the critical initial condition $p_{crit}=1-\xi(R_{0}d\xi)^{-\beta/\hat{\beta}}$, where $\xi=(R_{1}-1/d)/(R_{1}-R_{0})$. If $p_{ic}<p_{crit}$ solutions converge to a point in the IFE as $t\to\infty$. If $p_{ic}>p_{crit}$ solutions converge to the EE, $p^{I*}=(1-\delta/(\hat{\beta}d))\mathbf{1}$, as $t\to\infty$. If $p_{ic}=p_{crit}$, the solution flows along the stable manifold of the point $\xi\mathbf{1}\in\mathcal{M}_{0}$ and converges to $\xi\mathbf{1}$. These results suggest more generally that the stable manifold of $\mathcal{M}_{0}$ separates solutions that converge to the IFE from those that converge to the EE, as in [27].

Further, as in [27], if $p_{ic}>p_{crit}$, the solution exhibits a resurgent epidemic in which $p^{I}_{j}$ initially decreases to a minimum value $p^{I}_{min}$ and then increases to the EE. As $(p_{ic}-p_{crit})\to 0$, $p^{I}_{min}\to 0$ and the time it takes for the solution to resurge goes to infinity. That is, the infection may look like it is gone for a long time before resurging without warning. We note that the SIR and SIS models do not admit a bistable regime and therefore fail to account for the possibility of a resurgent epidemic. This implies the SIR and SIS models are not robust to variability in infection and reinfection rates.

Figure 2 illustrates the bistability and resurgent epidemic phenomena for an example network of $N=4$ agents for Case 5 (compromised immunity). $B=0.7A$, $\hat{B}=\textrm{diag}([1.5,0.7,0.7,0.7])A$, and $D=\mathbb{I}.$ That is, agents $2,3,4$ are stubborn and agent 1 acquires compromised immunity to all its infected neighbors (agents 2 and 4). So $R_{max}=R_{1}=1.28$ and $R_{min}=R_{0}=0.85$, placing the system in the bistable regime. In both panels of Figure 2, $p^{S}(0)=\mathbf{1}-p^{I}(0)$. In the left panel $p^{I}(0)=[0,0.05,0.1,0]^{T}$ and the solution can be observed to converge to the IFE, i.e., $p^{I}_{j}\to 0$ for all $j$. In the right panel $p^{I}(0)=[0,0.08,0.1,0]^{T}$ and there is a resurgent epidemic: each $p^{I}_{j}$ initially decays and then remains close to zero until $t\approx 100$, after which the $p^{I}_{j}$ increase rapidly to the EE, which is $p^{I*}=[0.29,0.11,0.17,0.17]^{T}$.

VIII Control Strategies

We apply our theory to design control strategies for technological, as well as biological and behavioral settings, that guarantee desired steady-state behavior, such as the eradication of an infection. We begin with an example network with mixed immunity in the endemic regime, which has an infected steady state. We show three strategies for changing parameters that modify the reproduction numbers $R_{0}$, $R_{1}$, $R_{min}$, and $R_{max}$, and control the dynamics to a behavioral regime that results in an infection-free steady state, according to Theorem 3. We then consider two example networks with mixed immunity in the bistable regime. We illustrate how vaccination of well-chosen agents increases the set of initial conditions that yield an infection-free steady state or at least delay a resurgent epidemic so that further control can be introduced.

VIII-A Control from Endemic to Infection-Free Steady State

Consider the network of four agents shown in the top left panel of Figure 3. Let $B=A$, the unweighted adjacency matrix for the digraph shown. Let $\hat{B}=\textrm{diag}([0.3,1,2,1])A$ and $D=\mathbb{I}$. The network has weak mixed immunity (Case 6a): agent 1 acquires partial immunity to reinfection, agent 3 acquires compromised immunity to reinfection, and agents 2 and 4 are stubborn. We compute the reproduction numbers: $R_{0}=1.32$, $R_{1}=1.22$, $R_{min}=1.13$, and $R_{max}=1.52$, which, by Theorem 3, imply dynamics in the endemic regime and an infected steady state for every initial condition. This is illustrated in the simulation of $p^{I}_{j}$ versus time $t$ for initial conditions $p^{S}(0)=\mathbf{1}-p^{I}(0)$, $p^{I}(0)=[0.01,0.01,0.01,0.2]^{T}$. By (5), the solution converges to the EE: $p^{I*}=[0.07,0.23,0.12,0.19]^{T}$.

VIII-A1 Modification of agent recovery rate

This strategy controls the network behavior by selecting one or more agents for treatment to increase its recovery rate. In the epidemiological setting, this could mean medication. In the behavioral setting, this could mean providing incentives or training. We make just one modification to the example network of four agents: $\delta_{2}=1$ becomes $\delta_{2}=3.5$. The corresponding reproduction numbers are $R_{0}=0.80$, $R_{1}=0.72$, $R_{min}=0.65$, and $R_{max}=0.94$, which, by Theorem 3, imply dynamics in the infection-free regime. Using the same initial conditions as in the top left panel, we simulate the modified system in the bottom left panel of Figure 3. The solution converges to an infection-free steady state as predicted by Theorem 3.

VIII-A2 Modification of agent reinfection rate

This strategy controls the network behavior by selecting one or more agents for treatment to decrease its reinfection rate. In the epidemiological setting, this is vaccination, and in the behavioral setting, it could be inoculation as in psychology research [2]. We make just one modification to the example network of four agents: $\hat{\beta}_{42}=1$ becomes $\hat{\beta}_{42}=0.3$. The corresponding reproduction numbers are $R_{0}=1.32$, $R_{1}=0.96$, $R_{min}=0.82$, and $R_{max}=1.52$, which, by Theorem 3, imply dynamics in the epidemic regime. Using the same initial conditions as in the top left panel, we simulate the modified system in the top right panel of Figure 3. After a small and short-lived epidemic, the solution converges to an infection-free steady state, as predicted by Theorem 3.

VIII-A3 Modification of network topology

This strategy controls the network behavior by selecting one or more edges in the network graph for re-wiring. In all settings, this means affecting who comes in contact with whom. For the example network, we move the connections between agents 4 and 1 to be between agents 4 and 2. The corresponding reproduction numbers are identical to those in the modification of reinfection rate example and therefore, by Theorem 3, also imply dynamics in the epidemic regime. Using the same initial conditions as in the top left panel, we simulate the modified system in the bottom right panel of Figure 3. The solution converges to an infection-free steady state, as predicted by Theorem 3, with an initial epidemic smaller than in the top right panel.

VIII-B Control in Bistable Regime

VIII-B1 Small network

A network of two agents with weak mixed immunity is shown in Figure 4. Agent 1 acquires compromised immunity while agent 2 acquires partial immunity: $\beta_{12}=\hat{\beta}_{21}=0.8$, $\hat{\beta}_{12}=\beta_{21}=1.3$, $\delta_{1}=\delta_{2}=1$. The corresponding reproduction numbers are $R_{0}=R_{1}=1.02$, $R_{min}=0.8$, and $R_{max}=1.3$. By Theorem 3, the system is in the bistable regime. Suppose initially there are no recovered agents, i.e., $p^{S}(0)=\mathbf{1}-p^{I}(0)$. Then, it can be shown that the solution will always converge to the EE. Now suppose we apply a control strategy in which we vaccinate the agent who acquires partial immunity, where vaccination is equivalent to exposing the agent to the infection. After the vaccination of agent 2 in our example, we have $p^{S}_{1}(0)=1-p^{I}_{1}(0)$ with $p^{I}_{1}(0)\in[0,1]$ and $p^{S}_{2}(0)=0$ with $p^{I}_{2}(0)\in[0,1]$.

We illustrate the results of the vaccination of agent 2 in the left panel of Figure 4. Initial conditions that lead to the EE are shown in magenta and to an infection-free steady state in yellow. We illustrate with simulations using the initial conditions $p^{I}_{1}(0)$ and $p^{I}_{2}(0)$ denoted by the star in the left panel of Figure 4. The top right simulation is of the system with no vaccination: there is an infected steady state as predicted. The bottom right simulation is of the system with agent 2 vaccinated: there is an infection-free steady state as predicted.

VIII-B2 Large network

A network of twenty agents with weak mixed immunity is shown in Figure 5. Agents 2, 3, 5, 7, 11, 13, 17, and 19 (dark gray) acquire partial immunity to reinfection: $\beta_{jk}=0.5$, $\hat{\beta}_{jk}=0.1$. Agents 1 and 20 (gray) are stubborn: $\beta_{jk}=\hat{\beta}_{jk}=0.5$. Agents 4, 6, 8, 9, 10, 12, 14, 15, 16, and 18 (light gray) acquire compromised immunity to reinfection: $\beta_{jk}=0.5$, $\hat{\beta}_{jk}=0.875$. $D=\mathbb{I}$. The corresponding reproduction numbers are $R_{0}=1.06$, $R_{1}=1.26$, $R_{min}=0.76$, and $R_{max}=1.59$. By Theorem 3, the system is in the bistable regime.

In Figure 5, we plot simulations of the average infected state $\bar{p}^{I}=\mathbf{1}^{T}p^{I}/N$ versus $t$ for no vaccinations (solid black) and for different sets of vaccinated agents: agents 2, 3, 5, and 7 (dotted blue), agent 11 (dashed orange), agents 7 and 11 (point-dashed green), and agents 7, 11, and 13 (dashed violet). Vaccinating agent 11 has the strongest effect. The infection appears to be eradicated by vaccinating agents 7, 11, and 13. The other vaccination cases delay the epidemic, which provides time for treatments or other control interventions.

IX Final Remarks

The network SIRI model generalizes the network SIS and SIR models by allowing agents to adapt their susceptibility to reinfection. We have proved new conditions on network structure and model parameters that distinguish four behavioral regimes in the network SIRI model. The conditions depend on four scalar quantities, which generalize the basic reproduction number used in the analysis and control of the SIS and SIR models. The SIRI model captures dynamic outcomes that the SIS and SIR models necessarily miss, most dramatic of which is the resurgent epidemic and the bistable regime.

The generality of the network SIRI model provides a means to assess robustness to uncertainty and changes in infection rates. Further, the model provides new flexibility in control design, including, as we have illustrated, selective vaccination of agents that acquire partial immunity or modification of recovery rates, reinfection rates, and wiring, to leverage what can happen when agents adapt their susceptibility.

Control strategies from the literature can also be extended to the SIRI setting, for example, optimal node removal, optimal link removal, and budget-constrained allocation [10, 46, 47, 48]. However, optimal node and link removal have been shown to be NP-complete and NP-hard problems, respectively [10, 49] and all of these strategies rely on knowledge of the network structure or system parameters. A practical alternative is to design feedback control strategies whereby agents actively adapt their susceptibility.

Our results were obtained using the IBMF approach [15, 9] to reduce the size of the state space. It has been shown that if solutions in an IBMF model converge to an infection-free equilibrium, then the expected time it takes the stochastic Markov model to reach the infection-free absorbing state is sublinear with respect to network size [10]. It is an open question if bistability in the stochastic Markov model will be sensitive to noise. Other questions to be explored include deriving centrality measures that facilitate optimal control design and extending to the case in which agents adapt their susceptibility to every reinfection.