Time-dependent solution of the NIMFA equations around the epidemic threshold

The majority of epidemic models are described by non-linear differential equations which do not have a closed-form solution. Due to the absence of a closed-form solution, the understanding of the precise dynamics of a virus is rather limited. We solve the differential equations of the N-intertwined mean-field approximation of the susceptible-infected-susceptible epidemic process with heterogeneous spreading parameters around the epidemic threshold for an arbitrary contact network, provided that the initial viral state vector is small or parallel to the steady-state vector. Numerical simulations demonstrate that the solution around the epidemic threshold is accurate, also above the epidemic threshold and for general initial viral states that are below the steady-state.


Introduction
Epidemiology originates from the study of infectious diseases such as gonorrhoea, cholera and the flu (Bailey 1975;Anderson and May 1992). Human beings do not only transmit infectious diseases from one individual to another, but also opinions, online social media content and innovations. Furthermore, man-made structures exhibit epidemic phenomena, such as the propagation of failures in power networks or the B Bastian Prasse b.prasse@tudelft.nl Piet Van Mieghem p.f.a.vanmieghem@tudelft.nl spread of a malicious computer virus. Modern epidemiology has evolved into the study of general spreading processes (Pastor-Satorras et al. 2015;Nowzari et al. 2016). Two properties are essential to a broad class of epidemic models. First, individuals are either infected with the disease (respectively, possess the information, opinion, etc.) or healthy. Second, individuals can infect one another only if they are in contact (e.g., by a friendship). In this work, we consider an epidemic model which describes the spread of a virus between groups of individuals.
We consider a contact network of N nodes, and every node i = 1, . . . , N corresponds to a group 1 of individuals. If the members of two groups i, j are in contact, then group i and group j can infect one another with the virus. We denote the symmetric N × N adjacency matrix by A and its elements by a i j . If there is a link between node i and node j, then a i j = 1, and a i j = 0 otherwise. Hence, the virus directly spreads between two nodes i and j only if a i j = 1. We stress that in most applications it holds that a ii = 0, since infected individuals in group i usually do infect susceptible individuals in the same group i. At any time t ≥ 0, we denote the viral state of node i by v i (t). The viral state v i (t) is in the interval [0, 1] and is interpreted as the fraction of infected individuals of group i. N -intertwined mean-field approximation (NIMFA) with heterogeneous spreading parameters (Lajmanovich and Yorke 1976;Van Mieghem and Omic 2014) assumes that the curing rates δ i and infection rates β i j depend on the nodes i and j.
Definition 1 (Heterogeneous NIMFA) At any time t ≥ 0, the NIMFA governing equation is for every group i = 1, . . . , N , where δ i > 0 is the curing rate of node i, andβ i j > 0 is the infection rate from node j to node i.
For a vector x ∈ R N , we denote the diagonal matrix with x on its diagonal by diag(x). We denote the N × N curing rate matrix S = diag(δ 1 , . . . , δ N ). Then, the matrix form of (1) is a vector differential equation where v(t) = (v 1 (t), . . . , v N (t)) T is the viral state vector at time t, the N × N infection rate matrix B is composed of the elements β i j =β i j a i j , and u is the N × 1 all-one vector. In this work, we assume that the matrix B is symmetric.
Definition 2 (Steady-State Vector) The N × 1 steady-state vector v ∞ is the non-zero equilibrium of NIMFA, which satisfies In its simplest form, NIMFA (Van Mieghem et al. 2009) assumes the same infection rate β and curing rate δ for all nodes. More precisely, for homogeneous NIMFA the governing equations (2) reduce to dv(t) dt = −δv(t) + βdiag (u − v(t)) Av(t).
For the vast majority of epidemiological, demographical, and ecological models, the basic reproduction number R 0 is an essential quantity (Hethcote 2000;Heesterbeek 2002). The basic reproduction number R 0 is defined (Diekmann et al. 1990) as "The expected number of secondary cases produced, in a completely susceptible population, by a typical infective individual during its entire period of infectiousness". Originally, the basic reproduction number R 0 was introduced for epidemiological models with only N = 1 group of individuals. Van den Driessche and Watmough (2002) proposed a definition of the basic reproduction number R 0 to epidemic models with N > 1 groups. For NIMFA (1), the basic reproduction number R 0 follows (Van den Driessche and Watmough 2002) as R 0 = ρ(S −1 B), where ρ(M) denotes the spectral radius of a square matrix M. For the stochastic Susceptible-Infected-Removed (SIR) epidemic process on data-driven contact networks, Liu et al. (2018) argue that the basic reproduction number R 0 is inadequate to characterise the behaviour of the viral dynamics, since the number of secondary cases produced by an infectious individual varies greatly with time t. In contrast to the stochastic SIR process, for the deterministic NIMFA equations (1), the basic reproduction number R 0 = ρ(S −1 B) is of crucial importance for the viral state dynamics. Lajmanovich and Yorke (1976) showed that there is a phase transition at the epidemic threshold criterion R 0 = 1: If R 0 ≤ 1, then the only equilibrium of NIMFA (1) is the origin, which is globally asymptotically stable. Else, if R 0 > 1, then there is a second equilibrium, the steady-state v ∞ , whose components are positive, and the steady-state v ∞ is globally asymptotically stable for every initial viral state v(0) = 0. For real-world epidemics, the regime around epidemic threshold criterion R 0 = 1 is of particular interest. In practice, the basic reproduction number R 0 cannot be arbitrarily great, since natural immunities and vaccinations lead to significant curing rates δ i and the frequency and intensity of human contacts constrain the infection rates β i j . Beyond the spread of infectious diseases, many real-world systems seem to operate in the critical regime around a phase transition (Kitzbichler et al. 2009;Nykter et al. 2008).
The basic reproduction number R 0 only provides a coarse description of the dynamics of NIMFA (1). Recently (Prasse and Van Mieghem 2019), we analysed the viral state dynamics for the discrete-time version of NIMFA (1), provided that the initial viral state v(0) is small (see also Assumption 2 in Sect. 3). Three results of Prasse and Van Mieghem (2019) are worth mentioning, since we believe that they could also apply to NIMFA (1) in continuous time. First, the steady-state v ∞ is exponentially stable. Second, the viral state is (almost always) monotonically increasing. Third, the viral state v(t) is bounded by linear time-invariant systems at any time t. In this work, we go a step further in analysing the dynamics of the viral state v(t), and we focus on the region around the threshold R 0 = 1. More precisely, we find the closed-form expression of the viral state v i (t) for every node i at every time t when R 0 ↓ 1, given that the initial state v(0) is small or parallel 2 to the steady-state vector v ∞ .
We introduce the assumptions in Sect. 3. Section 4 gives an explicit expression for the steady-state vector v ∞ when R 0 ↓ 1. In Sect. 5, we derive the closed-form expression for the viral state vector v(t) at any time t ≥ 0. The closed-form solution for R 0 ↓ 1 gives an accurate approximation also for R 0 > 1 as demonstrated by numerical evaluations in Sect. 6. Lajmanovich and Yorke (1976) originally proposed the differential equations (1) to model the spread of gonorrhoea and proved the existence and global asymptotic stability of the steady-state v ∞ for strongly connected directed graphs. In Lajmanovich and Yorke (1976), Fall et al. (2007), Wan et al. (2008), Rami et al. (2013), Prasse and Van Mieghem (2018) and Paré et al. (2018), the differential equations (1) are considered as the exact description of the virus spread between groups of individuals. Van Mieghem et al. (2009) derived the differential equations (1) as an approximation of the Markovian Susceptible-Infected-Susceptible (SIS) epidemic process (Pastor-Satorras et al. 2015;Nowzari et al. 2016), which lead to the acronym "NIMFA" for "N -Intertwined Mean-Field Approximation" (Van Mieghem 2011;Van Mieghem and Omic 2014;Devriendt and Van Mieghem 2017). The approximation of the SIS epidemic process by NIMFA is least accurate around the epidemic threshold (Van Mieghem et al. 2009;Van Mieghem and van de Bovenkamp 2015). Thus, the solution of NIMFA when R 0 ↓ 1, which is derived in this work, might be inaccurate for the description of the probabilistic SIS process. Fall et al. (2007) analysed the generalisation of the differential equations (1) of Lajmanovich and Yorke (1976) to a non-diagonal curing rate matrix S. Khanafer et al. (2016) showed that the steady-state v ∞ is globally asymptotically stable, also for weakly connected directed graphs. Furthermore, NIMFA (1) has been generalised to time-varying parameters. Paré et al. (2017) consider that the infection rates 3 β i j (t) depend continuously on time t. Rami et al. (2013) consider a switched model in which both the infection rates β i j (t) and the curing rates δ i (t) change with time t. NIMFA (1) in discrete time has been analysed in Ahn and Hassibi (2013), Paré et al. (2018), Prasse and Van Mieghem (2019) and Liu et al. (2020).

