A Low-Dimensional Network Model for an SIS Epidemic: Analysis of the Super Compact Pairwise Model

Network-based models of epidemic spread have become increasingly popular in recent decades. Despite a rich foundation of such models, few low-dimensional systems for modeling SIS-type diseases have been proposed that manage to capture the complex dynamics induced by the network structure. We analyze one recently introduced model and derive important epidemiological quantities for the system. We derive the epidemic threshold and analyze the bifurcation that occurs, and we use asymptotic techniques to derive an approximation for the endemic equilibrium when it exists. We consider the sensitivity of this approximation to network parameters, and the implications for disease control measures are found to be in line with the results of existing studies.


Introduction
In the past few decades, network-based models of epidemic spread have become a central topic (Kiss et al. 2017;Pastor-Satorras et al. 2015) in epidemiology. Their ability to capture mathematically the complex structure of transmission interactions makes them an invaluable theoretical paradigm. Mathematically, a network is modeled as a graph consisting of a set of nodes that are connected by a set of links (called edges). In the context of epidemiology, typically nodes represent individuals, and edges represent interactions that can transmit the infection. Used in conjunction with compartment models, the disease natural history determines the number of possible states an individual node might be in at any point in time. When disease spread is modeled as a continuous time Markov chain, the network size and disease natural history can lead to high dimensional state spaces. For example, in a network with N nodes where individual nodes can be in m possible states, the size of the state space for the network is m N . Efforts to describe this process with a system of ordinary differential equations are similarly hampered by size-the Kolmogorov equations governing this system are exact, but prohibitively large. Thus, an important goal in network-based modeling has been to find a (relatively) low-dimensional system that accurately approximates the underlying high-dimensional system.
Many approaches (Pastor-Satorras and Vespignani 2001;Pastor-Satorras et al. 2015;Miller et al. 2012; Barnard et al. 2019;Karrer and Newman 2010) in recent years have sought to introduce models with systems of a manageable size. Pairwise models (Keeling 1999;Eames and Keeling 2002;House and Keeling 2011) have been a popular and fruitful approach to this question. Derived from the Kolmogorov equations and exact in their unclosed form (Taylor et al. 2012), pairwise models consider the evolution of not just the expected number of nodes in a given state, but also pairs and triples of nodes. The dynamical variables are of the form [A] (the expected number of nodes in state A), [AB] (the expected number of pairs in state A − B), and [ABC] (the expected number of triples in state A − B − C). Higher-order groupings of nodes are also considered but rarely written, as dimension-reduction efforts often focus on approximating the expected number of triples in terms of pairs and individual nodes. Pairwise models have been successful with a variety of different network types, with models developed for networks with heterogeneous degree (Eames and Keeling 2002), weighted networks (Rattana et al. 2013), directed networks (Sharkey et al. 2006), and networks with motifs (House et al. 2009;Keeling et al. 2016) to name a few. Moreover, pairwise models have been developed for a variety of disease natural histories, with particular focus on SIR (susceptible-infectious-recovered) and SIS (susceptible-infectious-susceptible) models.
In this paper, we consider an SIS pairwise model for networks with heterogeneous degree. SIS dynamics are used to model diseases where no long-term immunity is conferred upon recovery, leading to their frequent application to sexually transmitted infections such as chlamydia or gonorrhea (Eames and Keeling 2002). Contact networks for diseases of this type frequently involve heterogeneity in the number of contacts for individuals, and thus node degree becomes an essential concept. The degree of a node in a network is the number of edges to which the node is connected, and thus the number of potential infectious contacts. In this way, heterogeneous networks can capture complex disease dynamics. An essential tool when working with such networks is the degree distribution, defined by p k which is the probability a randomly selected node has degree k. The degree distribution has played an important role in dimension reduction approximations for pairwise models.
For the SIR-type diseases, accurate low-dimensional models have been derived from the pairwise family using probability generating functions (Miller et al. 2012), complete with conditions for finding the final size of the epidemic. Despite the successes of the SIR case, the dimension reduction techniques in Miller et al. (2012) do not apply to the SIS case. Instead, the development of low-dimensional models of SIS-type disease spread on networks has relied on moment closure approximations. Under the assumption of a heterogeneous network with no clustering, House and Keeling (2011) introduced an approximation reducing the system size from O(N 2 ) to O(N ), where N is the number of nodes in the network. Termed the compact pairwise model (CPW), it has shown good agreement with stochastic simulations despite its considerably smaller size. However, the number of model equations still grows as the maximum degree of the network, making its application challenging for large networks with significant degree heterogeneity. Perhaps the most successful model in reducing the number of equations of the CPW for SIS-type diseases is the super compact pairwise model (SCPW) (Simon and Kiss 2016). The system consists of only four equations, with network structure being encoded to the model through the first three moments of the degree distribution. While Simon and Kiss demonstrated excellent agreement between the CPW and the SCPW, bifurcation analysis of the model and an explicit formula for the endemic steady state remains to be done.
This paper sets out on that analysis of the SCPW model. A common point of investigation among models of SIS-type diseases is the disease-free equilibrium (DFE) that loses stability as a relevant parameter passes a critical value known as the epidemic threshold Vespignani 2001, 2002;Boguñá and Pastor-Satorras 2002). The epidemic threshold serves as a dividing point between two qualitatively different types of outbreaks. Below the epidemic threshold, any outbreak will die out; above the epidemic threshold, the system converges asymptotically to a stable equilibrium where the disease remains endemic in the population. Many studies follow the "next generation matrix" approach for the basic reproduction number R 0 (van den Driessche and Watmough 2002) to characterize the epidemic threshold. We follow a more conventional bifurcation analysis to derive the epidemic threshold and offer a proof that the system undergoes a transcritical bifurcation, as one might expect. Perhaps more importantly, the SCPW's small fixed number of equations presents an excellent opportunity to investigate the endemic equilibrium for SIS models on heterogeneous networks, which has been heretofore inhibited by large system size. We present a novel asymptotic approach to approximating the endemic equilibrium, leveraging the low-dimensionality of the model. The results presented further our understanding of the SCPW model specifically and suggest potential new avenues in the challenging problem of analytically determining the nontrivial steady state of pairwise models of SIS-type diseases.
The paper is structured as follows: in Sect. 2, we nondimensionalize the model and reduce the number of equations to 3 to facilitate computations. In Sect. 3, we derive the epidemic threshold and show that the system undergoes a forward transcritical bifurcation. In Sect. 4, we tackle the endemic steady state that emerges through the bifurcation. We use asymptotic methods to approximate the size of the endemic steady state under two regimes-the system near the epidemic threshold and the system far away from the epidemic threshold-and give examples of the efficacy of these approximations on prototypical networks. Finally, we examine the implications of these two approximations. In line with existing studies (Eames and Keeling 2002), we find that control measures for reducing the prevalence at the endemic equilibrium may require different tactics depending on the regime.

