key: cord-0658056-9hey9722
authors: Tak'acs, B.; Farag'o, I.; Horv'ath, R.; Repovvs, D.
title: Qualitative properties of space-dependent SIR models with constant delay and their numerical solutions
date: 2021-12-13
journal: nan
DOI: nan
sha: 7bf0706bda5342998546e8d7112b2a900464ca71
doc_id: 658056
cord_uid: 9hey9722

In this article a space-dependent epidemic model equipped with a constant latency period is examined. We construct a delay partial integro-differential equation and show that its solution possesses some biologically reasonable features. We propose some numerical schemes and show that by choosing the time step to be sufficiently small the schemes preserve the qualitative properties of the original continuous model. Finally, some numerical experiments are presented that confirm the aforementioned theoretical results.

The increasing rate of globalization in recent decades led to an interconnected world in which diseases can spread faster than ever: this threat became reality in 2020 with the outbreak of the Covid-19 virus. Thus, the importance of adequate modeling of epidemics is apparent.

One of the most frequently used tools in mathematics to model the spread of diseases is the SIR model [5, 16] . These models can be used to describe the spread of some feature among a group of individuals. We split our population into three categories: class S contains those who do not have the property yet, class I includes those who have the property and they can also transmit it to others, and class R has those who had the property, but they do not have it any more. Note that these models can be used not only in epidemics but also in several other fields, see e.g. [1, 3, 19, 25] -however, in this article we are going to focus on epidemic models. The SIR models are generally systems of ordinary differential equations that do not take into the account the spatial location of the individuals of the above groups, although Kendall proposed a possible extension almost 60 years ago [14, 15] . This extension results in a system of partial integro-differential equations that can be solved only numerically. The numerical solution can be produced by some appropriate methods but it is a natural requirement that it should possess the characteristic qualitative properties of the disease propagation process. These properties were investigated in the recent papers [7, 11, 23, 24] . The qualitative properties of diffusion SIR models, which are also used to model spatial disease propagation, can be found e.g. in [5, 9, 26] .

Some diseases take some time to develop inside an infected individual, so they do not start to infect upon their infection but only after a short period of time (so-called latency period). For this, we generalize the previous models by using delay integro-differential equations. In this paper, we formulate the basic qualitative properties of these models, we show that the solution of the continuous problems possesses these features and give sufficient conditions that guarantee the properties to the numerical solutions.

In Section 2, we construct our mathematical models, and then we prove some properties of the analytical solution in Section 3. Section 4 deals with the spatial discretization of the model, while Section 5 discusses the application of the explicit Euler and Runge-Kutta methods. Then some numerical results are shown in Section 6.

The first, and possibly best known mathematical formulation of the SIR models comes from Kermack and McKendrick [16] , in which the authors propose the following model: 

where the functions S(t), I(t), R(t) : [0, T ] → R (T ∈ R + ) denote the number of individuals in class S (susceptible), class I (infected) and class R (recovered) at time t, respectively, while a, b ∈ R + are parameters describing the rate of infection and recovery, respectively. In this paper, we will consider the more realistic model

∂S(t, x, y) ∂t = −S(t, x, y)F I (t − σ, x, y) − cS(t, x, y),

that also contains a delay term besides the generalizations of papers [7, 23, 24] . In (2), X : (0, T ] × Ω → R, (t, x, y) → X(t, x, y) (X ∈ {S, I, R}) denotes the spatial density of the class X at the point (x, y) ∈ Ω ⊂ R 2 (Ω is a connected, bounded open set) and at the time instant t. The parameter b is the same as in (1) . The term cS(t, x, y) (c ∈ R + ) describes the effect of possible vaccination, which means that it is possible for some individuals to get from class S to class R without entering class I. The infection of a given individual is caused by the infected individuals in its surroundings, due to the assumed latency period at a previous time instant, and this effect depends on the spatial position of the individuals. Here we suppose that one can only be infected by others in a δ-radius (δ ∈ R + ) neighborhood around itself and the effect of the infection is described by a given non-negative, continuous and bounded weight function W . The constant delay σ > 0 is the length of the latency period of the disease (see e.g. [6, 18, 27] ). This phenomenon is included into the model by the delay term