Related work
In Van Mieghem (2014b), NIMFA (4) was solved for a special case: If the adjacency matrix A corresponds to a regular graph and the initial state v i (0) is the same 4 for every node i, then NIMFA with time-varying, homogeneous spreading parameters β(t), δ(t) has a closed-form solution. In this work, we focus on time-invariant but heterogeneous spreading parameters δ i , β i j . We solve NIMFA (1) for arbitrary graphs around the threshold criterion R 0 = 1 and for an initial viral state v(0) that is small or parallel to the steady-state vector v ∞ .

Notations and assumptions
The basic reproduction number R 0 = ρ(S −1 B) is determined by the infection rate matrix B and the curing rate matrix S. Thus, the notation R 0 ↓ 1 is imprecise, since there are infinitely many matrices B, S such that the basic reproduction number R 0 equals 1. To be more precise, we consider a sequence B (n) , S (n) n∈N of infection rate matrices B (n) and curing rate matrices S (n) that converges 5 to a limit (B * , S * ), such that ρ (S * ) −1 B * = 1 and For the ease of exposition, we drop the index n and replace B (n) and S (n) by the notation B and S, respectively. In particular, we emphasise that the assumptions below apply to every element B (n) , S (n) of the sequence. In Sects. 4 to 6, we formally abbreviated the limit process B (n) , S (n) → (B * , S * ) by the notation R 0 ↓ 1. For the proofs in the appendices, we use the lengthier but clearer notation (B, S) → (B * , S * ). Furthermore, we use the superscript notation Ξ * to denote the limit of any variable Ξ that depends on the infection rate matrix B and the curing rate matrix S. For instance, δ * i denotes the limit of the curing rate δ i of node i when (B, S) → (B * , S * ). The Landau-notation In the remainder of this work, we rely on three assumptions, which we state for clarity in this section.
Assumption 1 For every basic reproduction number R 0 > 1, the curing rates are positive and the infection rates are non-negative, i.e., δ i > 0 and β i j ≥ 0 for all nodes i, j. Furthermore, in the limit R 0 ↓ 1, it holds that δ i 0 and δ i ∞ for all nodes i.
We consider Assumption 1 a rather technical assumption, since only non-negative rates δ i and β i j have a physical meaning. Furthermore, if the curing rates δ i were zero, then the differential equations (1) would describe a Susceptible-Infected (SI) epidemic process. In this work, we focus on the SIS epidemic process, for which it holds that δ i > 0.
Assumption 2 For every basic reproduction number R 0 > 1, it holds that v i (0) ≥ 0 and v i (0) ≤ v ∞,i for every node i = 1, . . . , N . Furthermore, it holds that v i (0) > 0 for at least one node i.
For the description of most real-world epidemics, Assumption 2 is reasonable for two reasons. First, the total number of infected individuals often is small in the beginning of an epidemic outbreak. (Sometimes, there is even a single patient zero.) Second, a group i often contains many individuals. For instance, the viral state v i (t) could describe the prevalence of virus in municipality i. Thus, even if there is a considerable total number of infected individuals in group i, the initial fraction v i (0) would be small.

Assumption 3
For every basic reproduction number R 0 > 1, the infection rate matrix B is symmetric and irreducible. Furthermore, in the limit R 0 ↓ 1, the infection rate matrix B converges to a symmetric and irreducible matrix.
Assumption 3 holds if and only if the infection rate matrix B (and its limit) corresponds to a connected undirected graph (Van Mieghem 2014a).

The steady-state around the epidemic threshold
We define the N × N effective infection rate matrix W as In this section, we state an essential property that we apply to solve the NIMFA equations (1) when the basic reproduction number R 0 is close to 1: The steady-state vector v ∞ converges to a scaled version of the principal eigenvector x 1 of the effective infection rate matrix W when R 0 ↓ 1.
Under Assumptions 1 and 3, the effective infection rate matrix W is non-negative and irreducible. Hence, the Perron-Frobenius Theorem (Van Mieghem 2014a) implies that the matrix W has a unique eigenvalue λ 1 which equals the spectral radius ρ(W ). As we show in the beginning of Appendix B, the eigenvalues of the effective infection rate matrix W are real and satisfy λ 1 = ρ(W ) > λ 2 ≥ · · · ≥ λ N . In particular, under Assumptions 1 and 3, the largest eigenvalue λ 1 , the spectral radius ρ(W ) and the basic reproduction number R 0 are the same quantity, i.e., R 0 = ρ(W ) = λ 1 .
In Van Mieghem (2012, Lemma 4) it was shown that, for homogeneous NIMFA (4), the steady-state vector v ∞ converges to a scaled version of the principal eigenvector of the adjacency matrix A when R 0 ↓ 1. We generalise the results of Van Mieghem (2012) to heterogeneous NIMFA (1):

Theorem 1 Under Assumptions 1 and 3, the steady-state vector
where the scalar γ equals and the N × 1 vector η satisfies η 2 ≤ O (R 0 − 1) 2 when the basic reproduction number R 0 approaches 1 from above.
Proof Appendix B.

The viral state dynamics around the epidemic threshold
In Sect. 5.1, we give an intuitive motivation of our solution approach for the NIMFA equations (1) when R 0 ↓ 1. In Sect. 5.2, we state our main result.

