Positive Numerical Approximation of Integro-Differential Epidemic Model

: In this paper, we study a dynamically consistent numerical method for the approximation of a nonlinear integro-differential equation modeling an epidemic with age of infection. The discrete scheme is based on direct quadrature methods with Gregory convolution weights and preserves, with no restrictive conditions on the step-length of integration h , some of the essential properties of the continuous system. In particular, the numerical solution is positive and bounded and, in cases of interest in applications, it is monotone. We prove an order of convergence theorem and show by numerical experiments that the discrete final size tends to its continuous equivalent as h tends to zero.


Introduction
Epidemic models based on nonlinear integral and integro-differential equations [1][2][3][4][5][6] allow to include the contribution of time since becoming infected to the total infectivity. This property makes these (age-of-infection) models suitable to describe diseases such as smallpox [7], cholera [8], SARS [9], COVID-19 [10][11][12] and AIDS [13,14]. Numerical simulations support the comprehension of the phenomenon and provide a tool for its quantitative description. Therefore, attention should be paid to set up a numerical framework that allows to supply real-time and reliable answers. We consider the integro-differential equation representing the Kermack and McKendrick age-of-infection epidemic model: for which we refer to [5,15] and references therein. Here, S(t) is the number of susceptibles at time t, β represents the rate of effective contacts and A(s) ≥ 0 is the mean infectivity of members of the population with infection age s ( A ∈ L 1 (R + )). Equation (1) represents a special nonlinear Volterra integral equation and we refer to [16,17] for the theory about these problems. From the numerical point of view, there is a growing interest in the development and analysis of methods for the approximation of the solution to integral and integro-differential equations of Volterra type (see, for example, [16,18] and references cited therein). Recently, great attention is paid to the construction of methods that reflect and preserve the essential properties of the problem and that efficiently integrate the memory term (see, for example, [19,20]). In particular, for epidemic models, non-standard discretization methods allow the design of numerical models that qualitatively and quantitatively behave like the reference problem [21][22][23][24][25]. In [26], a non-standard discretization is used to integrate problem (1). Such a scheme is explicit, it preserves positivity of solutions, boundedness and the final size relation for the epidemic, and therefore it assumes the role of a discrete model. On the other hand, the method is only linearly convergent, and this can be a disadvantage when you have to integrate over long times, balancing accuracy and computational cost. Our aim here is to develop a higher order method. For numerical solutions obtained by direct discretizations of (1), results that guarantee the unconditional dynamic consistency properties of non-standard schemes are not straightforward. These considerations motivated us to use a different approach. The idea is to exploit the exponential form of the evolution operator in (1) and integrate it by direct quadrature methods of any order. The derived discrete equation has the advantage of automatically ensuring the positivity of the solution and, at the same time, allowing an asymptotic analysis that leads to results consistent with the continuous model. This paper is organized as follows: In Section 2, we report the main facts about the age-of-infection Model (1) as developed in [15,27]. Then, in Section 3, we formulate the numerical method, we study the convergence and prove the positivity and boundedness of the numerical solution, for any value of the stepsize. Furthermore, we analyze the asymptotic dynamics of the numerical model by deriving the numerical final size relation. In Section 4, we obtain additional properties of the numerical solution in cases of infectivity functions of interest in applications. Finally, numerical experiments are reported in Section 5, to show the theoretical results obtained. Some remarks, in Section 6, conclude the paper.

The Continuous Model
We refer to [15,27] for a complete description of Model (1). Here, we mention the main features that will be useful in the rest of the paper. In the following, we assume that the population has a constant size N and that where S 0 = S(0) ≤ N (see, for example, [28]). The basic reproduction number is defined as We want to point out that, for the solution S(t) of (1), it holds: S ∞ is the unique solution of the final size relation.
Our starting point is the following equation which is equivalent to (1) with (2). This equivalence is obtained by using the identity into Equation (1) with (2). It comes out that Then, integration in [0, t], and straightforward calculations lead to the form (5).

The Numerical Model
Let h > 0 and {t n } n≥0 be a uniform mesh such that t n = nh. The direct discretization of (5) by direct quadrature (DQ) methods with Gregory convolution weights (see, for example, [18]) reads where S n ≈ S(t n ), n ≥ n 0 ≥ 1, and w nj and ω j are the non-negative quadrature weights.
The starting values S 0 , S 1 , . . . , S n 0 −1 are given. When n 0 = 1, ω 0 = 0 and ω j = 1, the numerical method above corresponds to the discrete-time Kermack-McKendrick model introduced by Diekmann in [29]. This motivates us to deepen the analysis of the dynamic behavior of (6) since, as it will be clear later in this section, from a numerical point of view it possesses good convergence properties and, in many cases, is easy to implement.

