A theta-scheme approximation of basic reproduction number for an age-structured epidemic system in a finite horizon

1.   Introduction

Many classical SIS (Susceptible-Infective-Susceptible) and SIRS (Susceptible-Infective-Recovered-Susceptible) models have been developed to study disease outbreaks ^[1,2,3,4,5]. Since certain diseases (e.g., childhood diseases) are age dependent, age-structured epidemic models have attracted the attention of many scholars ^{[6,7,8,9,10,11]}. In ^[6], Busenberg found that a sharp threshold (defined by the spectral radius of a positive linear operator) exists and can determine the global behavior of an age-structured epidemic model. In ^[7], Cao investigated the existence and global stability of all equilibria for an age-structured epidemic model with imperfect vaccination and relapse. It was found that, if the threshold is less than 1, the disease-free equilibrium is globally and asymptotically stable; if the threshold is greater than 1, the endemic equilibrium is globally stable. In reference ^[9], by discretizing the multigroup model, the authors transformed a PDE (Partial Differential Equations) system into an ODE (Ordinary Differential Equations) system, and proved that the global asymptotic stability of each equilibrium of the discretized system is completely determined by threshold $R_0$. The threshold is defined as the basic reproduction number, which denotes the expected value of secondary cases produced by infective individuals during the entire infectious period when the entire population are susceptible ^[12].

As the threshold that controls disease outbreaks, $R_0$ plays an extremely important role in assessing disease transmission trend and in reducing disease burden. However, for most age-structured epidemic equations such as the system in ^[13], the basic reproduction number, $R_0 = \int_0^{a_\dagger}k(\sigma)e^{-\int_0^\sigma\mu(\eta)d\eta}\frac{1}{\gamma}(1-e^{-\gamma\sigma})d\sigma\int_\mathbb{R}\tilde{P}(\omega)d\omega$, is merely a theoretical expression of the next generation operator, and is always difficult to calculate. It is a common practice to use numerical approaches to approximate the threshold value ^[10,14]. Since many widely used epidemic models do not satisfy the global Lipschitz coefficients required for using the explicit Euler-Maruyama (EM) scheme, we propose the semi-implicit theta-scheme ^[15,16], which is known as the backward EM when $\theta = 1$, to approximate the exact basic reproduction number. We also estimate the approximate error of the exact basic reproduction number and the numerical threshold.

The novelty of this paper is that we use the theta scheme to discrete the linear operator produced by the infective population in a finite dimensional horizon, so that we can find out the spectral radius, which is the positive dominant eigenvalue of a nonnegative irreducible matrix defined by the next-generation operator. Subsequently, based on the spectral approximation theory ^[17], we obtain the threshold that converges to the exact basic reproduction number under a relatively weak condition (i.e., the compactness of the next-generation operator needs to be satisfied). These results are expected to be useful for studying infectious diseases.

The rest of this paper is organized as follows: in Section 2, the theta scheme is constructed based on the operator theory, and the scheme yields the numerical approximation of the basic reproduction number for a deterministic and a stochastic age-structure epidemic system. Section 3 presented several numerical simulations to illustrate the theoretical results. Concluding remarks are given in Section 4.

2.   Numerical approximation for the basic reproduction number

2.1.   Theta scheme approximation for the deterministic age-structured SIRS system

In this section, we first present the age-structured SIRS epidemic model developed by ^[11],

where $S(t, a)$, $I(t, a)$ and $R(t, a)$ denote the density of susceptible, infective and recovered individuals of age $a$ at time $t$, respectively. Define the force of infectious $\lambda(a, t)$ by

The condition $S(t, 0) = \Lambda$ means that the newborns are all susceptible, $\Lambda$ is the recruitment rate of the population. $S_0(a)$, $I_0(a)$ and $R_0(a)\in L^1(0, A)$ for $\forall a\in[0, A]$. All parameters are positive and their meanings are shown in Table 1.

Let us consider system (2.1) on the Banach space $X: = L^1(0, A)\times L^1(0, A)\times L^1(0, A)$. Let $T$ be a linear operator defined by

$\varphi(a) = (\varphi_1(a), \varphi_2(a), \varphi_3(a))^\top\in D(T)$, where the domain $D(T)$ is given as