where B δ (x, y) is the δ radius open ball centered at (x, y) and W is the weight function. The natural births and deaths are not taken into consideration in the model. Because (2) is a delay system, to obtain a properly posed problem we also need the values of the function I on the time interval [−σ, 0]. This is the so-called history function that will be denoted by I h : [−σ, 0] × Ω → R, (t, x, y) → I h (t, x, y). Later we will show that some assumptions are needed on the history function to assure that the solution is continuous in space and time and behaves as expected. The history functions of S ans R are denoted similarly but these functions do not appear in the model.

The model does not have boundary condition in a classical sense but due to the integral in (3) we assume that I is equal to zero outside the domain Ω.

In the next section we show that our problem has a unique solution, which behaves in a biologically reasonable way.

In this section we discuss the solvability of system (2) . The key tool will be the method of steps introduced by Bellman [2] , which involves the splitting of our time interval [0, T ] into smaller intervals with length σ with the grid points 0, σ, 2σ, 3σ, . . . , T . The (k + 1)th element of this list will be denoted by t k . Let us denote the solution on the time interval (t k−1 , t k ] by S k (t, x, y), I k (t, x, y) and R k (t, x, y). Then the delay differential equation (2) on a given interval (t k , t k+1 ] (k ≥ 1) becomes a classical partial integro-differential equation

with initial conditions

In the case of t ∈ (0, σ], (4) has the form

with initial conditions S 1 (0, x, y) = S h (0, x, y), I 1 (0, x, y) = I h (0, x, y), R 1 (0, x, y) = R h (0, x, y).

Because the structures of (4) and (5) are the same, we will only consider the case of (4) but all the statements also will be valid for (5).

If we proceed to solve the equation step by step (i.e. solving it on (0, σ] then on (σ, 2σ] and so on) then the term F I k (t − σ, x, y) is known, since it only depends on I k (t, x, y). It is also important to notice that the first equation of (4) does not depend on the later ones, thus can be solved independently of the others. After this S k+1 (t, x, y) will be known, which means that the only unknown is I k+1 (t, x, y) in the second equation, so it can also be solved directly. Finally, the last equation can also be solved by integration, resulting in the solution of system (4) .

For the sake of simplicity, we introduce the functionF k+1 (t, x, y) := F I k (t − σ, x, y). In other words, we shift the domain of

∂R k+1 (t, x, y) ∂t = bI k+1 (t, x, y) + cS k+1 (t, x, y).

Theorem 1. Assume that the history functions S h (t, x, y), I h (t, x, y) and R h (t, x, y) are continuous in time and also in spatial variables (x, y). Then system (2) has a unique solution, which is continuously differentiable in time on (0, T ] and is continuous in space.

Proof. We can express the solution of (6) as

For the solution to be continuous in time at point t k , we take

since in this case lim t→t k X k+1 (t, x, y) = X k (t k , x, y) for X ∈ {S, I, R}. (Note that the values of K 1 , K 2 and K 3 are defined differently for different choices of k.) Also, because of the form of the solutions, it can be proved by induction that if the history functions are continuous in time, then our solution will also be continuous in time on the whole solution domain. It can be shown similarly that if the history functions are continuous in space, then the solution is also continuous in space since the weighting function W was assumed to be continuous. We show that the solution is not only continuous in time but is also continuously differentiable in that variable. By the fundamental theorem of calculus, it is also easy to see that the solutions are continuously differentiable in t on each (t k , t k+1 ). At t k (k ≥ 1), we show that