Motivation of the solution approach
For simplicity, this subsection is confined to the homogeneous NIMFA equations (4). In numerical simulations (Prasse and Van Mieghem 2018), we observed that the N × N viral state matrix V = (v(t 1 ), . . . , v(t N )), for arbitrary observation times t 1 < · · · < t N , is severely ill-conditioned. Thus, the viral state v(t) at any time t ≥ 0 approximately equals the linear combination of m << N orthogonal vectors y 1 , . . . , y m , and we can write v(t) ≈ c 1 (t)y 1 + · · · + c m (t)y m , see also Prasse and Van Mieghem (2020). Here, the functions c 1 (t), . . . , c m (t) are scalar. We consider the most extreme case by representing the viral state v(t) by a scaled version of only m = 1 vector y 1 , which corresponds to v(t) ≈ c(t)y 1 for a scalar function c(t). The viral state v(t) converges to the steady-state vector v ∞ as t → ∞. Hence, a natural choice for the vector y 1 is y 1 = v ∞ , which implies that c(t) → 1 as t → ∞. If R 0 ≈ 1 and v(0) ≈ 0, then the approximation v(t) ≈ c(t)v ∞ is accurate at all times t ≥ 0 due to two intuitive reasons.
1. If v(t) ≈ 0 when t ≈ 0, then NIMFA (4) is approximated by the linearisation around zero. Hence, it holds that when t ≈ 0. The state v(t) of the linear system (8) converges rapidly to a scaled version of the principal eigenvector x 1 of the matrix (β A − δ I ). Furthermore, Theorem 1 states that v ∞ ≈ γ x 1 when R 0 ≈ 1. Thus, the viral state v(t) rapidly converges to a scaled version of the steady-state v ∞ : 2. Suppose that the viral state v(t) approximately equals to a scaled version of the steady-state vector v ∞ . (In other words, the viral state v(t) is "almost parallel" to the vector v ∞ .) Then, it holds that for some scalar c(t). We insert (9) into the NIMFA equations (4), which yields that For homogeneous NIMFA (4), the steady-state equation (3) becomes We substitute (11) in (10) and obtain that Since v ∞ ≈ γ x 1 around the epidemic threshold, it holds that Av ∞ ≈ ρ(A)v ∞ . Hence, we obtain that Left-multiplying (13) by The logistic differential equation (14) has been introduced by Verhulst (1838) as a population growth model and has a closed-form solution.
Due to the two intuitive steps above, NIMFA (4) reduces around the threshold R 0 ≈ 1 to the one-dimension differential equation (14). Solving (14) for the function c(t) gives an approximation of the viral state v(t) by (9). The solution approach is applicable to other dynamics on networks, see for instance (Devriendt and Lambiotte 2020). However, the reasoning above is not rigorous for two reasons. First, the viral state vector v(t) is not exactly parallel to the steady state v ∞ . To be more specific, instead of (9) it holds that for some N × 1 error vector ξ(t) which is orthogonal to the steady-state vector v ∞ . In Sect. 5.2, we use (15) as an ansatz for solving NIMFA (1). Second, the steady-state vector v ∞ is not exactly parallel to the principal eigenvector x 1 . More precisely, we must consider the vector η in (6). Since η = 0, the step from (12) to (13) is affected by an error.

The solution around the epidemic threshold
Based on the motivation in Sect. 5.1, we aim to solve the NIMFA differential equations (1) around the epidemic threshold criterion R 0 = 1. The ansatz (15) forms the basis for our solution approach. From the orthogonality of the error vector ξ(t) and the steady-state vector v ∞ , it follows that the function c(t) at time t equals The error vector ξ(t) at time t follows from (15) and (16) as Our solution approach is based on two steps. First, we show that 6 the error term ξ(t) satisfies ξ(t) = O((R 0 − 1) 2 ) at every time t when R 0 ↓ 1. Hence, the error term ξ(t) converges to zero uniformly in time t. Second, we find the solution of the scalar function c(t) at the limit R 0 ↓ 1. Assumption 2 implies that 7 the viral state v(t) does not overshoot the steadystate v ∞ : Theorem 2 states that the error term ξ(t) converges to zero in the order of (R 0 − 1) 2 when R 0 ↓ 1.
Theorem 2 Under Assumptions 1 to 3, there exist constants σ 1 , σ 2 > 0 such that the error term ξ(t) at any time t ≥ 0 is bounded by when the basic reproduction number R 0 approaches 1 from above.

Proof Appendix D.
Under Assumption 2, the steady-state v ∞ is exponentially stable for NIMFA in discrete time (Prasse and Van Mieghem 2019). If the steady-state v ∞ is exponentially stable, then the error vector ξ(t) goes to zero exponentially fast, since ξ(t) is orthogonal to v ∞ . Thus, the first addend on the right-hand side in (18) is rather expectable, under the conjecture that the steady-state v ∞ is exponentially stable also for continuous-time NIMFA (1). Regarding this work, the most important implication of Theorem 2 is that ξ(t) = O (R 0 − 1) 2 uniformly in time t when R 0 ↓ 1, provided the initial value ξ(0) of the error vector is negligibly small. We define the constant Υ (0), which depends on the initial viral state v(0), as 6 Theorem 1 implies that the steady- Thus, a linear convergence of the error term ξ(t) to zero, i.e., ξ(t) 2 = O (R 0 − 1), would not be sufficient to show that the viral state v(t) converges to c(t)v ∞ when R 0 ↓ 1. 7 In Prasse and Van Mieghem (2019), an analogous statement has been proved for the discrete-time version of the NIMFA equations (2). Furthermore, we define the viral slope w, which determines the speed of convergence to the steady-state v ∞ , as Then, building on Theorems 1 and 2, we obtain our main result: Theorem 3 Suppose that Assumptions 1 to 3 hold and that, for some constant p > 1, Then, there exists some constant σ > 0 such that where s = min{ p, 2}, when the basic reproduction number R 0 approaches 1 from above.
Proof Appendix E.
We emphasise that Theorem 3 holds for any connected graph corresponding to the infection rate matrix B. Theorem 3 is in agreement with the universality of the SIS prevalence (Van Mieghem 2016). The bound (21) states a convergence of the viral state v(t) to the approximation v apx (t) which is uniform in time t. Furthermore, since both the viral state v(t) and the approximation v apx (t) converge to the steady-state v ∞ , it holds that v(t) − v apx (t) 2 → 0 when t → ∞. At time t = 0, we obtain from Theorem 3 and (17) that Hence, for general t ≥ 0 the approximation error v(t) − v apx (t) 2 / v ∞ 2 does not converge to zero faster than O (R 0 − 1) p−1 , and the bound (21) is best possible (up to the constant σ ) when p ≤ 2. With (17), the term ξ(0) 2 in Theorem 2 can be expressed explicitly with respect to the initial viral state v(0) and the steady-state v ∞ . In particular, it holds that ξ(0) 2 ≤ v(0) 2 . Furthermore, if the initial viral state v(0) is parallel to the steady-state vector v ∞ , then it holds that ξ(0) = 0. Thus, if the initial viral state v(0) is small or parallel to the steady-state vector v ∞ , then it holds that ξ(0) = 0 and the bound (21) on the approximation error vector becomes The time-dependent solution to NIMFA (1) at the epidemic threshold criterion R 0 = 1 depends solely on the viral slope w, the steady-state vector v ∞ and the initial viral state v(0). The viral slope w converges to zero as R 0 ↓ 1. Thus, Theorem 3 implies that the convergence time to the steady-state v ∞ goes to infinity when R 0 ↓ 1, even though the steady-state v ∞ converges to zero. More precisely, it holds: Corollary 1 Suppose that Assumptions 1 and 3 hold and that the initial viral state v(0) equals v(0) = r 0 v ∞ for some scalar r 0 ∈ (0, 1). Then, for any scalar r 1 ∈ [r 0 , 1), the largest time t 01 at which the viral state satisfies v i (t 01 ) ≤ r 1 v ∞,i for every node i converges to when the basic reproduction number R 0 approaches 1 from above.
Proof Appendix F.
We combine Theorem 1 and Theorem 3 to obtain Corollary 2.
Corollary 2 Suppose that Assumptions 1 to 3 hold and that, for some constant p > 1, Then, there exists some constant σ > 0 such that where s = min{ p, 2}, when the basic reproduction number R 0 approaches 1 from above.
In contrast to Theorem 3, the approximation error v(t) −ṽ apx (t) 2 in Corollary 2 does not converge to zero when t → ∞, since we replaced the steady-state v ∞ by the first-order approximation of Theorem 1. Corollary 2 implies that at every time t when R 0 ↓ 1, provided that the initial viral state v(0) is small or parallel to the steady-state vector v ∞ . From (24) it follows that, around the epidemic threshold criterion R 0 = 1, the eigenvector centrality (Van Mieghem 2010) fully determines the "dynamical importance" of node i versus node j. For homogeneous NIMFA (4), the infection rate matrix B and the curing rate matrix S reduce to B = β A and S = δ I , respectively. Hence, the effective infection rate matrix becomes W = β δ A, and the principal eigenvector x 1 of the effective infection rate matrix W equals the principal eigenvector of the adjacency matrix A. Furthermore, the limit process R 0 ↓ 1 reduces to τ ↓ τ c , with the effective infection rate τ = β δ and the epidemic threshold τ c = 1/ρ(A). For homogeneous NIMFA (4), Theorem 3 reduces to:

