An SIR pairwise epidemic model with infection age and demography

ABSTRACT The demography and infection age play an important role in the spread of slowly progressive diseases. To investigate their effects on the disease spreading, we propose a pairwise epidemic model with infection age and demography on dynamic networks. The basic reproduction number of this model is derived. It is proved that there is a disease-free equilibrium which is globally asymptotically stable if the basic reproduction number is less that unity. Besides, sensitivity analysis is performed and shows that increasing the variance in recovery time and decreasing the variance in infection time can effectively control the diseases. The complex interaction between the death rate and equilibrium prevalence suggests that it is imperative to correctly estimate the parameters of demography in order to assess the disease transmission dynamics accurately. Moreover, numerical simulations show that the endemic equilibrium is globally asymptotically stable.


Introduction
Mathematical modelling plays an important role in the understanding of infectious disease spreading [1,19]. It provides theoretical support for effective control measures for infectious diseases [12,13]. Since the pioneering work of Kermack and McKendrick in [10], compartmental models have been widely used as one of the main methods to study the transmission mechanism of infectious diseases. One fundamental assumption of these models is that the population is homogeneous mixing, i.e. all individuals have the same probability to contact with others. However, the connectivity pattern of the population is usually heterogeneous in real society, which cannot be described by classical compartmental models. In order to capture the heterogeneous contacts among individuals in epidemic modelling, the idea of complex networks is introduced. In the network, individuals are represented as nodes and contacts between individuals are represented as links between nodes. The contact rate of a node is then proportional to the number of links it has. This representation of contact patterns can offer us a more reliable framework for uncovering the impact of heterogeneous contacts on epidemic dynamics. With complex networks, epidemic models provide accurate predictions on the future epidemic outbreaks. Therefore, network-based models have been studied extensively over the last few decades (to name a few, see [6,14,15,22,23,28]).
The demography of the population is one of the most important factors and needs to be included in mathematical modelling. However, most network-based models neglect the effect of demography and assume that the underlying network is static. This is valid for fast diseases, such as influenza and childhood diseases, that develop on a much shorter timescale than the lifespan of individuals. Whereas for slowly progressive diseases, such as tuberculosis, HIV/HIDS, and hepatitis C, the diseases last for a long time, sometimes even for the lifespan of the hosts. In this situation, the birth and death of individuals change the number and links of nodes in networks, leading to a dynamic network which in turn affects the epidemic spreading. It is therefore crucial to incorporate the effects of demography when modelling epidemic spreading on networks. Kamp [9] established the SID and SI 1 I 2 D models with HIV as a case study to investigate the effects of network dynamics on epidemic spreading. Jin et al. [8] considered an SIS model with demographics on a time-varying network and got the global dynamics of the model. Piccardi et al. [24] proposed an SIRS model to explore the effects of vital dynamics and age of individuals on epidemics spreading on networks. Luo et al. [16] established an SIS pairwise model with demographics and found that demographics induce the extinction of epidemics.
In addition to the demography of the population, infection age (i.e. the time since the moment at which the infectious individual has been infected) also has a significant impact on the epidemic spreading for slowly progressive diseases. Currently, most network-based epidemic models neglect the infection age by simply assuming that the transmission and recovery rates are constant. This implies that the epidemic process exhibits Markovian behaviour, that is, the time to infect a neighbour and the recovery time are exponentially distributed. However, evidence shows that the time to infect a neighbour and the recovery time for real diseases are usually not exponentially distributed (see, for example, [2-4, 19, 21]). The non-exponential distribution suggests that the epidemic process is non-Markovian, and constant transmission and recovery rates are not satisfactory in describing the real-life diseases. It is therefore natural and essential to incorporate the infection age into epidemic models in which the transmission and recovery rates of infectious individuals depend on the infection age. Van Mieghem et al. [31] investigated the effect of Weibullean infection times with the same mean but different power exponents on the epidemic threshold. Kiss et al. [11] proposed a system of delay differential equations to study the disease transmission with fixed infectious periods. Yang et al. [34] studied an SIS heterogeneous mean-field model with infection age on complex network and discussed the control of diseases. Röst et al. [26] proposed a generalized pairwise model on static networks for non-Markovian epidemic processes with exponential infection time distribution and arbitrary recovery time distributions. Yang et al. [33] established a susceptible-infectious-recovered (SIR) epidemic model with demography and infection age on complex networks and got a threshold of the model. Sherborne et al. [27] extended the edge-based compartmental model [20] from Markovian to non-Markovian epidemic dynamics where the transmission and recovery rates are driven by general independent distributions. In conclusion, much attention has been paid to the effects of demography and infection age in modelling epidemic spreading on networks. However, to the best of our knowledge, there are only a few significant results considering both the demography and infection age in pairwise epidemic models on complex networks.
Based on the above discussion, in this paper, we establish an SIR pairwise model to investigate the effects of demography and infection age on epidemic spreading on complex networks. The rest of this paper is organized as follows. Section 2 introduces the pairwise model with demography and infection age. In Section 3, we obtain the basic reproduction number of the model and analyse the stability of the disease-free equilibrium. In Section 4, numerical simulations are performed to investigate the effects of infection age and demography on disease spreading. Finally, conclusions are given in Section 5.