However, because of the previous arguments we know that the right-hand side of (4) is continuous in time on [t k , t k+1 ), so e.g. in the case of S k+1 , we have

A similar argument holds in the case of I and R. This shows the continuous differentiability of the solution in variable t. Another question is whether the solution (7) is the only one. It can easily be shown that taking account only the first equation of (6), it can be thought of as an ordinary differential equation, and its right-hand side fulfils the Lipschitz property. The same can also be shown for the second one (assuming that S k+1 (t, x, y) is a known function), and also for the last one. Consequently, solution (7) is unique.

In the next section we observe whether the continuous solution possesses some biological features.

A mathematical model is considered to be reasonable not only when the system has only one solution but the behavior of such solution should also possess the properties of the biological model. Here we are going to observe four of these properties.

A natural requirement is that the density of each species cannot be negative.

The functions S(t, x, y), I(t, x, y) and R(t, x, y) should be non-negative.

Since we assume that there are no births or natural deaths in our region (or their rate is assumed to be equal) and the individuals do not move, the sum of the densities of the three classes should remain constant at each point, and this should also hold for the total number of individuals in the observed domain (which is given by the integral of the sum of the densities).

The sum S(t, x, y) + I(t, x, y) + R(t, x, y) should be constant in t for all (x, y) points inside our domain. Also, a consequence of this is that the following integral is also constant in time:

Since there is no way one individual can get to class S, the density of it cannot increase.

Also, no individual can get out of class R, so its density cannot decrease.

The function R(t, x, y) cannot decrease in t at each point (x, y).

Our goal is to show that the properties C 1 -C 4 are satisfied by the solution of system (2).

Theorem 2. Suppose that properties C 1 -C 4 hold for our history functions, which are also continuous in space as well as in time. Then C 1 -C 4 also hold for the solution of system (2) on a time interval (0, T ].

Proof. By Theorem 1, system (2) has a unique solution which is continuously differentiable in variable t for t ∈ (0, T ].

We will proceed by induction: first we prove the properties on (t 0 , t 1 ], then by supposing that they hold on (t k−1 , t k ], we prove them on (t k , t k+1 ]. However, since the proof on the first interval does not differ that much from the proof on any arbitrary one, we present here only the proof of the latter one.

By adding up the equations of (6), it is clear that property C 2 holds. For property C 1 , let us consider system (6) . The solution of this system is (7) . From this form, it is evident that since K 1 was chosen to be S k (t k , x, y)e ct k ≥ 0, the values of the function S k+1 (t, x, y) are also non-negative. (If k = 0 then S 1 (0, x, y) = S h (0, x, y) ≥ 0.) From the second equality of (7), since K 2 ≥ 0, S k+1 (s, x, y) ≥ 0 andF k+1 (s, x, y) ≥ 0 for s ∈ [t k , t] (by the induction assumption), the relation I k+1 (t, x, y) ≥ 0 is satisfied. The non-negativity of R k+1 (t, x, y) can be easily seen from the third equality of (7), since both I k+1 (s, x, y) and S k+1 (s, x, y) are non-negative for all s ∈ [t k , t], and K 3 = R k (t k , x, y) is also non-negative (if k = 0, then R h (0, x, y) is non-negative because C 1 holds for the history function by assumption). Thus, we proved that C 1 holds.

Property C 3 can be seen if we consider the form of the function S k+1 (t, x, y) in (7): it is evident that since c > 0 andF k+1 (s, x, y) ≥ 0, the function is decreasing (here we also use that K 1 ≥ 0, as mentioned before).

Property C 4 is a consequence of the non-negativity of I k+1 (t, x, y) and S k+1 (t, x, y), since the right-hand side of the third equation of system (2) is non-negative.

In the next sections we examine the semi-discretized and the fully discretized versions of (2) and check whether their solutions satisfy the analogous versions of C 1 -C 4 .