Corollary 3 Suppose that Assumptions 1 to 3 hold and consider the viral state
where s = min{ p, 2}, when the effective infection rate τ approaches the epidemic threshold τ c from above.
Proof Appendix G.
From Corollary 3, we can obtain the analogue to Corollary 2 for NIMFA (4) with homogeneous spreading parameters β, δ. Furthermore, the approximation v apx (t) defined by (25) equals the exact solution (Van Mieghem 2014b) of homogeneous NIMFA (4) on a regular graph, provided that the initial state v i (0) is the same for every node i. In particular, the net dose (t), a crucial quantity in Van Mieghem (2014b); Kendall (1948), is related to the viral slope w via (t) = wt. Theorem 3 and Corollary 3 suggest that, around the epidemic threshold criterion R 0 = 1, the dynamics of heterogeneous NIMFA (1) closely resembles the dynamics of homogeneous NIMFA (4). In particular, we pose the question: Can heterogeneous NIMFA (1) be reduced to homogeneous NIMFA (4) around the epidemic threshold criterion R 0 = 1 by choosing the homogeneous spreading parameters β, δ and the adjacency matrix A accordingly?
Theorem 4 Consider heterogeneous NIMFA (1) with given spreading parameters β i j , δ i . Suppose that Assumptions 1 to 3 hold and that, for some constant p > 1, when the basic reproduction number R 0 approaches 1 from above. Define the homogeneous NIMFA system where the homogeneous curing rate δ hom equals the homogeneous infection rate β hom equals with the variable γ defined by (7), and the self-infection rates β ii,hom equal where s = min{ p, 2}, when the basic reproduction number R 0 approaches 1 from above.
Proof Appendix H.
In other words, when R 0 ↓ 1, for any contact network and any spreading parameters δ i , β i j , heterogeneous NIMFA (1) can be reduced to homogeneous NIMFA (4) on a complete graph plus self-infection rates β ii,hom . We emphasise that the sole influence of the topology on the viral spread is given by the self-infection rates β ii,hom . Thus, under Assumptions 1 to 3, the network topology has a surprisingly small impact on the viral spread around the epidemic threshold.

Numerical evaluation
We are interested in evaluating the accuracy of the closed-form expression v apx (t), given by (20), when the basic reproduction number R 0 is close, but not equal, to one. We generate an adjacency matrix A according to different random graph models. If a i j = 1, then we set the infection rates β i j to a uniformly distributed random number in [0.4, 0.6] and, if a i j = 0, then we set β i j = 0. We set the initial curing rates δ (0) l to a uniformly distributed random number in [0.4, 0.6]. To set the basic reproduction number R 0 , we set the curing rates δ l to a multiple of the initial curing rates δ (0) l , i.e. δ l = σ δ (0) l for every node l and some scalar σ such that ρ(W ) = R 0 . Thus, we realise the limit process R 0 ↓ 1 by changing the scalar σ . Only in Sect. 6.2, we consider homogeneous spreading parameters by setting β i j = 0.5 and δ (0) i = 0.5 for all nodes i, j. Numerically, we obtain the "exact" NIMFA viral state sequence v(t) by Euler's method for discretisation, i.e., for a small sampling time T and a discrete time slot k ∈ N. In Prasse and Van Mieghem (2019), we derived an upper bound T max on the sampling time T which ensures that the discretisation (29) of NIMFA (1) converges to the steady-state v ∞ . We set the sampling time T to T = T max /100. Except for Sect. 6.3, we set the initial viral state to v(0) = 0.01v ∞ . We define the convergence time t conv as the smallest time t at which holds for every node i. Thus, at the convergence time t conv the viral state v(t conv ) has practically converged to the steady-state v ∞ . We evaluate Theorem 3 with respect to the approximation error V , which we define as All results are averaged over 100 randomly generated networks.

Approximation accuracy around the epidemic threshold
We generate a Barabási-Albert random graph (Barabási and Albert 1999) with N = 500 nodes and the parameters m 0 = 5, m = 2. Figure 1 gives an impression of the accuracy of the approximation of Theorem 3 around the epidemic threshold criterion R 0 = 1. For a basic reproduction number R 0 ≤ 1.1, the difference of the closed-form expression of Theorem 3 to the exact NIMFA viral state trace is negligible. We aim for a better understanding of the accuracy of the closed-form expression of Theorem 3 when the basic reproduction number R 0 converges to one. We generate Barabási-Albert and Erdős-Rényi connected random graphs with N = 100, . . . , 1000 nodes. The link probability of the Erdős-Rényi graphs (Erdős and Rényi 1960) is set to p ER = 0.05. Figure 2 illustrates the convergence of the approximation of Theorem 3 to the exact solution of NIMFA (1). Around the threshold criterion R 0 = 1, the approximation error V converges linearly to zero with respect to the basic reproduction number R 0 , which is in agreement with Theorem 3. The greater the network size N , the greater is the approximation error V for Barabási-Albert networks. The greater the network size N , the lower is the approximation error V for Erdős-Rényi graphs.

Impact of degree heterogeneity on the approximation accuracy
For NIMFA (4) with homogeneous spreading parameters β, δ, the approximation v apx (t) defined by (4) is exact if the contact network is a regular graph. We are interested how the approximation accuracy changes with respect to the heterogeneity of the node degrees. We generate Watts-Strogatz (Watts and Strogatz 1998) Barabási-Albert random graphs. with N = 100 nodes and an average node degree of 4. We vary the link rewiring probability p WS from p WS = 0, which correspond to a regular graph, to p WS = 1, which corresponds to a "completely random" graph. Figure 3 depicts the approximation error V versus the rewiring probability p WS for homogeneous spreading parameters β, δ. Interestingly, the approximation error reaches a maximum and improves when the adjacency matrix A is more random.