Model
Pairwise models of SIS-type diseases provide a network-based analog of the classical SIS model (Diekmann and Heesterbeek 2000).The essential characteristics of pairwise models of SIS epidemics are dynamical equations for not just the expected number of nodes in each state, but also pairs and triples of nodes. At the node level, The CPW closes the system by approximating the expected number of triples as where S 1 and S 2 are the first and second moments of the distribution of susceptible nodes; that is  Simon and Kiss (2016) is given aṡ where k n is the nth moment of the degree distribution, τ is the transmission rate, and γ is the recovery rate. Here, the quantity Q serves as an approximation of (S 2 − S 1 )/S 2 1 . As well, the quantities With the goal of performing bifurcation and asymptotic analyses in mind, nondimensionalizing the SCPW model is a natural first step. To do so, we will rearrange the Eqs.
(3)-(5) so that the network parameters k , k 2 , k 3 are consolidated into more workable constants. First, we rewrite Q as where As well, a natural rescaling of time is T = t/γ , which prompts the defining of the transmission-recovery rate ratio δ = τ/γ. The introduction of δ consolidates the two epidemiological parameters τ and γ into a single nondimensional parameter, so any changes to epidemiology of the disease will be captured in δ alone. With these substitutions, the system (1)-(5) becomesv where the dot notation represents the derivative with respect to the nondimensional time variable d dT . The conservation Eqs. (6) and (7) become respectively. At this point, the conservation equations can be used to reduce the system to a mere 3 equations. However, the elimination of different equations for different analyses will be convenient. For characterizing the bifurcation undergone by the disease-free equilibrium (DFE), it is convenient to work with variables that are 0 at the DFE. For approximating the endemic steady state using asymptotic methods, the most parsimonious equations will make the algebraic manipulation required easier. Thus, we will work with slightly different (but equivalent) characterizations of (10)-(14) in the sections that follow.

