key: cord-0002713-3vkxrsbd
authors: Chen, Shanshan; Wang, Kaihua; Sun, Mengfeng; Fu, Xinchu
title: Spread of competing viruses on heterogeneous networks
date: 2017-06-28
journal: Philos Trans A Math Phys Eng Sci
DOI: 10.1098/rsta.2016.0284
sha: b04bbb12e2051520be42a62e23a90b3a80f75d6a
doc_id: 2713
cord_uid: 3vkxrsbd

In this paper, we propose a model where two strains compete with each other at the expense of common susceptible individuals on heterogeneous networks by using pair-wise approximation closed by the probability-generating function (PGF). All of the strains obey the susceptible–infected–recovered (SIR) mechanism. From a special perspective, we first study the dynamical behaviour of an SIR model closed by the PGF, and obtain the basic reproduction number via two methods. Then we build a model to study the spreading dynamics of competing viruses and discuss the conditions for the local stability of equilibria, which is different from the condition obtained by using the heterogeneous mean-field approach. Finally, we perform numerical simulations on Barabási–Albert networks to complement our theoretical research, and show some dynamical properties of the model with competing viruses. This article is part of the themed issue ‘Mathematical methods in medicine: neuroscience, cardiology and pathology’.

Despite centuries of efforts to improve public health and mitigate epidemic disease effects, the threat of infectious diseases remains. An important research method is epidemiology, with mathematical modelling as an analytical approach [1, 2] . Mathematical modelling is an effective tool in understanding the spread of infectious diseases: it provides the best way of translating our knowledge about the status of disease transmission [33] are insufficient and that the clustering coefficient also affects the final state of infection.

2. The SIR model with the probability-generating function and its epidemic threshold PA is the simplest kind of 'moment closure mode'. It is closed on the tuple level, and gives a set of differential equations about tuple changes. In a static network, research on regular and random networks, whose degree of volatility is not too large, has been carried out previously [10] [11] [12] [13] .

In [26, 27] , Volz gave the pair-wise SIR model with the PGF and worked out the basic dynamical behaviour by some approximate methods. Then, Miller and Volz further simplified the Volz model [28] [29] [30] . However, this model did not involve 'clustering' or the heterogeneous degree, which are important characteristics of real contact networks. However, it provided an idea to reduce the dimensionality of the transmission models by using the pair-wise method.

(a) The model with the probability-generating function First, we revisit some results from the Volz model [27] . Volz defined a PGF g(θ) = k p(k)θ k in an undirected random network with heterogeneous connectivity. Nodes can be in any of the three exclusive states: susceptible (S), infectious (I) or recovered (R It is desirable to determine the dynamics of the number of susceptible and infected individuals and to develop equations in terms of these quantities. We note S(t), I(t) as the fraction of nodes susceptible, infected at time t, where S(t) = [S]/N, I = [I]/N. S k (t) is the fraction of nodes susceptible with degree k at time t, where S k = [S k ]/N k . P I , P S is the probability that an arc with a susceptible ego and an infectious, susceptible alter, where P

A degree k susceptible node has an expected number kP I of contacts with infectious nodes, and an expected number βkP I dt of degree k susceptible node contacts will transmit disease to that node in a small time dt. It is clear that S k is the fraction of degree k nodes in the network remaining susceptible at time t, and S k = −βkp I (t)S k . We note S 1 = θ; therefore, the fraction of

An edge in class θ loses its status only when it transmits, i.e. when a transmission occurs along it. The dynamics of θ are given by

Unfortunately, this does not form a closed system of differential equations, as both θ (t) and I (t) depend further on the dynamical variable P I . From the definition of P I , we have

. To obtain the derivative of [SI] and k k[S k ], it requires careful consideration of the rearrangements of arcs from S to S and from S to I, as −S nodes become infected in a small time interval. Firstly, the derivative of k k[S k ] is easily placed in terms of S k and θ.

To calculate [SI] , Volz introduces the notation δ XY (Z) to represent the average excess Z-degree of nodes currently in disease state X selected by following a randomly chosen arc from node Y to X [27] . For example, δ SI (I) represents randomly choosing an edge from I to S, then following that arc to its destination (susceptible) node and finally counting all of the other edges emanating from that node to other infected nodes (ignoring the one along which it arrived). δ SI (S) give the average number of contacts to other susceptible nodes. The calculation of δ XY (Z) is straightforward and based on the current degree distribution of susceptible nodes in Volz [27] .

Accordingly, [SI] includes the increasing or decreasing rate of [SI] at time t. Thence, [SI] is reduced at rate −β[SI], which accounts for all arcs which have an infectious ego which transmits to a susceptible alter; −γ [SI], which is due to the recovery of the I nodes; −S δ SI (I)/g (1), which represents S infected by its other infectious contacts by following that arc from one infectious node to this susceptible node. Similarly, [SI] is increased at rate −S δ SI (S)/g (1), which represents S infected by its other infectious contacts following that arc from one susceptible node to this susceptible node. Similarly, the dynamics for P S can be derived analogously to the equation for P I .

The details of calculations are given in the Volz model [27] and result in the following model:

Assume that a small fraction of the nodes in the network are selected uniformly at random and initially infected, that is, P S (t) ≈ 1, θ ≈ 1, P I (t) ≈ 0, [I(t)] ≈ 0; Volz gave the epidemic threshold in terms of the transmissibility. We can obtained the corresponding basic reproductive number, which denotes the expected number of secondary infections caused by a single infected individual in a completely susceptible population. Accordingly, if R 0 < 1, the disease will eventually disappear, and as long as R 0 > 1 there will be a positive equilibrium representing an endemic status. The basic reproductive number of the Volz model is Then, based on the above work and the results in [25] [26] [27] , House and Keeling studied the full SIR pairwise equations of Eames & Keeling [25] , and gave a 'clustered PGF' model using a relatively small number of ODEs [31] , which is capable of capturing epidemic dynamics on clustered networks with a heterogeneous link distribution. The model consists of a system of six ODEs:

2)