Impact of general initial viral states on the approximation accuracy
Theorem 3 required that the initial error ξ(0) converges to zero, which means that the initial viral state v(0) must be parallel to the steady-state v ∞ or, since ξ(0) 2 ≤ v(0) , converge to zero. To investigate whether the approximation of Theorem 3 is accurate also when the initial error ξ(0) does not converge to zero, we set the initial viral state v i (0) of every node i to a uniformly distributed random number in (0, r 0 v ∞,i ] for some scalar r 0 ∈ (0, 1]. By increasing the scalar r 0 , the initial viral state v(0) is "more random". Figure 4 shows that the approximation error V is almost unaffected by an initial viral state v(0) that is neither parallel to the steady-state v ∞ nor small.
The approximation error V versus the scalar r 0 , which controls the variance of the randomly generated initial viral state v(0), for Barabási-Albert networks with N = 250 nodes Figure 5 shows that the viral state v(t) converges rapidly to the approximation v apx (t) as time t increases.
For general initial viral states v(0) with ξ(0) = 0, it holds that v apx (0) = v(0) since the approximation v apx (0) is parallel to the steady-state vector v ∞ . Hence, the approximation v apx (t) does not converge point-wise to the viral state v(t) when R 0 ↓ 1. However, based on the results shown in Figs. 4 and 5, we conjecture convergence with respect to the L 2 -norm for general initial viral states v(0) when R 0 ↓ 1.

Conjecture 1
Suppose that Assumptions 1 to 3 hold. Then, it holds for the approxima- when the basic reproduction number R 0 approaches 1 from above.

Directed infection rate matrix
The proof of Theorem 3 relies on a symmetric infection rate matrix B as stated by Assumption 3. We perform the same numerical evaluation as shown in Fig. 2 in Sect. 6.1 with the only difference that we generate strongly connected directed Erdős-Rényi random graphs. Figure 6 demonstrates the accuracy of the approximation v apx (t) for a directed infection rate matrix B, which leads us to: Conjecture 2 Suppose that Assumptions 1 and 2 hold and that the infection rate matrix B is irreducible but, in contrast to Assumption 3, not necessarily symmetric. Then, the viral state v(t) is "accurately described" by the approximation v apx (t) when the basic reproduction number R 0 approaches 1 from above.

Accuracy of the approximation of the convergence time
Corollary 1 gives the expression of the convergence time t 01 from the initial viral state v(0) = r 0 v ∞ to the viral state v(t 01 ) ≤ r 1 v ∞ for any scalars 0 < r 0 ≤ r 1 < 1 around the epidemic threshold criterion R 0 = 1. We set the scalars to r 0 = 0.01 and r 1 = 0.9   where t 01 denotes the exact convergence time andt 01 denotes the approximate expression of Corollary 1. We generate Barabási-Albert and Erdős-Rényi random graphs with N = 100, . . . , 1000 nodes. Figure 7 shows that Corollary 1 gives an accurate approximation of the convergence time t 01 when the basic reproduction number R 0 is reasonably close to one.

Reduction to a complete graph with homogeneous spreading parameters
Theorem 4 states that, around the epidemic threshold, heterogeneous NIMFA (1) on any graph can be reduced to homogeneous NIMFA (4) on a complete graph. Figures 8  and 9 show the approximation accuracy of Theorem 4 for Erdős-Rényi and Barabási-Albert random graphs, respectively. To accurately approximate heterogeneous NIMFA on Barabási-Albert graphs by homogeneous NIMFA on a complete graph, the basic reproduction number R 0 must be closer to 1 than for Erdős-Rényi graphs.

Conclusion
We solved the NIMFA governing equations (1) with heterogeneous spreading parameters around the epidemic threshold when the initial viral state v(0) is small or parallel to the steady-state v ∞ , provided that the infection rates are symmetric (β i j = β ji ). Numerical simulations demonstrate the accuracy of the solution when the basic reproduction number R 0 is close, but not equal, to one. Furthermore, the solution serves as an accurate approximation also when the initial viral state v(0) is neither small nor parallel to the steady-state v ∞ . We observe four important implications of the solution of NIMFA around the epidemic threshold. First, the viral state v(t) is almost parallel to the steady-state v ∞ for every time t ≥ 0. On the one hand, since the viral dynamics approximately remain in a one-dimensional subspace of R N , an accurate network reconstruction is numerically not viable around the epidemic threshold (Prasse and Van Mieghem 2018). Furthermore, when the basic reproduction number R 0 is large, then the viral state v(t) rapidly converges to the steady-state v ∞ , which, again, prevents an accurate network reconstruction. On the other hand, only the principal eigenvector x 1 of the effective infection rate matrix W and the viral slope w are required to predict the viral state dynamics around the epidemic threshold. Thus, around the epidemic threshold, the prediction of an epidemic does not require an accurate network reconstruction.
Second, the eigenvector centrality (with respect to the principal eigenvector x 1 of the effective infection rate matrix W ) gives a complete description of the dynamical importance of a node i around the epidemic threshold. In particular, the ratio v i (t)/v j (t) of the viral states of two nodes i, j does not change over time t.
Third, around the epidemic threshold, we gave an expression of the convergence time t 01 to approach the steady-state v ∞ . The viral state v(t) converges to the steadystate v ∞ exponentially fast. However, as the basic reproduction number R 0 approaches one, the convergence time t 01 goes to infinity.
Fourth, around the epidemic threshold, NIMFA with heterogeneous spreading parameter on any graph can be reduced to NIMFA with homogeneous spreading parameters on the complete graph plus self-infection rates.
Potential generalisations of the solution of NIMFA to non-symmetric infection rate matrices B or time-dependent spreading parameters β i j (t), δ l (t) stand on the agenda of future research.

Acknowledgements
We are grateful to Karel Devriendt for his help in proving Theorem 4.
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/.

A Nomenclature
The eigenvalues of the effective infection rate matrix W are denoted, in decreasing order, by |λ 1 | ≥ · · · ≥ |λ N |. The principal eigenvector of unit length of the matrix W is denoted by x 1 and satisfies W x 1 = λ 1 x 1 . The greatest and smallest curing rate in {δ 1 , . . . , δ N } are denoted by δ max and δ min , respectively. The numerical radius r (M) for an N × N matrix M is defined as (Horn and Johnson 1990) where z H is the conjugate transpose of a complex N × 1 vector z. For a square matrix M, we denote the 2-norm by M 2 , which equals the largest singular value of M. In particular, it holds that the 2-norm of the curing rate matrix S equals S 2 = δ max . Table 1 summarises the nomenclature.