Epidemic Threshold
To derive the epidemic threshold, we consider the stability of the DFE in terms of the epidemiological parameter δ. We will show that as δ increases through a critical value δ c , the DFE loses stability. Typically as the DFE loses stability, an asymptotically stable endemic equilibrium emerges. The SCPW is no exception, and here we derive the epidemic threshold, with a proof that the system undergoes a transcritical bifurcation (and thus an endemic equilibrium emerges) when δ = δ c included in Appendix A.
First, we use the conservation Eqs. (15) and (16) to eliminate Eqs. (10) and (13). The resulting system iṡ Though ostensibly a messier choice of equation reduction, we note that at the DFE, will be convenient moving forward. To determine the stability of the DFE, we compute the Jacobian at x = 0 : A straightforward computation shows that We can write DF as a block triangular matrix as where the dimensions A and B, respectively, are 1 × 2 and 2 × 2. The properties of determinants of block matrices tell us that the eigenvalues of DF are −1 and the eigenvalues of B, which will determine the stability of the DFE. We appeal here to the trace-determinant theorem, which tells us the eigenvalues ξ of the 2 × 2 matrix B are given by First, we observe that these eigenvalues are real, as which is clearly positive. As a consequence, for the DFE to be stable we must have Tr(B) < 0 and Det(B) > 0. The determinant can be written and is thus positive if and only if δ < 1/k. Moreover, if δ < 1/k, then Therefore, we conclude that the DFE is stable for δ < 1/k and unstable for δ > 1/k. Thus, the epidemic threshold is the critical value of the bifurcation parameter δ : Notably, this threshold value is identical to that of the CPW as shown in Kiss et al. (2017). However, it remains to be shown that a bifurcation actually does occur here, and that an asymptotically stable endemic steady state emerges. To prove this, we apply a theorem of Castillo-Chavez and Song (2004) in Appendix A. We note that both the CPW and SCPW models are approximations to the true SIS dynamics on a network, so while (25) is a good approximation of the true epidemic threshold, it may not be appropriate in some cases. For instance, (25) is greater than zero for networks with a power law degree distribution ( p k ∼ k −d ) with d > 3 in the large network limit (N → ∞). However, exact results show that the true epidemic threshold is zero in the large network limit (Chatterjee and Durrett 2009).