It is evident that Bellman's method can only be used in practice, when the integral of the history function I h at the time t = −σ is known, which is usually not the case. Because of this, in the following sections we propose a numerical approach to approximate these analytic solutions. Upon looking at system (2), it is evident that the most problematic part is the fact that it contains integrals on its right-hand side. In Section 4.1 we approximate the integral term F I (t − σ, x, y) by using some cubatures, and then in Section 4.2 we define a spatial grid on the domain, approximating the partial differential equation by a class of ordinary differential equations, by defining a separate equation for each grid point.

We define some two-dimensional cubature formula on the disc B δ (x, y) with positive weights to approximate the integral F I (t, x, y). Let us introduce a set of cubature points (x+η i , y+ξ i ) on the disc B δ (x, y) and the positive cubature weights w i > 0 (i = 1, . . . , p). Here the values η i , ξ i and w i might be also dependent on (x, y) in the most general setting but in this article for simplicity we use the same cubature formula

for every point of the domain.

Using approximation (8) we obtain the system of partial differential equations

It is clear that the arguments detailed in Section 3 can be used similarly, which results in the following theorem.

Theorem 3. Assume that the history functions are continuous in space and time and properties C 1 -C 4 hold for them. Then system (9) has a unique solution, which is continuously differentiable in time, continuous in space and also has properties C 1 -C 4 .

A natural question is the choice of the numerical approximation of the integral. In [24] two separate choices of cubatures were investigated. One of them, the Elhay-Kautsky cubature results in a uniform distribution of points on the unit disc, while the other, the Gauss-Legendre cubature (which involves a transformation of the integral to a unit square) results in a non-uniform distribution. Numerical experiments show that while the first one works well for polynomials, the second one is better for arbitrary nonlinear functions. Since we cannot guarantee that our function I(t, x, y) is a polynomial, we are going to use here the latter one. For further details of the different methods, see [24] .

To obtain the numerical solution, we assume that our domain Ω is a rectangle with one vertex at the origin, i.e. Ω = (0, A) × (0, B), A, B ∈ R + . Note that the following arguments also hold for domains in a more general form, but involve much more careful choice of spatial grids.

Let us discretize our rectangle shaped domain using the spatial grid

This grid consists of KL points with spatial step sizes h x and h y , and we approximate the continuous solutions by a matrix containing the values at these grid points. After this semi-discretisation, we get the following set of equations:

in which X k,l (t) (X ∈ {S, I, R}) denotes the approximation of the function at grid point (x k , y l ) and at time t and

Note that the points (x k +η i , y l +ξ i ) might not be part of G, so there are no I k,l values assigned to them. Because of this, we will approximate them by some interpolation method using the nearest known I k,l values and positive coefficients. (This is the reason for theÎ notation.) In order to satisfy the qualitative properties, it is important to choose such interpolations that preserve non-negativity in the sense that if the known values are non-negative, then the function we get at the end of our process should also be non-negative. Such interpolations include monotone interpolation that uses piecewise cubic Hermite interpolating polynomials [8, 12] ('pchip' for short), which will be used in the numerical experiments. As in the previous section, the methods described in Section 3 can be used again for system (10) , which results in the following theorem.

Theorem 4. Assume that the history functions are continuous and properties C 1 -C 4 hold for them. Then system (10) has a unique solution, which is continuously differentiable in time and also has properties C 1 -C 4 .

In Section 5 we present two different numerical methods for system (10): first we solve it using the explicit Euler method via the Elsgolts approach [10] , and later positivity-preserving Runge-Kutta methods [20, 21, 22] .

One of the key elements in the solution of delay differential equations is the fact that the discontinuities should be included in the mesh of the time discretisation. Since our history function is smooth, the only discontinuities in the higher order derivatives might appear at the points kσ, k ∈ N. Based on this, we define them as

where m is a positive integer. On this above mesh, we can define the scheme