B Proof of Theorem 1
The steady-state v ∞ solely depends on the effective infection rate matrix W : By leftmultiplication of (3) with the diagonal matrix S −1 , we obtain that In general, the effective infection rate matrix W , defined in (5) as W = S −1 B, is asymmetric, which prevents a straightforward adaptation of the proof in Van Mieghem (2012, Lemma 4). However, the matrix W is similar to the matrix Since the infection rate matrix B is symmetric under Assumption 3, the matrixW is symmetric. Hence, the matrixW , and also the effective infection rate matrix W , are diagonalisable. With (32), we write the steady-state (31) with respect to the symmetric matrixW as We decompose the matrixW as where the eigenvalues ofW are real and equal to λ 1 > λ 2 ≥ · · · ≥ λ N with the corresponding normalized eigenvectors denoted byx 1 , . . . ,x N . Then, the steady-state vector v ∞ can be expressed as linear combination where the coefficients equal ψ l = v T ∞x l . To prove Theorem 1, we would like to express the coefficients ψ 1 , . . . , ψ N as a power series around R 0 = 1. However, in the limit process (B, S) → (B * , S * ), the eigenvectorsx 1 , . . . ,x N of the matrixW are not necessarily constant. Hence, the coefficients ψ l depend on the full matrixW and not only on the basic reproduction number R 0 . To overcome the challenge of nonconstant eigenvectorsx 1 , . . . ,x N in the limit process (B, S) → (B * , S * ), we define the symmetric auxiliary matrix for a scalar z ≥ 1. Thus, the matrix M(z) is obtained from the matrixW by replacing the largest eigenvalue λ 1 ofW by z. In particular, the definition of the matrix M(z) in (35) and (34) illustrate that M(λ 1 ) =W . When the matrixW is formally replaced by the matrix M(z), the steady-state equation (33) becomes where the N × 1 vectorṽ(z) denotes the solution of (36). Since M(R 0 ) =W , the solution of (36) at z = R 0 and the solution to (33) coincide, i.e.,ṽ(R 0 ) = v ∞ . Lemma 2 expresses the solution of the equation (36) as a power series.

Proof
The proof is an adaptation of the proof (Van Mieghem 2012, Lemma 4). We express the solutionṽ(z) of (36) as linear combination of the vectors S − 1 2x 1 , . . . , S − 1 2x N , i.e.,ṽ Since the diagonal matrix S − 1 2 is full rank, the vectors S − 1 2x k , where k = 1, . . . , N , are linearly independent. Furthermore, we express the coefficients ψ k (z) as a power series where g 0 (k) = 0 for every k = 1, . . . , N , since (Lajmanovich and Yorke 1976) it holds thatṽ(z) = 0 when z = 1. We denote the eigenvalues of the matrix M(z) by By substituting (38) into (36), we obtain that and left-multiplying with the eigenvectorx T m , for any m = 1, . . . , N , yields We define Then, we rewrite (41) as First, we focus on the left-hand side of (42), which we denote by With the power series (39), we obtain that Further rewriting yields that Second, we rearrange the right-hand side of (42) as By the definition of λ k (z) in (40) it holds that λ 1 (z) = z, and we obtain that Introducing the power series (39) into (44) and executing the Cauchy product for We shift the index j in the first term and obtain Finally, we equate powers in (z − 1) j in (43) and (45), which yields for j = 1 that for every m = 1, . . . , N . The spectral radius of the limit W * of the effective infection rate matrix W equals 1. Furthermore, the limit W * is a non-negative and irreducible matrix. Thus, the eigenvalues of the limit W * obey λ * 1 = 1 > |λ * m | for every m ≥ 2, which implies that |λ m | < 1 for every m ≥ 2 provided that (B, S) is sufficiently close to (B * , S * ). With the definition of λ m (z) in (40), we obtain from (46) that g 1 (m) = 0 when m ≥ 2 provided that (B, S) is sufficiently close to (B * , S * ), since z ≥ 1.
We believe that, based on (47), a recursion for the coefficients g j (k) can be obtained for powers j ≥ 2, similar to the proof of Van Mieghem (2012, Lemma 4). The radius of convergence of the power series (49) is an open problem, see also He and Van Mieghem (2020). To express the solutionṽ(z) in (37) in terms of the principal eigenvector x 1 of the effective infection rate matrix W , we propose Lemma 3.

Lemma 3 Under Assumptions 1 and 3, it holds that
Proof From (32), it follows that the principal eigenvectorx 1 of the matrixW and the principal eigenvector x 1 of the effective infection rate matrix W are related viã or, component-wise, Then, we rewrite the left-hand side of (50) as Writing out the quadratic form in the numerator completes the proof.
The basic reproduction number R 0 converges to 1 when (B, S) → (B * , S * ). Hence, if (B, S) is sufficiently close to (B * , S * ), then the basic reproduction number R 0 is smaller than the radius of convergence of the power series (38). Thus, if (B, S) is sufficiently close to (B * , S * ), then the solutionṽ(R 0 ) to (36) at z = R 0 follows with Lemma 2 asṽ where the last equality follows from Lemma 3 and the definition of the scalar γ in (7). We emphasise that Lemma 2 implies that γ = O(R 0 − 1) and, hence, Since M(R 0 ) =W , the solution of (36) at z = R 0 and the solution to (33) coincide, i.e.,ṽ(R 0 ) = v ∞ . Thus, from the definition of the vector η in (6), we obtain that when (B, S) → (B * , S * ). Lemma 2 states that φ(z) 2 = O (z − 1) 2 as z ↓ 1. Hence, we obtain from (51) that

C Proof of Lemma 1
We divide Lemma 1 into two parts. In Sect. C.1, we prove that the viral state v(t) does not overshoot the steady-state v ∞ . In Sect. C.2, we show that the function c(t) lies in the interval [0, 1].

C.1 Absence of overshoot
The proof follows the same reasoning as Prasse and Van Mieghem (2019, Corollary 1). Assume that at some time Since v j (t 0 ) ≤ v ∞, j for every node j, we obtain that where the last equality follows from the steady-state equation (3).
dt t=t 0 ≤ 0, which means that, at time t 0 , the viral state v i (t 0 ) does not increase. Hence, the viral state v i (t 0 ) cannot exceed the steady-state v ∞,i at any time t ≥ 0.

C.2 Boundedness of the function c(t)
Relation (16) indicates that Section C.1 shows that Assumption 2 implies that v i (t) ≤ v ∞,i for all nodes i and every time t. Thus, we obtain from (53) that Analogously, since v i (t) ≥ 0 for all nodes i and every time t, we obtain from (53) that c(t) ≥ 0.

D Proof of Theorem 2
By inserting the ansatz (15) into the NIMFA equations (2), we obtain that Here, the function Λ 1 (t) is given by which simplifies, with the steady-state equation (3), to The function Λ 2 (t) is given by To show that the error term ξ(t) converges to zero at every time t when (B, S) → (B * , S * ), we consider the squared Euclidean norm ξ(t) 2 2 . The convergence of the squared norm ξ(t) 2 2 to zero implies the convergence of the error term ξ(t) to zero. The derivative of the squared norm ξ(t) 2 2 is given by Thus, we obtain from (54) that since ξ T (t)v ∞ = 0 by definition of ξ(t). We do not know how to solve (57) exactly, and we resort to bounding the two addends on the right-hand side of (57) in Sects. D.1 and D.2, respectively. In Sect. D.3 we complete the proof of Theorem 2 by deriving an upper bound on the squared norm ξ(t) 2 2 .