The Endemic Equilibrium
With the existence of an endemic steady state established, we turn to the question of finding an approximate analytic expression for it. In general, this is a difficult proposition with epidemic models on networks owing to the frequently high-dimensional nature of the dynamical systems. An exact closed-form expression for the endemic equilibrium of the SCPW model requires solving a system of polynomial equations in multiple variables, which we do not attempt here. However, with asymptotic techniques, a workable approximation can be derived for two cases of δ: near the epidemic threshold (δ ≈ δ c ), and far away from it (δ >> δ c ). We do not have a good approximation in the intermediate case. Two challenges are apparent. First, how to eliminate equations to facilitate asymptotic expansions of the equilibrium and second, the choice of small nondimensional parameter in each case. Unlike in Sect. 3, the most parsimonious characterization of (10)- (14) is desirable. So we eliminate (11) and (14) with the conservation equations. To promote the finding of a small nondimensional parameter, we rewrite the resulting system using δ = δ c · δ δ c and incorporate the constants σ = k δ c , λ = αδ c / k , μ = βδ c . With these alterations, the system becomeṡ At the endemic equilibrium,v =ẋ =ẏ = 0. We can solve (26) for v and substitute into (27) and (28). With some rearrangement of terms and a little algebra (and adding (28) to (27)), we arrive at the system of polynomial equations that determines the endemic steady state: Note that in (30), we have dropped a factor of x that corresponds to the DFE. For the endemic steady state, we are interested in knowing the prevalence when the system is at equilibrium: w * . We use the following procedure to approximate the solution.
1. Express δ c /δ in terms of a small parameter. 2. Use the Implicit Function Theorem to linearize P(x, y) = 0 as around a point (x,ỹ) that is mathematically and/or biologically justified for the given regime. 3. Expand x, y, and other relevant quantities in terms of the small parameter. 4. Substitute the expansions into Q(x, y) = 0 and obtain a regular perturbation problem and find an asymptotic solution for the equilibrium value x, which approximates x * . 5. Apply the relation w * = (δ c /δ) −1 σ x * to obtain an asymptotic series for the prevalence at the endemic equilibrium.
We describe the results of this procedure for each case in the remainder of this sectionthe details of the computations are included in Appendix B.

Case 1: Near the Epidemic Threshold (ı ≈ ı c )
For δ ≈ δ c , we choose η = 1 − δ c /δ as a small parameter. In terms of this small parameter, (29) and (30) become: When δ ≈ δ c , an endemic steady state has just emerged, so we can view this equilibrium as a small perturbation to the steady state x = 0, y = 1. Linearizing P(x, y) = 0 about this point gives we have Substituting into (32) and equating coefficients to 0, we find an η-order expansion of the approximate equilibrium value x * as Using the relation w To demonstrate the efficacy of this approximation, we compare the approximation (38) to the actual endemic equilibrium using bifurcation diagrams (Fig. 1). We consider two example configuration model random networks (Molloy and Reed 1995) with N = 10, 000. In Fig. 1a, a bimodal network is considered with 5000 degree 3 nodes and 5000 degree 5 nodes. In Fig. 1b, a network with a Poisson degree distribution (with average degree k = 10) is considered. As is clear in both examples, the agreement between the actual and approximate endemic equilibrium is quite good near the epidemic threshold. Interestingly, the approximate value of w * is greater than the exact value for the bimodal network and less than the exact value for the Poisson network. We suspect that this is due to network structure and higher order terms in the asymptotic expansion, which we have not computed. An analogous situation is found in the δ >> δ c case.

Case 2: Far Away from the Epidemic Threshold (ı >> ı c )
For δ >> δ c , our small parameter of choice is = δ c /δ. We can rewrite (29) and (30) in terms of this parameter: When δ >> δ c , the transmission rate τ is large relative to the recovery rate γ. Thus, we expect the disease to affect much of the population, and consequently there will be very few remaining [SS] links, and therefore y ≈ 0. Solving P(φ, 0) = 0 for φ yields φ( ) = 2 − λ and slope of the linearization is then Next, we seek to expand y in terms of only. The relevant expansions for φ, ψ, and x are: To ease the writing of coefficients, we let φ α and ψ α refer to the coefficients on α for the respective series. From this, it follows that Substituting into (40), and equating the coefficients to 0, we find that we need the coefficients up to order 4 in order to find a 2 order expansion of the approximate equilibrium value of x * . The result is Finally, as w * = σ −1 x * , we arrive at an −order approximation for size of the endemic steady state as As with the δ ≈ δ c case, we compare the approximation (49) to the actual endemic equilibrium in Fig. 2 for the same networks as previously described. Again, the agreement is quite good, even for relatively small values of δ. In this case, the approximation for the endemic equilibrium also provides an approximation to the epidemic threshold. Whether this approximation is an overestimate or underestimate of the exact threshold depends on network structure. If k 2 ≥ k 2 + k , the approximation is an overestimate. On the other hand, if k 2 < k 2 + k , the approximation being an overestimate or underestimate depends on the relationship between k 3 and the other two moments.