The disease-free equilibrium of model (2.1) is $E = \big(E^0(a), 0, E^r(a)\big)$, where $E^r(a) = e^{-\int_0^a(\mu(\eta)+\gamma(\eta))d\eta}$, and $E^0(a) = \gamma(a)E^r(a)\int_0^ae^{-\int_\varrho^a\mu(\eta)d\eta}d\varrho$ is the density of the susceptible population at age $a$ in the disease-free state. Then we define a nonlinear operator $F:X\rightarrow X$ by

Let $u(t) = (S(t, \cdot), I(t, \cdot).R(t, \cdot))^\top$, together with (2.2) and (2.3), system (2.1) has been rewritten as the following abstract Cauchy problem

Next, we mainly consider the second equation of (2.1). By simple calculation, the positive inverse $(-T_2)^{-1}$ is defined as follows

Then, according to ^[11], we can give the next generation operator $\mathcal{K}$ by

Based on the definition in ^[12], the basic reproduction number $\mathcal{\mathcal{R}}_0$ is defined as $r(\mathcal{K})$, where $r(\mathcal{K})$ is the spectral radius of the operator $\mathcal{K}$.

Since the form of $r(\mathcal{K})$ is abstract, we can not calculate $\mathcal{\mathcal{R}}_0$ explicitly. To avoid misunderstanding, we let $B = T_2, G = F_2$, $\varphi_2 = \hbar\in D(B)$,

Hence, we discretize the following system

into a system of ordinary differential equations in $Y_n: = \mathbb{R}^n$, $n\in\mathbb{N}$. Let $\Delta a = A/n$, $a_k: = k\Delta a, \beta_{kj}: = \beta(a_k, a_j), \mu_k: = \mu(a_k), \nu_k: = \nu(a_k)$ and $\delta_k: = \delta(a_k), k = 0, 1, \ldots, n, j = 1, 2, \ldots, n.$ Then the abstract Cauchy system (2.5) is discretized as

where $I(t)$ and $I_0$ are $n-$column vectors, $B_n$ and $G_n$ are $n-$square matrices with the following form

where $M_i = \mu_i+\nu_i+\delta_i(i = 1, \cdots, n)$, $N_i = (1-\theta) E^0_i+\theta E^0_{i+1}(i = 0, \cdots, n-1)$. The additional parameter $\theta\in[0, 1]$ allows us to control the implicitness of the numerical scheme ^[16], for technical reasons we always require $\theta\geq\frac{1}{2}$. Here we denote the next generation matrix $\mathcal{K}_n: = G_n(-B_n)^{-1}$, $\mathcal{R}_{0, n}: = r(\mathcal{K}_n)$ is the threshold corresponding to $\mathcal{R}_0$, and $\mathcal{R}_{0, n}$ can be analyzed in a finite horizon. Since $-B_n$ is a nonsingular M-matrix, and $(-B_n)^{-1}$ is positive. Hence, from the Perron-Frobenius theorem ^[18], we know that $r(\mathcal{K}_n)$ is the positive dominant eigenvalue with algebraic multiplicity 1.

We give two bounded linear operators $\mathcal{P}: Y\rightarrow Y_n$ and $\mathcal{J}: Y_n\rightarrow Y$ as follows

where $k$ is the $k$th entry of a vector, $\top$ is the transpose of matrix $\psi$, and $\chi_{(a_{k}, a_{k+1}]}(a)$ is the indicator function which implies that

From Section 4.1 in ^[19], we know that for all $n\in\mathbb{N}$, $\|\mathcal{P}_n\|\leq1$ and $\|\mathcal{J}_n\|\leq1$. We denote $\|\cdot\|_{Y_n}$ is the norm in $Y_n$, and

Next, we apply the spectral approximation theory to present the convergence theorem of the basic reproduction number.

Theorem 2.1. Assuming that $\mathcal{K}$ is compact, if for any $\hbar\in Y$, $\lim\limits_{n\rightarrow+\infty}\|\mathcal{J}_n\mathcal{K}_n\mathcal{P}_n\hbar-\mathcal{K}\hbar\|_Y = 0$, then $\mathcal{R}_{0, n}\rightarrow\mathcal{R}_0$ as $n\rightarrow+\infty$, preserving algebraic multiplicity 1.

Proof. It is easy to see $\mathcal{K}$ is strictly positive and irreducible, then from Theorem 3 in ^[20] and the Krein-Rutman theorem in ^[21], yield that $\mathcal{R}_0 = r(\mathcal{K}) > 0$ is the maximum eigenvalue of operator $\mathcal{K}$. By a simple calculation, the inverse matrix of $-B_n$ is shown as follows