The above system is determined by the following model:

House and Keeling closed the two triples [SSI] and [ISI] that appear in the above unclosed pairwise SIR equations using the standard pairwise closure, the triple closure

, so that the heterogeneity in both the link distribution and clustering can be analysed by using a low-dimensional model with a small number of dynamical variables. Then, the triples can be approximated by tuples and singles [31] ,

Although the dynamical model with respect to 'closure (Φ)' by the PGF method is low dimensional, there are few studies on its dynamical behaviour because of the complicated formula containing some variables, 'Y, θ ', etc., in the denominators.

In this section, we will undertake an analysis of the dynamical behaviour under some assumptions and approximations. We first discuss the stability of equilibria for the model (2.2). Obviously, system (2.2) has zero equilibrium E 0 = (1, 0, 0, nN, 0, 0) ; we will use the Jacobian matrix near the equilibrium to derive the basic reproduction number of system (2.3).

In the early spreading stage, if the population size is large enough, it is reasonable to assume that

Furthermore, [SI]/Y is assumed to be equal to 1 initially, and the evolution of [SI]/Y with time is depicted in figure 1b . Therefore, this shows that the assumption of [SI]/Y ≈ 1 is reasonable in initial time.

Then, under the above basic assumptions, we obtain the approximation of some modules in the Jacobian matrix near E 0 . Finally, we note that g (1)/g (1) = ξ , and the corresponding Jacobian matrix near E 0 is

In order to analyse the eigenvalues of J| E 0 , we give the characteristic equation:

.

Then we only need to calculate the eigenvalues ofJ; its characteristic equation is: By simply applying the 'Sheng Jin formula', which is the improvement of Cardano's formula for solving the cubic equation [44] , equation (2.4) has two conjugate imaginary roots and one real root. Then in order to judge the local stability of the disease-free equilibrium, we need to make sure that all eigenvalues of (2.4) have negative real parts. That is to say, the roots need to satisfy the following conditions:

From the above analysis, we have the following result.

, then the disease-free equilibrium E 0 of the system (2.3) is locally stable; otherwise, it is unstable.

Although the above method is convenient, it cannot derive the specific expression of the basic reproduction number. Let us consider the initial phase of an infection invading a total susceptible population. In line with requiring the infection states to increase in initial time, we can also study the issue according to the method presented in [45, 46] under some conditions. Because [S] is assumed to be equal to N initially,

Theorem 2.2. Define R 0 = (βn/γ )C SI , when γ − 2βΦξ ≥ max{−2β, ((2βφξ − β(ξ 2 − 1))/(ξ − 1))}. If R 0 < 1, the disease-free equilibrium E 0 of system (2.3) is locally asymptotically stable; otherwise, it is unstable.