Sensitivity Analysis
With any model of infectious disease, its implications in preventing or mitigating spread should be considered. For network models, some pharmaceutical and nonpharmaceutical interventions can alter the contact network structure in the effort to contain or mitigate outbreaks (Salathé and Jones 2010). For an SIS-type disease, particularly when containment is impossible, one such goal may be to decrease the size of the endemic equilibrium. To that end, we examine the sensitivity of our approximations of w * to network parameters in the SCPW model. One benefit of explicit asymptotic expressions for the endemic equilibrium is that sensitivity analyses are straightforward to implement. For a fixed δ, we have a three-dimensional parameter space. To visualize these parameter combinations, we use two-dimensional heat maps taken at slices of the third network parameter. In this case, we have decided to look at several fixed values of k 3 and draw sensitivity heat maps in the variables ( k , k 2 ). Further complicating matters is the fact that moments of a distribution are subject to many inequalities which restrict the domain of the sensitivity heat maps. Two natural restrictions to include are the results of Jensen's Inequality and the Cauchy-Schwarz Inequality, respectively: For a fixed value of k 3 , these restrictions give a wedge-shaped feasible region of ( k , k 2 ). We plot the sensitivities for k 3 = 20, 100, and 400 to display a range of possible parameter combinations.
In the δ ≈ δ c case, calculating the partial derivatives is straightforward. To compute the sensitivities, we evaluate the partial derivatives at the epidemic threshold: δ = δ c . Table 1 shows the expressions for these sensitivities, and Fig. 3 shows corresponding plots. Clearly ∂w * ∂ k ≤ 0 and ∂w * ∂ k 2 ≥ 0, with more extreme values near the upper-right corner of the feasible region.
For the δ >> δ c case, the partial derivatives (Table 2) all depend on a factor of 1/δ, so the choice of δ for computing sensitivities does not affect the relative magnitudes of the partial derivatives. For convenience, we select δ = 1.5. The sensitivity plots in (a) (b) Fig. 3 Sensitivities a ∂w * ∂ k and b ∂w * ∂ k 2 for the δ ≈ δ c approximation. White denotes regions of the ( k , k 2 ) plane outside of the feasible region. Sensitivities are evaluated at δ = δ c Fig. 4 show that ∂w * ∂ k ≥ 0, ∂w * ∂ k 2 ≤ 0, and ∂w * ∂ k 3 ≥ 0, with the greatest sensitivity near the curve k 2 2 = k 3 k , though the large magnitude appears to be due to the partial derivatives being undefined there. A significant observation from these sensitivities is that ∂w * ∂ k and ∂w * ∂ k 2 change signs depending on the regime considered. If the goal of an intervention is to reduce the size of the endemic equilibrium, near the epidemic threshold, this can be accomplished in principle by increasing k or decreasing k 2 , which will in effect increase δ c as well. This is intuitive, as an effort to push the system below the epidemic threshold would also decrease the endemic equilibrium for a fixed δ. However, in the δ >> δ c regime, the system is far from the epidemic threshold, and reducing the size of the endemic equilibrium can be accomplished by decreasing k or increasing k 2 . This suggests that containment and mitigation strategies that depend on altering the structure of the contact network may require different goals in terms of the moments of the degree distribution.