where 0 ≤ n ≤ m T σ and τ = σ m . The symbol • is the element-by-element or Hadamard product of the matrices and the (k, l) element of the matrix X n (X ∈ {S, I, R}) is the approximation of the value X k,l (t n/m ), moreover the (k, l) element of T n−m gives the approximation ofF(t n/m − σ, x k , y l ) by interpolating the elements of I n−m . Instead of analyzing the previous numerical method in terms of its convergence, we are going to observe how well the model describes the real-life processes: more precisely, whether our model preserves the discrete versions of qualitative properties C 1 -C 4 . From now on, we denote these by D 1 -D 4 . Now we prove that for a sufficiently small time-step (or in other words, a sufficiently large m) properties D 1 -D 4 hold for n > 0.

Theorem 5. Suppose that properties D 1 -D 4 hold for the history functions discretized on the grids G and G t . Property D 2 holds without restrictions. Furthermore, if we assume that the time step satisfies

in whichT

and M = max We proceed through induction in n: first we prove that the properties hold for the first few steps, then we prove that if the properties are true up to some step n then they also hold for step n + 1. Since the proof for some arbitrary step and the first few ones are very similar, we only present here the one for the latter. Note that the aforementioned properties should be proved for every element in the matrices S n+1 , I n+1 and R n+1 , but since they are all similar, here we present the proof for an arbitrary element.

By adding up all the equations of (11) we get property D 2 . The first equation in the scheme can also be rewritten as Properties C 1 and C 3 hold iff

The higher bound is true if T n−m k,l ≥ 0 but this holds because of either the assumption of the induction condition or the first assumption of the theorem. Also, the inequality of the lower bound can be rephrased as

It is easy to see that because of the construction ofT and properties D 1 and D 2 (which hold at t (n−m)/m by the induction condition),

which means that if (12) holds then property D 1 holds for S n+1 . Also, since S n , I n and T n−m all have non-negative elements, D 1 holds also for I n+1 if (12) holds. By similar considerations, D 1 and D 4 also hold for R n+1 , which proves the statement.

An important remark is that the step size must be in the form σ m , which means that by condition (12), the theoretically best step size is in the form σ m , wherẽ

To achieve higher order convergence in our numerical approximation, we need to use higher order methods. One of the most widely used ones are the Runge-Kutta methods. Since in this paper we consider a constant delay, these methods are easily applicable. In the sequel, the investigations will be based on the Shu-Osher form. Now we introduce the necessary notations used in the case of the classical ordinary differential equations, and then generalize them for the delay equations.

Consider a time-dependent problem in the form

and a Runge-Kutta method given in the Butcher form [4] with coefficients (a ij ) ∈ R s×s and b ∈ R s . Let B be the following matrix containing all the aforementioned coefficients of the method:

in which 0 is the zero vector with length s (the number of the stages in the method). Let us also denote the (s + 1)-dimensional identity matrix by I. If there is such a number r > 0 for which (I + rB) is invertible, then the explicit Runge-Kutta method for equation (14) can be expressed in the canonical Shu-Osher form [20, 21, 22] u n

, i = 1, . . . , s + 1,

in which α ij , v i are real constants, ∆t is the time step and u n approximates u(n∆t). Moreover, we have α r = (α ij ) = r(I + rB) −1 B ∈ R (s+1)×(s+1) and v r = (v i ) = (I + rB) −1ē ∈ R s+1 (herē e = (1, 1, . . . , 1) T ). The reason for using such representation is that on the right-hand side of (15) we have the linear combinations of the steps of the explicit Euler method with time step ∆t/r. This means that the qualitative properties of the explicit Euler method may be transmitted to the Runge-Kutta methods.

Depending on the values of the parameter r we might have different representations of the same scheme. The Shu-Osher representation having the largest possible value of r for which (I + rB) −1 exists and α r and v r have non-negative components is called optimal, and we define C := max{r ≥ 0 : ∃(I + rB) −1 and α r ≥ 0, v r ≥ 0}.