Proof. R 0 = (βn/γ )C SI relates to the parameters of the network structure average degree n and clustering coefficient Φ. We know that the correlation coefficient of state C SI has a critical relationship with the basic reproduction number. So we consider the correlation coefficient of states C SI = N

We also suppose that θ(0) ≈ 

It is clear that the change rate also depends on the value of [II]/Y. It is of order 1, even when the density of infectious individuals is small [43] Now we discuss the existence and uniqueness of the positive solution of system (2.6). Obviously, the positive equilibrium (x * , y * ) of system (2.6) satisfies the following equations:

and y * = 2β γ − 2βΦξ + β(1 + ξ ) . We get that, if γ − 2βΦξ > (2βφξ − β(ξ 2 − 1))/(ξ − 1), it can ensure that 0 < x * , y * < 1. Then, we consider G = {(x, y) | 0 < x, y < 1}. Therefore, system (2.6) has the unique positive solution (x * , y * ). Now it is necessary to prove the global stability of (x * , y * ). We calculate the Jacobian matrix in the positive equilibrium (x * , y * ) of (2.6),

Clearly,

In that way, all real parts of the eigenvalues of the Jacobian matrix are negative. Hence, the positive equilibrium (x * , y * ) is locally asymptotically stable. Similarly, we find that (0, y * ) is unstable. Now we further discuss the global stability of the positive equilibrium (x * , y * ). For system (2.6) ,

, so the condition γ − 2βΦξ ≥ −2β should be ensured for the sake of the fixed sign of D, and, also, it is not identically equal to zero in any subregion. Therefore, based on the 'Bendixson-Dulac criteria' [47, 48] , when γ − 2βΦξ ≥ −2β and γ − 2βΦξ > ((2βξ − β(ξ 2 − 1))/(ξ − 1)), system (2.6) does not have limit cycles in G. So, the positive equilibrium point is globally stable. So, we obtain R 0 = βnx * /γ .

The process of one epidemic spreading usually prevents many other viruses spreading, or one pathogen sometimes generates many strains with different spreading features [32, 49, 50] . A problem of some interest in multi-strain spreading is the behaviour of competing viruses. They will compete with populations with different infectious rates and final states.

In this section, we will construct an SIAR model depicting two kinds of viruses from the same pathogen or two distinct pathogens competing with each other at the expense of common susceptible individuals by using the pair-wise method and closed by PGF. The two kinds of viruses are named I and A. In our competitive spreading model, the two viruses are exclusive: a node cannot be infected by virus I and virus A simultaneously, and the 'virus' may refer to pathogens, products or two opposite views. Let the spreading rate of I be β and that of strain A be μ; the recovery rates are α 1 and α 2 , respectively. On a network, every node represents an individual, then the edges are the links between the individuals. Each individual can be in one of the three states: susceptible to the disease, infectious when they can spread the disease to the susceptible, and recovered when they have been infected but can no longer spread or catch the disease. However, the equations describe the behaviour of [AB] pairs, instead of just the behaviour of individuals, which is different from [33] . Out of many different types of PAs, we follow in our discussion the deterministic model described by House & Keeling in [31] and some studies in §2. 

According to the discussion in §2, our two-virus model with a competing mechanism using the moment closure PGF can be described and determined as follows: 

The triples will be approximated by tuples and singles, which can capture epidemic dynamics on clustered networks with a heterogeneous link distribution:

YQ .

In the following, we will derive the basic reproduction number for system (3.1) from the eigenvalues of the Jacobian matrix at disease-free equilibrium. 

Now we discuss the eigenvalues of the above matrix in the following two cases.

Case 3.1. Consider the special case first: the clustering coefficient Φ = 0. In this case,

The eigenvalues corresponding to this Jacobian satisfy the following identity:

which implies that, in order to ensure the stability of E 0 , if and only if βξ − β − α 1 < 0 and μξ − μ − α 2 < 0, we note R 01 = β(ξ − 1)/α 1 and R 02 = β(ξ − 1)/α 1 ; that is, R 01 < 1, R 02 < 1. According to the results in [27] , we know that they are the basic reproduction numbers of the SIR and SAR models, respectively, when there is only one virus spreading in the network. Therefore, we derive the basic reproduction number R 0 = max{R 01 , R 02 }. Obviously, the disease-free equilibrium E 0 is stable in the SIAR model if and only if the disease-free equilibria are stable in both the SIR and SAR models. 