Conclusion
In this paper, we have analyzed the super compact pairwise model presented in Simon and Kiss (2016). A non-dimensional version of the model was considered, and a bifurcation analysis was performed, demonstrating that the SCPW and CPW models share an epidemic threshold. Moreover, we derived approximate formulas for the endemic equilibrium in two regimes: when the transmission/recovery ratio is near the epidemic threshold, and far away from it. While the asymptotic techniques used here are ad hoc, similar techniques may prove fruitful in other low-dimensional models of infectious disease spread on networks. However, an exact expression for the endemic equilibrium remains elusive.
Before explaining the advantages of our approach, we acknowledge two limitations of our approximation. First, approximations of the endemic equilibrium for diseases between the two regimes are lacking. Second, while the examples of simulated networks show good agreement between the exact and approximate prevalence, we have not quantified the approximation error generally. As such, there may be types of networks for which our approximation of the endemic equilibrium is less accurate or inappropriate.
Our approximation of the endemic equilibrium is very useful in providing a more detailed look into the interactions of the moments of the degree distribution as they relate to the size of an outbreak. This has implications for disease control measures, particularly those that work by altering the contact network structure. Our results suggest that for SIS-type diseases, strategies to contain (near the epidemic threshold) or mitigate (far away from the epidemic threshold) an outbreak may require different goals. In the mitigation scenario where the prevalence is high, measures might be employed that decrease the first moment k of the degree distribution. In effect, this may mean initiatives aimed at reducing the number of contacts of individuals alone. On the other hand, in the containment scenario where the prevalence is low, decreasing the second moment k 2 may be efficient. When couched in degree distribution terms this goal is hard to conceptualize, but using probability generating functions (Newman et al. 2001) one can show that k 2 is the average number of first and second neighbors of nodes in the network. Thus, measures that reduce both the contacts of individuals and their partners are effective in this scenario. This suggests the importance of contact tracing. We note that the sensitivities also suggest that increasing k 2 in the high prevalence case and increasing k in the low prevalence case may lead to a reduction in the size of the endemic equilibrium, though it is not clear why from a biological perspective.
Our results complement the findings of Eames and Keeling (2002), who observed that the effectiveness of two common control measures, screening and contact tracing, depend on the prevalence at the endemic equilibrium. Screening, which targets and treats individuals, is efficient when the prevalence is high. Contact tracing, which targets and treats individuals and their partners, if efficient when the prevalence is low. Unlike this paper, Eames and Keeling implement these measures through epidemiological parameters (rather than through changing network structure). In this way, our results can be viewed as a network-structure analog for their conclusions and confirm that control measures appropriate in a network setting can be found. Further work in this area may include investigating this phenomenon with alternative models of SIS diseases on networks.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Appendix A Bifurcation and Endemic Steady State
We begin with Theorem 4.1 from Castillo-Chavez and Song (2004), referring to the specific conditions that will be relevant for this analysis. Consider a system of ODEs with a parameter φ : Assume that 0 is an equilibrium for all values of φ. Assume further that D is the linearization matrix of (A.1) around the equilibrium 0 and with φ = 0, and zero is a simple eigenvalue of this matrix with all other eigenvalues having negative real parts. Assume as well that this matrix has a nonnegative right eigenvector w and left eigenvector v corresponding to the zero eigenvalue. Let F k be the kth component of f and If a < 0 and b > 0, then when φ changes from negative to positive, 0 changes its stability from stable to unstable. Correspondingly, a negative unstable equilibrium becomes positive and locally asymptotically stable. We apply this theorem to (17)- (19), where the equilibrium occurs at w = x = z = 0. Moreover, we define a bifurcation parameter φ = δ − δ c , so φ = 0 corresponds to δ = δ c , and ∂ ∂φ = ∂ ∂δ . For consistency with previously established notation, we will treat δ as our parameter, with φ increasing through 0 as δ increases through δ c . The Jacobian given in (21) when w = 0, x = 0, z = 0, and δ = δ c is .4) and the characteristic polynomial is given by The left and right eigenvectors (v and w, respectively) corresponding to the eigenvalue To compute a and b, it is convenient to express (A.2) and (A.3) in matrix-vector form: where H 2 and H 3 are the Hessians of F 2 and F 3 , respectively, at 0. These Hessians are As k 3 > k , it follows that a < 0. The computation for b is simpler. We note that Finally, as a < 0 and b > 0, we conclude that as δ increases through δ c , a positive, asymptotically stable equilibrium emerges, which is the endemic equilibrium.

Appendix B Asymptotic Approximations of the Endemic Equilibrium
The full derivations of the approximations (38) and (49) are presented in this appendix.