C is called the SSP (strong stability preserving) coefficient and this constant will be used in the scheme (r = C). For further reading, see [13] . In our case, the problem (10) can be written in the form

Because of this, the explicit Runge-Kutta method applied to (16) takes the form

(0 ≤ n ≤ mT /σ) with the notations of Section 5.1. The next theorem states that for a sufficiently small time step, properties D 1 -D 4 hold for the above scheme. Theorem 6. Consider an explicit Runge-Kutta method in the form (17) with SSP coefficient C > 0 and applied to system (10) with non-negative history function. Then D 2 holds, and the properties D 1 , D 3 and D 4 also hold if the step-size satisfies

Proof. We prove the theorem by induction in n but since the proof for the first and an arbitrary step is similar, here we present only the latter one. Also, we are proving the required properties for a fixed k, l point, since if they hold for an arbitrary point then they also hold for all of them. To make notations simpler, in the next sections we are using the notation X n = X n k,l for X ∈ {S, I, R, T } and k, l are arbitrary fixed indices. By applying method (17) to our equation (10), we get the following scheme for an ith

Let Z n (i) = S n (i) + I n (i) + R n (i) . Then adding up equations (19) yields 2) , . . . , Z n (s+1) ) T and Z n = S n + I n + R n . It is well known that for our Runge-Kutta method to be consistent, we need that v i + s+1 j=1 α ij = 1 for all 1 ≤ i ≤ s + 1 and consequently (I − α C ) −1 v C =ē, which shows that Z n (s+1) = Z n . Then because of Z n (s+1) = Z n+1 , we get that Z n+1 = S n+1 + I n+1 + R n+1 = S n + I n + R n = Z n , ∀n, so property D 2 holds. The remaining three properties require that S n+1 , I n+1 and R n+1 are non-negative, and S n+1 is non-increasing while R n+1 is non-decreasing. Let us rewrite the first two equations in the internal stage (19) as

By the definition of T n−m and by the assumption of the induction, and since we are using non-negative cubature and positivity-preserving interpolation, T n−m is also non-negative. It is clear that S n (1) = v 1 S n , I n (1) = v 1 I n and R n (1) = v 1 R n (by the Shu-Osher form of explicit Runge-Kutta methods) and they are non-negative. Also, from the form of (20) it is clear that the non-negative property of S n (i) and I n (i) (and also of R n (i) ) holds for every 2 ≤ i ≤ s + 1 if the condition

holds. The first part of this condition is satisfied because of the second term of the righthand side of (18) . For the second part of (21) , observe that by the definition ofT in (13), we know that T n−m ≤T .

Therefore the second part of condition (21) is also satisfied because of (18), so D 1 holds. Moreover, since T n−m ≥ 0, we get that

With this, we get the following estimate for S n (i) from the first equation of (20)

As mentioned before, to achieve consistency we need

By rearranging, we have

that is S n (p) ≤ S n . Consequently, S n (i) ≤ S n for all 1 ≤ i ≤ s + 1 and also S n+1 = S n (s+1) ≤ S n which gives property D 3 .

By the third equation of (19) , it is evident that R n (i) ≥ R n , hence R n+1 = R n (s+1) ≥ R n , which gives property D 4 .

We should also note here that like in the case of the explicit Euler method, we cannot use arbitrary values for time steps, but they should be in the form σ/m. Therefore, the theoretically best time step is σ/m, wherẽ m = min m ∈ N + σ m < C min 1 T + c , 1 b .

In this section we present some numerical experiments to confirm our previous results. First we show that the bounds we got in Theorems 5 and 6 are sharp in the sense that the use of bigger time steps results in qualitatively bad behavior. Then, in the second part we present some graphs on which the solutions are compared for different values of σ and their qualitative properties are checked.

