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, 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 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 adjacency matrix by A and its elements by . If there is a link between node i and node j, then , and otherwise. Hence, the virus directly spreads between two nodes i and j only if . We stress that in most applications it holds that , since infected individuals in group i usually do infect susceptible individuals in the same group i. At any time , we denote the viral state of node i by . The viral state 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 and infection rates depend on the nodes i and j.

Definition 1

(Heterogeneous NIMFA) At any time , the NIMFA governing equation isfor every group , where is the curing rate of node i, and is the infection rate from node j to node i.

For a vector , we denote the diagonal matrix with x on its diagonal by . We denote the curing rate matrix . Then, the matrix form of (1) is a vector differential equationwhere is the viral state vector at time t, the infection rate matrix B is composed of the elements , and u is the all-one vector. In this work, we assume that the matrix B is symmetric.

Definition 2

(Steady-State Vector) The steady-state vector 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 toFor the vast majority of epidemiological, demographical, and ecological models, the basic reproduction number is an essential quantity (Hethcote 2000; Heesterbeek 2002). The basic reproduction number 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 was introduced for epidemiological models with only group of individuals. Van den Driessche and Watmough (2002) proposed a definition of the basic reproduction number to epidemic models with groups. For NIMFA (1), the basic reproduction number follows (Van den Driessche and Watmough 2002) as , where 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 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 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 : If , then the only equilibrium of NIMFA (1) is the origin, which is globally asymptotically stable. Else, if , then there is a second equilibrium, the steady-state , whose components are positive, and the steady-state is globally asymptotically stable for every initial viral state . For real-world epidemics, the regime around epidemic threshold criterion is of particular interest. In practice, the basic reproduction number cannot be arbitrarily great, since natural immunities and vaccinations lead to significant curing rates and the frequency and intensity of human contacts constrain the infection rates . 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 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 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 . More precisely, we find the closed-form expression of the viral state for every node i at every time t when , given that the initial state v(0) is small or parallel2 to the steady-state vector .

We introduce the assumptions in Sect. 3. Section 4 gives an explicit expression for the steady-state vector when . In Sect. 5, we derive the closed-form expression for the viral state vector v(t) at any time . The closed-form solution for gives an accurate approximation also for 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 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 , 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 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 depend continuously on time t. Rami et al. (2013) consider a switched model in which both the infection rates and the curing rates 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 is the same4 for every node i, then NIMFA with time-varying, homogeneous spreading parameters has a closed-form solution. In this work, we focus on time-invariant but heterogeneous spreading parameters . We solve NIMFA (1) for arbitrary graphs around the threshold criterion and for an initial viral state v(0) that is small or parallel to the steady-state vector .

Notations and assumptions

The basic reproduction number is determined by the infection rate matrix B and the curing rate matrix S. Thus, the notation is imprecise, since there are infinitely many matrices B, S such that the basic reproduction number equals 1. To be more precise, we consider a sequence of infection rate matrices and curing rate matrices that converges5 to a limit , such that andFor the ease of exposition, we drop the index n and replace and by the notation B and S, respectively. In particular, we emphasise that the assumptions below apply to every element of the sequence. In Sects. 4 to 6, we formally abbreviated the limit process by the notation . For the proofs in the appendices, we use the lengthier but clearer notation . 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, denotes the limit of the curing rate of node i when . The Landau-notation as denotes that for some constant as . For instance, it holds that as .

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 , the curing rates are positive and the infection rates are non-negative, i.e., and for all nodes i, j. Furthermore, in the limit , it holds that and for all nodes i.

We consider Assumption 1 a rather technical assumption, since only non-negative rates and have a physical meaning. Furthermore, if the curing rates 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 .

Assumption 2

For every basic reproduction number , it holds that and for every node . Furthermore, it holds that 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 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 would be small.

Assumption 3

For every basic reproduction number , the infection rate matrix B is symmetric and irreducible. Furthermore, in the limit , 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 effective infection rate matrix W asIn this section, we state an essential property that we apply to solve the NIMFA equations (1) when the basic reproduction number is close to 1: The steady-state vector converges to a scaled version of the principal eigenvector of the effective infection rate matrix W when .

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 which equals the spectral radius . As we show in the beginning of Appendix B, the eigenvalues of the effective infection rate matrix W are real and satisfy . In particular, under Assumptions 1 and 3, the largest eigenvalue , the spectral radius and the basic reproduction number are the same quantity, i.e., .

In Van Mieghem (2012, Lemma 4) it was shown that, for homogeneous NIMFA (4), the steady-state vector converges to a scaled version of the principal eigenvector of the adjacency matrix A when . We generalise the results of Van Mieghem (2012) to heterogeneous NIMFA (1):