then we have

where $\bar{E}^0$ and $\bar{\beta}$ denote the upper bounds of $E^0$ and $\beta$, respectively. $\underline{\mu}$, $\underline{\nu}$ and $\underline{\delta}$ denote the lower bounds of $\mu$, $\nu$ and $\delta$, respectively. They are both finite positive.

In addition, we give the following assumption to make that $\mathcal{K}$ is compact.

Assumption 2.1. For any $h > 0$,

where $E^0\beta$ is extended by $E^0(a)\beta(a, \varrho) = 0$ for any $a, \varrho\in(-\infty, 0)\cup(A, \infty)$.

The above assumption implies that the operator $\mathcal{K}$ keep the compactness [11, Assumption 4.4]. In order to prove $\mathcal{J}_n\mathcal{K}_n\mathcal{P}_n$ converges to $\mathcal{K}$ point by point, we provide the following lemma.

Lemma 2.1. For all $\hbar\in Y$, $\lim\limits_{n\rightarrow+\infty}\|\mathcal{J}_n\mathcal{K}_n\mathcal{P}_n\hbar-\mathcal{K}\hbar\|_Y = 0.$

Proof. For any $\hbar\in Y$, we obtain

Since $\|\mathcal{J}_n\|\leq1$, and for any $n\in\mathbb{N}$, $\|G_n\|\leq A\bar{E^0}\bar{\beta}$, we have $L = \|\mathcal{J}_n\|\; \|G_n\| = A\bar{E^0}\bar{\beta}$. Next we estimate the first term in the right-hand of (2.11), then

where $\phi: = (-B)^{-1}\hbar\in D(B)$, and for any $\psi = (\psi_1, \psi_2, \cdots, \psi_n)^\top\in Y_n$,

namely, $\|(-B_n)^{-1}\|\leq A$. From (2.7), we obtain

where $a_0 = a_{-1} = 0$. By the mean value theorem, we have

where $\omega(f, r)$ denotes the modulus of continuity. We know that $\omega(f, r)$ is defined by $\sup_{|x-y|\leq r}|f(x)-f(y)|$ with the following property

Hence, $\|(-B_n)^{-1}\mathcal{P}_n\hbar-(-B)^{-1}\mathcal{P}_n\hbar\|_{Y_n}\rightarrow0$ holds. Then we consider the second term of (2.11) as follows

where $\omega(E^0, \Delta a)\rightarrow0(\Delta a\rightarrow0)$ and $\omega(\beta, \Delta a)\rightarrow0(\Delta a\rightarrow0)$, respectively. Hence,

Combine the above discussion, we have $\lim\limits_{n\rightarrow+\infty}\|\mathcal{J}_n\mathcal{K}_n\mathcal{P}_n\hbar-\mathcal{K}\hbar\|_Y = 0$.

By virtue of Assumption 2.1 and Lemma 2.1, we know that Theorem 2.1 holds. Namely, $\mathcal{R}_{0, n}\rightarrow\mathcal{R}_0$ as $n\rightarrow+\infty$, preserving algebraic multiplicity 1.

2.2.   Theta scheme approximation for the stochastic age-structured SIRS system

In this section, we seem the natural mortality $\mu(a)$ as a random variable $\mu(a)-\sigma\frac{dB_t}{dt}$, where $B_t$ is a standard Brownian motion, $\sigma$ is the intensity of noise perturbation. Then, replace $\mu(a)$ with $\mu(a)-\sigma\frac{dB_t}{dt}$ in system (2.1), we can obtain a stochastic age-structured SIRS model

Next, we analysis the stochastic basic reproduction number. In the same way, we take the infective population of system (2.13) into account, and substitute $S(t, a) = E^0(a)$ into it, we derive

According to the general definition of the stochastic basic reproduction number, the following two operators are defined on $Y: = L^1(0, A)$

and

Using $\mathcal{A}$ and $\mathcal{F}$ to rewrite (2.14) as

Then we have

The next generation operator $\mathcal{T}$ is shown by

Samely, we define $r(\mathcal{T})$ as the basic reproduction number $\mathcal{R}{_0^s}$ of the stochastic system (2.13), and $\mathcal{R}_{0, n}^s: = r(\mathcal{T})$ is the threshold corresponding to $\mathcal{R}_0^s$.