where we denote λ(λ + α 1 ), (λ + α 2 ) and λ − [(β + μ)Φξ − (α 1 + α 2 )] as a, b and c, respectively, andĴ Through some similar transformation of the matrixĴ, it is simplified to the following form:

Obviously, according to the discussions in §2, matrices A and B correspond to the characteristic equations of the Jacobian matrices of the SIR and SAR models near their disease-free equilibria E 0 , respectively. Therefore, in the SIAR model, the characteristic equation of the Jacobian matrix near E 0 satisfies the following equation:

Therefore, the disease-free equilibrium E 0 is locally stable if and only if all of the eigenvalues are negative, equivalently, R 01 < 1,

where R 01 and R 02 are the basic reproduction numbers of the SIR and SAR models, respectively, which are the expected numbers of secondary infections caused by a single-infected individual in a completely susceptible population if there is only one virus spreading in the network. The expression of (β + μ)Φξ /(α 1 + α 2 ) shows that R 0 of the SIAR model also depends on the clustering coefficient Φ during the period of the transmission of virus A and I. Then, in the SIAR model, as long as we require the basic reproduction number

we can obtain the following result.

This shows clearly that, in the two-strain competing model, only controlling their own spreading thresholds is not enough. The value of (β + μ)Φξ /(α 1 + α 2 ) also has a vital influence on the state of infection, that is to say, the spreading threshold of the two-strain competing model not only depends on the degree distribution, but also is related to the clustering coefficient of the network.

Because of the complicated nature of the SIAR model, it is difficult to analyse its global dynamics by a rigorous theoretical method. So in this section, we will provide further analysis of the transmission dynamics by numerical simulations.

We perform the simulations on networks not only to verify our theoretical results obtained in the previous sections, but also to reveal some new phenomena which are difficult to obtain by 3) in §2b, R 0 = βn/γ x * , which depends on the clustering coefficient Φ. Figure 1a shows the effect of the clustering coefficient Φ on R 0 . It is clear that, with increasing Φ, R 0 is decreasing under the condition that the degree distribution remains unchanged, which is the same as the conclusion in [34] for the model with random and regular contacts by the pair-wise method.

As shown in figure 2, we consider how the value of ((β + μ)/(α 1 + α 2 )) influences the epidemic spreading of SIAR model in §3 with four sets of parameters: It is clear that the final state is critical to the relationship between β and μ. As long as β < μ (β > μ) and (β + μ)/(α 1 + α 2 ) > 1,

) occupying the network, which corresponds to the solution (0, A * ) ((I * , 0)) to system (2.8). However, if (β + μ)/(α 1 + α 2 ) < 1, the condition is the same, but the infection scale will decrease. Then, paying attention to the second and third sets of parameters, when β = μ, no matter whether (β + μ)/(α 1 + α 2 ) > 1 or < 1, the two competing viruses will coexist. However, it is noteworthy that when (β + μ)/(α 1 + α 2 ) < 1 the infection scale is lower than the case of (β + μ)/(α 1 + α 2 ) > 1. Figure 3 depicts the relationship between the ultimate sizes of R 1 and R 2 for the competing diseases I and A under different cases in logarithmic coordinates. We note that, according to §2, R 01 = βnx * 1 /α 1 , R 02 = μnx * 2 /α 2 , and R 03 = (β + μ)Φξ /(α 1 + α 2 ). The lines C 1 , C 2 show that two competing viruses will lead to large-scale outbreaks under the conditions R 01 > 1, R 02 > 1, R 03 < 1, R 01 > 1, R 02 > 1, R 03 > 1. From C 3 and C 4 , we can see that only I(t) has large-scale outbreaks on the network when R 01 > 1, R 02 < 1. While we see similar phenomena in figure 2 , the outbreak scale is lower when R 03 < 1 than when R 03 > 1. Therefore, from figures 2 and 3, we see that, as long as there is a value R 01 or R 02 greater than 1, disease will be prevalent on the network; whereas if R 03 > 1, disease will have a massive outbreak on the network.