Theorem 1

Under Assumptions 1 and 3, the steady-state vector obeyswhere the scalar equalsand the vector satisfies when the basic reproduction number 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 . 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 viral state matrix , for arbitrary observation times , is severely ill-conditioned. Thus, the viral state v(t) at any time approximately equals the linear combination of orthogonal vectors , and we can write , see also Prasse and Van Mieghem (2020). Here, the functions are scalar. We consider the most extreme case by representing the viral state v(t) by a scaled version of only vector , which corresponds to for a scalar function c(t). The viral state v(t) converges to the steady-state vector as . Hence, a natural choice for the vector is , which implies that as . If and , then the approximation is accurate at all times due to two intuitive reasons. Due to the two intuitive steps above, NIMFA (4) reduces around the threshold 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).

If when , then NIMFA (4) is approximated by the linearisation around zero. Hence, it holds that when . The state v(t) of the linear system (8) converges rapidly to a scaled version of the principal eigenvector of the matrix . Furthermore, Theorem 1 states that when . Thus, the viral state v(t) rapidly converges to a scaled version of the steady-state :

Suppose that the viral state v(t) approximately equals to a scaled version of the steady-state vector . (In other words, the viral state v(t) is “almost parallel” to the vector .) 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 around the epidemic threshold, it holds that . Hence, we obtain that Left-multiplying (13) by and dividing by yields that The logistic differential equation (14) has been introduced by Verhulst (1838) as a population growth model and has a closed-form solution.

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 . To be more specific, instead of (9) it holds thatfor some error vector which is orthogonal to the steady-state vector . In Sect. 5.2, we use (15) as an ansatz for solving NIMFA (1).

Second, the steady-state vector is not exactly parallel to the principal eigenvector . More precisely, we must consider the vector in (6). Since , 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 . The ansatz (15) forms the basis for our solution approach. From the orthogonality of the error vector and the steady-state vector , it follows that the function c(t) at time t equalsThe error vector at time t follows from (15) and (16) asOur solution approach is based on two steps. First, we show that6 the error term satisfies at every time t when . Hence, the error term converges to zero uniformly in time t. Second, we find the solution of the scalar function c(t) at the limit .

Assumption 2 implies that7 the viral state v(t) does not overshoot the steady-state :

Lemma 1

Under Assumptions 1 to 3, it holds that for all nodes i at every time . Furthermore, it holds that at every time .

Proof

Appendix C.

Theorem 2 states that the error term converges to zero in the order of when .

Theorem 2

Under Assumptions 1 to 3, there exist constants such that the error term at any time is bounded bywhen the basic reproduction number approaches 1 from above.

Appendix D.

Under Assumption 2, the steady-state is exponentially stable for NIMFA in discrete time (Prasse and Van Mieghem 2019). If the steady-state is exponentially stable, then the error vector goes to zero exponentially fast, since is orthogonal to . Thus, the first addend on the right-hand side in (18) is rather expectable, under the conjecture that the steady-state is exponentially stable also for continuous-time NIMFA (1). Regarding this work, the most important implication of Theorem 2 is that uniformly in time t when , provided the initial value of the error vector is negligibly small.

We define the constant , which depends on the initial viral state v(0), asFurthermore, we define the viral slope w, which determines the speed of convergence to the steady-state , asThen, 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 , when . Furthermore, defineThen, there exists some constant such thatwhere , when the basic reproduction number approaches 1 from above.

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 which is uniform in time t. Furthermore, since both the viral state v(t) and the approximation converge to the steady-state , it holds that when . At time , we obtain from Theorem 3 and (17) thatSince and, by Theorem 1, , we obtain thatHence, for general the approximation error does not converge to zero faster than , and the bound (21) is best possible (up to the constant ) when . With (17), the term in Theorem 2 can be expressed explicitly with respect to the initial viral state v(0) and the steady-state . In particular, it holds that . Furthermore, if the initial viral state v(0) is parallel to the steady-state vector , then it holds that . Thus, if the initial viral state v(0) is small or parallel to the steady-state vector , then it holds that and the bound (21) on the approximation error vector becomesThe time-dependent solution to NIMFA (1) at the epidemic threshold criterion depends solely on the viral slope w, the steady-state vector and the initial viral state v(0). The viral slope w converges to zero as . Thus, Theorem 3 implies that the convergence time to the steady-state goes to infinity when , even though the steady-state 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 for some scalar . Then, for any scalar , the largest time at which the viral state satisfies for every node i converges towhen the basic reproduction number approaches 1 from above.

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 , when . Furthermore, defineThen, there exists some constant such thatwhere , when the basic reproduction number approaches 1 from above.

