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.

KEYWORDS:

• Demography
• infection age
• complex network
• non-Markovian transmission and recovery

2010 MSC SUBJECT CLASSIFICATIONS:

• 92D30
• 37N25

1. 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.

2. The model formulation

We consider an SIR epidemic spreading on a dynamic network with $N\left(t\right)$ 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\left(t\right)$ 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\left(t\right)$ the number of recovered nodes of degree k at time t. Let $S\left(t\right)$, $I(t,a)$, and $R\left(t\right)$ 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+\mathrm{\Delta }a)$ at time t is approximately $I(t,a)\mathrm{\Delta }a$, where $\mathrm{\Delta }a$ is a small increment. Thus, $I\left(t\right)={\int }_0^{\mathrm{\infty }}I(t,a)\mathrm{d}a$ gives the total number of infectious nodes at time t. It follows that $S\left(t\right)=\sum _k{S}_k\left(t\right),I(t,a)=\sum _k{I}_k(t,a),R\left(t\right)=\sum _k{R}_k\left(t\right),N\left(t\right)=S\left(t\right)+I\left(t\right)+R\left(t\right).$ Besides, let $\left[SI\right(t,a\left)\right]$ represent the density of edges with one node being susceptible and the other one being infectious with infection age a at time t. For undirected networks, $\left[SI\right(t,a\left)\right]=\left[I\right(t,a\left)S\right]$. $\left[SS\right]\left(t\right)$ doubles the number of edges with two edge nodes being susceptible at time t. $\left[SSI\right(t,a\left)\right]$ is the density of triplets with the middle node being susceptible and the other two edge nodes being susceptible and infectious with infection age a at time t. $\left[I\right(t,b\left)SI\right(t,a\left)\right]$ denotes the number of triples with the middle node being susceptible and the other two edge nodes being infectious with infection ages a and b, respectively, at time t.

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 ${\theta }_k$ and $\varsigma =\sum _{k}k{\theta }_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.

3. An infectious individual with infection age a infects a susceptible neighbour with rate $\lambda \left(a\right)$.

4. An infectious individual with infection age a recovers with rate $\gamma \left(a\right)$. The recovered individuals are assumed to have permanent immunity.

5. It is naturally assumed that $\lambda \left(a\right)$ and $\gamma \left(a\right)$ are nonnegative and bounded over the interval $[0,\mathrm{\infty })$.