The equation (2) is solved on the rectangle domain Ω = (0, A) × (0, B) (A, B ∈ R + ). We set the parameters as A = B = 1 and c = 0.01. In order to be able to define the F I (t − σ, x, y) function, we choose the weight function

In the numerical experiments, the choice a = 100 is used. The radius δ, delay-parameter σ and the rate of recovery b will be set later. The history functions are chosen as

for t ∈ [−σ, 0], where I h (t, x, y) is a scaled Gaussian distribution with standard deviation s = 1/10 concentrated at the middle of the domain (1/2, 1/2). Note that since I h is monotone increasing and continuous in time, continuous in space, is non-negative and 1/(2πs 2 ) ≈ 15.92 < 20, functions (23) fulfill properties C 1 -C 4 . The graphs of the history functions at time t = 0 can be seen in Figure 1 . The semi-discretization is carried out on a standard rectangular mesh with step sizes h x = h y = 1/19. As mentioned before, we can use different cubatures to approximate the integral F I (t, x, y) (see [24] ) -we define the cubature as follows. First, we transform the disklike infection domain with radius δ to the rectangle [0, δ] × [0, 2π) on the plane (r, ϑ) using a polar transformation with x = x + r cos ϑ and y = y + r sin ϑ. Then, we transform this rectangle to the square [0, 1] × [0, 1) on the plane (r , θ ) by using the linear transformation r = δr and ϑ = 2πϑ with the Jacobian determinant 2πδ. Using the transformations above, the integral F I (t − σ, x, y) has the form 1 0 1 0 a(−r δ + δ)I(t − σ, x + r δ cos(2πϑ ), y + r δ sin(2πϑ )) r 2πδ 2 dr dϑ .

For the integration over the interior of the aforementioned square, we take the 40-point generalised Gaussian quadrature rule [17] with the quadrature points by µ 1 , µ 2 , . . . , µ 40 and corresponding weights ω 1 , ω 2 , . . . , ω 40 . Therefore, the cubature has the form

where i = 40(j − 1) + l, η i = µ j δ cos(2πµ l ), ξ i = µ j δ sin(2πµ l ) and w i = ω j ω l 2πδ 2 µ j . Based on this, at the given spatial grid point (x k , y l ), the approximation

is used, whereÎ(t−σ, x k +η i , y l +ξ i ) is computed using piecewise cubic Hermite interpolation.

With the help of the previous constructions, the time discretization methods discussed in Section 5 now can be applied (see the next subsections).

In the previous sections, namely in Theorems 5 and 6 we gave sufficient conditions for the qualitatively good behavior of the numerical solution, i.e. if we use a smaller time step than the bounds, then our numerical solution possesses properties D 1 -D 4 . A natural question which might arise in the case of such sufficient conditions is the effect of the use of bigger time steps.

In Table 1 , we can see the theoretical bound (12) τ ≤ min 

denoted by "theor. b.", the actual time step in the form σ/m (see the end of Section 5.1) denoted by "time step", and the bound calculated by experiments, i.e. the time step in the form σ/m exp for which the method works as expected but for σ/(m exp − 1) it gives qualitatively inaccurate results -this value is denoted by "real b." in the table. Also, the difference m exp −m is denoted by "diff." in the table. The bound can be considered sharp when this value is zero. The last column shows the ratio of the "time step" and the "real bound", i.e. the sharpness of the bound we got from our theorem. As it can be seen, this ratio is 1 for several parameter values. and b = 0.05 (the other parameters are given before).

We show an example for the qualitatively bad behaviour of the method and show the sharpness of the obtained time step bound (24) in the second row of Table 1 . On the left panel of Figure 2 , the numerical solution S n can be seen at the time level t = 3 obtained with time step τ = σ/4 = 1/4 = 0.25. We can see that the solution is qualitatively correct, namely, the values are non-negative. But when we use the next possible time step τ = σ/3 = 1/3 = 0.3333, we also get negative values. On the right panel of Figure 2 the white area corresponds to those grid points at which the solution is negative. Thus, the obtained bound is sharp.