The model formulation
We consider an SIR epidemic spreading on a dynamic network with N(t) nodes at time t, where nodes represent individuals in the population and edges are contacts among them. During the epidemic spreading process, each node may have only one of the three possible states: susceptible (S), infectious (I), and recovered (R). Let S k (t) be the number of susceptible nodes of degree k at time t, I k (t, a) the density of infectious nodes of degree k and with infection age a at time t, R k (t) the number of recovered nodes of degree k at time t. Let S(t), I(t, a), and R(t) denote the number of susceptible nodes at time t, the density of infectious nodes with infection age a at time t, and the number of recovered nodes at time t in the network, respectively. The fact that I(t, a) gives 'density' rather than 'number' means that the number of infectious nodes with infection age in the interval (a, a + a) at time t is approximately I(t, a) a, where a is a small increment. Thus, I(t) = ∞ 0 I(t, a) da gives the total number of infectious nodes at time t. It follows that Besides To build the SIR pairwise model with demography and infection age, we make the following hypotheses: (1) Assume that the population grows at a constant rate independent of population size, and all individuals are born susceptible. A newborn individual randomly attaches to k existing individuals, where k is extracted from the newcomer degree distribution θ k and ς = k kθ k is the average degree of newborns. (2) Individuals die at a per capita death rate η. We assume that the disease is not fatal, thus the disease-induced death rate can be neglected. (4) An infectious individual with infection age a recovers with rate γ (a). The recovered individuals are assumed to have permanent immunity. (5) It is naturally assumed that λ(a) and γ (a) are nonnegative and bounded over the interval [0, ∞).
Then the SIR pairwise model on networks with demography and infection age takes the formṠ with the boundary conditions and the initial conditions (S k (0), I k (0, a), R k (0), [SI(0, a)], [SS](0)). Equations (1a)-(1c) describe the dynamics of nodes and Equations (4d)-(1e) describe the dynamics of links. On the right-hand side of Equation (1a), the first (or second) term describes the recruitment (or loss) of S k due to natural birth (or death). The third term represents the loss of S k due to infection. The fourth term considers the fact that a node in S k (or S k−1 ) flows into S k+1 (or S k ) when it receives a link emanating from the newborns. The last term corresponds to that a node in S k (or S k+1 ) flows into S k−1 (or S k ) due to the natural death of its neighbours.  (1). To break the dependency on higher order moments and obtain a closed model, the triple closure approximation in [25] [ is used. Focusing on the total number of susceptible nodes and the density of infectious nodes with infection age a, and replacing the dynamics of R(t) by the evolution of the population size N(t), model (1) is further reduced tȯ with the boundary conditions and the initial conditions (6) where N 0 , S 0 , [SS] 0 ∈ R + = [0, ∞), and ϕ(a), ψ(a) ∈ L 1 + , and the initial conditions are determined by the initial network. Here L 1 + : [0, ∞) → R + is the space of functions that are nonnegative and Lebesgue integrable over the specified interval. All the parameters are positive. Moreover, we assume that model (4) satisfies the compatibility conditions Then model (4) is well posed [7,32]. In the following sections, we mainly focus on the low-dimensional model (4). Preliminary results. We first consider the existence of equilibria to model (4). For simplicity of notation, define π : Model (4) can be rewritten as the following Volterra-type equations: where β(t) is the unique solution of the following nonlinear Volterra equation: This is derived by using We then apply the approach introduced by Thieme [30] to reformulate model (4) with the boundary conditions (5) and initial conditions (6) as a semilinear Cauchy problem.
Rearrange the variables of model (4)

as (S(t), [SS](t), N(t), I(t, a), [SI(t, a)]).
In order to take care of the boundary conditions, we enlarge the state space and consider endowed the usual product norm, and the set We consider the linear operator A : where W 1,1 is a Sobolev space, and we define F : Then by defining we can reformulate the partial differential equation problem (4) as the following abstract Cauchy problem: By using the results in [17,18,30], we derive the existence and the uniqueness of the semiflow , it can be proved that the semiflow coincides with the one generated by using the Volterra integral formulation. Moreover, by the positivity of variables, we deduce the following differential inequalities for model (4): Thus, if are satisfied for some t = t 0 , they are satisfied for all t t 0 . Hence, model (4) leaves the set This means that {U(t)} t 0 is bounded dissipative on X 0+ (see [5]).

Equilibria and the basic reproduction number
In this section, we will derive the basic reproduction number R 0 by analysing the equilibria. To find the equilibria, we look for time-independent solutions (S,Ī(a), [SI(a)], [SS],N) that satisfy model (4) with the time derivatives equal to zero. Then the equilibria satisfy the following equations: It is easy to see that 0 Ī (0) , and E 0 = ( /η, 0, 0, ( /η)ς, /η) is one solution of Equation (11). This solution produces the disease-free equilibrium, and it always exists. Next, an endemic equilibrium will be given by Since we are looking for a nonzero solution, I * (0) = 0. Cancelling I * (0) from both sides of Equation (14) and denoting by F(I * (0)) the left-hand side, we have and F (I * (0)) < 0 for I * (0) ∈ (0, ).
Thus, Equation (14) has a unique positive solution in (0, ) if and only if Given the solution I * (0), according to Equation (12) In epidemiology, R 0 is referred to as the basic reproduction number of model (4). Furthermore, for model (4), we have the following results.

Theorem 3.1: The disease-free equilibrium E 0 of model (4) is locally asymptotically stable if
Proof: First, we introduce the change of variables as follows: After dropping the nonlinear quadratic terms, we obtain the linearized system at the disease-free equilibrium E 0 as follows: [ss](t) = −2ς  (18) and cancelling e zt , we obtain the following equations: Plugging Equation (20)  Thus, Equation (21) has a unique real root z * . Noting that P(0) = R 0 , we have z * > 0 if R 0 > 1, and z * < 0 if R 0 < 1. Set z = x+yi to be an arbitrary complex root to Equation (21). Then 1 = |P(z)| = |P(x + yi)| |P(x)|, which implies that z * x. Thus, all roots of Equation (21) have negative part if and only if R 0 < 1. This implies that the disease-free equilibrium E 0 is locally asymptotically stable if R 0 < 1 and unstable if R 0 > 1. This completes the proof of Theorem 3.1.
Theorem 3.1 provides the local stability of the disease-free equilibrium E 0 of model (4). Furthermore, we can also prove that the equilibrium E 0 is globally asymptotically stable.

Theorem 3.2:
The disease-free equilibrium E 0 of model (4) is globally asymptotically stable if R 0 < 1. SI(t, a)]da. Plugging Equation (8c) into the expression of β(t) gives Let x(t) = [SS](t)/S(t). By Equations (4a) and (4d), we havė Here, we have used the inequality lim t→∞ S(t) /η. Combining Equation (23) with lim t→∞ N(t) = /η, we arrive at lim t→∞ x(t) ς . Using this result and Fatou's lemma, and taking the limit supremum when t → ∞ on both sides of Equation (22), we have As it is assumed that R 0 < 1, the only condition that inequality (24)  Accordingly, the equilibrium E 0 is globally asymptotically stable if R 0 < 1. This completes the proof of Theorem 3.2.