In contrast to Theorem 3, the approximation error in Corollary 2 does not converge to zero when , since we replaced the steady-state by the first-order approximation of Theorem 1. Corollary 2 implies thatat every time t when , provided that the initial viral state v(0) is small or parallel to the steady-state vector . From (24) it follows that, around the epidemic threshold criterion , 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 and , respectively. Hence, the effective infection rate matrix becomes , and the principal eigenvector of the effective infection rate matrix W equals the principal eigenvector of the adjacency matrix A. Furthermore, the limit process reduces to , with the effective infection rate and the epidemic threshold . 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 for some constant when and defineThen, there exists some constant such thatwhere , when the effective infection rate approaches the epidemic threshold from above.

Appendix G.

From Corollary 3, we can obtain the analogue to Corollary 2 for NIMFA (4) with homogeneous spreading parameters . Furthermore, the approximation defined by (25) equals the exact solution (Van Mieghem 2014b) of homogeneous NIMFA (4) on a regular graph, provided that the initial state is the same for every node i. In particular, the net dose , a crucial quantity in Van Mieghem (2014b); Kendall (1948), is related to the viral slope w via .

Theorem 3 and Corollary 3 suggest that, around the epidemic threshold criterion , 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 by choosing the homogeneous spreading parameters and the adjacency matrix A accordingly?

Theorem 4

Consider heterogeneous NIMFA (1) with given spreading parameters . Suppose that Assumptions 1 to 3 hold and that, for some constant , when the basic reproduction number approaches 1 from above. Define the homogeneous NIMFA systemwhere the homogeneous curing rate equalsthe homogeneous infection rate equalswith the variable defined by (7), and the self-infection rates equalThen, if , there exists some constant such thatwhere , when the basic reproduction number approaches 1 from above.

Appendix H.

In other words, when , for any contact network and any spreading parameters , heterogeneous NIMFA (1) can be reduced to homogeneous NIMFA (4) on a complete graph plus self-infection rates . We emphasise that the sole influence of the topology on the viral spread is given by the self-infection rates . 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 , given by (20), when the basic reproduction number is close, but not equal, to one. We generate an adjacency matrix A according to different random graph models. If , then we set the infection rates to a uniformly distributed random number in [0.4, 0.6] and, if , then we set . We set the initial curing rates to a uniformly distributed random number in [0.4, 0.6]. To set the basic reproduction number , we set the curing rates to a multiple of the initial curing rates , i.e. for every node l and some scalar such that . Thus, we realise the limit process by changing the scalar . Only in Sect. 6.2, we consider homogeneous spreading parameters by setting and 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 . In Prasse and Van Mieghem (2019), we derived an upper bound on the sampling time T which ensures that the discretisation (29) of NIMFA (1) converges to the steady-state . We set the sampling time T to . Except for Sect. 6.3, we set the initial viral state to . We define the convergence time as the smallest time t at whichholds for every node i. Thus, at the convergence time the viral state has practically converged to the steady-state . We evaluate Theorem 3 with respect to the approximation error , which we define asAll 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 nodes and the parameters , . Figure 1 gives an impression of the accuracy of the approximation of Theorem 3 around the epidemic threshold criterion . For a basic reproduction number , 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 nodes, the approximation accuracy of Theorem 3 is depicted. Each of the sub-plots shows the viral state traces of seven different nodes i, including the node i with the greatest steady-state

We aim for a better understanding of the accuracy of the closed-form expression of Theorem 3 when the basic reproduction number converges to one. We generate Barabási–Albert and Erdős–Rényi connected random graphs with nodes. The link probability of the Erdős–Rényi graphs (Erdős and Rényi 1960) is set to . Figure 2 illustrates the convergence of the approximation of Theorem 3 to the exact solution of NIMFA (1). Around the threshold criterion , the approximation error converges linearly to zero with respect to the basic reproduction number , which is in agreement with Theorem 3. The greater the network size N, the greater is the approximation error for Barabási–Albert networks. The greater the network size N, the lower is the approximation error for Erdős–Rényi graphs.

Fig. 2

The approximation error of the NIMFA solution versus the basic reproduction number for different network sizes N

Impact of degree heterogeneity on the approximation accuracy

For NIMFA (4) with homogeneous spreading parameters , the approximation 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 nodes and an average node degree of 4. We vary the link rewiring probability from , which correspond to a regular graph, to , which corresponds to a “completely random” graph. Figure 3 depicts the approximation error versus the rewiring probability 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 versus the link rewiring probability for Watts–Strogatz random graphs with nodes and homogeneous spreading parameters