The results for the second order Runge-Kutta method are presented in Table 2 . Here we can see that the theoretical bound is not that sharp in all of the cases, so our condition is only sufficient, but not necessary. However, in some cases the theory gives time steps that are not far from the best possible one (when we have small numbers in the 'diff.' column). For different choices of the parameters and initial conditions, we might get even better results, resulting in a sharp bound even in the higher order case.

In this section, we present some graphs of the numerical solutions of equation (2) with different values of σ. We are plotting the solution at final time T = 7 with parameters b = 0.1 and δ = 0.1, computed with the possible largest time step below the theoretical bound 0.4752 (computed similarly as (24) with C = 1) and the second order Runge-Kutta method is used. As we can see in Figure 3 , the increase of parameter σ results in a slower spread of the infection, which corresponds to the biological requirements (a longer latent period results in a slower pandemic).

In this article we extended our previous works concerning space-dependent SIR models with the addition of constant delay to the equation. We could show using the method of steps that this new equation has a unique solution with biologically reasonable properties. Then, it was shown that numerical methods with carefully chosen step sizes will preserve the discrete versions of the aforementioned features. One remaining question concerns the use of time steps not in the form σ/m, (m ∈ N + ). In these cases, an additional interpolation is needed to compute the values of our functions at times t n/m − σ. However, if we use positivity-preserving interpolations (like the method 'pchip' mentioned in Section 4.2) then similar statements can be formulated as Theorems 5 or 6, although the computational time increases considerably.

Another possible extension of the previous methods include the introduction of nonconstant delay, i.e. given by a time-or space-dependent function. In these cases the addition of an interpolation step to our algorithm is needed, and a more careful choice of time steps is also required.

A convection-diffusion model for gang territoriality

On the computational solution of differential-difference equations

Epidemiological Modeling of the 2005 French Riots: A Spreading Wave and the Role of Contagion

Numerical methods for ordinary differential equations

Mathematical Structures of Epidemic Systems

Stability analysis for a vector disease model

Operator splitting for space-dependent epidemic model

Nonnegativity-, monotonicity-, or convexity-preserving cubic and quintic Hermite interpolation

Qualitative properties of spatial epidemiological models

Qualitative Methods in Mathematical Analysis

Qualitative properties of some discrete models of disease propagation

Monotone piecewise cubic interpolation

Strong stability preserving Runge-Kutta and multistep time discretizations

Measles periodicity and community size

Mathematical models of the spread of infection

A contribution to the mathematical theory of epidemics

Generalized Gaussian cubature rules for systems of arbitrary functions

Global stability of an SIR epidemic model with time delay

Simulation and separation by principal components of multiple demic expansions in Europe

Total-variation-diminishing time discretizations

Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws

Efficient implementation of essentially nonoscillatory shockcapturing schemes

Space dependent models for studying the spread of some diseases

High Order discretization methods for spatial dependent SIR models

Forecasting elections using compartmental models of infections

Global attractivity, spreading speeds and traveling waves of delayed nonlocal reaction-diffusion systems

Global stability of a SIR epidemic model with nonlinear incidence rate and time delay

The research by the authors B.M.T., I.F. and R.H. reported in this paper and carried out at BME has been supported by the NRDI Fund (TKP2020 NC, Grant No. BME-NC) based on the charter of bolster issued by the NRDI Office under the auspices of the Ministry for Innovation and Technology, and the Hungarian Ministry of Human Capacities OTKA grant SNN125119.The work of the author I.F. was completed in the ELTE Institutional Excellence Program (TKP2020-IKA-05) financed by the Hungarian Ministry of Human Capacities.The research of the author D.R. reported in this paper was supported by the Slovenian Research Agency grants P1-0292, N1-0114, N1-0083, N1-0064, and J1-8131.