Numerical simulations
In this section, we will perform sensitivity analysis to investigate the effects of demography and infection age on disease spreading and propose effective methods to control diseases. Here, we assume that [SS] 0 ≈ (ς/N 0 )S 2 0 and ψ(a) ≈ (ς/N 0 )S 0 ϕ(a), which mean that the initial network is homogeneous and newborns have the same average degree with the initial network.

Markovian transmission and non-Markovian recovery
We consider a Markovian transmission process with rate λ 0 (i.e. λ(a) ≡ λ 0 ) and a general recovery process. In the simulations, the recovery rate is determined by the recovery time distribution f R (t). Specifically, the instantaneous recovery rate at infection age a is given by the hazard function of f R (a), i.e.
where F R (a) = a 0 f R (t) dt is the cumulative distribution function of the recovery time [26,27]. Particularly, if the recovery time is exponentially distributed with rate γ 0 , it corresponds to a Markovian recovery process. In this case, γ (a) ≡ γ 0 , and the basic reproduction number R 0,1 of model (4) is Next, we study three other recovery time distributions: Weibull distribution, Gamma distribution, and uniform distribution [29,35].
the recovery rate is When α 1 = 1, the Weibull distribution reduces to an exponential distribution, corresponding to a Markovian recovery process. (2) If the recovery time is gamma distributed with shape α 2 and scale β 2 (denoted as G(α 2 , β 2 )), i.e.
When α 2 = 1, the Gamma distribution reduces to an exponential distribution, corresponding to a Markovian recovery process.
In order to compare these distributions, we fix the average recovery time to 1/γ 0 . Then Figure 1 shows the density of infectious nodes (I(t, a)) as a function of time and infection age from model (4) with different recovery time distributions. From Figure 1, we see that the recovery time distributions have a significant impact on epidemic spreading in initial transmission stage. Moreover, Figure 2 demonstrates that non-exponential recovery time distributions can lead to the change of multiple peaks of infection, and a larger shape value of α 1 (or α 2 , α 3 ) results in a larger value of the first peak. Besides, Figure 3 shows that    where V X denotes the variance of distribution X. Figure 3 illustrates that higher variance in recovery time leads to larger transmission threshold and smaller equilibrium prevalence for Weibull, Gamma, or Uniform distributions with the same average recovery time. It indicates that the disease can be effectively controlled by increasing the variance in recovery time.
In Figure 4, we compare the epidemic curves with different types of recovery time distributions. For the sake of comparison, we select the Weibull distribution and Gamma distribution. The parameters of these two distributions are set to have the same mean and variance. It can be seen from Figure 4(a) that the Gamma recovery distribution leads to lower equilibrium prevalence than the Weibull distribution. However, the opposite is true in Figure 4(b). Thus, it cannot identify which distribution is better for disease control. Besides, Figure 4 illustrates that the mean and variance of the recovery time distribution alone are not sufficient to predict the epidemic dynamics. Therefore, for a real disease, it is crucial to estimate the recovery time distribution as accurately as possible in order to get a better prediction of the epidemic dynamics.
Finally, Figure 5 suggests that the endemic equilibrium is globally asymptotically stable.