Next, we discretize (2.16) in $Y_n: = \mathbb{R}^n$, $n\in\mathbb{N}$. Then the system (2.16) is discretized into the following equation

where $\mathcal{A}_n$ is defined as the same as $B_n$($\mathcal{A}_n: = B_n$), and

where $\theta\in[\dfrac{1}{2}, 1]$, $Q_i = (1-\dfrac{\sigma^2}{2})\Big[(1-\theta) E^0_i+\theta E^0_{i+1}\Big](i = 0, 1, \cdots, n-1)$.

Theorem 2.2. From Theorem 2.1, we know that $\mathcal{T}$ is irreducible, compact and strictly positive. If

for any $\hbar\in Y$, then

where $\mathcal{J}_n$ and $\mathcal{P}_n$ are defined as (2.7).

Proof. Obviously, $\mathcal{R}_0^s = r(\mathcal{T}) > 0$, and $r(\mathcal{T})$ is the spectral radius of operator $\mathcal{T}$. We know that $(-\mathcal{A}_n)^{-1} = (-B_n)^{-1}$, and $(-B_n)^{-1}$ is given by (2.9). Then we have

where $\underline{E}^0$ is the lower bound of $E^0$.

Next, we verify that $\lim\limits_{\Delta\rightarrow0}\|\mathcal{J}_n\mathcal{T}_n\mathcal{P}_n\hbar-\mathcal{T}\hbar\|_Y = 0$. For any $\hbar\in Y$, we have

where the first term of (2.18)

is similar to the first term in the right-hand of (2.11), so it is easy to see that

Next we estimated the second term of (2.18). Let $\varpi: = (-\mathcal{A})^{-1}\hbar\in D(\mathcal{A})$, we obtain

Thus, $\|\mathcal{J}_n\mathcal{F}_n\mathcal{P}_n(-\mathcal{A})^{-1}\hbar-\mathcal{F}(-\mathcal{A})^{-1}\hbar\|_Y\rightarrow0$ holds. Hence, we obtain the desired assertion.

In conclusion, Theorem 2.2 holds, which implies that $\mathcal{R}_{0, n}^s\rightarrow\mathcal{R}_0^s$ as $\Delta\rightarrow0$, preserving algebraic multiplicity 1.

Remark 2.1. Compared with ^[10], our paper has two advantages:

● ^[10] employed a backward Euler method to approximate $R_0$, and obtain the numerical threshold $R_0^n\rightarrow R_0$ as $n\rightarrow\infty$. In present paper, we propose a $\theta$ method is know as the backward EM when $\theta = 1$, and the explicit Euler-Maruyama(EM) scheme when $\theta = 0$. The $\theta$ scheme has the parameter $\theta$, and different $\theta$ values give different convergence rates. Therefore, we can use the $\theta$ method to find the optimal convergence rate. And the backward Euler method is a special case when $\theta = 1$ of our method. Our work provides an extension of ^[10].

● A deterministic age-structured epidemic model is discussed in ^[10], but in present paper, we studied not only the deterministic system but also the stochastic age-structured epidemic model, and the stochastic system is more practical.

3.   Numerical simulations

In this section, numerical examples are shown to verify our Theorems. In what follows, let $A = 100$, $\mu(a) = 0.2(1+\frac{a^3}{10^3})$ (^[13], see Fig. 1(a)), $\gamma(a) = \gamma = 0.25$, $\nu(a) = \nu = 0.1$ and $\delta(a) = \delta = 0.05$ (see ^[22]). Thus, $E^0(a) = \gamma(a)E^r(a)\int_0^ae^{-\int_\varrho^a\mu(\eta)d\eta}d\varrho = 0.25e^{(-0.45a-\frac{a^4}{2\times10^4})}\int_0^ae^{\varrho(\frac{\varrho^3}{2\times10^4}+0.2)-a(\frac{a^3}{2\times10^4}+0.2)}d\varrho$. Based on numerical integration for $E^0(a)$, we obtain Fig. 1(b).