Then the SIR pairwise model on networks with demography and infection age takes the form (1a) $\begin{array}{rl}{\dot{S}}_k\left(t\right)& =\mathrm{\Lambda }{\theta }_k-\eta S_k\left(t\right)-\frac{k{S}_k\left(t\right)}{\sum _{k}k{S}_k\left(t\right)}{\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[SI\right(t,a\left)\right]\mathrm{d}a\\ & +\mathrm{\Lambda }\varsigma \left(\frac{S_{k-1}\left(t\right)}{N\left(t\right)}-\frac{S_k\left(t\right)}{N\left(t\right)}\right)-\eta \left(k{S}_k\right(t)-(k+1\left)S_{k+1}\right(t\left)\right),\end{array}$(1a) (1b) $\begin{array}{rl}\left(\frac{\mathrm{\partial }}{\mathrm{\partial }t}+\frac{\mathrm{\partial }}{\mathrm{\partial }a}\right)I_k(t,a)& =-\gamma \left(a\right)I_k(t,a)+\mathrm{\Lambda }\varsigma \left(\frac{I_{k-1}(t,a)}{N\left(t\right)}-\frac{I_k(t,a)}{N\left(t\right)}\right)\\ & -\eta I_k(t,a)-\eta \left(k{I}_k\right(t,a)-(k+1\left)I_{k+1}\right(t,a\left)\right),\end{array}$(1b) (1c) $\begin{array}{rl}{\dot{R}}_k\left(t\right)& ={\int }_0^{\mathrm{\infty }}\gamma \left(a\right)I_k(t,a)\mathrm{d}a+\mathrm{\Lambda }\varsigma \left(\frac{R_{k-1}\left(t\right)}{N\left(t\right)}-\frac{R_k\left(t\right)}{N\left(t\right)}\right)\\ & -\eta R_k\left(t\right)-\eta \left(k{R}_k\right(t)-(k+1\left)R_{k+1}\right(t\left)\right),\end{array}$(1c) (1d) $\begin{array}{rl}\left(\frac{\mathrm{\partial }}{\mathrm{\partial }t}+\frac{\mathrm{\partial }}{\mathrm{\partial }a}\right)\left[SI\right(t,a\left)\right]& =-\left(\lambda \right(a)+\gamma (a)+2\eta )\left[SI\right(t,a\left)\right]\\ & -{\int }_0^{\mathrm{\infty }}\lambda \left(b\right)\left[I\right(t,b\left)SI\right(t,a\left)\right]\mathrm{d}b+\mathrm{\Lambda }\varsigma \frac{I(t,a)}{N\left(t\right)},\end{array}$(1d) (1e) $\begin{array}{rl}\dot{\left[SS\right]}\left(t\right)& =-2{\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[SSI\right(t,a\left)\right]\mathrm{d}a-2\eta \left[SS\right]\left(t\right)+2\mathrm{\Lambda }\varsigma \frac{S\left(t\right)}{N\left(t\right)},\end{array}$(1e) with the boundary conditions (2a) $\begin{array}{rl}{I}_k(t,0)& =\frac{k{S}_k\left(t\right)}{\sum _{k}k{S}_k\left(t\right)}{\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[SI\right(t,a\left)\right]\mathrm{d}a,\end{array}$(2a) (2b) $\begin{array}{rl}\left[SI\right(t,0\left)\right]& ={\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[SSI\right(t,a\left)\right]\mathrm{d}a,\end{array}$(2b) and the initial conditions $\left(S_k\right(0),I_k(0,a),R_k(0),[SI(0,a)],[SS\left]\right(0\left)\right)$. Equations (1a(1a) $\begin{array}{rl}{\dot{S}}_k\left(t\right)& =\mathrm{\Lambda }{\theta }_k-\eta S_k\left(t\right)-\frac{k{S}_k\left(t\right)}{\sum _{k}k{S}_k\left(t\right)}{\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[SI\right(t,a\left)\right]\mathrm{d}a\\ & +\mathrm{\Lambda }\varsigma \left(\frac{S_{k-1}\left(t\right)}{N\left(t\right)}-\frac{S_k\left(t\right)}{N\left(t\right)}\right)-\eta \left(k{S}_k\right(t)-(k+1\left)S_{k+1}\right(t\left)\right),\end{array}$(1a) )–(1c(1c) $\begin{array}{rl}{\dot{R}}_k\left(t\right)& ={\int }_0^{\mathrm{\infty }}\gamma \left(a\right)I_k(t,a)\mathrm{d}a+\mathrm{\Lambda }\varsigma \left(\frac{R_{k-1}\left(t\right)}{N\left(t\right)}-\frac{R_k\left(t\right)}{N\left(t\right)}\right)\\ & -\eta R_k\left(t\right)-\eta \left(k{R}_k\right(t)-(k+1\left)R_{k+1}\right(t\left)\right),\end{array}$(1c) ) describe the dynamics of nodes and Equations (4d(4d) $\begin{array}{rl}& \dot{\left[SS\right]}\left(t\right)=-2\frac{\left[SS\right]\left(t\right)}{S\left(t\right)}{\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[SI\right(t,a\left)\right]\mathrm{d}a-2\eta \left[SS\right]\left(t\right)+2\mathrm{\Lambda }\varsigma \frac{S\left(t\right)}{N\left(t\right)},\end{array}$(4d) )–(1e(1e) $\begin{array}{rl}\dot{\left[SS\right]}\left(t\right)& =-2{\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[SSI\right(t,a\left)\right]\mathrm{d}a-2\eta \left[SS\right]\left(t\right)+2\mathrm{\Lambda }\varsigma \frac{S\left(t\right)}{N\left(t\right)},\end{array}$(1e) ) describe the dynamics of links. On the right-hand side of Equation (1a(1a) $\begin{array}{rl}{\dot{S}}_k\left(t\right)& =\mathrm{\Lambda }{\theta }_k-\eta S_k\left(t\right)-\frac{k{S}_k\left(t\right)}{\sum _{k}k{S}_k\left(t\right)}{\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[SI\right(t,a\left)\right]\mathrm{d}a\\ & +\mathrm{\Lambda }\varsigma \left(\frac{S_{k-1}\left(t\right)}{N\left(t\right)}-\frac{S_k\left(t\right)}{N\left(t\right)}\right)-\eta \left(k{S}_k\right(t)-(k+1\left)S_{k+1}\right(t\left)\right),\end{array}$(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. Equations (1b(1b) $\begin{array}{rl}\left(\frac{\mathrm{\partial }}{\mathrm{\partial }t}+\frac{\mathrm{\partial }}{\mathrm{\partial }a}\right)I_k(t,a)& =-\gamma \left(a\right)I_k(t,a)+\mathrm{\Lambda }\varsigma \left(\frac{I_{k-1}(t,a)}{N\left(t\right)}-\frac{I_k(t,a)}{N\left(t\right)}\right)\\ & -\eta I_k(t,a)-\eta \left(k{I}_k\right(t,a)-(k+1\left)I_{k+1}\right(t,a\left)\right),\end{array}$(1b) ) and (1c(1c) $\begin{array}{rl}{\dot{R}}_k\left(t\right)& ={\int }_0^{\mathrm{\infty }}\gamma \left(a\right)I_k(t,a)\mathrm{d}a+\mathrm{\Lambda }\varsigma \left(\frac{R_{k-1}\left(t\right)}{N\left(t\right)}-\frac{R_k\left(t\right)}{N\left(t\right)}\right)\\ & -\eta R_k\left(t\right)-\eta \left(k{R}_k\right(t)-(k+1\left)R_{k+1}\right(t\left)\right),\end{array}$(1c) ) can be explained in a similar manner. For Equation (1d(1d) $\begin{array}{rl}\left(\frac{\mathrm{\partial }}{\mathrm{\partial }t}+\frac{\mathrm{\partial }}{\mathrm{\partial }a}\right)\left[SI\right(t,a\left)\right]& =-\left(\lambda \right(a)+\gamma (a)+2\eta )\left[SI\right(t,a\left)\right]\\ & -{\int }_0^{\mathrm{\infty }}\lambda \left(b\right)\left[I\right(t,b\left)SI\right(t,a\left)\right]\mathrm{d}b+\mathrm{\Lambda }\varsigma \frac{I(t,a)}{N\left(t\right)},\end{array}$(1d) ), the first term represents the decrease of $\left[SI\right(t,a\left)\right]$ due to infection of susceptible nodes, recovery of infectious nodes and natural death of nodes along the links. The second term corresponds to the decrease of $\left[SI\right(t,a\left)\right]$ due to infection of susceptible nodes by another infectious neighbour with arbitrary infection age. The last term describes the fact that the number of $\left[SI\right(t,a\left)\right]$ increases when newborns are connected to infectious nodes with infection age a in the network at time t. Equation (1e(1e) $\begin{array}{rl}\dot{\left[SS\right]}\left(t\right)& =-2{\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[SSI\right(t,a\left)\right]\mathrm{d}a-2\eta \left[SS\right]\left(t\right)+2\mathrm{\Lambda }\varsigma \frac{S\left(t\right)}{N\left(t\right)},\end{array}$(1e) ) is analogous. Equation (2a(2a) $\begin{array}{rl}{I}_k(t,0)& =\frac{k{S}_k\left(t\right)}{\sum _{k}k{S}_k\left(t\right)}{\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[SI\right(t,a\left)\right]\mathrm{d}a,\end{array}$(2a) ) represents that susceptible nodes with degree k become new infectious nodes when they are infected by one infectious neighbour. Equation (2b(2b) $\begin{array}{rl}\left[SI\right(t,0\left)\right]& ={\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[SSI\right(t,a\left)\right]\mathrm{d}a,\end{array}$(2b) ) denotes the number of S−S links becoming S−I links when one of the susceptible nodes is infected by its infectious neighbour.

It is observed that the number of pairs depends on the number of triples in Equation (1(1a) $\begin{array}{rl}{\dot{S}}_k\left(t\right)& =\mathrm{\Lambda }{\theta }_k-\eta S_k\left(t\right)-\frac{k{S}_k\left(t\right)}{\sum _{k}k{S}_k\left(t\right)}{\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[SI\right(t,a\left)\right]\mathrm{d}a\\ & +\mathrm{\Lambda }\varsigma \left(\frac{S_{k-1}\left(t\right)}{N\left(t\right)}-\frac{S_k\left(t\right)}{N\left(t\right)}\right)-\eta \left(k{S}_k\right(t)-(k+1\left)S_{k+1}\right(t\left)\right),\end{array}$(1a) ). To break the dependency on higher order moments and obtain a closed model, the triple closure approximation in [25] (3) $\left[XSY\right]\left(t\right)=\frac{\left[SX\right]\left(t\right)\left[SY\right]\left(t\right)}{S\left(t\right)}$(3) 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\left(t\right)$ by the evolution of the population size $N\left(t\right)$, model (1(1a) $\begin{array}{rl}{\dot{S}}_k\left(t\right)& =\mathrm{\Lambda }{\theta }_k-\eta S_k\left(t\right)-\frac{k{S}_k\left(t\right)}{\sum _{k}k{S}_k\left(t\right)}{\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[SI\right(t,a\left)\right]\mathrm{d}a\\ & +\mathrm{\Lambda }\varsigma \left(\frac{S_{k-1}\left(t\right)}{N\left(t\right)}-\frac{S_k\left(t\right)}{N\left(t\right)}\right)-\eta \left(k{S}_k\right(t)-(k+1\left)S_{k+1}\right(t\left)\right),\end{array}$(1a) ) is further reduced to (4a) $\begin{array}{rl}& \dot{S}\left(t\right)=\mathrm{\Lambda }-{\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[SI\right(t,a\left)\right]\mathrm{d}a-\eta S\left(t\right),\end{array}$(4a) (4b) $\begin{array}{rl}& \left(\frac{\mathrm{\partial }}{\mathrm{\partial }t}+\frac{\mathrm{\partial }}{\mathrm{\partial }a}\right)I(t,a)=-\left(\gamma \right(a)+\eta )I(t,a),\end{array}$(4b) (4c) $\begin{array}{rl}& \left(\frac{\mathrm{\partial }}{\mathrm{\partial }t}+\frac{\mathrm{\partial }}{\mathrm{\partial }a}\right)\left[SI\right(t,a\left)\right]=-\left(\lambda \right(a)+\gamma (a)+2\eta )\left[SI\right(t,a\left)\right]\\ & -\frac{\left[SI\right(t,a\left)\right]}{S\left(t\right)}{\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[SI\right(t,a\left)\right]\mathrm{d}a+\mathrm{\Lambda }\varsigma \frac{I(t,a)}{N\left(t\right)},\end{array}$(4c) (4d) $\begin{array}{rl}& \dot{\left[SS\right]}\left(t\right)=-2\frac{\left[SS\right]\left(t\right)}{S\left(t\right)}{\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[SI\right(t,a\left)\right]\mathrm{d}a-2\eta \left[SS\right]\left(t\right)+2\mathrm{\Lambda }\varsigma \frac{S\left(t\right)}{N\left(t\right)},\end{array}$(4d) (4e) $\begin{array}{rl}& \dot{N}\left(t\right)=\mathrm{\Lambda }-\eta N\left(t\right),\end{array}$(4e) with the boundary conditions (5a) $\begin{array}{rl}I(t,0)& ={\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[SI\right(t,a\left)\right]\mathrm{d}a,\end{array}$(5a) (5b) $\begin{array}{rl}\left[SI\right(t,0\left)\right]& =\frac{\left[SS\right]\left(t\right)}{S\left(t\right)}{\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[SI\right(t,a\left)\right]\mathrm{d}a,\end{array}$(5b) and the initial conditions (6) $S\left(0\right)=S_0,I(0,a)=\varphi \left(a\right),\left[SI\right(0,a\left)\right]=\psi \left(a\right),\left[SS\right]\left(0\right)=[SS{]}_0,N(0)=N_0,$(6) where $N_0,S_0,[SS{]}_0\in {\mathbb{R}}_{+}=[0,\mathrm{\infty })$, and $\varphi \left(a\right),\psi \left(a\right)\in L_{+}^1$, and the initial conditions are determined by the initial network. Here $L_{+}^1:[0,\mathrm{\infty })\to {\mathbb{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(4a) $\begin{array}{rl}& \dot{S}\left(t\right)=\mathrm{\Lambda }-{\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[SI\right(t,a\left)\right]\mathrm{d}a-\eta S\left(t\right),\end{array}$(4a) ) satisfies the compatibility conditions (7) $\varphi \left(0\right)={\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\psi \left(a\right)\mathrm{d}a,\psi \left(0\right)=\frac{[SS{]}_0}{S_0}{\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\psi \left(a\right)\mathrm{d}a.$(7) Then model (4(4a) $\begin{array}{rl}& \dot{S}\left(t\right)=\mathrm{\Lambda }-{\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[SI\right(t,a\left)\right]\mathrm{d}a-\eta S\left(t\right),\end{array}$(4a) ) is well posed [7, 32]. In the following sections, we mainly focus on the low-dimensional model (4(4a) $\begin{array}{rl}& \dot{S}\left(t\right)=\mathrm{\Lambda }-{\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[SI\right(t,a\left)\right]\mathrm{d}a-\eta S\left(t\right),\end{array}$(4a) ).

Preliminary results. We first consider the existence of equilibria to model (4(4a) $\begin{array}{rl}& \dot{S}\left(t\right)=\mathrm{\Lambda }-{\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[SI\right(t,a\left)\right]\mathrm{d}a-\eta S\left(t\right),\end{array}$(4a) ). For simplicity of notation, define $\pi :{\mathbb{R}}_{+}\to {\mathbb{R}}_{+}$ by $\pi \left(a\right)={\mathrm{e}}^{-{\int }_0^a\left(\gamma \right(s)+\eta )\mathrm{d}s}$ and $\nu :{\mathbb{R}}_{+}\to {\mathbb{R}}_{+}$ by $\nu \left(a\right)={\mathrm{e}}^{-{\int }_0^a\left(\lambda \right(s)+\gamma (s)+2\eta )\mathrm{d}s}.$ Model (4(4a) $\begin{array}{rl}& \dot{S}\left(t\right)=\mathrm{\Lambda }-{\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[SI\right(t,a\left)\right]\mathrm{d}a-\eta S\left(t\right),\end{array}$(4a) ) can be rewritten as the following Volterra-type equations: (8a) $\begin{array}{rl}\dot{S}\left(t\right)& =\mathrm{\Lambda }-{\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[SI\right(t,a\left)\right]\mathrm{d}a-\eta S\left(t\right),\end{array}$(8a) (8b) $\begin{array}{rl}I(t,a)& =\left\{\begin{array}{ll}\beta (t-a)\pi \left(a\right),& t>a,\\ \varphi (a-t)\frac{\pi \left(a\right)}{\pi (a-t)},& t⩽a,\end{array}\right.\end{array}$(8b) (8c) $\begin{array}{rl}\left[SI\right(t,a\left)\right]& =\left\{\begin{array}{ll}\beta (t-a)\nu \left(a\right)\left\{\frac{\left[SS\right](t-a)}{S(t-a)}{\mathrm{e}}^{-{\int }_{t-a}^t\left(\beta \left(s\right)/S\left(s\right)\right)\mathrm{d}s}\right.& \\ \left.+{\int }_0^a\frac{\mathrm{\Lambda }\varsigma \pi \left(s\right)}{N(t-a+s)\nu \left(s\right)}{\mathrm{e}}^{-{\int }_{t-a+s}^t\left(\beta \left(w\right)/S\left(w\right)\right)\mathrm{d}w}\mathrm{d}s\right\},& t>a,\\ \frac{\nu \left(a\right)\psi (a-t)}{\nu (a-t)}{\mathrm{e}}^{-{\int }_0^t\left(\beta \left(s\right)/S\left(s\right)\right)\mathrm{d}s}& \\ +{\int }_0^t\frac{\mathrm{\Lambda }\varsigma \varphi (a-t)\pi (a-t+s)\nu \left(a\right)}{N\left(s\right)\pi (a-t)\nu (a-t+s)}{\mathrm{e}}^{-{\int }_s^t\left(\beta \left(w\right)/S\left(w\right)\right)\mathrm{d}w}\mathrm{d}s,& t⩽a,\end{array}\right.\end{array}$(8c) (8d) $\begin{array}{rl}\dot{\left[SS\right]}\left(t\right)& =-2\frac{\left[SS\right]\left(t\right)}{S\left(t\right)}{\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[SI\right(t,a\left)\right]\mathrm{d}a-2\eta \left[SS\right]\left(t\right)+2\mathrm{\Lambda }\varsigma \frac{S\left(t\right)}{N\left(t\right)},\end{array}$(8d) (8e) $\begin{array}{rl}\dot{N}\left(t\right)& =\mathrm{\Lambda }-\eta N\left(t\right),\end{array}$(8e) where $\beta \left(t\right)$ is the unique solution of the following nonlinear Volterra equation: $\begin{array}{rl}\beta \left(t\right)& ={\int }_0^t\lambda \left(a\right)\nu \left(a\right)\beta (t-a)\left(\frac{\left[SS\right](t-a)}{S(t-a)}{\mathrm{e}}^{-{\int }_{t-a}^t\left(\beta \left(s\right)/S\left(s\right)\right)\mathrm{d}s}\right.\\ & \left.+{\int }_0^a\frac{\mathrm{\Lambda }\varsigma \pi \left(s\right)}{N(t-a+s)\nu \left(s\right)}{\mathrm{e}}^{-{\int }_{t-a+s}^t\left(\beta \left(\omega \right)/S\left(\omega \right)\right)\mathrm{d}\omega }\mathrm{d}s\right)\mathrm{d}a\\ & +{\int }_t^{\mathrm{\infty }}\lambda \left(a\right)\nu \left(a\right)\left(\frac{\psi (a-t)}{\nu (a-t)}{\mathrm{e}}^{-{\int }_0^t\left(\beta \left(s\right)/S\left(s\right)\right)\mathrm{d}s}\right.\\ & \left.+{\int }_0^t\frac{\mathrm{\Lambda }\varsigma \varphi (a-t)\pi (a-t+s)}{N\left(s\right)\pi (a-t)\nu (a-t+s)}{\mathrm{e}}^{-{\int }_s^t\left(\beta \left(\omega \right)/S\left(\omega \right)\right)\mathrm{d}\omega }\mathrm{d}s\right)\mathrm{d}a.\end{array}$ This is derived by using (9) $\beta \left(t\right)={\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[SI\right(t,a\left)\right]\mathrm{d}a=I(t,0).$(9) We then apply the approach introduced by Thieme [30] to reformulate model (4(4a) $\begin{array}{rl}& \dot{S}\left(t\right)=\mathrm{\Lambda }-{\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[SI\right(t,a\left)\right]\mathrm{d}a-\eta S\left(t\right),\end{array}$(4a) ) with the boundary conditions (5(5a) $\begin{array}{rl}I(t,0)& ={\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[SI\right(t,a\left)\right]\mathrm{d}a,\end{array}$(5a) ) and initial conditions (6(6) $S\left(0\right)=S_0,I(0,a)=\varphi \left(a\right),\left[SI\right(0,a\left)\right]=\psi \left(a\right),\left[SS\right]\left(0\right)=[SS{]}_0,N(0)=N_0,$(6) ) as a semilinear Cauchy problem.

Rearrange the variables of model (4(4a) $\begin{array}{rl}& \dot{S}\left(t\right)=\mathrm{\Lambda }-{\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[SI\right(t,a\left)\right]\mathrm{d}a-\eta S\left(t\right),\end{array}$(4a) ) as $\left(S\right(t),[SS\left]\right(t),N(t),I(t,a),[SI(t,a)\left]\right)$. In order to take care of the boundary conditions, we enlarge the state space and consider $\mathcal{X}=\mathbb{R}\times \mathbb{R}\times \mathbb{R}\times \mathbb{R}\times \mathbb{R}\times L^1\left(\right(0,\mathrm{\infty }),\mathbb{R})\times L^1\left(\right(0,\mathrm{\infty }),\mathbb{R}),$ endowed the usual product norm, and the set $\begin{array}{rl}{\mathcal{X}}_0& =\mathbb{R}\times \mathbb{R}\times \mathbb{R}\times \left\{0\right\}\times \left\{0\right\}\times L^1\left(\right(0,\mathrm{\infty }),\mathbb{R})\times L^1\left(\right(0,\mathrm{\infty }),\mathbb{R}),\\ {\mathcal{X}}_{+}& ={\mathbb{R}}_{+}\times {\mathbb{R}}_{+}\times {\mathbb{R}}_{+}\times {\mathbb{R}}_{+}\times {\mathbb{R}}_{+}\times L_{+}^1\left(\right(0,\mathrm{\infty }),\mathbb{R})\times L_{+}^1\left(\right(0,\mathrm{\infty }),\mathbb{R}),\end{array}$ and ${\mathcal{X}}_{0+}={\mathcal{X}}_0\cap {\mathcal{X}}_{+}.$ We consider the linear operator $\mathcal{A}:\mathrm{Dom}\left(\mathcal{A}\right)\subset \mathcal{X}\to \mathcal{X}$ defined by $\mathcal{A}\left(\begin{array}{c}S\\ \left[SS\right]\\ N\\ \left(\begin{array}{c}0\\ 0\\ I\\ \left[SI\right]\end{array}\right)\end{array}\right)=\left(\begin{array}{c}-\eta S\\ -2\eta \left[SS\right]\\ -\eta N\\ \left(\begin{array}{c}-I\left(0\right)\\ \left[SI\right(0\left)\right]\\ -\left(\gamma \right(a)+\eta )I-I^{\prime }\\ -\left(\lambda \right(a)+\gamma (a)+2\eta )\left[SI\right]-{\left[SI\right]}^{\prime }\end{array}\right)\end{array}\right),$ with $\text{Dom}\left(\mathcal{A}\right)=\mathbb{R}\times \mathbb{R}\times \mathbb{R}\times \left\{0\right\}\times \left\{0\right\}\times W^{1,1}\left(\right(0,\mathrm{\infty }),\mathbb{R})\times W^{1,1}\left(\right(0,\mathrm{\infty }),\mathbb{R}),$ where $W^{1,1}$ is a Sobolev space, and we define $\mathcal{F}:{\mathcal{X}}_0\to \mathcal{X}$ by $\mathcal{F}\left(\begin{array}{c}S\\ \left[SS\right]\\ N\\ \left(\begin{array}{c}0\\ 0\\ I\\ \left[SI\right]\end{array}\right)\end{array}\right)=\left(\begin{array}{c}\mathrm{\Lambda }-{\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[SI\right]\left(a\right)\mathrm{d}a\\ -2\frac{\left[SS\right]}{S}{\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[SI\right]\left(a\right)\mathrm{d}a+2\mathrm{\Lambda }\varsigma \frac{S}{N}\\ \mathrm{\Lambda }\\ \left(\begin{array}{c}{\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[SI\right]\left(a\right)\mathrm{d}a\\ \frac{\left[SS\right]}{S}{\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[SI\right]\left(a\right)\mathrm{d}a\\ 0_{L^1}\\ -\frac{\left[SI\right]\left(a\right)}{S}{\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[SI\right]\left(a\right)\mathrm{d}a+\mathrm{\Lambda }\varsigma \frac{I\left(a\right)}{N}+0_{L^1}\end{array}\right)\end{array}\right).$ Then by defining $\mathit{u}\left(t\right)=\left(\begin{array}{c}S\left(t\right)\\ \left[SS\right]\left(t\right)\\ N\left(t\right)\\ \left(\begin{array}{c}0\\ 0\\ I(t,\cdot )\\ \left[SI\right(t,\cdot \left)\right]\end{array}\right)\end{array}\right),$ we can reformulate the partial differential equation problem (4(4a) $\begin{array}{rl}& \dot{S}\left(t\right)=\mathrm{\Lambda }-{\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[SI\right(t,a\left)\right]\mathrm{d}a-\eta S\left(t\right),\end{array}$(4a) ) as the following abstract Cauchy problem: (10) $\frac{\mathrm{d}\mathit{u}\left(t\right)}{\mathrm{d}t}=\mathcal{A}\mathit{u}\left(t\right)+\mathcal{F}\left(\mathit{u}\right(t\left)\right)\mathrm{for}t⩾0\mathrm{and}\mathit{u}\left(0\right)={\mathit{u}}_0\in {\mathcal{X}}_{0+}.$(10) By using the results in [17, 18, 30], we derive the existence and the uniqueness of the semiflow $\left\{U\right(t){\}}_{t⩾0}$ on ${\mathcal{X}}_{0+}$. By identifying $\left(S\right(t),[SS\left]\right(t),N(t),0,0,$ $I(t,\cdot ),\left[SI\right(t,\cdot \left)\right])$ with $\left(S\right(t),[SS\left]\right(t),N(t),I(t,\cdot ),[SI(t,\cdot )\left]\right)$, 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(4a) $\begin{array}{rl}& \dot{S}\left(t\right)=\mathrm{\Lambda }-{\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[SI\right(t,a\left)\right]\mathrm{d}a-\eta S\left(t\right),\end{array}$(4a) ): $\begin{array}{rl}& \frac{\mathrm{d}N\left(t\right)}{\mathrm{d}t}=\mathrm{\Lambda }-\eta N\left(t\right),\\ & \frac{\mathrm{d}\left(S\right(t)+{\int }_0^{\mathrm{\infty }}I(t,a)\mathrm{d}a-N(t\left)\right)}{\mathrm{d}t}⩽-\eta \left(S\right(t)+{\int }_0^{\mathrm{\infty }}I(t,a)\mathrm{d}a-N(t\left)\right),\end{array}$ and $\frac{\mathrm{d}\left(\right[SS\left]\right(t)+2{\int }_0^{\mathrm{\infty }}[SI(t,a)\left]\mathrm{d}a\right)}{\mathrm{d}t}⩽-2\eta \left(\right[SS\left]\right(t)+2{\int }_0^{\mathrm{\infty }}[SI(t,a)\left]\right)+2\mathrm{\Lambda }\varsigma .$ Thus, if $S\left(t\right)+{\int }_0^{\mathrm{\infty }}I(t,a)\mathrm{d}a⩽N\left(t\right),\left[SS\right]\left(t\right)+2{\int }_0^{\mathrm{\infty }}\left[SI\right(t,a\left)\right]\mathrm{d}a⩽\frac{\mathrm{\Lambda }\varsigma }{\eta },N\left(t\right)⩽\frac{\mathrm{\Lambda }}{\eta }$ are satisfied for some $t=t_0$, they are satisfied for all $t⩾t_0$. Hence, model (4(4a) $\begin{array}{rl}& \dot{S}\left(t\right)=\mathrm{\Lambda }-{\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[SI\right(t,a\left)\right]\mathrm{d}a-\eta S\left(t\right),\end{array}$(4a) ) leaves the set $\begin{array}{rl}\mathrm{\Omega }& =\left\{(S,[SS],N,I(\cdot ),[SI(\cdot )\left]\right)\in {\mathbb{R}}_{+}\times {\mathbb{R}}_{+}\times {\mathbb{R}}_{+}\times L_{+}^1\left(\right(0,\mathrm{\infty }),\mathbb{R})\times L_{+}^1\left(\right(0,\mathrm{\infty }),\mathbb{R})|\phantom{\frac{\mathrm{\Lambda }\varsigma }{\eta }}\right.\\ & \left.S+{\int }_0^{\mathrm{\infty }}I\left(a\right)\mathrm{d}a⩽N,\left[SS\right]+2{\int }_0^{\mathrm{\infty }}\left[SI\right(a\left)\right]\mathrm{d}a⩽\frac{\mathrm{\Lambda }\varsigma }{\eta },N⩽\frac{\mathrm{\Lambda }}{\eta }\right\}\end{array}$ positively invariant. Furthermore, it follows that the set $\begin{array}{rl}B& =\left\{\left(\begin{array}{c}S\\ \left[SS\right]\\ N\\ \left(\begin{array}{c}0\\ 0\\ I\\ \left[SI\right]\end{array}\right)\end{array}\right)\in {\mathcal{X}}_{0+}:S+{\int }_0^{\mathrm{\infty }}I\left(a\right)\mathrm{d}a⩽N,\right.\\ & \left.\phantom{\left(\begin{array}{c}S\\ \left[SS\right]\\ N\\ \left(\begin{array}{c}0\\ 0\\ I\\ \left[SI\right]\end{array}\right)\end{array}\right)}\left[SS\right]+2{\int }_0^{\mathrm{\infty }}\left[SI\right(a\left)\right]\mathrm{d}a⩽\frac{\mathrm{\Lambda }\varsigma }{\eta },N⩽\frac{\mathrm{\Lambda }}{\eta }\right\}\end{array}$ is a positively invariant absorbing set under the semiflow $\left\{U\right(t){\}}_{t⩾0}$ on ${\mathcal{X}}_{0+}$, i.e. $U\left(t\right)B\subseteq B$ and for each $x=(S_0,[SS{]}_0,N_0,\varphi (\cdot ),\psi (\cdot ))\in {\mathcal{X}}_{0+}$, $\mathrm{d}\left(U\right(t)x,B):=\underset{y\in B}{inf}\parallel U\left(t\right)x-y\parallel \to 0\text{as}t\to \mathrm{\infty }.$ This means that $\left\{U\right(t){\}}_{t⩾0}$ is bounded dissipative on ${\mathcal{X}}_{0+}$ (see [5]).

3. Equilibria and the basic reproduction number

In this section, we will derive the basic reproduction number ${\mathcal{R}}_0$ by analysing the equilibria. To find the equilibria, we look for time-independent solutions $(\overline{S},\overline{I}(a),[\overline{SI}\left(a\right)],[\overline{SS}],\overline{N})$ that satisfy model (4(4a) $\begin{array}{rl}& \dot{S}\left(t\right)=\mathrm{\Lambda }-{\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[SI\right(t,a\left)\right]\mathrm{d}a-\eta S\left(t\right),\end{array}$(4a) ) with the time derivatives equal to zero. Then the equilibria satisfy the following equations: (11a) $\begin{array}{rl}& \mathrm{\Lambda }-{\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[\overline{SI}\right(a\left)\right]\mathrm{d}a-\eta \overline{S}=0,\end{array}$(11a) (11b) $\begin{array}{rl}& \frac{\mathrm{d}\overline{I}\left(a\right)}{\mathrm{d}a}=-\left(\gamma \right(a)+\eta )\overline{I}\left(a\right),\end{array}$(11b) (11c) $\begin{array}{rl}& \frac{\mathrm{d}\left[\overline{SI}\right(a\left)\right]}{\mathrm{d}a}=-\left(\lambda \right(a)+\gamma (a)+2\eta )\left[\overline{SI}\right(a\left)\right]\\ & -\frac{{\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[\overline{SI}\right(a\left)\right]\mathrm{d}a}{\overline{S}}\left[\overline{SI}\right(a\left)\right]+\mathrm{\Lambda }\varsigma \frac{\overline{I}\left(a\right)}{\overline{N}},\end{array}$(11c) (11d) $\begin{array}{rl}& -2\frac{{\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[\overline{SI}\right(a\left)\right]\mathrm{d}a}{\overline{S}}\left[\overline{SS}\right]-2\eta \left[\overline{SS}\right]+2\mathrm{\Lambda }\varsigma \frac{\overline{S}}{\overline{N}}=0,\end{array}$(11d) (11e) $\begin{array}{rl}& \mathrm{\Lambda }-\eta \overline{N}=0,\end{array}$(11e) (11f) $\begin{array}{rl}& \left[\overline{SI}\right(0\left)\right]=\frac{{\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[\overline{SI}\right(a\left)\right]\mathrm{d}a}{\overline{S}}\left[\overline{SS}\right],\end{array}$(11f) (11g) $\begin{array}{rl}& \overline{I}\left(0\right)={\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[\overline{SI}\right(a\left)\right]\mathrm{d}a.\end{array}$(11g) It is easy to see that $0⩽\overline{I}\left(0\right)⩽\mathrm{\Lambda }$, and $E_0=(\mathrm{\Lambda }/\eta ,0,0,(\mathrm{\Lambda }/\eta )\varsigma ,\mathrm{\Lambda }/\eta )$ is one solution of Equation (11(11a) $\begin{array}{rl}& \mathrm{\Lambda }-{\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[\overline{SI}\right(a\left)\right]\mathrm{d}a-\eta \overline{S}=0,\end{array}$(11a) ). This solution produces the disease-free equilibrium, and it always exists. Next, an endemic equilibrium will be given by a nontrivial positive solution $E^{\ast }=(S^{\ast },I^{\ast }(a),[SI\left(a\right){]}^{\ast },[SS{]}^{\ast },N^{\ast })$, in which $I^{\ast }\left(a\right)\ne 0$. By Equations (11a(11a) $\begin{array}{rl}& \mathrm{\Lambda }-{\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[\overline{SI}\right(a\left)\right]\mathrm{d}a-\eta \overline{S}=0,\end{array}$(11a) ), (11b(11b) $\begin{array}{rl}& \frac{\mathrm{d}\overline{I}\left(a\right)}{\mathrm{d}a}=-\left(\gamma \right(a)+\eta )\overline{I}\left(a\right),\end{array}$(11b) ), (11d(11d) $\begin{array}{rl}& -2\frac{{\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[\overline{SI}\right(a\left)\right]\mathrm{d}a}{\overline{S}}\left[\overline{SS}\right]-2\eta \left[\overline{SS}\right]+2\mathrm{\Lambda }\varsigma \frac{\overline{S}}{\overline{N}}=0,\end{array}$(11d) ), (11e(11e) $\begin{array}{rl}& \mathrm{\Lambda }-\eta \overline{N}=0,\end{array}$(11e) ), and (11g(11g) $\begin{array}{rl}& \overline{I}\left(0\right)={\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[\overline{SI}\right(a\left)\right]\mathrm{d}a.\end{array}$(11g) ), we can express $S^{\ast },I^{\ast }\left(a\right),N^{\ast },[SS{]}^{\ast }$, and $\left[SI\right(0){]}^{\ast }$ in terms of $I^{\ast }\left(0\right)$: (12a) $\begin{array}{rl}& S^{\ast }=\frac{\mathrm{\Lambda }-I^{\ast }\left(0\right)}{\eta },I^{\ast }\left(a\right)=I^{\ast }\left(0\right)\pi \left(a\right),N^{\ast }=\frac{\mathrm{\Lambda }}{\eta },\end{array}$(12a) (12b) $\begin{array}{rl}& [SS{]}^{\ast }=\frac{\varsigma (\mathrm{\Lambda }-I^{\ast }(0){)}^2}{\eta \mathrm{\Lambda }},[SI\left(0\right){]}^{\ast }=\frac{\varsigma (\mathrm{\Lambda }-I^{\ast }(0\left)\right)I^{\ast }\left(0\right)}{\mathrm{\Lambda }}.\end{array}$(12b) Substituting them into Equation (11c(11c) $\begin{array}{rl}& \frac{\mathrm{d}\left[\overline{SI}\right(a\left)\right]}{\mathrm{d}a}=-\left(\lambda \right(a)+\gamma (a)+2\eta )\left[\overline{SI}\right(a\left)\right]\\ & -\frac{{\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[\overline{SI}\right(a\left)\right]\mathrm{d}a}{\overline{S}}\left[\overline{SI}\right(a\left)\right]+\mathrm{\Lambda }\varsigma \frac{\overline{I}\left(a\right)}{\overline{N}},\end{array}$(11c) ) yields a differential equation for $\left[SI\right(a){]}^{\ast }$: $\frac{\mathrm{d}\left[SI\right(a){]}^{\ast }}{\mathrm{d}a}=-\left(\lambda \left(a\right)+\gamma \left(a\right)+2\eta +\frac{\eta I^{\ast }\left(0\right)}{\mathrm{\Lambda }-I^{\ast }\left(0\right)}\right)\left[SI\right(a){]}^{\ast }+\eta \varsigma \pi (a\left)I^{\ast }\right(0),$ whose solution can be given as follows: (13) $\left[SI\right(a){]}^{\ast }=\nu (a){\mathrm{e}}^{-\left(\eta I^{\ast }\left(0\right)/(\mathrm{\Lambda }-I^{\ast }(0\left)\right)\right)a}\left(\left[SI\right(0){]}^{\ast }+\eta \varsigma I^{\ast }(0){\int }_0^a\frac{\pi \left(s\right)}{\nu \left(s\right)}{e}^{\left(\eta I^{\ast }\left(0\right)/(\mathrm{\Lambda }-I^{\ast }(0\left)\right)\right)s}\mathrm{d}s\right).$(13) To find $I^{\ast }\left(0\right)$, substituting Equation (13(13) $\left[SI\right(a){]}^{\ast }=\nu (a){\mathrm{e}}^{-\left(\eta I^{\ast }\left(0\right)/(\mathrm{\Lambda }-I^{\ast }(0\left)\right)\right)a}\left(\left[SI\right(0){]}^{\ast }+\eta \varsigma I^{\ast }(0){\int }_0^a\frac{\pi \left(s\right)}{\nu \left(s\right)}{e}^{\left(\eta I^{\ast }\left(0\right)/(\mathrm{\Lambda }-I^{\ast }(0\left)\right)\right)s}\mathrm{d}s\right).$(13) ) into Equation (11g(11g) $\begin{array}{rl}& \overline{I}\left(0\right)={\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[\overline{SI}\right(a\left)\right]\mathrm{d}a.\end{array}$(11g) ), we obtain (14) $\begin{array}{rl}& I^{\ast }\left(0\right)\varsigma {\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\nu \left(a\right)\left({\mathrm{e}}^{-\left(\eta I^{\ast }\left(0\right)/(\mathrm{\Lambda }-I^{\ast }(0\left)\right)\right)a}\frac{\mathrm{\Lambda }-I^{\ast }\left(0\right)}{\mathrm{\Lambda }}\right.\\ & \left.+\eta {\int }_0^a\frac{\pi \left(s\right)}{\nu \left(s\right)}{\mathrm{e}}^{-\left(\eta I^{\ast }\left(0\right)/(\mathrm{\Lambda }-I^{\ast }(0\left)\right)\right)(a-s)}\mathrm{d}s\right)\mathrm{d}a-I^{\ast }\left(0\right)=0.\end{array}$(14) Since we are looking for a nonzero solution, $I^{\ast }\left(0\right)\ne 0$. Cancelling $I^{\ast }\left(0\right)$ from both sides of Equation (14(14) $\begin{array}{rl}& I^{\ast }\left(0\right)\varsigma {\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\nu \left(a\right)\left({\mathrm{e}}^{-\left(\eta I^{\ast }\left(0\right)/(\mathrm{\Lambda }-I^{\ast }(0\left)\right)\right)a}\frac{\mathrm{\Lambda }-I^{\ast }\left(0\right)}{\mathrm{\Lambda }}\right.\\ & \left.+\eta {\int }_0^a\frac{\pi \left(s\right)}{\nu \left(s\right)}{\mathrm{e}}^{-\left(\eta I^{\ast }\left(0\right)/(\mathrm{\Lambda }-I^{\ast }(0\left)\right)\right)(a-s)}\mathrm{d}s\right)\mathrm{d}a-I^{\ast }\left(0\right)=0.\end{array}$(14) ) and denoting by $F\left(I^{\ast }\right(0\left)\right)$ the left-hand side, we have $\begin{array}{rl}& F\left(I^{\ast }\right(0\left)\right)=\varsigma {\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\nu \left(a\right)\left({\mathrm{e}}^{-\left(\eta I^{\ast }\left(0\right)/(\mathrm{\Lambda }-I^{\ast }(0\left)\right)\right)a}\frac{\mathrm{\Lambda }-I^{\ast }\left(0\right)}{\mathrm{\Lambda }}\right.\\ & \left.+\eta {\int }_0^a\frac{\pi \left(s\right)}{\nu \left(s\right)}{\mathrm{e}}^{-\left(\eta I^{\ast }\left(0\right)/(\mathrm{\Lambda }-I^{\ast }(0\left)\right)\right)(a-s)}\mathrm{d}s\right)\mathrm{d}a-1.\end{array}$ It follows that (15) $\underset{I^{\ast }\left(0\right)\to {\mathrm{\Lambda }}^{-}}{lim}F\left(I^{\ast }\right(0\left)\right)=-1$(15) and (16) $F^{\prime }\left(I^{\ast }\right(0\left)\right)<0\text{for}{I}^{\ast }\left(0\right)\in (0,\mathrm{\Lambda }).$(16) Thus, Equation (14(14) $\begin{array}{rl}& I^{\ast }\left(0\right)\varsigma {\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\nu \left(a\right)\left({\mathrm{e}}^{-\left(\eta I^{\ast }\left(0\right)/(\mathrm{\Lambda }-I^{\ast }(0\left)\right)\right)a}\frac{\mathrm{\Lambda }-I^{\ast }\left(0\right)}{\mathrm{\Lambda }}\right.\\ & \left.+\eta {\int }_0^a\frac{\pi \left(s\right)}{\nu \left(s\right)}{\mathrm{e}}^{-\left(\eta I^{\ast }\left(0\right)/(\mathrm{\Lambda }-I^{\ast }(0\left)\right)\right)(a-s)}\mathrm{d}s\right)\mathrm{d}a-I^{\ast }\left(0\right)=0.\end{array}$(14) ) has a unique positive solution in $(0,\mathrm{\Lambda })$ if and only if $F\left(0\right)={\mathcal{R}}_0-1>0$, where (17) ${\mathcal{R}}_0=\varsigma {\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\nu \left(a\right)\mathrm{d}a+\eta \varsigma {\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\nu \left(a\right)\left({\int }_0^a\frac{\pi \left(s\right)}{\nu \left(s\right)}\mathrm{d}s\right)\mathrm{d}a.$(17) Given the solution $I^{\ast }\left(0\right)$, according to Equation (12(12a) $\begin{array}{rl}& S^{\ast }=\frac{\mathrm{\Lambda }-I^{\ast }\left(0\right)}{\eta },I^{\ast }\left(a\right)=I^{\ast }\left(0\right)\pi \left(a\right),N^{\ast }=\frac{\mathrm{\Lambda }}{\eta },\end{array}$(12a) ), there correspond unique nonzero values $S^{\ast }$, $I^{\ast }\left(a\right)$, $\left[SI\right(a){]}^{\ast }$, $[SS{]}^{\ast }$, and $N^{\ast }$. We summarize this result in the following proposition.

Proposition 3.1

If ${\mathcal{R}}_0<1,$ model (4(4a) $\begin{array}{rl}& \dot{S}\left(t\right)=\mathrm{\Lambda }-{\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[SI\right(t,a\left)\right]\mathrm{d}a-\eta S\left(t\right),\end{array}$(4a) ) has only the disease-free equilibrium $E_0$. If ${\mathcal{R}}_0>1,$ model (4(4a) $\begin{array}{rl}& \dot{S}\left(t\right)=\mathrm{\Lambda }-{\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[SI\right(t,a\left)\right]\mathrm{d}a-\eta S\left(t\right),\end{array}$(4a) ) has a disease-free equilibrium $E_0$ and a unique endemic equilibrium $E^{\ast }$.

In epidemiology, ${\mathcal{R}}_0$ is referred to as the basic reproduction number of model (4(4a) $\begin{array}{rl}& \dot{S}\left(t\right)=\mathrm{\Lambda }-{\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[SI\right(t,a\left)\right]\mathrm{d}a-\eta S\left(t\right),\end{array}$(4a) ). Furthermore, for model (4(4a) $\begin{array}{rl}& \dot{S}\left(t\right)=\mathrm{\Lambda }-{\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[SI\right(t,a\left)\right]\mathrm{d}a-\eta S\left(t\right),\end{array}$(4a) ), we have the following results.

Theorem 3.1

The disease-free equilibrium $E_0$ of model (4(4a) $\begin{array}{rl}& \dot{S}\left(t\right)=\mathrm{\Lambda }-{\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[SI\right(t,a\left)\right]\mathrm{d}a-\eta S\left(t\right),\end{array}$(4a) ) is locally asymptotically stable if ${\mathcal{R}}_0<1$ and it is unstable if ${\mathcal{R}}_0>1$.

Proof.

First, we introduce the change of variables as follows: $\begin{array}{rl}s\left(t\right)& =S\left(t\right)-\frac{\mathrm{\Lambda }}{\eta },i(t,a)=I(t,a),\left[si\right(t,a\left)\right]=\left[SI\right(t,a\left)\right],\\ \left[ss\right]\left(t\right)& =\left[SS\right]\left(t\right)-\frac{\mathrm{\Lambda }}{\eta }\varsigma ,n\left(t\right)=N\left(t\right)-\frac{\mathrm{\Lambda }}{\eta }.\end{array}$ After dropping the nonlinear quadratic terms, we obtain the linearized system at the disease-free equilibrium $E_0$ as follows: (18a) $\dot{s}\left(t\right)=-{\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[si\right(t,a\left)\right]\mathrm{d}a-\eta s\left(t\right),$(18a) (18b) $\left(\frac{\mathrm{\partial }}{\mathrm{\partial }t}+\frac{\mathrm{\partial }}{\mathrm{\partial }a}\right)i(t,a)=-\left(\gamma \right(a)+\eta )i(t,a),$(18b) (18c) $\left(\frac{\mathrm{\partial }}{\mathrm{\partial }t}+\frac{\mathrm{\partial }}{\mathrm{\partial }a}\right)\left[si\right(t,a\left)\right]=-\left(\lambda \right(a)+\gamma (a)+2\eta )\left[si\right(t,a\left)\right]+\eta \varsigma i(t,a),$(18c) (18d) $\dot{\left[ss\right]}\left(t\right)=-2\varsigma {\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[si\right(t,a\left)\right]\mathrm{d}a-2\eta \left[ss\right]\left(t\right)+2\eta \varsigma s\left(t\right)-2\eta \varsigma n\left(t\right),$(18d) (18e) $\dot{n}\left(t\right)=-\eta n\left(t\right),$(18e) (18f) $i(t,0)={\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[si\right(t,a\left)\right]\mathrm{d}a,$(18f) (18g) $\left[si\right(t,0\left)\right]=\varsigma {\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[si\right(t,a\left)\right]\mathrm{d}a.$(18g) Since the system is linear, it has separable solutions in which the time-dependent function is an exponential. Hence, it is sensible to look for the solutions in the form $s\left(t\right)=s_0{\mathrm{e}}^{zt},$ $i(t,a)=i\left(a\right){\mathrm{e}}^{zt}$, $\left[si\right(t,a\left)\right]=\left[si\right(a\left)\right]{\mathrm{e}}^{zt}$, $\left[ss\right]\left(t\right)=[ss{]}_0{\mathrm{e}}^{zt},$ $n\left(t\right)=n_0{\mathrm{e}}^{zt}$, where $s_0$, $i\left(a\right)$, $\left[si\right(a\left)\right]$, $[ss{]}_0$, $n_0$ will be determined later. Substituting the solution into Equation (18(18a) $\dot{s}\left(t\right)=-{\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[si\right(t,a\left)\right]\mathrm{d}a-\eta s\left(t\right),$(18a) ) and cancelling ${\mathrm{e}}^{zt}$, we obtain the following equations: (19a) $\begin{array}{rl}& (z+\eta )s_0+{\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[si\right(a\left)\right]\mathrm{d}a=0,\end{array}$(19a) (19b) $\begin{array}{rl}& \frac{\mathrm{d}i\left(a\right)}{\mathrm{d}a}=-(z+\gamma (a)+\eta )i\left(a\right),\end{array}$(19b) (19c) $\begin{array}{rl}& \frac{\mathrm{d}\left[si\right(a\left)\right]}{\mathrm{d}a}=-(z+\lambda (a)+\gamma (a)+2\eta )\left[si\right(a\left)\right]+\eta \varsigma i\left(a\right),\end{array}$(19c) (19d) $\begin{array}{rl}& -(z+2\eta )[ss{]}_0-2\varsigma {\int }_0^{\mathrm{\infty }}\lambda (a\left)\right[si\left(a\right)]\mathrm{d}a+2\eta \varsigma s_0-2\eta \varsigma n_0=0,\end{array}$(19d) (19e) $\begin{array}{rl}& -(z+\eta )n_0=0,\end{array}$(19e) (19f) $\begin{array}{rl}& i\left(0\right)={\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[si\right(a\left)\right]\mathrm{d}a,\end{array}$(19f) (19g) $\begin{array}{rl}& \left[si\right(0\left)\right]=\varsigma {\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[si\right(a\left)\right]\mathrm{d}a.\end{array}$(19g) Integrating Equation (19b(19b) $\begin{array}{rl}& \frac{\mathrm{d}i\left(a\right)}{\mathrm{d}a}=-(z+\gamma (a)+\eta )i\left(a\right),\end{array}$(19b) ) from 0 to a yields $i\left(a\right)=i\left(0\right)\pi \left(a\right){\mathrm{e}}^{-za}.$ Substituting it into Equation (19c(19c) $\begin{array}{rl}& \frac{\mathrm{d}\left[si\right(a\left)\right]}{\mathrm{d}a}=-(z+\lambda (a)+\gamma (a)+2\eta )\left[si\right(a\left)\right]+\eta \varsigma i\left(a\right),\end{array}$(19c) ) and solving the differential equation yield (20) $\left[si\right(a\left)\right]={\mathrm{e}}^{-za}\nu \left(a\right)\left[si\right(0\left)\right]+\left(\eta \varsigma \nu \left(a\right){\mathrm{e}}^{-za}{\int }_0^a\frac{\pi \left(s\right)}{\nu \left(s\right)}\mathrm{d}s\right)i\left(0\right).$(20) Plugging Equation (20(20) $\left[si\right(a\left)\right]={\mathrm{e}}^{-za}\nu \left(a\right)\left[si\right(0\left)\right]+\left(\eta \varsigma \nu \left(a\right){\mathrm{e}}^{-za}{\int }_0^a\frac{\pi \left(s\right)}{\nu \left(s\right)}\mathrm{d}s\right)i\left(0\right).$(20) ) into Equation (19f(19f) $\begin{array}{rl}& i\left(0\right)={\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[si\right(a\left)\right]\mathrm{d}a,\end{array}$(19f) ), and combining with Equation (19g(19g) $\begin{array}{rl}& \left[si\right(0\left)\right]=\varsigma {\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[si\right(a\left)\right]\mathrm{d}a.\end{array}$(19g) ), we obtain the characteristic equation for z: (21) $P\left(z\right):=\varsigma {\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\nu \left(a\right){\mathrm{e}}^{-za}\mathrm{d}a+\eta \varsigma {\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\nu \left(a\right){\mathrm{e}}^{-za}\left({\int }_0^a\frac{\pi \left(s\right)}{\nu \left(s\right)}\mathrm{d}s\right)\mathrm{d}a=1.$(21) Obviously, $P\left(z\right)$ is a continuous function satisfying $\underset{z\to +\mathrm{\infty }}{lim}P\left(z\right)=0,\underset{z\to -\mathrm{\infty }}{lim}P\left(z\right)=+\mathrm{\infty },\text{and}{P}^{\prime }\left(z\right)<0.$ Thus, Equation (21(21) $P\left(z\right):=\varsigma {\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\nu \left(a\right){\mathrm{e}}^{-za}\mathrm{d}a+\eta \varsigma {\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\nu \left(a\right){\mathrm{e}}^{-za}\left({\int }_0^a\frac{\pi \left(s\right)}{\nu \left(s\right)}\mathrm{d}s\right)\mathrm{d}a=1.$(21) ) has a unique real root $z^{\ast }$. Noting that $P\left(0\right)={\mathcal{R}}_0$, we have $z^{\ast }>0$ if ${\mathcal{R}}_0>1$, and $z^{\ast }<0$ if ${\mathcal{R}}_0<1$. Set z=x+yi to be an arbitrary complex root to Equation (21(21) $P\left(z\right):=\varsigma {\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\nu \left(a\right){\mathrm{e}}^{-za}\mathrm{d}a+\eta \varsigma {\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\nu \left(a\right){\mathrm{e}}^{-za}\left({\int }_0^a\frac{\pi \left(s\right)}{\nu \left(s\right)}\mathrm{d}s\right)\mathrm{d}a=1.$(21) ). Then $1=|P\left(z\right)|=|P(x+yi)|⩽|P\left(x\right)|$, which implies that $z^{\ast }⩾x$. Thus, all roots of Equation (21(21) $P\left(z\right):=\varsigma {\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\nu \left(a\right){\mathrm{e}}^{-za}\mathrm{d}a+\eta \varsigma {\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\nu \left(a\right){\mathrm{e}}^{-za}\left({\int }_0^a\frac{\pi \left(s\right)}{\nu \left(s\right)}\mathrm{d}s\right)\mathrm{d}a=1.$(21) ) have negative part if and only if ${\mathcal{R}}_0<1$. This implies that the disease-free equilibrium $E_0$ is locally asymptotically stable if ${\mathcal{R}}_0<1$ and unstable if ${\mathcal{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(4a) $\begin{array}{rl}& \dot{S}\left(t\right)=\mathrm{\Lambda }-{\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[SI\right(t,a\left)\right]\mathrm{d}a-\eta S\left(t\right),\end{array}$(4a) ). 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(4a) $\begin{array}{rl}& \dot{S}\left(t\right)=\mathrm{\Lambda }-{\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[SI\right(t,a\left)\right]\mathrm{d}a-\eta S\left(t\right),\end{array}$(4a) ) is globally asymptotically stable if ${\mathcal{R}}_0<1$.

Proof.

Let $\beta \left(t\right)={\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[SI\right(t,a\left)\right]\mathrm{d}a$. Plugging Equation (8c(8c) $\begin{array}{rl}\left[SI\right(t,a\left)\right]& =\left\{\begin{array}{ll}\beta (t-a)\nu \left(a\right)\left\{\frac{\left[SS\right](t-a)}{S(t-a)}{\mathrm{e}}^{-{\int }_{t-a}^t\left(\beta \left(s\right)/S\left(s\right)\right)\mathrm{d}s}\right.& \\ \left.+{\int }_0^a\frac{\mathrm{\Lambda }\varsigma \pi \left(s\right)}{N(t-a+s)\nu \left(s\right)}{\mathrm{e}}^{-{\int }_{t-a+s}^t\left(\beta \left(w\right)/S\left(w\right)\right)\mathrm{d}w}\mathrm{d}s\right\},& t>a,\\ \frac{\nu \left(a\right)\psi (a-t)}{\nu (a-t)}{\mathrm{e}}^{-{\int }_0^t\left(\beta \left(s\right)/S\left(s\right)\right)\mathrm{d}s}& \\ +{\int }_0^t\frac{\mathrm{\Lambda }\varsigma \varphi (a-t)\pi (a-t+s)\nu \left(a\right)}{N\left(s\right)\pi (a-t)\nu (a-t+s)}{\mathrm{e}}^{-{\int }_s^t\left(\beta \left(w\right)/S\left(w\right)\right)\mathrm{d}w}\mathrm{d}s,& t⩽a,\end{array}\right.\end{array}$(8c) ) into the expression of $\beta \left(t\right)$ gives (22) $\begin{array}{rl}\beta \left(t\right)& ={\int }_0^t\lambda \left(a\right)\nu \left(a\right)\beta (t-a)\left(\frac{\left[SS\right](t-a)}{S(t-a)}{\mathrm{e}}^{-{\int }_{t-a}^t\left(\beta \left(s\right)/S\left(s\right)\right)\mathrm{d}s}\right.\\ & \left.+{\int }_0^a\frac{\mathrm{\Lambda }\varsigma \pi \left(s\right)}{N(t-a+s)\nu \left(s\right)}{\mathrm{e}}^{-{\int }_{t-a+s}^t\left(\beta \left(\omega \right)/S\left(\omega \right)\right)\mathrm{d}\omega }\mathrm{d}s\right)\mathrm{d}a\\ & +{\int }_t^{\mathrm{\infty }}\lambda \left(a\right)\nu \left(a\right)\left(\frac{\psi (a-t)}{\nu (a-t)}{\mathrm{e}}^{-{\int }_0^t\left(\beta \left(s\right)/S\left(s\right)\right)\mathrm{d}s}\right.\\ & \left.+{\int }_0^t\frac{\mathrm{\Lambda }\varsigma \varphi (a-t)\pi (a-t+s)}{N\left(s\right)\pi (a-t)\nu (a-t+s)}{\mathrm{e}}^{-{\int }_s^t\left(\beta \left(\omega \right)/S\left(\omega \right)\right)\mathrm{d}\omega }\mathrm{d}s\right)\mathrm{d}a\\ & ⩽{\int }_0^t\lambda \left(a\right)\nu \left(a\right)\beta (t-a)\left(\frac{\left[SS\right](t-a)}{S(t-a)}+{\int }_0^a\frac{\mathrm{\Lambda }\varsigma \pi \left(s\right)}{N(t-a+s)\nu \left(s\right)}\mathrm{d}s\right)\mathrm{d}a\\ & +{\int }_t^{\mathrm{\infty }}\lambda \left(a\right)\nu \left(a\right)\left(\frac{\psi (a-t)}{\nu (a-t)}+{\int }_0^t\frac{\mathrm{\Lambda }\varsigma \varphi (a-t)\pi (a-t+s)}{N\left(s\right)\pi (a-t)\nu (a-t+s)}\mathrm{d}s\right)\mathrm{d}a.\end{array}$(22) Let $x\left(t\right)=\left[SS\right]\left(t\right)/S\left(t\right)$. By Equations (4a(4a) $\begin{array}{rl}& \dot{S}\left(t\right)=\mathrm{\Lambda }-{\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[SI\right(t,a\left)\right]\mathrm{d}a-\eta S\left(t\right),\end{array}$(4a) ) and (4d(4d) $\begin{array}{rl}& \dot{\left[SS\right]}\left(t\right)=-2\frac{\left[SS\right]\left(t\right)}{S\left(t\right)}{\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[SI\right(t,a\left)\right]\mathrm{d}a-2\eta \left[SS\right]\left(t\right)+2\mathrm{\Lambda }\varsigma \frac{S\left(t\right)}{N\left(t\right)},\end{array}$(4d) ), we have (23) $\begin{array}{rl}\dot{x}\left(t\right)& =-\frac{\beta \left(t\right)}{S\left(t\right)}x\left(t\right)-\eta x\left(t\right)+2\mathrm{\Lambda }\frac{\varsigma }{N\left(t\right)}-\frac{\mathrm{\Lambda }}{S\left(t\right)}x\left(t\right)\\ & ⩽-\eta x\left(t\right)+2\mathrm{\Lambda }\frac{\varsigma }{N\left(t\right)}-\frac{\mathrm{\Lambda }}{S\left(t\right)}x\left(t\right)\\ & ⩽-2\eta x\left(t\right)+2\mathrm{\Lambda }\frac{\varsigma }{N\left(t\right)}.\end{array}$(23) Here, we have used the inequality $\underset{t\to \mathrm{\infty }}{lim}S\left(t\right)⩽\mathrm{\Lambda }/\eta .$ Combining Equation (23(23) $\begin{array}{rl}\dot{x}\left(t\right)& =-\frac{\beta \left(t\right)}{S\left(t\right)}x\left(t\right)-\eta x\left(t\right)+2\mathrm{\Lambda }\frac{\varsigma }{N\left(t\right)}-\frac{\mathrm{\Lambda }}{S\left(t\right)}x\left(t\right)\\ & ⩽-\eta x\left(t\right)+2\mathrm{\Lambda }\frac{\varsigma }{N\left(t\right)}-\frac{\mathrm{\Lambda }}{S\left(t\right)}x\left(t\right)\\ & ⩽-2\eta x\left(t\right)+2\mathrm{\Lambda }\frac{\varsigma }{N\left(t\right)}.\end{array}$(23) ) with $\underset{t\to \mathrm{\infty }}{lim}N\left(t\right)=\mathrm{\Lambda }/\eta$, we arrive at $\underset{t\to \mathrm{\infty }}{lim}x\left(t\right)⩽\varsigma .$ Using this result and Fatou's lemma, and taking the limit supremum when $t\to \mathrm{\infty }$ on both sides of Equation (22(22) $\begin{array}{rl}\beta \left(t\right)& ={\int }_0^t\lambda \left(a\right)\nu \left(a\right)\beta (t-a)\left(\frac{\left[SS\right](t-a)}{S(t-a)}{\mathrm{e}}^{-{\int }_{t-a}^t\left(\beta \left(s\right)/S\left(s\right)\right)\mathrm{d}s}\right.\\ & \left.+{\int }_0^a\frac{\mathrm{\Lambda }\varsigma \pi \left(s\right)}{N(t-a+s)\nu \left(s\right)}{\mathrm{e}}^{-{\int }_{t-a+s}^t\left(\beta \left(\omega \right)/S\left(\omega \right)\right)\mathrm{d}\omega }\mathrm{d}s\right)\mathrm{d}a\\ & +{\int }_t^{\mathrm{\infty }}\lambda \left(a\right)\nu \left(a\right)\left(\frac{\psi (a-t)}{\nu (a-t)}{\mathrm{e}}^{-{\int }_0^t\left(\beta \left(s\right)/S\left(s\right)\right)\mathrm{d}s}\right.\\ & \left.+{\int }_0^t\frac{\mathrm{\Lambda }\varsigma \varphi (a-t)\pi (a-t+s)}{N\left(s\right)\pi (a-t)\nu (a-t+s)}{\mathrm{e}}^{-{\int }_s^t\left(\beta \left(\omega \right)/S\left(\omega \right)\right)\mathrm{d}\omega }\mathrm{d}s\right)\mathrm{d}a\\ & ⩽{\int }_0^t\lambda \left(a\right)\nu \left(a\right)\beta (t-a)\left(\frac{\left[SS\right](t-a)}{S(t-a)}+{\int }_0^a\frac{\mathrm{\Lambda }\varsigma \pi \left(s\right)}{N(t-a+s)\nu \left(s\right)}\mathrm{d}s\right)\mathrm{d}a\\ & +{\int }_t^{\mathrm{\infty }}\lambda \left(a\right)\nu \left(a\right)\left(\frac{\psi (a-t)}{\nu (a-t)}+{\int }_0^t\frac{\mathrm{\Lambda }\varsigma \varphi (a-t)\pi (a-t+s)}{N\left(s\right)\pi (a-t)\nu (a-t+s)}\mathrm{d}s\right)\mathrm{d}a.\end{array}$(22) ), we have (24) $\underset{t\to \mathrm{\infty }}{lim sup}\beta \left(t\right)⩽{\mathcal{R}}_0\underset{t\to \mathrm{\infty }}{lim sup}\beta \left(t\right).$(24) As it is assumed that ${\mathcal{R}}_0<1$, the only condition that inequality (24(24) $\underset{t\to \mathrm{\infty }}{lim sup}\beta \left(t\right)⩽{\mathcal{R}}_0\underset{t\to \mathrm{\infty }}{lim sup}\beta \left(t\right).$(24) ) holds is $\underset{t\to \mathrm{\infty }}{lim sup}\beta \left(t\right)=0.$ From Equations (8b(8b) $\begin{array}{rl}I(t,a)& =\left\{\begin{array}{ll}\beta (t-a)\pi \left(a\right),& t>a,\\ \varphi (a-t)\frac{\pi \left(a\right)}{\pi (a-t)},& t⩽a,\end{array}\right.\end{array}$(8b) ) and (8c(8c) $\begin{array}{rl}\left[SI\right(t,a\left)\right]& =\left\{\begin{array}{ll}\beta (t-a)\nu \left(a\right)\left\{\frac{\left[SS\right](t-a)}{S(t-a)}{\mathrm{e}}^{-{\int }_{t-a}^t\left(\beta \left(s\right)/S\left(s\right)\right)\mathrm{d}s}\right.& \\ \left.+{\int }_0^a\frac{\mathrm{\Lambda }\varsigma \pi \left(s\right)}{N(t-a+s)\nu \left(s\right)}{\mathrm{e}}^{-{\int }_{t-a+s}^t\left(\beta \left(w\right)/S\left(w\right)\right)\mathrm{d}w}\mathrm{d}s\right\},& t>a,\\ \frac{\nu \left(a\right)\psi (a-t)}{\nu (a-t)}{\mathrm{e}}^{-{\int }_0^t\left(\beta \left(s\right)/S\left(s\right)\right)\mathrm{d}s}& \\ +{\int }_0^t\frac{\mathrm{\Lambda }\varsigma \varphi (a-t)\pi (a-t+s)\nu \left(a\right)}{N\left(s\right)\pi (a-t)\nu (a-t+s)}{\mathrm{e}}^{-{\int }_s^t\left(\beta \left(w\right)/S\left(w\right)\right)\mathrm{d}w}\mathrm{d}s,& t⩽a,\end{array}\right.\end{array}$(8c) ), we obtain that $\underset{t\to \mathrm{\infty }}{lim sup}I(t,a)=0\text{and}\underset{t\to \mathrm{\infty }}{lim sup}\left[SI\right(t,a\left)\right]=0.$ Furthermore, solving Equation (4a(4a) $\begin{array}{rl}& \dot{S}\left(t\right)=\mathrm{\Lambda }-{\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[SI\right(t,a\left)\right]\mathrm{d}a-\eta S\left(t\right),\end{array}$(4a) ) gives $S\left(t\right)={\mathrm{e}}^{-\eta t}\left(S_0+{\int }_0^t(\mathrm{\Lambda }-\beta (s\left)\right){\mathrm{e}}^{\eta s}\mathrm{d}s\right).$ Thus $\underset{t\to \mathrm{\infty }}{lim sup}S\left(t\right)=\frac{\mathrm{\Lambda }}{\eta }.$ Similarly, it has $\underset{t\to \mathrm{\infty }}{lim sup}\left[SS\right]\left(t\right)=\frac{\mathrm{\Lambda }}{\eta }\varsigma .$

Accordingly, the equilibrium $E_0$ is globally asymptotically stable if ${\mathcal{R}}_0<1$. This completes the proof of Theorem 3.2.

4. 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\approx (\varsigma /N_0)S_0^2$ and $\psi \left(a\right)\approx \left(\varsigma /N_0\right)S_0\varphi \left(a\right)$, which mean that the initial network is homogeneous and newborns have the same average degree with the initial network.

4.1. Markovian transmission and non-Markovian recovery

We consider a Markovian transmission process with rate ${\lambda }_0$ (i.e. $\lambda \left(a\right)\equiv {\lambda }_0$) and a general recovery process. In the simulations, the recovery rate is determined by the recovery time distribution $f_R\left(t\right)$. Specifically, the instantaneous recovery rate at infection age a is given by the hazard function of $f_R\left(a\right)$, i.e. $\gamma \left(a\right)=\frac{f_R\left(a\right)}{1-F_R\left(a\right)},$ where $F_R\left(a\right)={\int }_0^a{f}_R\left(t\right)\mathrm{d}t$ is the cumulative distribution function of the recovery time [26, 27]. Particularly, if the recovery time is exponentially distributed with rate ${\gamma }_0$, it corresponds to a Markovian recovery process. In this case, $\gamma \left(a\right)\equiv {\gamma }_0$, and the basic reproduction number ${\mathcal{R}}_{0,1}$ of model (4(4a) $\begin{array}{rl}& \dot{S}\left(t\right)=\mathrm{\Lambda }-{\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[SI\right(t,a\left)\right]\mathrm{d}a-\eta S\left(t\right),\end{array}$(4a) ) is (25) ${\mathcal{R}}_{0,1}=\frac{{\lambda }_0^2\varsigma }{({\lambda }_0+\eta )({\lambda }_0+{\gamma }_0+2\eta )}+\frac{\eta {\lambda }_0\varsigma }{({\lambda }_0+\eta )({\gamma }_0+\eta )}.$(25) Next, we study three other recovery time distributions: Weibull distribution, Gamma distribution, and uniform distribution [29, 35].

1. If the recovery time distribution is Weibullean with shape ${\alpha }_1$ and scale ${\beta }_1$ (denoted as $W({\alpha }_1,{\beta }_1)$), i.e. $f_{R,1}\left(t\right)=\left\{\begin{array}{ll}\frac{{\alpha }_1}{{\beta }_1}{\left(\frac{t}{{\beta }_1}\right)}^{{\alpha }_1-1}{\mathrm{e}}^{-(t/{\beta }_1{)}^{{\alpha }_1}},& t⩾0,\\ 0,& t<0,\end{array}\right.$ the recovery rate is (26) $\gamma \left(a\right)=\frac{{\alpha }_1}{{\beta }_1}{\left(\frac{a}{{\beta }_1}\right)}^{{\alpha }_1-1}\text{for }a⩾0.$(26) When ${\alpha }_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 ${\alpha }_2$ and scale ${\beta }_2$ (denoted as $G({\alpha }_2,{\beta }_2)$), i.e. $f_{R,2}\left(t\right)=\left\{\begin{array}{ll}\frac{t^{{\alpha }_2-1}{\mathrm{e}}^{-t/{\beta }_2}}{{\beta }_2^{{\alpha }_2}\mathrm{\Gamma }\left({\alpha }_2\right)},& t⩾0,\\ 0,& t<0,\end{array}\right.$ where $\mathrm{\Gamma }\left({\alpha }_2\right)={\int }_0^{\mathrm{\infty }}{t}^{{\alpha }_2-1}{\mathrm{e}}^{-t}\mathrm{d}t$, the recovery rate is (27) $\gamma \left(a\right)=\frac{a^{{\alpha }_2-1}{\mathrm{e}}^{-a/{\beta }_2}}{{\beta }_2^{{\alpha }_2}\left(\mathrm{\Gamma }\right({\alpha }_2)-{\int }_0^{a/{\beta }_2}{t}^{{\alpha }_2-1}{\mathrm{e}}^{-t}\mathrm{d}t)}\text{for }a⩾0.$(27) When ${\alpha }_2=1$, the Gamma distribution reduces to an exponential distribution, corresponding to a Markovian recovery process.

3. If the recovery time has a uniform distribution on the interval $[{\alpha }_3,{\beta }_3]$ (denoted as $U({\alpha }_3,{\beta }_3)$), i.e. $f_{R,3}\left(t\right)=\left\{\begin{array}{ll}\frac{1}{{\beta }_3-{\alpha }_3}& \text{if}{\alpha }_3⩽t⩽{\beta }_3,\\ 0& \text{otherwise},\end{array}\right.$ the recovery rate is (28) $\gamma \left(a\right)=\frac{1}{{\beta }_3-a}\text{for }a\in [{\alpha }_3,{\beta }_3].$(28)

In order to compare these distributions, we fix the average recovery time to $1/{\gamma }_0$. Then ${\beta }_1={\left[\mathrm{\Gamma }\left(1+\frac{1}{{\alpha }_1}\right){\gamma }_0\right]}^{-1},{\beta }_2=\frac{1}{{\alpha }_2{\gamma }_0},{\beta }_3=\frac{2}{{\gamma }_0}-{\alpha }_3.$

Figure 1 shows the density of infectious nodes ($I(t,a)$) as a function of time and infection age from model (4(4a) $\begin{array}{rl}& \dot{S}\left(t\right)=\mathrm{\Lambda }-{\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[SI\right(t,a\left)\right]\mathrm{d}a-\eta S\left(t\right),\end{array}$(4a) ) 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 ${\alpha }_1$ (or ${\alpha }_2$, ${\alpha }_3$) results in a larger value of the first peak. Besides, Figure 3 shows that decreasing ${\alpha }_1$ (or ${\alpha }_2$, ${\alpha }_3$) leads to a larger transmission threshold ${\lambda }_c$ (above which the epidemic outbreaks) and a smaller equilibrium prevalence $I^{\ast }/N^{\ast }$, where $I^{\ast }={\int }_0^{\mathrm{\infty }}{I}^{\ast }\left(a\right)\mathrm{d}a$. Therefore, decreasing ${\alpha }_1$ (or ${\alpha }_2$, ${\alpha }_3$) is beneficial to the control of diseases. It is worth noting that the variances of the recovery time distributions satisfy the following inequalities: $\begin{array}{rl}& V_{W(1/2,1/\mathrm{\Gamma }(3\left)\right)}>V_{W(3/4,1/\mathrm{\Gamma }(7/3\left)\right)}>V_{W(1,1/\mathrm{\Gamma }(2\left)\right)}>V_{W(5/4,1/\mathrm{\Gamma }(9/5\left)\right)}>V_{W(3/2,1/\mathrm{\Gamma }(5/3\left)\right)},\\ & V_{G(1/2,2)}>V_{G(3/4,4/3)}>V_{G(1,1)}>V_{G(5/4,4/5)}>V_{G(3/2,2/3)}>V_{G(7/4,4/7)}>V_{G(2,1/2)},\end{array}$ and $V_{U(0,2)}>V_{U(0.25,1.75)}>V_{U(0.5,1.5)}>V_{U(0.75,1.25)},$ 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.

Figure 1. The phase diagrams showing the density of infectious nodes as a function of time and infection age ($I(t,a)$) from model (4(4a) $\begin{array}{rl}& \dot{S}\left(t\right)=\mathrm{\Lambda }-{\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[SI\right(t,a\left)\right]\mathrm{d}a-\eta S\left(t\right),\end{array}$(4a) ) with different recovery time distributions: (a) Weibull distribution $W(\frac{3}{2},1/\mathrm{\Gamma }(\frac{5}{3}\left)\right)$, (b) Gamma distribution $G(\frac{1}{2},2)$, and (c) Uniform distribution $U(0.75,1.25)$. The parameters are: $\eta =0.2$, ${\lambda }_0=0.3$, ${\gamma }_0=1$, $\mathrm{\Lambda }=2000$, $\varsigma =8$. The initial values are: $N\left(0\right)=10,000$, $S\left(0\right)=9900$, $\left[SS\right]\left(0\right)=\varsigma S\left(0\right)$, $I\left(0\right)=100$, $I(0,a)=I\left(0\right)\overline{\varphi }\left(a\right)$, $\left[SI\right(0,a\left)\right]=\varsigma I(0,a)$, where $\overline{\varphi }\left(a\right)$ is the value of $\overline{\varphi }\left(x\right)$ at a, and $\overline{\varphi }\left(x\right)$ denotes the uniform distribution on interval $[0,K]$. K is the maximal infection age.

Figure 2. The time dependence of the total number of infectious nodes from model (4(4a) $\begin{array}{rl}& \dot{S}\left(t\right)=\mathrm{\Lambda }-{\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[SI\right(t,a\left)\right]\mathrm{d}a-\eta S\left(t\right),\end{array}$(4a) ) with different recovery time distributions: (a) Weibull distribution $W({\alpha }_1,{\beta }_1)$, (b) Gamma distribution $G({\alpha }_2,{\beta }_2),$ and (c) Uniform distribution $U({\alpha }_3,{\beta }_3)$. Other parameters and initial values are the same as those in Figure 1.

Figure 3. Equilibrium prevalence as a function of the transmission rate from model (4(4a) $\begin{array}{rl}& \dot{S}\left(t\right)=\mathrm{\Lambda }-{\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[SI\right(t,a\left)\right]\mathrm{d}a-\eta S\left(t\right),\end{array}$(4a) ) with different recovery time distributions: (a) Weibull distribution $W({\alpha }_1,{\beta }_1)$, (b) Gamma distribution $G({\alpha }_2,{\beta }_2)$, and (c) Uniform distribution $U({\alpha }_3,{\beta }_3)$. The initial values are: $N\left(0\right)=\mathrm{10,000}$, $S\left(0\right)=9900$, $\left[SS\right]\left(0\right)=\varsigma S\left(0\right)$, $I\left(0\right)=100$, $I(0,a)=I\left(0\right)\overline{\varphi }\left(a\right)$, $\left[SI\right(0,a\left)\right]=\varsigma I(0,a)$, where $\overline{\varphi }\left(a\right)$ is the value of $\overline{\varphi }\left(x\right)$ at a, and $\overline{\varphi }\left(x\right)$ denotes the uniform distribution on interval $[0,K]$. K is the maximal infection age. The parameters are: $\eta =0.2$, ${\lambda }_0=0.3$, ${\gamma }_0=1$, $\mathrm{\Lambda }=2000$, $\varsigma =8$.

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.

Figure 4. The time dependence of the total number of infectious nodes from model (4(4a) $\begin{array}{rl}& \dot{S}\left(t\right)=\mathrm{\Lambda }-{\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[SI\right(t,a\left)\right]\mathrm{d}a-\eta S\left(t\right),\end{array}$(4a) ) with different recovery time distributions. The parameters of the distributions are set to have the same mean and variance. The mean of these distributions is $1/{\gamma }_0=0.9$. The variances are: (a) 0.3734 (solid), 0.2213 (dashed), 0.1483 (dashdot), and (b) 1.7323 (solid), 1.2870 (dashed), 1.0035 (dashdot). The parameters are: $\eta =0.2$, ${\lambda }_0=0.2$, $\mathrm{\Lambda }=2000$, $\varsigma =8$. The initial values are: $N\left(0\right)=\mathrm{10,000}$, $S\left(0\right)=N\left(0\right)-I\left(0\right)$, $\left[SS\right]\left(0\right)=\varsigma S\left(0\right)$, $I(0,a)=I\left(0\right)\overline{\varphi }\left(a\right)$, $\left[SI\right(0,a\left)\right]=\varsigma I(0,a)$, where $\overline{\varphi }\left(a\right)$ is the value of $\overline{\varphi }\left(x\right)$ at a, and $\overline{\varphi }\left(x\right)$ denotes the uniform distribution on interval $[0,K]$. K is the maximal infection age.

Finally, Figure 5 suggests that the endemic equilibrium is globally asymptotically stable.

Figure 5. The time dependence of the total number of infectious nodes from model (4(4a) $\begin{array}{rl}& \dot{S}\left(t\right)=\mathrm{\Lambda }-{\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[SI\right(t,a\left)\right]\mathrm{d}a-\eta S\left(t\right),\end{array}$(4a) ) with different recovery time distributions and initial values: (a) Weibull distribution $W(\frac{1}{2},\frac{1}{\mathrm{\Gamma }\left(3\right)})$, (b) Gamma distribution $G(\frac{1}{2},2)$, and (c) Uniform distribution $U(0,2)$. The parameters are: $\eta =0.2$, ${\lambda }_0=0.3$, ${\gamma }_0=1$, $\mathrm{\Lambda }=2000$, $\varsigma =8$. The initial values are: $N\left(0\right)=\mathrm{10,000}$, $S\left(0\right)=N\left(0\right)-I\left(0\right)$, $\left[SS\right]\left(0\right)=\varsigma S\left(0\right)$, $I(0,a)=I\left(0\right)\overline{\varphi }\left(a\right)$, $\left[SI\right(0,a\left)\right]=\varsigma I(0,a)$, where $\overline{\varphi }\left(a\right)$ is the value of $\overline{\varphi }\left(x\right)$ at a, and $\overline{\varphi }\left(x\right)$ denotes the uniform distribution on interval $[0,K]$. K is the maximal infection age.

4.2. Non-Markovian transmission and Markovian recovery

We consider a Markovian recovery process with rate ${\gamma }_0$ (i.e. $\gamma \left(a\right)\equiv {\gamma }_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\left(t\right)$, the transmission rate $\lambda \left(a\right)$ is given by the hazard function of $f_I\left(t\right)$, i.e. $\lambda \left(a\right)=\frac{f_I\left(a\right)}{1-F_I\left(a\right)},$ where $F_I\left(a\right)={\int }_0^a{f}_I\left(t\right)\mathrm{d}t$ is the cumulative distribution function of the infection time [11, 27]. Here, we study three different infection time distributions: Weibull distribution $W({\alpha }_4,{\beta }_4)$, Gamma distribution $G({\alpha }_5,{\beta }_5)$, and Uniform distribution $U({\alpha }_6,{\beta }_6)$. In order to compare these distributions, we fix the average infection time to $1/{\lambda }_0$. Then it has ${\beta }_4={\left[\mathrm{\Gamma }\left(1+\frac{1}{{\alpha }_4}\right){\lambda }_0\right]}^{-1},{\beta }_5=\frac{1}{{\alpha }_5{\lambda }_0},{\beta }_6=\frac{2}{{\lambda }_0}-{\alpha }_6.$

Figure 6 plots the time dependence of the total number of infectious nodes from model (4(4a) $\begin{array}{rl}& \dot{S}\left(t\right)=\mathrm{\Lambda }-{\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[SI\right(t,a\left)\right]\mathrm{d}a-\eta S\left(t\right),\end{array}$(4a) ) with different infection time distributions. It shows that the prevalence is quite sensitive to the variation in the infection time. Moreover, the variances of the infection time distributions satisfy the following inequalities: $\begin{array}{rl}& V_{W(1/2,10/(3\mathrm{\Gamma }\left(3\right)\left)\right)}>V_{W(3/4,10/(3\mathrm{\Gamma }\left(7/3\right)\left)\right)}>V_{W(1,10/(3\mathrm{\Gamma }\left(2\right)\left)\right)}\\ & >V_{W(5/4,10/(3\mathrm{\Gamma }\left(9/5\right)\left)\right)}>V_{W(3/2,10/3\mathrm{\Gamma }(5/3\left)\right)},\\ & V_{G(1/2,20/3)}>V_{G(3/4,40/9)}>V_{G(1,10/3)}>V_{G(5/4,8/3)}\\ & >V_{G(3/2,20/9)}>V_{G(7/4,40/21)}>V_{G(2,5/3)},\end{array}$ and $V_{U(0,20/3)}>V_{U(1/4,77/12)}>V_{U(1/2,37/6)}>V_{U(3/4,71/12)}.$ 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 6. The time dependence of the total number of infectious nodes from model (4(4a) $\begin{array}{rl}& \dot{S}\left(t\right)=\mathrm{\Lambda }-{\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[SI\right(t,a\left)\right]\mathrm{d}a-\eta S\left(t\right),\end{array}$(4a) ) with different infection time distributions: (a) Weibull distribution $W({\alpha }_4,{\beta }_4)$, (b) Gamma distribution $G({\alpha }_5,{\beta }_5)$, and (c) Uniform distribution $U({\alpha }_6,{\beta }_6)$. The parameters are: $\eta =0.2$, ${\lambda }_0=0.3$, ${\gamma }_0=1$, $\mathrm{\Lambda }=2000$, $\varsigma =8$. The initial values are: $N\left(0\right)=\mathrm{10,000}$, $S\left(0\right)=9900$, $\left[SS\right]\left(0\right)=\varsigma S\left(0\right)$, $I\left(0\right)=100$, $I(0,a)=I\left(0\right)\overline{\varphi }\left(a\right)$, $\left[SI\right(0,a\left)\right]=\varsigma I(0,a)$, where $\overline{\varphi }\left(a\right)$ is the value of $\overline{\varphi }\left(x\right)$ at a, and $\overline{\varphi }\left(x\right)$ denotes the uniform distribution on interval $[0,K]$. K is the maximal infection age.

4.3. The effects of demography on the epidemic dynamics

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.

Figure 7. The time dependence of the total number of infectious nodes from model (4(4a) $\begin{array}{rl}& \dot{S}\left(t\right)=\mathrm{\Lambda }-{\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[SI\right(t,a\left)\right]\mathrm{d}a-\eta S\left(t\right),\end{array}$(4a) ) with different values of η. The initial values are: $N\left(0\right)=\mathrm{10,000}$, $S\left(0\right)=9900$, $\left[SS\right]\left(0\right)=\varsigma S\left(0\right)$, $I\left(0\right)=100$, $I(0,a)=I\left(0\right)\overline{\varphi }\left(a\right)$, $\left[SI\right(0,a\left)\right]=\varsigma I(0,a)$, where $\overline{\varphi }\left(a\right)$ is the value of $\overline{\varphi }\left(x\right)$ at a, and $\overline{\varphi }\left(x\right)$ denotes the uniform distribution on interval $[0,K]$. K is the maximal infection age. The parameters are: $\gamma \left(a\right)=1$ for $a\in [0,10]$ and 0 otherwise, $\lambda \left(a\right)=0.2$ for $a\in [0,10]$ and 0 otherwise, $\mathrm{\Lambda }=N\left(0\right)/\eta$, $\varsigma =8$.

5. 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 ${\mathcal{R}}_0>1$. Moreover, the disease-free equilibrium is globally asymptotically stable if ${\mathcal{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(4a) $\begin{array}{rl}& \dot{S}\left(t\right)=\mathrm{\Lambda }-{\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[SI\right(t,a\left)\right]\mathrm{d}a-\eta S\left(t\right),\end{array}$(4a) ) 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(4a) $\begin{array}{rl}& \dot{S}\left(t\right)=\mathrm{\Lambda }-{\int }_0^{\mathrm{\infty }}\lambda \left(a\right)\left[SI\right(t,a\left)\right]\mathrm{d}a-\eta S\left(t\right),\end{array}$(4a) ). 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.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This work was supported by National Natural Science Foundation of China [grant numbers 11331009, 11501340, 11601294, and 11701348].

Appendix

In this appendix, we give the derivation of $\left[SI\right(t,a\left)\right]$. Let us consider the following first-order PDE: (A1) $\left(\frac{\mathrm{\partial }}{\mathrm{\partial }t}+\frac{\mathrm{\partial }}{\mathrm{\partial }a}\right)x(t,a)=-f\left(t\right)x(t,a)-g\left(a\right)x(t,a)+h(t,a)$(A1) with boundary condition $x(t,0)=\phi \left(t\right)$ and initial condition $x(0,a)$. In the method of characteristics, it is assumed that the partial differential equation can be expressed as a system of ordinary differential equations along the characteristic curves. The characteristic curves are found by solving the following ordinary differential equations: $\frac{\mathrm{d}t}{\mathrm{d}\tau }=1,\frac{\mathrm{d}a}{\mathrm{d}\tau }=1,\frac{\mathrm{d}y\left(\tau \right)}{\mathrm{d}\tau }=-f\left(t\right(\tau \left)\right)y\left(\tau \right)-g\left(a\right(\tau \left)\right)y\left(\tau \right)+h\left(t\right(\tau ),a(\tau \left)\right),$ with initial conditions $t\left(0\right)=t_0,a\left(0\right)=a_0,y\left(0\right)=x(t_0,a_0),$ where $y\left(\tau \right)=x\left(t\right(\tau ),a(\tau \left)\right)$. The solution of the system of ordinary differential equations along the characteristic curves is $\begin{array}{rl}& t\left(\tau \right)=t_0+\tau ,a\left(\tau \right)=a_0+\tau ,\\ & y\left(\tau \right)=y\left(0\right){\mathrm{e}}^{-{\int }_0^{\tau }\left[f\right(t\left(s\right))+g(a\left(s\right)\left)\right]\mathrm{d}s}+{\int }_0^{\tau }h\left(t\right(s),a(s\left)\right){\mathrm{e}}^{-{\int }_s^{\tau }\left[f\right(t\left(w\right))+g(a\left(w\right)\left)\right]\mathrm{d}w}\mathrm{d}s.\end{array}$ Furthermore, it has (A2) $\begin{array}{rl}x(t,a)& =x(t_0,a_0){\mathrm{e}}^{-{\int }_0^{\tau }\left[f\right(t\left(s\right))+g(a\left(s\right)\left)\right]\mathrm{d}s}+{\int }_0^{\tau }h\left(t\right(s),a(s\left)\right){\mathrm{e}}^{-{\int }_s^{\tau }\left[f\right(t\left(w\right))+g(a\left(w\right)\left)\right]\mathrm{d}w}\mathrm{d}s\\ & =x(t_0,a_0){\mathrm{e}}^{-{\int }_0^{\tau }\left[f\right(t_0+s)+g(a_0+s\left)\right]\mathrm{d}s}+{\int }_0^{\tau }h(t_0+s,a_0+s){\mathrm{e}}^{-{\int }_s^{\tau }\left[f\right(t_0+w)+g(a_0+w\left)\right]\mathrm{d}w}\mathrm{d}s.\end{array}$(A2) If t>a, set $a_0=0$. Then $\tau =a$ and $t_0=t-a$. Substituting these expressions into Equation (A.2(A2) $\begin{array}{rl}x(t,a)& =x(t_0,a_0){\mathrm{e}}^{-{\int }_0^{\tau }\left[f\right(t\left(s\right))+g(a\left(s\right)\left)\right]\mathrm{d}s}+{\int }_0^{\tau }h\left(t\right(s),a(s\left)\right){\mathrm{e}}^{-{\int }_s^{\tau }\left[f\right(t\left(w\right))+g(a\left(w\right)\left)\right]\mathrm{d}w}\mathrm{d}s\\ & =x(t_0,a_0){\mathrm{e}}^{-{\int }_0^{\tau }\left[f\right(t_0+s)+g(a_0+s\left)\right]\mathrm{d}s}+{\int }_0^{\tau }h(t_0+s,a_0+s){\mathrm{e}}^{-{\int }_s^{\tau }\left[f\right(t_0+w)+g(a_0+w\left)\right]\mathrm{d}w}\mathrm{d}s.\end{array}$(A2) ), we obtain the solution when t>a (Ae) $\begin{array}{rl}x(t,a)& =x(t-a,0){\mathrm{e}}^{-{\int }_0^a\left[f\right(t-a+s)+g(s\left)\right]\mathrm{d}s}+{\int }_0^{a}h(t-a+s,s){\mathrm{e}}^{-{\int }_s^a\left[f\right(t-a+w)+g(w\left)\right]\mathrm{d}w}\mathrm{d}s\\ & =x(t-a,0){\mathrm{e}}^{-{\int }_{t-a}^{t}f\left(s\right)\mathrm{d}s-{\int }_0^{a}g\left(s\right)\mathrm{d}s}+{\int }_0^{a}h(t-a+s,s){\mathrm{e}}^{-{\int }_s^{a}g\left(w\right)\mathrm{d}w-{\int }_{t-a+s}^{t}f\left(w\right)\mathrm{d}w}\mathrm{d}s.\end{array}$(Ae) If $t⩽a$, set $t_0=0$. Then $\tau =t$ and $a_0=a-t$. Substituting these expressions into Equation (A.2(A2) $\begin{array}{rl}x(t,a)& =x(t_0,a_0){\mathrm{e}}^{-{\int }_0^{\tau }\left[f\right(t\left(s\right))+g(a\left(s\right)\left)\right]\mathrm{d}s}+{\int }_0^{\tau }h\left(t\right(s),a(s\left)\right){\mathrm{e}}^{-{\int }_s^{\tau }\left[f\right(t\left(w\right))+g(a\left(w\right)\left)\right]\mathrm{d}w}\mathrm{d}s\\ & =x(t_0,a_0){\mathrm{e}}^{-{\int }_0^{\tau }\left[f\right(t_0+s)+g(a_0+s\left)\right]\mathrm{d}s}+{\int }_0^{\tau }h(t_0+s,a_0+s){\mathrm{e}}^{-{\int }_s^{\tau }\left[f\right(t_0+w)+g(a_0+w\left)\right]\mathrm{d}w}\mathrm{d}s.\end{array}$(A2) ), we obtain the solution when $t⩽a$ (A4) $\begin{array}{rl}x(t,a)& =x(0,a-t)e^{-{\int }_0^t\left[f\right(s)+g(a-t+s\left)\right]\mathrm{d}s}+{\int }_0^{t}h(s,a-t+s){\mathrm{e}}^{-{\int }_s^t\left[f\right(w)+g(a-t+w\left)\right]\mathrm{d}w}\mathrm{d}s\\ & =x(0,a-t){\mathrm{e}}^{-{\int }_0^{t}f\left(s\right)\mathrm{d}s-{\int }_{a-t}^{a}g\left(s\right)\mathrm{d}s}+{\int }_0^{t}h(s,a-t+s){\mathrm{e}}^{-{\int }_s^{t}f\left(w\right)\mathrm{d}w-{\int }_{a-t+s}^{a}g\left(w\right)\mathrm{d}w}\mathrm{d}s.\end{array}$(A4) So finally, the solution of Equation (A.1(A1) $\left(\frac{\mathrm{\partial }}{\mathrm{\partial }t}+\frac{\mathrm{\partial }}{\mathrm{\partial }a}\right)x(t,a)=-f\left(t\right)x(t,a)-g\left(a\right)x(t,a)+h(t,a)$(A1) ) can be written in the form (A5) $x(t,a)=\left\{\begin{array}{ll}x(t-a,0){\mathrm{e}}^{-{\int }_{t-a}^{t}f\left(s\right)\mathrm{d}s-{\int }_0^{a}g\left(s\right)\mathrm{d}s}& \\ +{\int }_0^{a}h(t-a+s,s){\mathrm{e}}^{-{\int }_s^{a}g\left(w\right)\mathrm{d}w-{\int }_{t-a+s}^{t}f\left(w\right)\mathrm{d}w}\mathrm{d}s,& t>a,\\ x(0,a-t){\mathrm{e}}^{-{\int }_0^{t}f\left(s\right)\mathrm{d}s-{\int }_{a-t}^{a}g\left(s\right)\mathrm{d}s}& \\ +{\int }_0^{t}h(s,a-t+s){\mathrm{e}}^{-{\int }_s^{t}f\left(w\right)\mathrm{d}w-{\int }_{a-t+s}^{a}g\left(w\right)\mathrm{d}w}\mathrm{d}s,& t⩽a.\end{array}\right.$(A5)

In our model, $x(t,a)=\left[SI\right(t,a\left)\right]$, $f\left(t\right)=\beta \left(t\right)/S\left(t\right)$, $g\left(a\right)=\lambda \left(a\right)+\gamma \left(a\right)+2\eta$, $x(t,0)=\left[SS\right]\left(t\right)/S\left(t\right)\beta \left(t\right)$, $x(0,a)=\psi \left(a\right)$ and $h(t,a)=\left\{\begin{array}{ll}\frac{\mathrm{\Lambda }\varsigma }{N\left(t\right)}\beta (t-a)\pi \left(a\right),& t>a,\\ \frac{\mathrm{\Lambda }\varsigma \varphi (a-t)}{N\left(t\right)}\frac{\pi \left(a\right)}{\pi (a-t)},& t⩽a.\end{array}\right.$ Substituting these expressions into Equation (A.5(A5) $x(t,a)=\left\{\begin{array}{ll}x(t-a,0){\mathrm{e}}^{-{\int }_{t-a}^{t}f\left(s\right)\mathrm{d}s-{\int }_0^{a}g\left(s\right)\mathrm{d}s}& \\ +{\int }_0^{a}h(t-a+s,s){\mathrm{e}}^{-{\int }_s^{a}g\left(w\right)\mathrm{d}w-{\int }_{t-a+s}^{t}f\left(w\right)\mathrm{d}w}\mathrm{d}s,& t>a,\\ x(0,a-t){\mathrm{e}}^{-{\int }_0^{t}f\left(s\right)\mathrm{d}s-{\int }_{a-t}^{a}g\left(s\right)\mathrm{d}s}& \\ +{\int }_0^{t}h(s,a-t+s){\mathrm{e}}^{-{\int }_s^{t}f\left(w\right)\mathrm{d}w-{\int }_{a-t+s}^{a}g\left(w\right)\mathrm{d}w}\mathrm{d}s,& t⩽a.\end{array}\right.$(A5) ), we obtain the solution of $\left[SI\right(t,a\left)\right]$ as follows: (A6) $\left[SI\right(t,a\left)\right]=\left\{\begin{array}{ll}\frac{\left[SS\right](t-a)}{S(t-a)}\beta (t-a)\nu \left(a\right){\mathrm{e}}^{-{\int }_{t-a}^t\left(\beta \left(s\right)/S\left(s\right)\right)\mathrm{d}s}& \\ +{\int }_0^a\frac{\mathrm{\Lambda }\varsigma \beta (t-a)\pi \left(s\right)\nu \left(a\right)}{N(t-a+s)\nu \left(s\right)}{\mathrm{e}}^{-{\int }_{t-a+s}^t\left(\beta \left(w\right)/S\left(w\right)\right)\mathrm{d}w}\mathrm{d}s,& t>a,\\ \frac{\nu \left(a\right)\psi (a-t)}{\nu (a-t)}{\mathrm{e}}^{-{\int }_0^t\left(\beta \left(s\right)/S\left(s\right)\right)\mathrm{d}s}\\ +{\int }_0^t\frac{\mathrm{\Lambda }\varsigma \varphi (a-t)\pi (a-t+s)\nu \left(a\right)}{N\left(s\right)\pi (a-t)\nu (a-t+s)}{\mathrm{e}}^{-{\int }_s^t\left(\beta \left(w\right)/S\left(w\right)\right)\mathrm{d}w}\mathrm{d}s,& t⩽a.\end{array}\right.$(A6)