Non-Markovian transmission and Markovian recovery
We consider a Markovian recovery process with rate γ 0 (i.e. γ (a) ≡ γ 0 ) and a general transmission process. In the simulations, the transmission rate is determined by the infection time (the time that an infectious node takes to infect a susceptible neighbour). Specifically, we assume that the infection time is distributed as f I (t), the transmission rate λ(a) is given by the hazard function of f I (t), i.e.
where F I (a) = a 0 f I (t) dt is the cumulative distribution function of the infection time [11,27]. Here, we study three different infection time distributions: Weibull distribution W(α 4 , β 4 ), Gamma distribution G(α 5 , β 5 ), and Uniform distribution U(α 6 , β 6 ). In order to compare these distributions, we fix the average infection time to 1/λ 0 . Then it has   In contrast with the non-Markovian recovery processes, when the infection time follows Weibull, Gamma, or Uniform distributions with the same average infection time, a higher variance in infection time leads to a larger equilibrium prevalence. It indicates that decreasing the variance in infection time can efficiently control the diseases. Figure 7 illustrates the effects of demography on the spread of an epidemic. The values of η in Figure 7(a) are smaller than those in Figure 7(b). We can see that smaller death rate η results in multiple peaks of infection in Figure 7(a), while this phenomenon disappears for larger η in Figure 7(b). Moreover, a larger value of η leads to a larger equilibrium prevalence in Figure 7(a), whereas the opposite result is true in Figure 7(b). This implies that the death rates have a significant impact on the epidemic dynamics, and there is a complex interaction between the death rate and equilibrium prevalence, enhancing the difficulty of disease control. Thus, it is imperative to correctly estimate the parameters of demography in order to assess the disease transmission dynamics accurately.

Conclusion and discussion
Demography and infection age are important factors in the spread of slowly progressive diseases. In order to investigate their impacts on the spread of infectious diseases, in this paper, we established a pairwise model on dynamic networks. In this model, the transmission and recovery rates depend on the infection age, and the demography of the population is incorporated. The basic reproduction number of the model was derived. It was proved that there is always a disease-free equilibrium and there is also a unique endemic equilibrium when R 0 > 1. Moreover, the disease-free equilibrium is globally asymptotically stable if R 0 < 1. Furthermore, we performed sensitivity analysis to seek for effective measures to control the diseases. We found that increasing the variance in recovery time and decreasing the variance in infection time are beneficial to the control of diseases. However, interaction between the demography and epidemic dynamics is complex, and it enhances the difficulty of disease control. This shows that it is important to correctly estimate the parameters of demography in order to assess the disease transmission dynamics accurately. We note that model (4) only includes the average degree of newborns, but ignores the degree heterogeneity of newcomers. And it is assumed that newborn nodes link to existing nodes randomly in model (4). However, in some situations, newborns may choose contacts adaptively. For example, they may prefer to link to nodes with large degree and avoid to contact with infectious nodes. Thus, it is important to seek for more accurate network model to describe the epidemic transmission processes under different network evolving mechanisms. We leave this for our future work.