Impact of general initial viral states on the approximation accuracy

Theorem 3 required that the initial error converges to zero, which means that the initial viral state v(0) must be parallel to the steady-state or, since , converge to zero. To investigate whether the approximation of Theorem 3 is accurate also when the initial error does not converge to zero, we set the initial viral state of every node i to a uniformly distributed random number in for some scalar . By increasing the scalar , the initial viral state v(0) is “more random”. Figure 4 shows that the approximation error is almost unaffected by an initial viral state v(0) that is neither parallel to the steady-state nor small. Figure 5 shows that the viral state v(t) converges rapidly to the approximation as time t increases.

Fig. 4

The approximation error versus the scalar , which controls the variance of the randomly generated initial viral state v(0), for Barabási–Albert networks with nodes

Fig. 5

For a Barabási–Albert random graph with nodes, a basic reproduction number and a randomly generated initial viral state v(0), the approximation accuracy of Theorem 3 is depicted. The viral state traces of seven different nodes i are depicted

For general initial viral states v(0) with , it holds that since the approximation is parallel to the steady-state vector . Hence, the approximation does not converge point-wise to the viral state v(t) when . However, based on the results shown in Figs. 4 and 5, we conjecture convergence with respect to the -norm for general initial viral states v(0) when .

Conjecture 1

Suppose that Assumptions 1 to 3 hold. Then, it holds for the approximation defined by (20) thatwhen the basic reproduction number 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 for a directed infection rate matrix B, which leads us to:

Fig. 6

The approximation error of the NIMFA solution versus the basic reproduction number for directed Erdős–Rényi graphs for different network sizes N

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 when the basic reproduction number approaches 1 from above.

Accuracy of the approximation of the convergence time

Corollary 1 gives the expression of the convergence time from the initial viral state to the viral state for any scalars around the epidemic threshold criterion . We set the scalars to and and define the approximation errorwhere denotes the exact convergence time and denotes the approximate expression of Corollary 1. We generate Barabási–Albert and Erdős–Rényi random graphs with nodes. Figure 7 shows that Corollary 1 gives an accurate approximation of the convergence time when the basic reproduction number is reasonably close to one.

Fig. 7

The approximation error of the convergence time versus the basic reproduction number for different network sizes N

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 must be closer to 1 than for Erdős–Rényi graphs.

Fig. 8

The approximation accuracy of Theorem 4 on a Erdős–Rényi random graph with nodes. Each of the sub-plots shows the viral state traces of seven different nodes i, including the node i with the greatest steady-state

Fig. 9

The approximation accuracy of Theorem 4 on a Barabási–Albert random graph with nodes. Each of the sub-plots shows the viral state traces of seven different nodes i, including the node i with the greatest steady-state

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 , provided that the infection rates are symmetric (). Numerical simulations demonstrate the accuracy of the solution when the basic reproduction number 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 . 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 for every time . On the one hand, since the viral dynamics approximately remain in a one-dimensional subspace of , an accurate network reconstruction is numerically not viable around the epidemic threshold (Prasse and Van Mieghem 2018). Furthermore, when the basic reproduction number is large, then the viral state v(t) rapidly converges to the steady-state , which, again, prevents an accurate network reconstruction. On the other hand, only the principal eigenvector 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 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 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 to approach the steady-state . The viral state v(t) converges to the steady-state exponentially fast. However, as the basic reproduction number approaches one, the convergence time 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 stand on the agenda of future research.

Table 1

Nomenclature

β_ij	Infection rate from node j to node i
B	Infection rate matrix; B_ij=β_ij
c(t)	Projection of the viral state v(t) on the steady-state v_∞; see (16)
δ_i	Curing rate of node i
diag(x)	N×N diagonal matrix with x∈R^N on its diagonal
I	N×N identity matrix
λ_k	k-th eigenvalue of the matrix W; λ_1>λ_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)
R_0	Basic reproduction number; R_0=ρ(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; v_i(t)∈[0,1] for i=1,…,N
v_∞	Steady-state vector, see Definition 2
w	The viral slope; w=(R_0-1)∑l=1Nδ_lx_1l2
W	Effective infection rate matrix W=S^-1B; ρ(W)>1
W~	Symmetric N×N matrix W~=S^-12BS^-12
x_k	k-th eigenvector of the matrix W; Wx_k=λ_kx_k
ξ(t)	Viral state component that is orthogonal to v_∞; ξ(t)=v(t)-c(t)v_∞