D.1 Upper bound on T (t)3 1 (t)
We obtain an upper bound on the projection of the function Λ 1 (t) onto the error vector ξ(t), which is linear with respect to the norm ξ(t) 2 : Lemma 4 Under Assumptions 1 to 3, it holds at every time t ≥ 0 that Proof From (55) and the definition of the matrix W in (5) it follows that With Theorem 1, we obtain The triangle inequality yields that With the Cauchy-Schwarz inequality, the first addend in (58) is upper-bounded by since x 1 2 = 1 and the matrix S is symmetric. The matrix 2-norm is submultiplicative, which yields that Thus, (58) gives that since γ > 0 and R 0 > 1. We consider the second addend in (59), which we write with (32) as From the Cauchy-Schwarz inequality and the sub-multiplicativity of the matrix norm we obtain The triangle inequality and the symmetry of the matrixW imply that Thus, we can upper bound the second added in (59) by Finally, Lemma 1 states that 0 ≤ c(t) ≤ 1, which implies that and completes the proof.

D.2 Upper bound on T (t)3 2 (t)
Lemma 5 states an intermediate result, which we will use to bound the projection of the function Λ 2 (t) onto the error vector ξ(t).
Lemma 5 Suppose that Assumptions 1 to 3 hold. Then, at every time t ≥ 0, it holds that Proof From (56) it follows that To simplify (60), we aim to bound the last addend of (60) by an expression that is quadratic in the error vector ξ(t). The last addend equals Since v(t) = c(t)v ∞ + ξ(t) and v i (t) ≥ 0 for every node i at every time t, it holds that By inserting (62) in (61), the last addend of (60) is upper bounded by By applying the upper bound (63) to (60), we obtain that With the definition of the matrixW in (32), we obtain and further rearranging completes the proof.
For any scalar ς ∈ [0, 1] and any vector υ ∈ R N , we define Then, we obtain from Lemma 5 that To upper-bound the term Θ(c(t), ξ(t), B, S), we make use of (parts of) the results of Issos (1966), which are analogues of the Perron-Frobenius Theorem for the numerical radius of a non-negative, irreducible matrix: Theorem 5 (Issos 1966) Let M be a real irreducible and non-negative N × N matrix. Then, there is a positive vector z ∈ R N of length z T z = 1 such that z T Mz = r (M). Furthermore, ifz T Mz = r (M) holds for a vectorz ∈ R N of lengthz Tz = 1, then eitherz = z orz = −z.
We refer the reader to Issos (1966), Maroulas et al. (2002) and Li et al. (2002) for further results on the numerical radius of non-negative matrices. We apply Theorem 5 to obtain: Lemma 6 Denote the set of N × 1 vectors with at least one positive and at least one negative component as Then, it holds Θ(ς, υ, B, S) < R 0 for every scalar ς ∈ [0, 1] and for every vector υ ∈ S.
Proof By introducing the N × 1 vector υ = S 1 2 υ and by using (32), we rewrite the term Θ(ς, υ, B, S) as For every scalar ς ∈ [0, 1] the matrix (diag(u − ςv ∞ )W ) is irreducible and nonnegative. Since υ ∈ S and the matrix S is a diagonal matrix with non-negative entries, it holds thatυ i < 0 andυ j > 0 for some i, j. Hence, at least two components of the vectorυ have different signs, and Theorem 5 implies that (65) is upper-bounded by Since the matrixW is irreducible and diag(u − ςv ∞ )W ≤W for every ς ∈ [0, 1], where the inequality holds element-wise, it holds (Li et al. 2002, Corollary 3.6.) that The matrixW is symmetric, and, hence, the numerical radius r W equals the spectral radius ρ W = R 0 , which yields that Finally, we obtain a bound on the projection of the function Λ 2 (t) onto the error vector ξ(t): Lemma 7 Under Assumptions 1 to 3, there is some constant ω > 0 such that As a first step, we consider the value of Θ max (B * , S * ) at the limit (B * , S * ). Since the steady-state v ∞ equals to zero at the limit (B * , S * ), we obtain from (65) that where we denoteW * = (S * ) − 1 2 B * (S * ) − 1 2 . Since it holds R 0 = 1 at the limit (B * , S * ), Lemma 6 implies that As a second step, we consider that the infection rate matrix B and the curing rate matrix S do not equal the respective limit B * and S * . Thus, there are non-zero N × N matrices ΔB, ΔS and ΔW such that B = B * + ΔB, S = S * + ΔS, and W =W * + ΔW . Then, we obtain from (65) that which is upper-bounded by Maximising every addend in (69) where the last equality follows from the definition of Θ max (B * , S * ) in (66). Regarding the second addend in (70) where the last equality follows from the definition the numerical radius. Hence, the second addend in (70) for some ς (2) opt ∈ [0, 1]. With (71), (72) and (73), we obtain from (70) that The numerical radius r (M) is a vector 8 norm (Horn and Johnson 1990)  Hence, for every scalar ω > 0 there is a ϑ(ω) such that B − B * 2 < ϑ(ω) and S − S * 2 < ϑ(ω) implies that We choose the scalar ω = (1 − Θ max (B * , S * ))/2, which is positive due to (68). Then, the right-hand side of (75) becomes Thus, we obtain from (75) that holds for all (B, S) which are sufficiently close to the limit (B * , S * ). By definition, the error vector ξ(t) at any time t ≥ 0 is orthogonal to the steady-state vector v ∞ . Since the steady-state v ∞ is positive, the error vector ξ(t) has at least one positive and one negative element, and, hence, it holds that ξ(t) ∈ S. Thus, we obtain from the definition of the term Θ max (B, S) in (66) that With (76), we obtain from (64) that From the sub-multiplicativity of the matrix norm, we obtain which completes the proof, since S 1 2 2 2 = δ max .

D.3 Bound on the error vector (t)
With Lemma 4 and Lemma 7, we upper-bound (57) by We denote and we obtain that The upper bound (78) is a linear first-order ordinary differential inequality, which is solved by (Arfken and Weber 1999) which simplifies to The triangle inequality yields that Furthermore, since e −ωδ max t ≤ 1 at every time t ≥ 0, we obtain from (79) that The maximum δ max of the curing rates converges to some limit δ * max when (B, S) → (B * , S * ). Hence, for any > 0 it holds that δ * max − < δ max when (B, S) approaches (B * , S * ). For some ∈ (0, δ * max ), we set the constant Then, it holds that σ 1 < ωδ max when (B, S) approaches (B * , S * ), and we obtain from (80) that

E Proof of Theorem 3
By projecting the differential equation (54) onto the steady-state vector v ∞ , we obtain that since v T ∞ ξ(t) = 0 by definition of the error term ξ(t). We divide by v ∞ 2 2 and obtain with (55) that The first addend in the differential equation (82) can be expressed in a simpler manner when (B, S) approaches (B * , S * ): where ζ = O (R 0 − 1) 2 when (B, S) approaches (B * , S * ).
Proof With Theorem 1 and the definition of the matrix W in (5), the numerator of the left-hand side of (83) becomes where the last equality follows from W x 1 = R 0 x 1 . Thus, it holds that Under Assumption 3, both matrices B and S are symmetric, which implies that Hence, we obtain from (84) that Since γ = O(R 0 − 1) and η 2 = O (R 0 − 1) 2 , we finally rewrite the numerator of the left-hand side of (83) as With Theorem 1, the denominator of the left-hand side of (83) equals Combining the approximate expressions for the numerator (85) and the denominator (86) completes the proof.
We define the viral slope w as and the function n(t) as Then, we obtain from (82) that The function n(t) is complicated and depends on the error vector ξ(t). Hence, we cannot solve the differential equation (89) for the function c(t) without knowing the solution for the error vector ξ(t). However, as (B, S) → (B * , S * ), the function n(t) converges to zero uniformly in time t as stated by the bound in Lemma 9.
Proof Regarding the first addend in the definition of the function n(t) in (88), Lemma 1 implies that 0 ≤ c(t) − c 2 (t) ≤ 1/4 at every time t. Hence, Lemma 8 yields that there is a constantσ 0 such that at every time t when (B, S) approaches (B * , S * ). Regarding the second addend of the function n(t) defined in (88), it follows from the definition of the function Λ 2 (t) in With the definition of the matrixW in (32), we obtain that

