==== Front J Math Biol J Math Biol Journal of Mathematical Biology 0303-6812 1432-1416 Springer Berlin Heidelberg Berlin/Heidelberg 32959068 1542 10.1007/s00285-020-01542-6 Article Time-dependent solution of the NIMFA equations around the epidemic threshold http://orcid.org/0000-0002-7935-9109Prasse Bastian b.prasse@tudelft.nl Van Mieghem Piet p.f.a.vanmieghem@tudelft.nl Faculty of Electrical Engineering, Mathematics and Computer Science, P.O Box 5031, 2600 GA Delft, The Netherlands 22 9 2020 22 9 2020 2020 81 6 1299 1355 26 6 2019 3 7 2020 13 9 2020 © The Author(s) 2020Open AccessThis 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/.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. Keywords NIMFA differential equationsSIS processEpidemic modelsViral state dynamicsMathematics Subject Classification 92D3092D2534A34Delft University of Technologyissue-copyright-statement© Springer-Verlag GmbH Germany, part of Springer Nature 2020 ==== Body 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, on-line social media content and innovations. Furthermore, man-made structures exhibit epidemic phenomena, such as the propagation of failures in power networks or the 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 group1 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 aij. If there is a link between node i and node j, then aij=1, and aij=0 otherwise. Hence, the virus directly spreads between two nodes i and j only if aij=1. We stress that in most applications it holds that aii≠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 vi(t). The viral state vi(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 βij depend on the nodes i and j. Definition 1 (Heterogeneous NIMFA) At any time t≥0, the NIMFA governing equation is 1 dvi(t)dt=-δivi(t)+1-vi(t)∑j=1Nβ~ijaijvj(t) for every group i=1,…,N, where δi>0 is the curing rate of node i, and β~ij>0 is the infection rate from node j to node i. For a vector x∈RN, 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 2 dv(t)dt=-Sv(t)+diagu-v(t)Bv(t), where v(t)=(v1(t),…,vN(t))T is the viral state vector at time t, the N×N infection rate matrix B is composed of the elements βij=β~ijaij, 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 3 B-Sv∞=diagv∞Bv∞. 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 4 dv(t)dt=-δv(t)+βdiagu-v(t)Av(t). For the vast majority of epidemiological, demographical, and ecological models, the basic reproduction number R0 is an essential quantity (Hethcote 2000; Heesterbeek 2002). The basic reproduction number R0 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 R0 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 R0 to epidemic models with N>1 groups. For NIMFA (1), the basic reproduction number R0 follows (Van den Driessche and Watmough 2002) as R0=ρ(S-1B), 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 R0 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 R0=ρ(S-1B) 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 R0=1: If R0≤1, then the only equilibrium of NIMFA (1) is the origin, which is globally asymptotically stable. Else, if R0>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 R0=1 is of particular interest. In practice, the basic reproduction number R0 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 βij. 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 R0 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 R0=1. More precisely, we find the closed-form expression of the viral state vi(t) for every node i at every time t when R0↓1, given that the initial state v(0) is small or parallel2 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 R0↓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 R0↓1 gives an accurate approximation also for R0>1 as demonstrated by numerical evaluations in Sect. 6. Related work 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 R0↓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 rates3βij(t) depend continuously on time t. Rami et al. (2013) consider a switched model in which both the infection rates βij(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). 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 vi(0) is the same4 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,βij. We solve NIMFA (1) for arbitrary graphs around the threshold criterion R0=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 R0=ρ(S-1B) is determined by the infection rate matrix B and the curing rate matrix S. Thus, the notation R0↓1 is imprecise, since there are infinitely many matrices B, S such that the basic reproduction number R0 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 converges5 to a limit (B∗,S∗), such that ρS∗-1B∗=1 and ρS(n)-1B(n)>1∀n∈N. 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 R0↓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 f(R0)=O(g(R0)) as R0↓1 denotes that |f(R0)|≤σ|g(R0)| for some constant σ as R0↓1. For instance, it holds that (R0-1)2=O(R0-1) as R0↓1. 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 R0>1, the curing rates are positive and the infection rates are non-negative, i.e., δi>0 and βij≥0 for all nodes i, j. Furthermore, in the limit R0↓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 βij 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 R0>1, it holds that vi(0)≥0 and vi(0)≤v∞,i for every node i=1,…,N. Furthermore, it holds that vi(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 vi(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 vi(0) would be small. Assumption 3 For every basic reproduction number R0>1, the infection rate matrix B is symmetric and irreducible. Furthermore, in the limit R0↓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 5 W=S-1B. In this section, we state an essential property that we apply to solve the NIMFA equations (1) when the basic reproduction number R0 is close to 1: The steady-state vector v∞ converges to a scaled version of the principal eigenvector x1 of the effective infection rate matrix W when R0↓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 R0 are the same quantity, i.e., R0=ρ(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 R0↓1. We generalise the results of Van Mieghem (2012) to heterogeneous NIMFA (1): Theorem 1 Under Assumptions 1 and 3, the steady-state vector v∞ obeys 6 v∞=γx1+η, where the scalar γ equals 7 γ=R0-1∑l=1Nδlx1l2∑l=1Nδlx1l3, and the N×1 vector η satisfies ‖η‖2≤OR0-12 when the basic reproduction number R0 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 R0↓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(t1),…,v(tN)), for arbitrary observation times t1<⋯0 such that the error term ξ(t) at any time t≥0 is bounded by 18 ‖ξ(t)‖2≤‖ξ(0)‖2e-σ1t+σ2(R0-1)2 when the basic reproduction number R0 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(R0-1)2 uniformly in time t when R0↓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 19 Υ(0)=artanh2v∞Tv(0)‖v∞‖22-1. Furthermore, we define the viral slope w, which determines the speed of convergence to the steady-state v∞, as w=(R0-1)∑l=1Nδlx1l2. 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, ‖ξ(0)‖2=O(R0-1)p when R0↓1. Furthermore, define 20 vapx(t)=121+tanhw2t+Υ(0)v∞. Then, there exists some constant σ>0 such that 21 ‖v(t)-vapx(t)‖2‖v∞‖2≤σ(R0-1)s-1∀t≥0, where s=min{p,2}, when the basic reproduction number R0 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 vapx(t) which is uniform in time t. Furthermore, since both the viral state v(t) and the approximation vapx(t) converge to the steady-state v∞, it holds that ‖v(t)-vapx(t)‖2→0 when t→∞. At time t=0, we obtain from Theorem 3 and (17) that ‖v(0)-vapx(0)‖2=‖ξ(0)‖2. Since ‖ξ(0)‖2=O(R0-1)p and, by Theorem 1, ‖v∞‖2=OR0-1, we obtain that ‖v(0)-vapx(0)‖2‖v∞‖2=O(R0-1)p-1. Hence, for general t≥0 the approximation error ‖v(t)-vapx(t)‖2/‖v∞‖2 does not converge to zero faster than O(R0-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 22 ‖v(t)-vapx(t)‖2‖v∞‖2≤σ(R0-1)∀t≥0. The time-dependent solution to NIMFA (1) at the epidemic threshold criterion R0=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 R0↓1. Thus, Theorem 3 implies that the convergence time to the steady-state v∞ goes to infinity when R0↓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)=r0v∞ for some scalar r0∈(0,1). Then, for any scalar r1∈[r0,1), the largest time t01 at which the viral state satisfies vi(t01)≤r1v∞,i for every node i converges to t01=1wlogr1r01-r01-r1 when the basic reproduction number R0 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, ‖ξ(0)‖2=O(R0-1)p when R0↓1. Furthermore, define 23 v~apx(t)=1+tanhw2t+Υ(0)γ2x1. Then, there exists some constant σ>0 such that ‖v(t)-v~apx(t)‖2‖v∞‖2≤σ(R0-1)s-1∀t≥0, where s=min{p,2}, when the basic reproduction number R0 approaches 1 from above. In contrast to Theorem 3, the approximation error ‖v(t)-v~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 24 vi(t)vj(t)→v~apx,i(t)v~apx,j(t)=(x1)i(x1)j at every time t when R0↓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 R0=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 x1 of the effective infection rate matrix W equals the principal eigenvector of the adjacency matrix A. Furthermore, the limit process R0↓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 v(t) of homogeneous NIMFA (4). Furthermore, suppose that ‖ξ(0)‖2=O(τ-τc)p for some constant p>1 when τ↓τc and define 25 vapx(t)=121+tanh(τ-τc)δ2τct+Υ(0)v∞. Then, there exists some constant σ>0 such that ‖v(t)-vapx(t)‖2‖v∞‖2≤σ(τ-τc)s-1∀t≥0, 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 vapx(t) defined by (25) equals the exact solution (Van Mieghem 2014b) of homogeneous NIMFA (4) on a regular graph, provided that the initial state vi(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 R0=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 R0=1 by choosing the homogeneous spreading parameters β,δ and the adjacency matrix A accordingly? Theorem 4 Consider heterogeneous NIMFA (1) with given spreading parameters βij,δi. Suppose that Assumptions 1 to 3 hold and that, for some constant p>1, ‖ξ(0)‖2=O(R0-1)p when the basic reproduction number R0 approaches 1 from above. Define the homogeneous NIMFA system 26 dvi,hom(t)dt=-δhomvi,hom(t)+βii,hom1-vi,hom(t)vi,hom(t)+1-vi,hom(t)βhom∑j=1,j≠iNvj,hom(t), where the homogeneous curing rate δhom equals 27 δhom=∑l=1Nδlx1l3∑l=1Nx1l3, the homogeneous infection rate βhom equals 28 βhom=δhom∑l=1Nx1l1+γ∑l=1Nx1l3minl=1,…,Nx1l with the variable γ defined by (7), and the self-infection rates βii,hom equal βii,hom=βhom1minl=1,…,Nx1l-1x1i∑j=1Nx1j+βhom. Then, if vhom(0)=v(0), there exists some constant σ>0 such that ‖v(t)-vhom(t)‖2‖v∞‖2≤σ(R0-1)s-1∀t≥0, where s=min{p,2}, when the basic reproduction number R0 approaches 1 from above. Proof Appendix H. □ In other words, when R0↓1, for any contact network and any spreading parameters δi,βij, 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 1to 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 vapx(t), given by (20), when the basic reproduction number R0 is close, but not equal, to one. We generate an adjacency matrix A according to different random graph models. If aij=1, then we set the infection rates βij to a uniformly distributed random number in [0.4, 0.6] and, if aij=0, then we set βij=0. We set the initial curing rates δl(0) to a uniformly distributed random number in [0.4, 0.6]. To set the basic reproduction number R0, we set the curing rates δl to a multiple of the initial curing rates δl(0), i.e. δl=σδl(0) for every node l and some scalar σ such that ρ(W)=R0. Thus, we realise the limit process R0↓1 by changing the scalar σ. Only in Sect. 6.2, we consider homogeneous spreading parameters by setting βij=0.5 and δi(0)=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., 29 dvi(t)dtt=Tk≈vi(Tk)-viT(k-1)T for a small sampling time T and a discrete time slot k∈N. In Prasse and Van Mieghem (2019), we derived an upper bound Tmax 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=Tmax/100. Except for Sect. 6.3, we set the initial viral state to v(0)=0.01v∞. We define the convergence time tconv as the smallest time t at which vi(tconv)-v∞,i≤0.01 holds for every node i. Thus, at the convergence time tconv the viral state v(tconv) has practically converged to the steady-state v∞. We evaluate Theorem 3 with respect to the approximation error ϵV, which we define as ϵV=1Ntconv∑i=1N∫0tconvvit~-vapx,it~v∞,idt~. 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 m0=5, m=2. Figure 1 gives an impression of the accuracy of the approximation of Theorem 3 around the epidemic threshold criterion R0=1. For a basic reproduction number R0≤1.1, the difference of the closed-form expression of Theorem 3 to the exact NIMFA viral state trace is negligible.Fig. 1 For a Barabási–Albert random graph with N=500 nodes, the approximation accuracy of Theorem 3 is depicted. Each of the sub-plots shows the viral state traces vi(t) of seven different nodes i, including the node i with the greatest steady-state v∞,i We aim for a better understanding of the accuracy of the closed-form expression of Theorem 3 when the basic reproduction number R0 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 pER=0.05. Figure 2 illustrates the convergence of the approximation of Theorem 3 to the exact solution of NIMFA (1). Around the threshold criterion R0=1, the approximation error ϵV converges linearly to zero with respect to the basic reproduction number R0, 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.Fig. 2 The approximation error ϵV of the NIMFA solution versus the basic reproduction number R0 for different network sizes N Impact of degree heterogeneity on the approximation accuracy For NIMFA (4) with homogeneous spreading parameters β,δ, the approximation vapx(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) random graphs with N=100 nodes and an average node degree of 4. We vary the link rewiring probability pWS from pWS=0, which correspond to a regular graph, to pWS=1, which corresponds to a “completely random” graph. Figure 3 depicts the approximation error ϵV versus the rewiring probability pWS for homogeneous spreading parameters β,δ. Interestingly, the approximation error reaches a maximum and improves when the adjacency matrix A is more random.Fig. 3 The approximation error ϵV versus the link rewiring probability pWS for Watts–Strogatz random graphs with N=100 nodes and homogeneous spreading parameters β,δ 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 vi(0) of every node i to a uniformly distributed random number in (0,r0v∞,i] for some scalar r0∈(0,1]. By increasing the scalar r0, 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. Figure 5 shows that the viral state v(t) converges rapidly to the approximation vapx(t) as time t increases. For general initial viral states v(0) with ξ(0)≠0, it holds that vapx(0)≠v(0) since the approximation vapx(0) is parallel to the steady-state vector v∞. Hence, the approximation vapx(t) does not converge point-wise to the viral state v(t) when R0↓1. However, based on the results shown in Figs. 4 and 5, we conjecture convergence with respect to the L2-norm for general initial viral states v(0) when R0↓1. Conjecture 1 Suppose that Assumptions 1 to 3 hold. Then, it holds for the approximation vapx(t) defined by (20) that 1‖v∞‖2∫0∞‖v(t)-vapx(t)‖2dt→0 when the basic reproduction number R0 approaches 1 from above. Fig. 4 The approximation error ϵV versus the scalar r0, which controls the variance of the randomly generated initial viral state v(0), for Barabási–Albert networks with N=250 nodes Fig. 5 For a Barabási–Albert random graph with N=500 nodes, a basic reproduction number R0=1.01 and a randomly generated initial viral state v(0), the approximation accuracy of Theorem 3 is depicted. The viral state traces vi(t) of seven different nodes i are depicted 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 vapx(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 vapx(t) when the basic reproduction number R0 approaches 1 from above. Fig. 6 The approximation error ϵV of the NIMFA solution versus the basic reproduction number R0 for directed Erdős–Rényi graphs for different network sizes N Fig. 7 The approximation error ϵt of the convergence time t01 versus the basic reproduction number R0 for different network sizes N Accuracy of the approximation of the convergence time Corollary 1 gives the expression of the convergence time t01 from the initial viral state v(0)=r0v∞ to the viral state v(t01)≤r1v∞ for any scalars 0λ2≥⋯≥λN N Number of nodes ‖M‖2 2-norm of a matrix M: largest singular value of M r(M) Numerical radius of a square matrix M; see (30) R0 Basic reproduction number; R0=ρ(W)=λ1 ρ(M) Spectral radius of a square matrix M S Curing rate matrix; S=diag(δ1,…,δN) u N×1 all-one vector u=(1,…,1)T v(t) N×1 viral state vector v(t) at time t≥0; vi(t)∈[0,1] for i=1,…,N v∞ Steady-state vector, see Definition 2 w The viral slope; w=(R0-1)∑l=1Nδlx1l2 W Effective infection rate matrix W=S-1B; ρ(W)>1 W~ Symmetric N×N matrix W~=S-12BS-12 xk k-th eigenvector of the matrix W; Wxk=λkxk ξ(t) Viral state component that is orthogonal to v∞; ξ(t)=v(t)-c(t)v∞ Proof of Theorem 1 The steady-state v∞ solely depends on the effective infection rate matrix W: By left-multiplication of (3) with the diagonal matrix S-1, we obtain that 31 W-Iv∞=diagv∞Wv∞. In general, the effective infection rate matrix W, defined in (5) as W=S-1B, 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 32 W~=S-12BS-12=S12WS-12. Since the infection rate matrix B is symmetric under Assumption 3, the matrix W~ is symmetric. Hence, the matrix W~, and also the effective infection rate matrix W, are diagonalisable. With (32), we write the steady-state (31) with respect to the symmetric matrix W~ as 33 W~-IS12v∞=diagv∞W~S12v∞. We decompose the matrix W~ as 34 W~=λ1x~1Tx~1+∑k=2Nλkx~kTx~k, where the eigenvalues of W~ are real and equal to λ1>λ2≥⋯≥λN with the corresponding normalized eigenvectors denoted by x~1,…,x~N. Then, the steady-state vector v∞ can be expressed as linear combination v∞=∑l=1Nψlx~l, where the coefficients equal ψl=v∞Tx~l. To prove Theorem 1, we would like to express the coefficients ψ1,…,ψN as a power series around R0=1. However, in the limit process B,S→B∗,S∗, the eigenvectors x~1,…,x~N of the matrix W~ are not necessarily constant. Hence, the coefficients ψl depend on the full matrix W~ and not only on the basic reproduction number R0. To overcome the challenge of non-constant eigenvectors x~1,…,x~N in the limit process B,S→B∗,S∗, we define the symmetric auxiliary matrix 35 M(z)=zx~1Tx~1+∑k=2Nλkx~kTx~k∈RN×N for a scalar z≥1. Thus, the matrix M(z) is obtained from the matrix W~ by replacing the largest eigenvalue λ1 of W~ by z. In particular, the definition of the matrix M(z) in (35) and (34) illustrate that M(λ1)=W~. When the matrix W~ is formally replaced by the matrix M(z), the steady-state equation (33) becomes 36 M(z)-IS12v~(z)=diagv~(z)M(z)S12v~(z) where the N×1 vector v~(z) denotes the solution of (36). Since M(R0)=W~, the solution of (36) at z=R0 and the solution to (33) coincide, i.e., v~(R0)=v∞. Lemma 2 expresses the solution of the equation (36) as a power series. Lemma 2 Suppose that Assumptions 1 and 3 hold. If (B, S) is sufficiently close to (B∗,S∗), then the N×1 vector v~(z) which satisfies (36) equals 37 v~(z)=(z-1)∑l=1N1δlx~1l3-1S-12x~1+ϕ(z), where the N×1 vector ϕ(z) satisfies ‖ϕ(z)‖2≤σ(B,S)(z-1)2 for some scalar σ(B,S) when z approaches 1 from above. Proof The proof is an adaptation of the proof (Van Mieghem 2012, Lemma 4). We express the solution v~(z) of (36) as linear combination of the vectors S-12x~1,…,S-12x~N, i.e., 38 v~(z)=∑k=1Nψk(z)S-12x~k. Since the diagonal matrix S-12 is full rank, the vectors S-12x~k, where k=1,…,N, are linearly independent. Furthermore, we express the coefficients ψk(z) as a power series 39 ψk(z)=∑j=0∞gj(k)z-1j, where g0(k)=0 for every k=1,…,N, since (Lajmanovich and Yorke 1976) it holds that v~(z)=0 when z=1. We denote the eigenvalues of the matrix M(z) by 40 λk(z)=zifk=1,λkifk≥2. By substituting (38) into (36), we obtain that ∑k=1Nλk(z)-1ψk(z)x~k=diag∑l=1Nψl(z)x~lS-12∑k=1Nλk(z)ψk(z)x~k and left-multiplying with the eigenvector x~mT, for any m=1,…,N, yields 41 λm(z)-1ψm(z)=∑n=1Nx~mn∑l=1Nψl(z)x~ln1δn∑k=1Nψk(z)λk(z)x~kn. We define X(m,l,k)=∑n=1N1δnx~mnx~lnx~kn. Then, we rewrite (41) as 42 λm(z)-1ψm(z)=∑l=1N∑k=1Nψl(z)ψk(z)λk(z)X(m,l,k). First, we focus on the left-hand side of (42), which we denote by θm(z)=λm(z)-1ψm(z). With the power series (39), we obtain that θm(z)=λm(z)-1∑j=1∞gj(m)z-1j. Further rewriting yields that 43 θm(z)=λm(z)-z+z-1∑j=1∞gj(m)z-1j=λm(z)-z∑j=1∞gj(m)z-1j+∑j=1∞gj(m)z-1j+1=∑j=1∞λm(z)-zgj(m)z-1j+∑j=2∞gj-1(m)z-1j. Second, we rearrange the right-hand side of (42) as θm(z)=λ1(z)∑l=1Nψl(z)ψ1(z)X(m,l,1)+∑l=1N∑k=2Nψl(z)ψk(z)λk(z)X(m,l,k). By the definition of λk(z) in (40) it holds that λ1(z)=z, and we obtain that 44 θm(z)=(z-1)∑l=1Nψl(z)ψ1(z)X(m,l,1)+∑l=1N∑k=1Nψl(z)ψk(z)λ~kX(m,l,k), where λ~k=1ifk=1,λkifk≥2. Introducing the power series (39) into (44) and executing the Cauchy product for ψl(z)ψk(z) yields θm(z)=∑j=1∞∑n=1j-1∑l=1Ngj-n(1)gn(l)X(m,l,1)z-1j+1+∑j=1∞∑n=1j-1∑l=1N∑k=1Ngj-n(l)gn(k)λ~kX(m,l,k)z-1j. We shift the index j in the first term and obtain 45 θm(z)=∑j=2∞∑n=1j-2∑l=1Ngj-1-n(1)gn(l)X(m,l,1)z-1j+∑j=1∞∑n=1j-1∑l=1N∑k=1Ngj-n(l)gn(k)λ~kX(m,l,k)z-1j. Finally, we equate powers in (z-1)j in (43) and (45), which yields for j=1 that 46 λm(z)-zg1(m)=0 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 g1(m)=0 when m≥2 provided that (B, S) is sufficiently close to (B∗,S∗), since z≥1. For j≥2, equating powers in (45) yields that 47 λm(z)-zgj(m)+gj-1(m)=∑n=1j-2∑l=1Ngj-1-n(1)gn(l)X(m,l,1)+∑n=1j-1∑l=1N∑k=1Ngj-n(l)gn(k)λ~kX(m,l,k). In particular, for the case j=2, we obtain 48 λm(z)-zg2(m)+g1(m)=∑l=1N∑k=1Ng1(l)g1(k)λ~kX(m,l,k)=g1(1)g1(1)X(m,1,1), since g1(l)=0 for all l≥2 and λ~1=1. Since λ1(z)=z, we obtain for m=1 from (48) that g1(1)=g1(1)2X(1,1,1) and, hence, g1(1)=1X(1,1,1)=∑l=1N1δlx~1l3-1. Since g1(m)=0 for m≥2, we obtain that the power series (38) for the solution v~(z) of (36) becomes 49 v~(z)=(z-1)g1(1)S-12x~1+ϕ(z), where the N×1 vector ϕ(z) equals ϕ(z)=∑k=1N∑j=2∞gj(k)z-1jS-12x~k. Thus, it holds ‖ϕ(z)‖2=O(z-1)2 when z approaches 1 from above, which proves Lemma 2. □ We believe that, based on (47), a recursion for the coefficients gj(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 v~(z) in (37) in terms of the principal eigenvector x1 of the effective infection rate matrix W, we propose Lemma 3. Lemma 3 Under Assumptions 1 and 3, it holds that 50 ∑l=1N1δlx~1l3-1S-12x~1=∑l=1Nδlx1l2∑l=1Nδlx1l3x1. Proof From (32), it follows that the principal eigenvector x~1 of the matrix W~ and the principal eigenvector x1 of the effective infection rate matrix W are related via x~1=1S12x12S12x1, or, component-wise, (x~1)l=1S12x12δl(x1)l. Then, we rewrite the left-hand side of (50) as ∑l=1N1δlx~1l3-1S-12x~1=∑l=1Nδlx1l3-1S12x122x1, which simplifies to ∑l=1N1δlx~1l3-1S-12x~1=x1TSx1∑l=1Nδlx1l3x1. Writing out the quadratic form in the numerator completes the proof. □ The basic reproduction number R0 converges to 1 when (B,S)→(B∗,S∗). Hence, if (B, S) is sufficiently close to (B∗,S∗), then the basic reproduction number R0 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 v~(R0) to (36) at z=R0 follows with Lemma 2 as v~(R0)=(R0-1)∑l=1N1δlx~1l3-1S-12x~1+ϕ(R0)=γx1+ϕ(R0), where the last equality follows from Lemma 3 and the definition of the scalar γ in (7). We emphasise that Lemma 2 implies that γ=O(R0-1) and, hence, ‖v~(R0)‖2=O(R0-1) as (B,S)→(B∗,S∗). Since M(R0)=W~, the solution of (36) at z=R0 and the solution to (33) coincide, i.e., v~(R0)=v∞. Thus, from the definition of the vector η in (6), we obtain that 51 η2=v∞-γx12=ϕ(R0)2 when (B,S)→(B∗,S∗). Lemma 2 states that ϕ(z)2=O(z-1)2 as z↓1. Hence, we obtain from (51) that 52 η2≤σ(B,S)(R0-1)2 for some scalar σ(B,S) when (B,S)→(B∗,S∗). Furthermore, when (B, S) converge to the limit (B∗,S∗), the scalar σ(B,S) converges to some limit σ(B∗,S∗). Hence, by defining the constant σ=σ(B∗,S∗)+ϵσ for some ϵσ>0, it holds that σ(B,S)<σ, for all (B, S) which are sufficiently close to (B∗,S∗). Finally, we obtain from (52) that η2≤σ(R0-1)2 when (B, S) approaches (B∗,S∗). 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]. Absence of overshoot The proof follows the same reasoning as Prasse and Van Mieghem (2019, Corollary 1). Assume that at some time t0 it holds vi(t0)=v∞,i for some node i and that vj(t0)≤v∞,j for every node j. Since vi(t0)=v∞,i, the NIMFA equation (1) yields that dvi(t)dtt=t0=-δiv∞,i+(1-v∞,i)∑j=1Nβijvj(t0). Since vj(t0)≤v∞,j for every node j, we obtain that dvi(t)dtt=t0≤-δiv∞,i+(1-v∞,i)∑j=1Nβijv∞,j=0, where the last equality follows from the steady-state equation (3). Thus, vi(t0)=v∞,i implies that dvi(t)dtt=t0≤0, which means that, at time t0, the viral state vi(t0) does not increase. Hence, the viral state vi(t0) cannot exceed the steady-state v∞,i at any time t≥0. Boundedness of the function c(t) Relation (16) indicates that 53 c(t)=1‖v∞‖22v∞Tv(t)=1‖v∞‖22v∞,1v1(t)+⋯+v∞,NvN(t) Section C.1 shows that Assumption 2 implies that vi(t)≤v∞,i for all nodes i and every time t. Thus, we obtain from (53) that c(t)≤1‖v∞‖22v∞,1v∞,1+⋯+v∞,Nv∞,N=1 Analogously, since vi(t)≥0 for all nodes i and every time t, we obtain from (53) that c(t)≥0. Proof of Theorem 2 By inserting the ansatz (15) into the NIMFA equations (2), we obtain that 54 dc(t)dtv∞+dξ(t)dt=Λ1(t)+Λ2(t). Here, the function Λ1(t) is given by Λ1(t)=B-Sc(t)v∞-c2(t)diag(v∞)Bv∞, which simplifies, with the steady-state equation (3), to 55 Λ1(t)=c(t)-c2(t)B-Sv∞. The function Λ2(t) is given by Λ2(t)=B-Sξ(t)-c(t)diag(ξ(t))Bv∞-c(t)diag(v∞)Bξ(t)-diag(ξ(t))Bξ(t). With diag(ξ(t))Bv∞=diag(Bv∞)ξ(t), we obtain that 56 Λ2(t)=B-S-c(t)diag(Bv∞)-c(t)diag(v∞)Bξ(t)-diag(ξ(t))Bξ(t). 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)‖22. The convergence of the squared norm ‖ξ(t)‖22 to zero implies the convergence of the error term ξ(t) to zero. The derivative of the squared norm ‖ξ(t)‖22 is given by d‖ξ(t)‖22dt=2ξT(t)dξ(t)dt. Thus, we obtain from (54) that 57 12d‖ξ(t)‖22dt=ξT(t)Λ1(t)+ξT(t)Λ2(t), 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)‖22. Upper bound on ξT(t)Λ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 ξT(t)Λ1(t)≤14δmaxγ(R0-1)+(R0+1)‖η‖2‖ξ(t)‖2. Proof From (55) and the definition of the matrix W in (5) it follows that ξT(t)Λ1(t)=c(t)-c2(t)ξT(t)SW-Iv∞. With Theorem 1, we obtain ξT(t)Λ1(t)=c(t)-c2(t)γ(R0-1)ξT(t)Sx1+ξT(t)S(W-I)η. The triangle inequality yields that 58 ξT(t)Λ1(t)≤c(t)-c2(t)|γ(R0-1)|ξT(t)Sx1+ξT(t)S(W-I)η. With the Cauchy–Schwarz inequality, the first addend in (58) is upper-bounded by ξT(t)Sx1≤‖STξ(t)‖2‖x1‖2=‖Sξ(t)‖2, since ‖x1‖2=1 and the matrix S is symmetric. The matrix 2-norm is sub-multiplicative, which yields that ξT(t)Sx1≤‖S‖2‖ξ(t)‖2=δmax‖ξ(t)‖2. Thus, (58) gives that 59 ξT(t)Λ1(t)≤c(t)-c2(t)γ(R0-1)δmax‖ξ(t)‖2+ξT(t)S(W-I)η, since γ>0 and R0>1. We consider the second addend in (59), which we write with (32) as ξT(t)S(W-I)η=ξT(t)S12(W~-I)S12η. From the Cauchy–Schwarz inequality and the sub-multiplicativity of the matrix norm we obtain ξT(t)S(W-I)η≤‖ξ(t)‖2‖S12‖2‖W~-I‖2‖S12‖2‖η‖2. The triangle inequality and the symmetry of the matrix W~ imply that ‖W~-I‖2≤‖W~‖2+‖I‖2=R0+1. Thus, we can upper bound the second added in (59) by ξT(t)S(W-I)η≤δmax(R0+1)‖ξ(t)‖2‖η‖2, since ‖S12‖2=δmax. Hence, (59) yields the upper bound ξT(t)Λ1(t)≤c(t)-c2(t)δmaxγ(R0-1)+(R0+1)‖η‖2‖ξ(t)‖2. Finally, Lemma 1 states that 0≤c(t)≤1, which implies that c(t)-c2(t)≤1/4 and completes the proof. □ Upper bound on ξT(t)Λ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 ξT(t)Λ2(t)≤-‖S12ξ(t)‖22+ξT(t)diag(u-c(t)v∞)Bξ(t). Proof From (56) it follows that 60 ξT(t)Λ2(t)=ξT(t)B-S-c(t)diag(Bv∞)-c(t)diag(v∞)Bξ(t)-ξT(t)diag(ξ(t))Bξ(t). 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 61 -ξT(t)diag(ξ(t))Bξ(t)=∑l=1Nξl2(t)∑j=1Nβlj-ξj(t). Since v(t)=c(t)v∞+ξ(t) and vi(t)≥0 for every node i at every time t, it holds that 62 -ξi(t)≤c(t)v∞,i,i=1,…,N. By inserting (62) in (61), the last addend of (60) is upper bounded by -ξT(t)diag(ξ(t))Bξ(t)≤∑l=1Nξl2(t)∑j=1Nβljc(t)v∞,j, which simplifies to 63 -ξT(t)diag(ξ(t))Bξ(t)≤c(t)ξT(t)diag(Bv∞)ξ(t). By applying the upper bound (63) to (60), we obtain that ξT(t)Λ2(t)≤ξT(t)B-S-c(t)diag(v∞)Bξ(t). With the definition of the matrix W~ in (32), we obtain ξT(t)Λ2(t)≤ξT(t)S12W~-I-c(t)diag(v∞)W~S12ξ(t), and further rearranging completes the proof. □ For any scalar ς∈[0,1] and any vector υ∈RN, we define Θ(ς,υ,B,S)=υTdiag(u-ςv∞)Bυ‖S12υ‖22. Then, we obtain from Lemma 5 that 64 ξT(t)Λ2(t)≤Θ(c(t),ξ(t),B,S)-1‖S12ξ(t)‖22. 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∈RN of length zTz=1 such that zTMz=r(M). Furthermore, if z~TMz~=r(M) holds for a vector z~∈RN of length z~Tz~=1, then either z~=z or z~=-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 S=υ∈RN|∃i,j:υj>0>υi. Then, it holds Θ(ς,υ,B,S)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 Θ(ς,υ,B,S)0 such that ξT(t)Λ2(t)≤-ωδmax‖ξ(t)‖22 holds at every time t≥0 when B,S approaches B∗,S∗. Proof We denote the maximum of the function Θ(ς,υ,B,S) with respect to ς∈[0,1] and υ∈S by 66 Θmax(B,S)=maxς∈[0,1],υ∈SΘ(ς,υ,B,S). 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 67 Θ(ς,υ,B∗,S∗)=1‖υ~‖22υ~TW~∗υ~, where we denote W~∗=S∗-12B∗S∗-12. Since it holds R0=1 at the limit (B∗,S∗), Lemma 6 implies that 68 Θmax(B∗,S∗)<1. 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 Θ(ς,υ,B,S)=1‖υ~‖22υ~TW~∗-ςdiag(v∞)W~∗+diag(u-ςv∞)ΔW~υ~, which is upper-bounded by 69 Θ(ς,υ,B,S)≤1‖υ~‖22υ~TW~∗υ~+1‖υ~‖22υ~Tςdiag(v∞)W~∗υ~+1‖υ~‖22υ~Tdiag(u-ςv∞)ΔW~υ~. Maximising every addend in (69) independently yields an upper bound on Θmax(B,S) as 70 Θmax(B,S)≤maxς∈[0,1],υ∈S1‖υ~‖22υ~TW~∗υ~+maxς∈[0,1],υ∈S1‖υ~‖22υ~Tςdiag(v∞)W~∗υ~+maxς∈[0,1],υ∈S1‖υ~‖22υ~Tdiag(u-ςv∞)ΔW~υ~. In the following, we state upper bounds for each of the three addends in (67) separately. With (67), we write the first addend in (70) as 71 maxς∈[0,1],υ∈S1‖υ~‖22υ~TW~∗υ~=maxς∈[0,1],υ∈SΘ(ς,υ,B∗,S∗)=Θmax(B∗,S∗), where the last equality follows from the definition of Θmax(B∗,S∗) in (66). Regarding the second addend in (70), it holds that maxς∈[0,1],υ∈S1‖υ~‖22υ~Tςdiag(v∞)W~∗υ~≤maxς∈[0,1]maxυ∈RN1‖υ~‖22υ~Tςdiag(v∞)W~∗υ~=maxς∈[0,1]rςdiag(v∞)W~∗, where the last equality follows from the definition the numerical radius. Hence, the second addend in (70) is upper-bounded by 72 maxς∈[0,1],υ∈S1‖υ~‖22υ~Tςdiag(v∞)W~∗υ~≤rςoptdiag(1)(v∞)W~∗ for some ςopt(1)∈[0,1]. Similarly, we obtain an upper bound on the third addend in (70) as 73 maxς∈[0,1],υ∈S1‖υ~‖22υ~Tdiag(u-ςv∞)ΔW~υ~≤rdiag(u-ςopt(2)v∞)ΔW~ for some ςopt(2)∈[0,1]. With (71), (72) and (73), we obtain from (70) that 74 Θmax(B,S)≤Θmax(B∗,S∗)+rςoptdiag(1)(v∞)W~∗+rdiag(u-ςopt(2)v∞)ΔW~. The numerical radius r(M) is a vector8 norm (Horn and Johnson 1990) on the space of N×N matrices M. Thus, the numerical radius r(M) converges to zero if the matrix M converges to zero. Since v∞→0 and ΔW~→0 as (B,S)→(B∗,S∗) and ςopt(1),ςopt(2) are bounded, the last two addends in (74) converge to zero as (B,S)→(B∗,S∗). Hence, for every scalar ω>0 there is a ϑ(ω) such that ‖B-B∗‖2<ϑ(ω) and ‖S-S∗‖2<ϑ(ω) implies that 75 Θmax(B,S)≤Θmax(B∗,S∗)+ω. We choose the scalar ω=(1-Θmax(B∗,S∗))/2, which is positive due to (68). Then, the right-hand side of (75) becomes Θmax(B∗,S∗)+ω=12+12Θmax(B∗,S∗)=1-ω. Thus, we obtain from (75) that 76 Θmax(B,S)≤1-ω 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 Θ(c(t),ξ(t),B,S)≤Θmax(B,S). With (76), we obtain from (64) that ξT(t)Λ2(t)≤-ω‖S12ξ(t)‖22. From the sub-multiplicativity of the matrix norm, we obtain ξT(t)Λ2(t)≤-ω‖S12‖22‖ξ(t)‖22, which completes the proof, since ‖S12‖22=δmax. □ Bound on the error vector ξ(t) With Lemma 4 and Lemma 7, we upper-bound (57) by 12d‖ξ(t)‖22dt≤14δmaxγ(R0-1)+(R0+1)‖η‖2‖ξ(t)‖2-ωδmax‖ξ(t)‖22. From d‖ξ(t)‖2dt=12‖ξ(t)‖2d‖ξ(t)‖22dt, it follows that d‖ξ(t)‖2dt≤14δmaxγ(R0-1)+(R0+1)‖η‖2-ωδmax‖ξ(t)‖2. We denote 77 φB,S=14γ(R0-1)+(R0+1)‖η‖2, and we obtain that 78 d‖ξ(t)‖2dt≤φB,Sδmax-ωδmax‖ξ(t)‖2. The upper bound (78) is a linear first-order ordinary differential inequality, which is solved by (Arfken and Weber 1999) ‖ξ(t)‖2≤e-ωδmaxt‖ξ(0)‖2+∫0tφB,Sδmaxeωδmaxt~dt~, which simplifies to ‖ξ(t)‖2≤‖ξ(0)‖2-φB,Sωe-ωδmaxt+φB,Sω. The triangle inequality yields that 79 ‖ξ(t)‖2≤‖ξ(0)‖2e-ωδmaxt+φB,Sω1+e-ωδmaxt. Furthermore, since e-ωδmaxt≤1 at every time t≥0, we obtain from (79) that 80 ‖ξ(t)‖2≤‖ξ(0)‖2e-ωδmaxt+2φB,Sω. 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 σ1=ω(δmax∗-ϵ). Then, it holds that σ1<ωδmax when B,S approaches B∗,S∗, and we obtain from (80) that 81 ‖ξ(t)‖2≤‖ξ(0)‖2e-σ1t+2φB,Sω. Theorem 1 states that γ=O(R0-1) and ‖η‖2=O(R0-1)2 when B,S approaches B∗,S∗. Thus, it follows from the definition of the term φB,S in (77) that φB,S=O(R0-1)2. Hence, there is a constant σ2>0 such that (81) yields ‖ξ(t)‖2≤‖ξ(0)‖2e-σ1t+σ2R0-12 when B,S approaches B∗,S∗. Proof of Theorem 3 By projecting the differential equation (54) onto the steady-state vector v∞, we obtain that dc(t)dtv∞Tv∞=v∞TΛ1(t)+v∞TΛ2(t), since v∞Tξ(t)=0 by definition of the error term ξ(t). We divide by ‖v∞‖22 and obtain with (55) that 82 dc(t)dt=c(t)-c2(t)v∞TB-Sv∞‖v∞‖22+v∞TΛ2(t)‖v∞‖22. The first addend in the differential equation (82) can be expressed in a simpler manner when B,S approaches B∗,S∗: Lemma 8 Under Assumptions 1 and 3, it holds 83 v∞TB-Sv∞v∞Tv∞=(R0-1)x1TSx1+ζ, where ζ=O(R0-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 v∞TB-Sv∞=(γx1+η)TS(W-I)γx1+(B-S)η=(γx1+η)T(γ(R0-1)Sx1+(B-S)η), where the last equality follows from Wx1=R0x1. Thus, it holds that 84 v∞TB-Sv∞=γ2(R0-1)x1TSx1+γx1T(B-S)η+γ(R0-1)ηTSx1+ηT(B-S)η. Under Assumption 3, both matrices B and S are symmetric, which implies that x1T(B-S)T=(B-S)x1=S(R0-1)x1. Hence, we obtain from (84) that v∞TB-Sv∞=γ2(R0-1)x1TSx1+γ(R0-1)x1TSη+γ(R0-1)ηTSx1+ηT(B-S)η. Since γ=O(R0-1) and ‖η‖2=O(R0-1)2, we finally rewrite the numerator of the left-hand side of (83) as 85 v∞TB-Sv∞=γ2(R0-1)x1TSx1+O(R0-1)4. With Theorem 1, the denominator of the left-hand side of (83) equals 86 v∞Tv∞=γ2+2γηTx1+‖η‖22=γ2+O(R0-1)3. Combining the approximate expressions for the numerator (85) and the denominator (86) completes the proof. □ We define the viral slope w as 87 w=(R0-1)x1TSx1 and the function n(t) as 88 n(t)=c(t)-c2(t)ζ+v∞TΛ2(t)‖v∞‖22. Then, we obtain from (82) that 89 dc(t)dt=c(t)-c2(t)w+n(t). 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. Lemma 9 Under Assumptions 1 to 3, it holds at every time t≥0 that n(t)≤σ1‖ξ(0)‖2e-σ2t+σ3(R0-1)2 for some constants σ1,σ2,σ3>0 when (B, S) approaches (B∗,S∗). Proof Regarding the first addend in the definition of the function n(t) in (88), Lemma 1 implies that 0≤c(t)-c2(t)≤1/4 at every time t. Hence, Lemma 8 yields that there is a constant σ~0 such that c(t)-c2(t)ζ≤σ~0(R0-1)2 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 (56) that v∞TΛ2(t)‖v∞‖22=1‖v∞‖22v∞TB-S-c(t)diag(Bv∞)-diag(v(t))Bξ(t), since v(t)=c(t)v∞+ξ(t). Thus, it holds that v∞TΛ2(t)‖v∞‖22=1‖v∞‖22v∞T-S+diag(u-v(t))B-c(t)diag(Bv∞)ξ(t). With the definition of the matrix W~ in (32), we obtain that v∞TΛ2(t)‖v∞‖22=1‖v∞‖22v∞TS12(-I+diag(u-v(t))W~-c(t)S-12diag(Bv∞)S-12)S12ξ(t). The Cauchy–Schwarz inequality yields an upper bound as v∞TΛ2(t)‖v∞‖22≤1‖v∞‖22S12ξ(t)2··-I+diag(u-v(t))W~-c(t)S-12diag(Bv∞)S-12S12v∞2 With S12ξ(t)2≤δmaxξ(t)2 and the triangle inequality, we obtain 90 v∞TΛ2(t)‖v∞‖22≤δmaxξ(t)2‖v∞‖22W~-IS12v∞2+δmaxξ(t)2‖v∞‖22diag(v(t))W~2S12v∞2+δmaxξ(t)2‖v∞‖22c(t)S-12diag(Bv∞)S-122S12v∞2. In the following, we consider the three addends in (90) separately. Regarding the first addend, we obtain with the definition of the matrix W~ in (32) that W~-IS12v∞=S12W-Iv∞=γ(R0-1)S12x1+S12W-Iη, where the last equality follows from Theorem 1. Thus, the triangle inequality yields W~-IS12v∞2≤γ(R0-1)S12x12+S12W-Iη2. With the sub-multiplicativity of the matrix 2-norm, we obtain W~-IS12v∞2≤δmaxγ(R0-1)+(R0+1)η2, since W-I2≤R0+1. Since γ=O(R0-1) and η2=O((R0-1)2) when (B,S)→(B∗,S∗), there is a constant σ~1 such that 91 W~-IS12v∞2≤σ~1(R0-1)2 when (B, S) approaches (B∗,S∗). Regarding the second addend in (90), it holds that diag(v(t))W~2≤diag(v(t))2W~2=R0maxl=1,…,Nv∞,l. Since v∞2=O(R0-1) when (B,S)→(B∗,S∗), it follows that there is a constant σ~2 such that 92 diag(v(t))W~2S12v∞2≤σ~2(R0-1)2 when (B, S) approaches (B∗,S∗). Regarding the third addend in (90), it holds per definition of the matrix 2-norm that c(t)S-12diag(Bv∞)S-122=c(t)maxl=1,…,N∑j=1Nβjlδlv∞,j≤maxl=1,…,NWv∞l, 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 upper-bound 93 c(t)S-12diag(Bv∞)S-122S12v∞2≤σ~3(R0-1)2 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 v∞TΛ2(t)‖v∞‖22≤δmaxσ~1+σ~2+σ~3(R0-1)2‖v∞‖22ξ(t)2 when (B, S) approaches (B∗,S∗). Since v∞22=O((R0-1)2) when (B,S)→(B∗,S∗), there is a constant σ~4 such that, as (B, S) approaches (B∗,S∗), it holds v∞TΛ2(t)‖v∞‖22≤σ~4ξ(t)2 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 94 dc(t)dt≈c(t)-c2(t)w. To make the statement (94) precise, we define the function cb(t,x), for any scalar x with |x|0 and p>1 when B,S approaches B∗,S∗. Then, the function c(t) is bounded by cb(t,-κ)≤c(t)≤cb(t,κ)∀t≥0, where the scalar κ equals κ=σ2(R0-1)s with s=min{p,2} and some constant σ2>0 as B,S approaches B∗,S∗. Proof With (96), Lemma 9 implies that it holds n(t)≤σ~1(R0-1)pe-σ~2t+σ~3(R0-1)2 for some constants σ~1,σ~2,σ~3>0. Since e-σ~2t≤1, we obtain that n(t)≤κ at every time t, where we define the scalar κ=σ~4(R0-1)s with the constants s=min{p,2} and σ~4=σ~1+σ~3. With n(t)≤κ, we obtain from the differential equation (89) for the function c(t) that 97 c(t)-c2(t)w-κ≤dc(t)dt≤c(t)-c2(t)w+κ∀t≥0. The upper and lower bound (97) give rise to a Riccati differential equation, which can be solved exactly, and we obtain that the function c(t) is bounded by c(t)≥12+121-κwtanhw(w-κ)2t+Υ(-κ) and c(t)≤12+121+κwtanhw(w+κ)2t+Υ(κ). at every time t≥0. □ 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 cb(t,0). Thus, based on the ansatz (15), we approximate the viral state v(t) by vapx(t)=cb(t,0)v∞. With the definition of the function cb(t,x) in (95), it holds that vapx(t)=121+tanhw2+Υ(0)v∞. Then, it follows from the ansatz (15) that the difference of the exact viral state v(t) to the approximation vapx(t) equals 98 ‖v(t)-vapx(t)‖2=c(t)-cb(t,0)‖v∞‖2+‖ξ(t)‖2. 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 cb(t,0) is bounded by 99 c(t)-cb(t,0)≤cb(t,κ)-cb(t,-κ). 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 cb(t,κ) converges to the lower bound cb(t,-κ) then (99) implies that the function c(t) converges to cb(t,0). Furthermore, we must show that the upper bound cb(t,κ) converges to the lower bound cb(t,-κ) uniformly in time t, since the upper bound on the approximation error ‖v(t)-vapx(t)‖2 in Theorem 3 does not depend on time t. From the definition of the function cb(t,x) in (95) we obtain that 100 c(t)-cb(t,0)≤121+κwg(t,κ)-121-κwg(t,-κ), where we denote 101 g(t,κ)=tanhww+κ2t+Υκ. Lemma 10 states that κ=O((R0-1)s) for some s>1 when B,S approaches B∗,S∗. Furthermore, Lemma 8 states that w=O(R0-1). Hence, it holds that κ/w=O((R0-1)s-1) when B,S approaches B∗,S∗. For small x, the series expansion of the square root yields that 121+x=12+14x+O(x2). Thus, for small values of κ/w, we obtain from (100) that c(t)-cb(t,0)≤12g(t,κ)-g(t,-κ)+14wκg(t,κ)+g(t,-κ)+g(t,κ)-g(t,-κ)·Oκ2w2. Since the magnitude of the hyperbolic tangent is bounded by 1, it follows from the definition of the function g(t,κ) in (101) that g(t,κ)-g(t,-κ)≤g(t,κ)+g(t,-κ)≤2, which yields that 102 c(t)-cb(t,0)≤12g(t,κ)-g(t,-κ)+12wκ+O(R0-1)2(s-1), since κ/w=O((R0-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. Lemma 11 Under Assumptions 1 to 3, there is some constant σ1>0 such that g(t,κ)-g(t,-κ)≤2σ1κ at every time t≥0 when the scalar κ approaches zero from above. Proof The mean value theorem gives that g(t,κ)=g(t,0)+∂κg(t,κ)κ=z(t)κ for some z(t)∈(0,κ). Thus, it holds that g(t,κ)-g(t,-κ)=∂κg(t,κ)κ=z1(t)κ+∂κg(t,κ)κ=z2(t)κ for some z1(t)∈(0,κ) and z2(t)∈(-κ,0), which yields that 103 g(t,κ)-g(t,-κ)=∂κg(t,κ)κ=z1(t)κ+∂κg(t,κ)κ=z2(t)κ. To express the derivative of the function g(t,κ), we write the function g(t, x) as g(t,κ)=tanh(h(t,κ)), where we define the function h(t,κ) as h(t,κ)=ww+κ2t+Υκ. Then, the derivative of the function g(t,κ) with respect to the scalar κ is given by ∂κg(t,κ)=4e-h(t,κ)+eh(t,κ)2∂κh(t,κ), which is upper-bounded by 104 ∂κg(t,κ)≤4e-2h(t,κ)∂κh(t,κ). With the derivative of the function h(t,κ), i.e. ∂κh(t,κ)=w4ww+κt+∂κΥκ, we obtain from (104) that ∂κg(t,κ)≤4e-ww+κt-2Υκw4ww+κt+∂κΥκ. The right-hand side of (104) is finite at every time t≥0. Furthermore, for every scalar κ, the right-hand side of (104) converges to zero when t→∞. Hence, we can upper-bound the derivative ∂κg(t,κ) by some constant σ1>0 for every time t. Thus, we obtain from (103) that g(t,κ)-g(t,-κ)=2σ1κ∀t≥0. □ With Lemma 11, we obtain from (102) that there is a constant σ1>0 such that c(t)-cb(t,0)≤σ1κ+12κw+O(R0-1)2(s-1)∀t≥0. Since κ=O((R0-1)s) and w=O(R0-1) when B,S approaches B∗,S∗, we obtain that there exists some constant σ2>0 such that c(t)-cb(t,0)≤σ2(R0-1)s-1. Thus, it follows from (98) that ‖v(t)-vapx(t)‖2≤σ2(R0-1)s-1‖v∞‖2+‖ξ(t)‖2,∀t≥0. Hence, we obtain an upper bound as ‖v(t)-vapx(t)‖2‖v∞‖2≤σ2(R0-1)s-1+‖ξ(t)‖2‖v∞‖2. Then, the upper bound on the error vector ξ(t) in Theorem 2 implies that there are constants σ3,σ4 such that ‖v(t)-vapx(t)‖2‖v∞‖2≤σ2(R0-1)s-1+‖ξ(0)‖2‖v∞‖2e-σ3t+σ4(R0-1)2‖v∞‖2. By assumption it holds that ‖ξ(0)‖2=O(R0-1)p for some constant p>1, and it holds that ‖v∞‖2=O(R0-1) as stated by Theorem 1. Thus, we obtain that ‖v(t)-vapx(t)‖2‖v∞‖2≤σ2(R0-1)s-1+σ5(R0-1)p-1+σ6(R0-1) for some constants σ5,σ6>0, since e-σ3t≤1. By using the definition s=min{p,2} of the scalar s, we complete the proof. Proof of Corollary 1 By assumption, it holds that v(0)=r0v∞, which implies that ξ(0)=0. Thus, we obtain from (22) that ‖v(t)-vapx(t)‖2‖v∞‖2≤σ1(R0-1)∀t≥0 when B,S approaches B∗,S∗. From the definition of the approximation vapx(t) in (20), we obtain that vapx,i(t01)=r1v∞,i for every node i is equivalent to tanhw2t01+Υ(0)=2r1-1 With the definition of the term Υ(0) in (19), it follows that w2t01=artanh(2r1-1)-artanh2v∞Tv(0)‖v∞‖22-1. From v(0)=r0v∞, we obtain that t01=2wartanh(2r1-1)-artanh2r0-1. The inverse hyperbolic tangent equals artanh(x)=12log(1+x)-log(1-x), which completes the proof. Proof of Corollary 3 For NIMFA (4) with homogeneous spreading parameters β,δ, the effective infection rate matrix reduces to W=βδA. Hence, the basic reproduction number reproduction becomes R0=βδρ(A)=ττc, where the last equation follows from the definition of the effective infection rate τ=β/δ and the epidemic threshold τc=1/ρ(A). Furthermore, it holds that δl=δ for every node l and ∑l=1N(x1)l2=1, since the principal eigenvector x1 is of unit length. Thus, the definition of the approximation vapx(t) in (20) yields that vapx(t)=121+tanh(τ-τc)δ2τct+Υ(0)v∞. 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 x1. The idea of proving Theorem 4 is based on Corollary 2: When R0↓1, the viral state dynamics of heterogeneous NIMFA (1) are determined by the four variables x1,w,γ,Υ(0). Thus, we aim to show that the corresponding four variables of the homogeneous NIMFA system (26), which we denote by x1,hom,whom,γhom and Υhom(0), are the same as the variables x1,w,γ,Υ(0) of heterogeneous NIMFA (1). Lemma 12 The homogeneous NIMFA system (26) and heterogeneous NIMFA (1) have the same principal eigenvector x1,hom=x1, the variable γhom=γ and viral slope whom=w. Proof First, we consider the principal eigenvector x1. The effective infection rate matrix of the homogeneous NIMFA system (26) equals 105 Whom=βhomδhomuuT+βhomδhom1minl=1,…,Nx1l∑j=1Nx1jI-βhomδhom∑j=1Nx1jdiag1x11,…,1x1N. We show that the principal eigenvector x1 of heterogeneous NIMFA (1) is also the principal eigenvector x1,hom of the matrix Whom. Indeed, Whomx1=βhomδhom∑j=1Nx1ju+βhomδhom1minl=1,…,Nx1l∑j=1Nx1jx1-βhomδhom∑j=1Nx1ju=βhomδhom1minl=1,…,Nx1l∑j=1Nx1jx1. Thus, x1 is an eigenvector of the effective infection rate matrix Whom of the homogeneous NIMFA system (26). The corresponding eigenvalue equals 106 λ1,hom=βhomδhom1minl=1,…,Nx1l∑j=1Nx1j. The effective infection rate matrix Whom 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 x1 equals the spectral radius ρWhom=λ1,hom and that x1,hom=x1. Second, we consider the variables γ, γhom in Theorem 1. By definition (7) and since x1 is a vector of length 1, it holds that γhom=λ1,hom-11∑l=1Nx1l3=βhomδhom1minl=1,…,Nx1l∑j=1Nx1j-11∑p=1Nx1p3, where the last equality follows from (106). With (28), we obtain further that γhom=1+γ∑l=1Nx1l3-11∑p=1Nx1p3=γ. Thus, the variable γhom of the homogeneous NIMFA (26) equals the variable γ of heterogeneous NIMFA (1). Third, we show that the viral slope whom of the homogeneous NIMFA (26) equals the viral slope w of heterogeneous NIMFA (1). From the definition (87), the variable whom of the homogeneous NIMFA system (26) follows as whom=λ1,hom-1δhom. With (106), we obtain that whom=βhom1minl=1,…,Nx1l∑j=1Nx1j-δhom. Then, the definition of the infection rate βhom in (28) yields that whom=δhom1+γ∑l=1Nx1l3-δhom=δhomγ∑l=1Nx1l3, which simplifies with the definition of δhom in (27) to whom=γ∑l=1Nδlx1l3. Then, the definition of γ in (7) yields that whom=R0-1∑l=1Nδlx1l2. Thus, the viral slope whom of the homogeneous NIMFA system (26) equals the viral slope w of heterogeneous NIMFA (1), which completes the proof. □ In contrast to the variables x1,γ,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 by v~apx(t) and v~apx,hom(t), respectively. The difference of the viral state vectors v(t) and vhom(t) can be written as v(t)-vhom(t)=v~apx(t)-v~apx,hom(t)+v(t)-v~apx(t)-vhom(t)-v~apx,hom(t). With the triangle inequality, we obtain that 107 v(t)-vhom(t)2≤v~apx(t)-v~apx,hom(t)2+v(t)-v~apx(t)2+vhom(t)-v~apx,hom(t)2. Corollary 2 states that there is some constant σ, such that, at every time t≥0, it holds that ‖v(t)-v~apx(t)‖2≤σ‖v∞‖2(R0-1)s-1=O(R0-1)s as R0↓1, since ‖v∞‖2=OR0-1 by Theorem 1. Similarly, Corollary 2 implies that ‖vhom(t)-v~apx, hom(t)‖2=O(R0-1)s as R0↓1. Thus, (107) yields that 108 v(t)-vhom(t)2≤vapx(t)-vapx,hom(t)2+O(R0-1)s. In the following, we bound the first addend on the right side of (108). We insert the expression (23) for the approximations vapx(t) and vapx,hom(t) to obtain that vapx(t)-vapx,hom(t)=1+tanhw2t+Υ(0)γ2x1-1+tanhwhom2t+Υhom(0)γhom2x1,hom=tanhw2t+Υ(0)-tanhw2t+Υhom(0)γ2x1, where the second equality follows from Lemma 12. From Abramowitz and Stegun (1965, 4.5.45), it follows that vapx(t)-vapx,hom(t)=sechw2t+Υ(0)sechw2t+Υhom(0)·sinhΥ(0)-Υhom(0)γ2x1. Since 00, there exists an n0(ϵ)∈N such that both ‖B(n)-B∗‖2<ϵ and ‖S(n)-S∗‖2<ϵ holds for all n≥n0(ϵ). 6 Theorem 1 implies that the steady-state v∞ satisfies ‖v∞‖2=OR0-1 when R0↓1. Thus, also ‖c(t)v∞‖2=OR0-1 at every time t. Thus, a linear convergence of the error term ξ(t) to zero, i.e., ‖ξ(t)‖2=OR0-1, would not be sufficient to show that the viral state v(t) converges to c(t)v∞ when R0↓1. 7 In Prasse and Van Mieghem (2019), an analogous statement has been proved for the discrete-time version of the NIMFA equations (2). 8 The numerical radius is not a matrix norm, since the numerical radius is not sub-multiplicative. Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. Acknowledgements We are grateful to Karel Devriendt for his help in proving Theorem 4. ==== Refs References Abramowitz M Stegun IA Handbook of mathematical functions: with formulas, graphs, and mathematical tables 1965 North Chelmsford Courier Corporation Ahn HJ, Hassibi B (2013) Global dynamics of epidemic spread over complex networks. In: Proceedings of 52nd IEEE conference on decision and control, CDC, IEEE, pp 4579–4585 Anderson RM May RM Infectious diseases of humans: dynamics and control 1992 Oxford Oxford University Press Arfken GB Weber HJ Mathematical methods for physicists Am J Phys 1999 67 165 Bailey NTJ The mathematical theory of infectious diseases and its applications 1975 2 London Charles Griffin & Company Barabási AL Albert R Emergence of scaling in random networks Science 1999 286 5439 509 512 10521342 Devriendt K, Lambiotte R (2020) Non-linear network dynamics with consensus-dissensus bifurcation. arXiv preprint arXiv:2002.08408 Devriendt K Van Mieghem P Unified mean-field framework for susceptible-infected-susceptible epidemics on networks, based on graph partitioning and the isoperimetric inequality Phys Rev E 2017 96 5 052314 29347672 Diekmann O Heesterbeek JAP Metz JA On the definition and the computation of the basic reproduction ratio \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_0$$\end{document} R 0 in models for infectious diseases in heterogeneous populations J Math Biol 1990 28 4 365 382 2117040 Van den Driessche P Watmough J Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission Math Biosci 2002 180 1–2 29 48 12387915 Erdős P Rényi A On the evolution of random graphs Publ Math Inst Hung Acad Sci 1960 5 1 17 60 Fall A Iggidr A Sallet G Tewa JJ Epidemiological models and Lyapunov functions Math Model Nat Phenom 2007 2 1 62 83 He Z Van Mieghem P Prevalence expansion in NIMFA Physica A 2020 540 123220 Heesterbeek JAP A brief history of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_0$$\end{document} R 0 and a recipe for its calculation Acta Biotheor 2002 50 3 189 204 12211331 Hethcote HW The mathematics of infectious diseases SIAM Rev 2000 42 4 599 653 Horn RA Johnson CR Matrix analysis 1990 Cambridge Cambridge University Press Issos JN (1966) The field of values of non-negative irreducible matrices. PhD thesis, Auburn University Kendall DG On the generalized birth-and-death process Ann Math Stat 1948 19 1 1 15 Khanafer A Başar T Gharesifard B Stability of epidemic models over directed graphs: a positive systems approach Automatica 2016 74 126 134 Kitzbichler MG Smith ML Christensen SR Bullmore E Broadband criticality of human brain network synchronization PLoS Comput Biol 2009 5 3 e1000314 19300473 Lajmanovich A Yorke JA A deterministic model for gonorrhea in a nonhomogeneous population Math Biosci 1976 28 3–4 221 236 Li CK Tam BS Wu PY The numerical range of a nonnegative matrix Linear Algebra Appl 2002 350 1–3 1 23 Liu QH Ajelli M Aleta A Merler S Moreno Y Vespignani A Measurability of the epidemic reproduction number in data-driven contact networks Proc Natl Acad Sci 2018 115 50 12680 12685 30463945 Liu F, Cui S, Li X, Buss M (2020) On the stability of the endemic equilibrium of a discrete-time networked epidemic mode. arXiv preprint arXiv:2001.07451 Maroulas J Psarrakos P Tsatsomeros M Perron-Frobenius type results on the numerical range Linear Algebra Appl 2002 348 1–3 49 62 Nowzari C Preciado VM Pappas GJ Analysis and control of epidemics: a survey of spreading processes on complex networks IEEE Control Syst Mag 2016 36 1 26 46 Nykter M Price ND Aldana M Ramsey SA Kauffman SA Hood LE Yli-Harja O Shmulevich I Gene expression dynamics in the macrophage exhibit criticality Proc Natl Acad Sci 2008 105 6 1897 1900 18250330 Paré PE Beck CL Nedić A Epidemic processes over time-varying networks IEEE Trans Control Netw Syst 2017 5 3 1322 1334 Paré PE Liu J Beck CL Kirwan BE Başar T Analysis, estimation, and validation of discrete-time epidemic processes IEEE Trans Control Syst Technol 2018 28 1 79 93 Pastor-Satorras R Castellano C Van Mieghem P Vespignani A Epidemic processes in complex networks Rev Mod Phys 2015 87 3 925 Prasse B, Van Mieghem P (2018) Network reconstruction and prediction of epidemic outbreaks for NIMFA processes. arXiv preprint arXiv:1811.06741 Prasse B Van Mieghem P The viral state dynamics of the discrete-time NIMFA epidemic model IEEE Trans. Netw. Sci. Eng. 2019 7 3 1667 1674 Prasse B, Van Mieghem P (2020) Predicting dynamics on networks hardly depends on the topology. arXiv preprint arXiv:2005.14575 Rami MA, Bokharaie VS, Mason O, Wirth FR (2013) Stability criteria for SIS epidemiological models under switching policies. arXiv preprint arXiv:1306.0135 Van Mieghem P Graph spectra for complex networks 2010 Cambridge Cambridge University Press Van Mieghem P The N-intertwined SIS epidemic network model Computing 2011 93 2–4 147 169 Van Mieghem P The viral conductance of a network Comput Commun 2012 35 12 1494 1506 Van Mieghem P Performance analysis of complex networks and systems 2014 Cambridge Cambridge University Press Van Mieghem P SIS epidemics with time-dependent rates describing ageing of information spread and mutation of pathogens Delft Univ Technol 2014 1 15 1 11 Van Mieghem P (2016) Universality of the SIS prevalence in networks. arXiv preprint arXiv:1612.01386 Van Mieghem P van de Bovenkamp R Accuracy criterion for the mean-field approximation in susceptible-infected-susceptible epidemics on networks Phys Rev E 2015 91 3 032812 Van Mieghem P, Omic J (2014) In-homogeneous virus spread in networks. arXiv preprint arXiv:1306.2588 Van Mieghem P Omic J Kooij R Virus spread in networks IEEE/ACM Trans Netw 2009 17 1 1 14 Verhulst PF Notice sur la loi que la population suit dans son accroissement Corresp Math Phys 1838 10 113 126 Wan Y Roy S Saberi A Designing spatially heterogeneous strategies for control of virus spread IET Syst Biol 2008 2 4 184 201 18681748 Watts DJ Strogatz SH Collective dynamics of “small-world” networks Nature 1998 393 6684 440 9623998