Basic Properties
At each step and for any fixed h > 0, Equation (6) implicitly defines S n as the solution of the nonlinear equation The nonlinear function (7) is concave and F n (0; h) < 0. If we assume that S 0 , . . . , S n−1 ∈ [0, N], then F n (N; h) > 0. Thus, F has exactly one root in [0, N]. Furthermore, this root is less than S 0 . The application of this result leads to state the following theorem. Theorem 1. Let {S n } n≥0 be the solution to the discrete Equation (6). Assume that 0 < S j ≤ N for j = 0, . . . , n 0 − 1, then for each h > 0: 1.
The sequence {S n } n≥0 is positive; 2.
The sequence {S n } n≥0 is bounded from above by S 0 .

Convergence
From now on, we assume that the quadrature weights in (6) satisfy (see, for example, [18] (p. 79, Theorem 2.6.10)) sup n≥0 max 0≤j<n 0 Consider δ(h; t n ) is the local error for (6). It is easy to prove that p = n 0 + 1 is the largest integer so that max 0≤n≤M |δ(h; t n )| ≤ ch p .
Proof. For n = n 0 , . . . , M, where δ(h; t n ) is the local truncation error defined in (9). Let e(h; t n ) = S(t n ) − S n , be the global error. Using standard manipulations and choosing h < The Gronwall discrete inequality (see, for example, [30] (p. 101)) yields Since the local error δ(h; t n ) and the starting errors η j (h) satisfy (10) and (11), respectively, it follows that with C positive constant not depending on h.

The Numerical Final Size
In this section, we assume that Moreover, we also assume that is ultimately monotonic, see [32] (p. 208)). We define the numerical basic reproduction number as the direct discretization (see [18]) of (3). In the Theorem below, we obtain a relation between S ∞ (h) and R 0 (h) which represents the discrete version of (4). (5) be sufficiently smooth and let {S n } n≥0 be the solution to the discrete Equation (6). If the starting values S j ≤ S 0 , for j = 0, . . . , n 0 − 1, are positive and (8) and (13) hold, then for each h > 0

Trapezoidal Discretization in Some Realistic Cases
In this section, we consider particular choices for A(t) in (5) which appear in applications (see, for example, [7], [27] (sec. 2.1), [11,33]). Indeed, we focus on cases where the infectivity function A(t) is identically zero on an initial time interval, The condition A(0) = 0 makes the numerical scheme (6) explicit, simplifying its analysis and reducing the computational cost at each step. We choose to approximate this type of model using the trapezoidal DQ method. In this case, the numerical solution not only agrees with the basic properties described in Theorem 1, without restrictions on the stepsize, but it is also monotone, and therefore retains the asymptotic dynamic of the analytical solution. The resulting explicit scheme has the form, for n > 0, with n 0 = 1, for which we state the following theorem.
Theorem 4. Let {S n } n≥0 be the solution to the discrete Equation (18). Let S 0 > 0 then, for each h > 0, the sequence {S n } n≥0 is non-increasing.
Proof. Theorem 1 assures non-negativity and boundedness of the sequence. Here, we prove that it is also non-increasing by induction, starting from Consider n > 1 and assume that the assertion is true for each j ≤ n − 1. Since we have Theorem 4 assures that the solution {S n (h)} n≥0 to the discrete Equation (18) is nonincreasing for each positive value of h and then it admits a finite limit S ∞ (h), as n goes to +∞. Hence, condition (13) holds for each h > 0, and the solution of Equation (14) is indeed the numerical final size S ∞ (h).
For the sake of comparison, we consider the reformulation of problem (1) given in [15] S(t) = S 0 − β t 0 S(s)ϕ(s) ds, and compare the behavior of the numerical solutions obtained by (18) and the following trapezoidal discretization of (19) for n > 0. Both methods are convergent of order 2. However, in (20), positiveness and monotonicity of the solution can be proved only under the following constraint on the discretization stepsize, h < 2/(βN · sup A(t)), leading to a severe limitation in case of large values of the population size N.