The Cauchy-Schwarz inequality yields an upper bound as
In the following, we consider the three addends in (90) separately. Regarding the first addend, we obtain with the definition of the matrixW in (32) that where the last equality follows from Theorem 1. Thus, the triangle inequality yields With the sub-multiplicativity of the matrix 2-norm, we obtain there is a constantσ 1 such that when (B, S) approaches (B * , S * ). Regarding the second addend in (90), it holds that Since when (B, S) approaches (B * , S * ). Regarding the third addend in (90), it holds per definition of the matrix 2-norm that where the last inequality follows from c(t) ≤ 1, as stated by Lemma 1, and the definition of the effective infection rate matrix W in (5). Hence, we obtain the upperbound for some constantσ 3 when (B, S) approaches (B * , S * ). We apply the three upper bounds (91), (92) and (93) to (90) and obtain that at every time t. Thus, we have obtained an upper bound, which is proportional to the norm of the error vector ξ(t). Finally, we apply Theorem 2 to bound the norm ξ(t) 2 , which completes the proof.
Lemma 9 suggests that, since n(t) → 0 when (B, S) → (B * , S * ), the differential equation (89) for the function c(t) is approximated by the logistic differential equation To make the statement (94) precise, we define the function c b (t, x), for any scalar x with |x| < w, as where the constant Υ (x) is set such that c b (0, x) = c(0), i.e., Lemma 10 states an upper and a lower bound on the function c(t).
When (B, S) approaches (B * , S * ), Theorem 2 states that the error term ξ(t) is negligible and, furthermore, Lemma 10 states that the function c(t) converges to c b (t, 0). Thus, based on the ansatz (15), we approximate the viral state v(t) by With the definition of the function c b (t, x) in (95) Then, it follows from the ansatz (15) that the difference of the exact viral state v(t) to the approximation v apx (t) equals The norm ξ(t) 2 of the error term ξ(t) is bounded by Theorem 2. Thus, it remains to bound the first addend of (98). With Lemma 10, the difference of the function c(t) to c b (t, 0) is bounded by Furthermore, the scalar κ converges to zero when (B, S) approaches (B * , S * ). Hence, if we show that, as the scalar κ converges to zero, the upper bound c b (t, κ) converges to the lower bound c b (t, −κ) then (99) implies that the function c(t) converges to c b (t, 0). Furthermore, we must show that the upper bound c b (t, κ) converges to the lower bound c b (t, −κ) uniformly in time t, since the upper bound on the approximation error v(t) − v apx (t) 2 in Theorem 3 does not depend on time t. From the definition of the function c b (t, x) in (95) we obtain that where we denote Lemma 10 states that κ = O((R 0 − 1) s ) for some s > 1 when (B, S) approaches (B * , S * ). Furthermore, Lemma 8 states that w = O(R 0 − 1). Hence, it holds that For small x, the series expansion of the square root yields that 1 2 Thus, for small values of κ/w, we obtain from (100) that Since the magnitude of the hyperbolic tangent is bounded by 1, it follows from the definition of the function g(t, κ) in (101) that which yields that since κ/w = O((R 0 − 1) s−1 ). The last two addends of (102) are independent of time t. Thus, it remains to show that first addend, i.e., the difference (g(t, κ) − g(t, −κ)), converges to zero uniformly in time t as κ → 0.

Proof
The mean value theorem gives that for some z(t) ∈ (0, κ). Thus, it holds that for some z 1 (t) ∈ (0, κ) and z 2 (t) ∈ (−κ, 0), which yields that To express the derivative of the function g(t, κ), we write the function g(t, x) as where we define the function h(t, κ) as Then, the derivative of the function g(t, κ) with respect to the scalar κ is given by which is upper-bounded by With the derivative of the function h(t, κ), i.e.

H Proof of Theorem 4
We acknowledge the help of Karel Devriendt, who constructed an effective infection rate matrix of homogeneous NIMFA with a given principal eigenvector x 1 . The idea of proving Theorem 4 is based on Corollary 2: When R 0 ↓ 1, the viral state dynamics of heterogeneous NIMFA (1) are determined by the four variables x 1 , w, γ, Υ (0). Thus, we aim to show that the corresponding four variables of the homogeneous NIMFA system (26), which we denote by x 1,hom , w hom , γ hom and Υ hom (0), are the same as the variables x 1 , w, γ, Υ (0) of heterogeneous NIMFA (1).
Proof First, we consider the principal eigenvector x 1 . The effective infection rate matrix of the homogeneous NIMFA system (26) We show that the principal eigenvector x 1 of heterogeneous NIMFA (1) is also the principal eigenvector x 1,hom of the matrix W hom . Indeed, Thus, x 1 is an eigenvector of the effective infection rate matrix W hom of the homogeneous NIMFA system (26). The corresponding eigenvalue equals λ 1,hom = β hom δ hom 1 min l=1,...,N The effective infection rate matrix W hom is non-negative and irreducible, by definition (105). Thus, the Perron-Frobenius Theorem (Van Mieghem 2010) yields that the eigenvalue λ 1,hom to the positive eigenvector x 1 equals the spectral radius ρ (W hom ) = λ 1,hom and that x 1,hom = x 1 . Second, we consider the variables γ , γ hom in Theorem 1. By definition (7) and since x 1 is a vector of length 1, it holds that where the last equality follows from (106). With (28), we obtain further that Thus, the variable γ hom of the homogeneous NIMFA (26) equals the variable γ of heterogeneous NIMFA (1). Third, we show that the viral slope w hom of the homogeneous NIMFA (26) equals the viral slope w of heterogeneous NIMFA (1). From the definition (87), the variable w hom of the homogeneous NIMFA system (26) follows as w hom = λ 1,hom − 1 δ hom .
With (106), we obtain that w hom = β hom 1 min l=1,...,N Then, the definition of the infection rate β hom in (28) yields that which simplifies with the definition of δ hom in (27) to Then, the definition of γ in (7) yields that Thus, the viral slope w hom of the homogeneous NIMFA system (26) equals the viral slope w of heterogeneous NIMFA (1), which completes the proof.
In contrast to the variables x 1 , γ, w in Lemma 12, the two variables Υ hom (0) and Υ (0), given by definition (19), are not necessarily equal, since the steady states v ∞ and v ∞,hom might be different. For the homogeneous NIMFA system (26) and heterogeneous NIMFA (1), we denote the viral state approximations of Corollary 2 bỹ v apx (t) andṽ apx,hom (t), respectively. The difference of the viral state vectors v(t) and v hom (t) can be written as − v hom (t) −ṽ apx,hom (t) .