When we consider the spread of an epidemic, it is the contact structure between individuals that determines the progress of the disease through the population. Correlation models, and, in particular, pair-wise models, have been primarily used to describe the behaviour of simple spatial models. In addition, generally speaking, the forms of interaction between two epidemic Where under the following cases: the red line with stars C 1 satisfies R 01 > 1, R 02 > 1, R 03 < 1; the green line with squares C 2 satisfies R 01 > 1, R 02 > 1, R 03 > 1; the blue line with circles C 3 satisfies R 01 > 1, R 02 < 1, R 03 < 1; and the black line with triangles C 4 satisfies R 01 > 1, R 02 < 1, R 03 > 1. (Online version in colour.) particles in a multi-strain epidemic model may contain many types, but the competing strains are relatively simple. However, because of the complicated nature of the approximated forms by pair-wise modelling, the further study of epidemic dynamics with this model is limited using the present methods. Now we give a brief summary of this paper. We have derived the basic reproduction number of the basic SIR model closed by PGF. Then, we proposed an SIAR model on heterogeneous networks with pair-wise modelling closed by PGF. By theoretical analysis and numerical simulation, we described the effect of a pair-wise model on the spreading dynamics of two competing viruses, and found the conditions to ensure the local stability of a disease-free equilibrium. The stable condition not only is related to the basic reproduction number of SIR and SAR models, but also connects with (β + μ)Φξ /(α 1 + α 2 ).

In this paper, we have only studied the simplest case of competing pathogens. The concept of competition between the two strains of infection (or pathogens) means that two strains from one pathogen cannot co-infect in a single host at the same time. There are also other forms of interaction between two pathogens except competing virus, such as super-infection, which means that one more virulent pathogen can outcompete the other less virulent pathogen; or co-infection, which means that two pathogens can be hosted in one individual. These processed may also be modelled by using pair-wise modelling closed by PGF, which will be further analysed in our future research.

The prevention of malaria

Infectious diseases of humans

A contribution to the mathematical theory of epidemics

Some evolutionary stochastic processes

A short history of mathematical population dynamics

Dynamical modelling and analysis of epidemics in networks

Velocity and hierarchical spread of epidemic outbreaks in scale-free networks

Thresholds for epidemic spreading in networks

Epidemic spreading on weighted complex networks

Susceptible-infected-removed epidemic models with dynamic partnerships

Correlation models for childhood epidemics

Representing spatial interactions in simple ecological models

Networks and the epidemiology of infectious disease

Pari-edge approximation for heterogeneous lattice population models

Pairlevel approximations to the spatio-temporal dynamics of epidemics on asymmetric contact networks

Reproduction numbers for epidemics on networks using pair approximation

Modelling disease spread through random and regular contacts in clustered populations

A motif-based approach to network epidemics

A closure approximation technique for epidemic models

Number of sexual encounters involving intercourse and the transmission of sexually transmitted infections

Contact tracing and disease control

The impact of contact tracing in clustered populations

Epidemic dynamics on adaptive networks

Untangling the interplay between epidemic spread and transmission network dynamics

Modeling dynamic and network heterogeneities in the spread of sexually transmitted diseases

Susceptible-infected-recovered epidemics in dynamic contact networks

SIR dynamics in random networks with heterogeneous connectivity

A note on a paper by Erik Volz: SIR dynamics in random networks

Incorporating disease and population structure into models of SIR disease in contact networks

Edge-based compartmental modelling for infectious disease spread

Insights from unifying modern approximations to infections on networks

Threshold effects for two pathogens spreading on a network

Epidemic thresholds in a heterogenous population with competing strains

Modelling disease spread through random and regular contacts in clustered populations

Modelling the spread of diseases in clustered networks

A model of multi-information dissemination with suppressed action

Dynamics of competing ideas in complex social systems

Marketing new products: bass models on random graphs

Competitive epidemic spreading over arbitrary multilayer networks

Clustering in complex networks. II. Percolation properties

Effects of heterogeneous and clustered contact patterns on infectious disease dynamics

Properties of highly clustered networks

The spread of epidemic disease on networks

A new extracting formula and a new distinguishing means on the one variable cubic equation

The spread of infectious diseases in spatially structured populations: an invasory pair approximation

The effects of local spatial structure on epidemiological invasions

Journal of dynamics and differential equations

Nonlinear oscillations, dynamical systems, and bifurcations of vector fields

Superinfection behaviors on scale-free networks with competing strains

Epidemic dynamics of two species of interacting particles on scale-free networks

Authors' contributions. S.C. conducted the literature survey, undertook the study and drafted the manuscript. K.W. discussed the content of the draft manuscript and performed the numerical simulations. M.S. discussed the content of the draft manuscript. X.F. designed the study, discussed the content and finalized the manuscript. All authors read and approved the manuscript.