Numerical Experiments
In this section, we report some numerical examples in order to show experimentally the theoretical results proved in Sections 3 and 4. For our experiments, we choose illustrative test equations and we use method (6). Our first example consists of problem (5) with for t ∈ [0, 1]. Here, we have integrated problem (21) by DQ method with trapezoidal, first and second Gregory weights (the orders of convergence are 2, 3 and 4, respectively). Since A(t) does not meet condition (17), the methods are implicit and require at each step the resolution of a nonlinear equation for which a fixed-point iteration process is used. Table 1 and Figure 1 show the numerical error behavior for different values of the stepsize h. The numerical errors E n = S n −S n , are computed against the reference solutionS n , obtained by second Gregory rule with h = 10 −6 . For the sake of comparison, we also report the errors obtained by integrating (21) with the first order non-standard finite difference (NSFD) method used in [26]. It is clear that the experimental order of convergence of DQ methods is coherent with the theoretical one predicted in Theorem 2. Figure 2 represents the work precision diagram of the number of function evaluations with respect to the accuracy of the numerical solutions by the four methods considered here. It shows better performances in higher order methods. Indeed, even though we are comparing an explicit method with implicit ones, we see that, at the same computational cost, the accuracy of DQ methods is higher than NSFD. In addition, for given accuracies, even low ones, the NSFD scheme requires much more effort. In order to show the long time behavior of the numerical solution, we consider problem (5) with: A(t) = 2t e −t , T = 80, N = 10 5 , S 0 = 9 × 10 4 , β = 10 −5 .
In Figure 3, the behavior of the numerical solution, computed by (6) and second Gregory weights, with h = 10 −3 , is reported. In this case, the theoretical value for the final state, obtained by solving the nonlinear Equation (4) by the MATLAB routine fzero (see [34]), is S ∞ = 1717.    We introduce the truncated numerical basic reproduction numberR 0 (h) and the truncated numerical final sizeS ∞ (h; T) as follows and define the residual on the numerical final size relation (14) as In this case, since (17) holds, Theorem 4 guarantees the existence of lim n→+∞ S n when the DQ method has trapezoidal weights. What is more, from our tests we observe that, for fixed h > 0, the value ofS ∞ (h) converges to a finite limit as T grows for any choice of the quadrature weights. Here, we have used T = 80 for computingS ∞ (h), and the effectiveness of this choice is confirmed by the values of r(h; T) in Table 2. In particular, Table 2 shows, for different values of h, the relative errors for the approximation of R 0 and S ∞ by (23), as well as the corresponding numerical final size residuals.
Numerical results represent a validation of Theorem 3 but also give evidence of unproven theoretical properties. Indeed, it is clear from Table 2 that, for h → 0, the numerical final size converges to the continuous one, with the same order p ≥ 1 of the quadrature rule employed in (6).  Here, we consider, referring to [33], a kernel function suitable for modeling influenza t 1 = 0.4 days, t 2 = 2.8 days, t 3 = 3.8 days, t 4 = 7.4 days, based on the assumptions that no member of the population is infectious before t 1 days or after t 4 days post-exposure and maximum infectivity occurs between t 2 and t 3 days after exposure. The latent period T L and the infectious period T I are We integrate problems (5)-(24) with T = 10, N = 10 5 , S 0 = 9 × 10 4 , β = 10 −2 , by methods (18) and (20), respectively. The results, reported in Figure 4, have been obtained with two different values of the stepsize. We remind that, according to Theorem 4, the numerical solution computed by the trapezoidal scheme (18) is positive and non-increasing, regardless of the stepsize h. Conversely, we observe in Figure 4 that the numerical solution computed by means of the direct trapezoidal scheme (20), with h = 0.5, is increasing or even negative in some points.

Conclusions
We have introduced a numerical method specially designed for age-of-infection epidemic models represented by integro-differential equations of the type (1), reformulated as an integral equation in exponential form. The proposed method discretizes the integral term by DQ methods with Gregory convolution weights, from which it inherits the order of convergence. We show that the numerical solution is positive and bounded, which are crucial properties taking into account the biophysical meaning of variables. In some cases of realistic infectivity functions (for example, incorporating a non-zero finite incubation period and a finite infectious period), we prove that the proposed scheme, based on the trapezoidal method, also retains the monotonicity of the solution. In general situations, we provide numerical evidence of this behavior. A comprehensive analysis allows us to state that the numerical solution is dynamically consistent with the continuous one. We introduce a discretization of the basic reproduction number and prove that the numerical final size satisfies a nonlinear equation which is the discrete counterpart of the continuous final size relation. Numerical simulations show convergence, with respect to h, of the numerical final size to the continuous one, at a rate that reflects the order of the employed quadrature rule. A numerical comparison clearly shows that the numerical approach described here outperforms the first order dynamic-preserving non-standard explicit scheme introduced in [26]. This fact is advantageous in long-time integrations.
Author Contributions: The authors E.M., M.P. and A.V. contributed equally to this work. All authors have read and agreed to the published version of the manuscript.