In this example, we do not specify what kinds of influenza-like disease it is, and the value of $\mathcal{R}_0$ is in the range of 2-3 ^[23]. We assumption that the disease is more likely to transmission between individuals with similar ages ^[10], then we let $\beta(a, \varrho) = kJ(a-\varrho)$, where $k = 200$ and $J(x) = 0.6(-x^2+100^2)\times10^{-6}+0.001$ is a normalized distance function. Thus, we can easily verify that Assumption 2.1 is true. Hence, Theorem 2.1 and 2.2 hold, which implies that $\mathcal{R}_{0, n}\rightarrow\mathcal{R}_0(\mathcal{R}_{0, n}^s\rightarrow\mathcal{R}_0^s)$ as $n\rightarrow+\infty$.

3.1.   Numerical approximation of $\mathcal{R}_{0, n}$ for the deterministic system

Let $\theta = 0.5$, and choose $\mathcal{R}_{0, 1000}\approx2.57673470573749 = :\mathcal{R}^*$ as a reference value for $\mathcal{R}_0$. From Fig. 2(a), we see that the threshold $\mathcal{R}_{0, n}$ for the discretized system (2.6) respect to the reference value $\mathcal{R}^*$ as $n$ increases. Further more, the error $\mathcal{R}^*-\mathcal{R}_{0, n}$ converges to zero as $n$ increases (see Fig. 2(b)). In Fig. 3, we show the numerical simulations of $\mathcal{R}_{0, n}$ at $\theta = 0.5$, $\theta = 0.7$ and $\theta = 0.9$, respectively. It is obvious to see that the value of $\theta$ has a certain impact on the convergence rate of $\mathcal{R}_{0, n}$. The bigger value of $\theta$, the faster rate of convergence. This implies that the backward EM method would make the convergence faster. Our paper verified the work of ^[10].

3.2.   Numerical approximation of $\mathcal{R}^s_{0, n}$ for the stochastic system

In this example, let $\sigma = 0.1$, and we also choose $\mathcal{R}^s_{0, 1000}\approx2.56385103220880 = :\mathcal{R}_s^*$ as a reference for $\mathcal{R}_0^s$. Similarly, the threshold value $\mathcal{R}^s_{0, n}$ for the discretized system (2.17) respect to the reference value $\mathcal{R}_s^*$ (see Fig. 4(a)) and the error $\mathcal{R}_s^*-\mathcal{R}^s_{0, n}$ converges to zero as $n$ increases (see Fig. 4(b)). Fig. 5(a) give a comparation for $\mathcal{R}^s_{0, n}$ at $\sigma = 0.1$, $\sigma = 0.5$ and $\sigma = 0.8$, respectively. We can see that the intensity of environmental disturbance has a great influence on the threshold $\mathcal{R}^s_{0, n}$. The higher value of $\sigma$, the smaller value of $\mathcal{R}^s_{0, n}$. This means that the intensity of environmental fluctuation can reduce the threshold of disease outbreak, which may be a better measure to control disease outbreak. We show a $3D$ simulation of $R_{0, n}^s$ corresponding to $\theta\in[0.5, 1]$ and $\sigma\in[0, 1]$ in Fig. 5(b), the effect of $\sigma$ on the threshold $R_{0, n}^s$ with the change of $\theta$ is further explained.

4.   Concluding remarks

For the age-structure epidemic model, the basic reproduction number is defined as an integral and difficult to be estimated. Hence, it is necessary to approximate it using numerical methods. This paper investigates the numerical approximation of two basic reproduction numbers for deterministic and stochastic age-structured epidemic systems, respectively. We use the theta scheme to discrete the infective population in a finite space, so that the two abstract basic reproduction numbers can be calculated explicitly. Afterward, using the spectral approximation theory, we obtain the numerical threshold that converges to the exact value as $n$ increases. We also estimate the approximation error between the exact basic reproduction number and its numerical approximation. Finally, several numerical simulations are shown to illustrate our theoretical results. The numerical results show that, for the deterministic system, the convergence rate of $\mathcal{R}_{0, n}$ is faster when $\theta$ is bigger under the condition of $\theta\in[\frac{1}{2}, 1]$. For $\theta\in[0, \frac{1}{2}]$, the proof of the pointwise convergence in Lemma 2.2 remains challenging, and is warranted to be investigated in a future study. For the stochastic system, the appropriate noise intensity can reduce the threshold of disease outbreak.

Acknowledgments

The research is supported by the Natural Science Foundation of China (Grant number 11661064).

Conflict of interest

The authors declare there is no conflict of interest.