A non-standard finite-difference-method for a non-autonomous epidemiological model: analysis, parameter identification and applications

1.   Introduction

1.1.   Introductory motivation and literature overview

Differential equations, especially those possessing non-negative solutions, play an important role in many areas of sciences such as biology ^[1,2,3,4,5], chemistry ^[6,7,8], epidemiology ^[9,10,11] or population dynamics ^[12,13,14] to name a few of these fields. Since we focus on examining a time-discrete, mathematical model of epidemiology in our work, we limit ourselves to a concise review of literature in this area.

In this article, we want to consider the time-continuous non-autonomous susceptible-infected-recovered (SIR) model

in epidemiology thoroughly investigated by Wacker and Schlüter in 2020 ^[15]. We provide further assumptions regarding this time-continuous model in Section 2. In 1927, one of the most famous papers in mathematical epidemiology, written by Kermack and McKendrick, was published ^[16]. Previously, Ross and Hudson proposed important foundations for epidemiological models ^[17,18]. Additional well-written books and review articles on mathematical epidemiology are authored by Hethcote ^[19] and Brauer and Castillo-Chávez ^[20]. Since the establishment of the classical SIR model, many modifications have been suggested and, due to the on-going COVID-19 pandemic, many articles have been published with respect to this topic and we can only summarize a brief selection of these works.

To this day, there are works which provide closed solution formulas for simple models in epidemiology ^[21,22,23]. However, more articles consider numerical approximations due to appearing nonlinearities. Especially during the COVID-19 pandemic, different approaches have been published ^[15,24]. In 2022, Chen and co-authors suggested a susceptible-exposed-infected-unreported-recovered (SEIUR) model for the spread of COVID-19 with piecewise constant dynamical parameters ^[25] as one possible extension of classical SIR models. Taking age- and sex-structure of a population into account, Wacker and Schlüter proposed an extension of classical SIR models by multiple groups which can be considered as, for example, age groups or locally related groups ^[26]. Using fractional differential operators, Xu and co-authors analyzed the spread of COVID-19 in 2020 with respect to the United States by means of a fractional-order susceptible-exposed-infected-recovered (SEIR) model ^[27]. Regarding parameter identification problems in epidemiological models, Marinov and Marinova designed adaptive SIR models to solve these inverse problems ^[28,29]. In context of inverse problems, Comunian, Gaburro and Giudici concluded that high-quality data are needed to calibrate constant parameters of epidemiological models ^[30]. Other interesting lines of research concern the consideration of epidemics on networks ^{[31,32,33,34,35]} or application of deep-learning approaches to, for example, forecast time-series data of epidemics by these methods ^[36,37,38].

1.2.   Article's motivation and contributions

One important aspect in epidemiological models is preservation of positivity. It plays a huge role in different scientific problems such as wave equations ^[39], simulations of chemical reactions ^[7], population dynamics ^[12,13,14] or other biological systems ^[3]. However, many classical time discretization schemes suffer from conserving non-negativity for arbitraty time step sizes ^[40]. Hence, designing non-negativity preserving numerical solution schemes for such problems can be regarded as an art as stated by Mickens ^[41]. For these reasons, Wacker and Schlüter examined a numerical solution scheme

for (1.1) based on implicit Eulerian time discretization in ^[15] whose algorithm is presented in Algorithm 1.

Such non-standard finite-difference-methods for time discretization are also important in other application areas such as magnetohydrodynamics ^[42]. Here, we use $f \left(t_{j} \right) : = f_{j}$ as an abbreviation for time-continuous functions at time points $t_{j}$ and $f_{h} \left(t_{j} \right) : = f_{h, j}$ as an abbreviation for time-discrete functions at time time points $t_{j}$ on possibly non-equidistant meshes $t_{1} < t_{2} < \ldots < t_{M - 1} < t_{M}$. Although Algorithm 1 preserves all desired properties of the time-continuous model (1.1) as shown in ^[15], there are still some disadvantages:

$1)$ From its formulation (1.1), there is to solve a quadratic equation

where the a-priori knowledge of non-negative solutions from the time-continuous model is used to rule out one of both possible solutions for Algorithm 1 to obtain unique solvability;

$2)$ In each time step, we must recalculate the auxiliary quantities $A_{j + 1}$ and $B_{j + 1}$;

$3)$ By both aforementioned arguments, we conclude that Algorithm 1 does not intrinsically calculate biologically correct solutions and, from an algorithm point of view, recalculation of these auxiliary quantities is undesirable.

We want to propose a different discretization by non-standard finite-difference-methods based on non-local approximations. To explain the idea of non-standard finite-difference-methods based on non-local approximations, we consider a time-continuous dynamical system

Standard time-discretizations normally approximate the right-hand-side function by function values of $\textbf{y} \left(t \right)$ at one previous or the current time point. Non-standard finite-difference-methods based on non-local approximations, on the contrary, approximate $\textbf{f} \left(t, \textbf{y} \left(t \right) \right)$ by function values from more then solely one previous or the current time point. For further details, we refer to the book ^[41] by Mickens. For the aforementioned three reasons, we suggest a non-standard finite-difference-method (NSFDM) given by

in Section 3 as our first main contribution, after summarizing analytical properties of the time-continuous problem formulation (1.1) in Section 2. In Section 3, we additionally prove unconditional preservation of non-negativity and linear convergence of (1.3) towards the time-continuous solution of (1.1). Additionally, (1.3) transfers to a uniquely solvable algorithm which preserves non-negativity intrinsically. Based on (1.3), we design a numerical algorithm to solve the inverse problem of finding the time-varying parameter functions in Section 4 as our next main contribution. Finally, we strengthen our theoretical findings in Section 5 by providing some numerical examples, especially one numerical example where linear convergence is stressed, in contrast to ^[15], where linear convergence was not shown numerically.

2.   Time-continuous problem formulation

Following ^[15], we take the following assumptions into consideration with regard to our time-continuous model (1.1):

$(1).$ The total population size $N$ is fixed for all time $t \geq 0$, i.e., it holds

for all $t \geq 0$ and we assume $I_{0} > 0$ to really get an evolution in time of the infected subgroup of the population;

$(2).$ We divide the total population into three homogeneous subgroups, namely susceptible people (S), infected people (I) and recovered people (R). The equation

is valid for all $t \geq 0$;

$(3).$ The time-varying transmission rate $\alpha \colon \left[0, \infty \right) \longrightarrow \left[0, \infty \right)$ is at least continuously differentiable once and there exist positive constants $\alpha_{ \text{min}} > 0$ and $\alpha_{ \text{max}} > 0$ such that

holds for all $t \geq 0$;

$(4).$ The time-varying recovery rate $\beta \colon \left[0, \infty \right) \longrightarrow \left[0, \infty \right)$ is at least continuously differentiable once and there exist positive constants $\beta_{ \text{min}} > 0$ and $\beta_{ \text{max}} > 0$ such that

holds for all $t \geq 0$.

We summarize the following analytical properties of the time-continuous problem formulation (1.1) from article ^[15]:

Theorem 1. Let all assumptions from (2.1)–(2.4) for our time-continuous model (1.1) be fulfilled. The following statements hold true:

(1). Possible solutions of (1.1) remain non-negative and bounded above, i.e.,

hold for all $t \geq 0$;

(2). The time-continuous problem formulation (1.1) possesses exactly one globally unique solution for all $t \geq 0$;

(3). First, $S$ is monotonically decreasing and there exists $S^{\star} > 0$ such that $\lim\limits_{t \to \infty} S \left(t \right) = S^{\star}$. Secondly, $R$ is monotonically increasing and there exists $R^{\star} > 0$ such that $\lim\limits_{t \to \infty} R \left(t \right) = R^{\star}$. Finally, it holds $\lim\limits_{t \to \infty} I \left(t \right) = 0$.

Proof. $(1).$ This statement's proof can be found in [15, Theorem 4].

$(2).$ Globally unique existence in time is proven in [15, Theorem 6].

$(3).$ Proof of the solution components' monotonicity and its long-time behavior is given in [15, Theorem 8].

We want to note that $S^{\star}$ and $R^{\star}$ depend on the initial conditions.

3.   A non-standard finite-difference-method

Let $0 = t_{1} < t_{2} < \ldots < t_{M} = T$ be an arbitrary and possibly non-equidistant partition of the simulation interval $\left[0, T \right]$ with final simulation time $T$. As abbreviations, we write $f \left(t_{j} \right) : = f_{j}$ for time-continuous functions and $f_{h} \left(t_{j} \right) : = f_{h, j}$ for time-discrete functions at time points $t_{j}$ for all $j \in \left\{ 1, 2, \ldots, M \right\}$. The time-discrete problem formulation of our time-continuous model (1.1) is given by (1.3) for all $j \in \left\{ 1, 2, \ldots, M - 1 \right\}$. Here, $h_{j + 1} : = t_{j + 1} - t_{j}$ denote non-equidistant time step sizes for all $j \in \left\{ 1, 2, \ldots, M - 1 \right\}$. Additionally, the supremum norm on $\left[0, T \right]$ for time-continuous functions is denoted by

while the same supremum norm on $\left[t_{p}, t_{p + 1} \right]$ reads

3.1.   Conservation of total population size $N$

As our first result, our time-discrete problem formulation (1.3) conserves total population size $N$ for all time points $t_{j} \in \left[0, T \right]$ with $j \in \left\{ 1, 2, \ldots, M \right\}$.

Theorem 2. Possible solutions of (1.3) fulfill

for all $j \in \left\{ 1, 2, \ldots, M \right\}$.

Proof. This statement can be proven by induction. For $j = 1$, it is clear by our initial conditions that

holds. Let us assume that

is valid for an arbitrary $j \in \left\{ 1, 2, \ldots, M - 1 \right\}$. Adding all three equations of (1.3), we obtain

and this implies

which proves our assertion (3.1).

3.2.   Boundedness and non-negativity

As our next result, we show that approximate solutions of our time-continuous model (1.1) from our discretization model (1.3) remain non-negative and bounded.

Theorem 3. Possible solutions of (1.3) remain non-negative and bounded for all $t_{j} \in \left[0, T \right]$ with arbitrary $j \in \left\{ 1, 2, \ldots, M \right\}$, i.e., these solultions fulfill

for all $j \in \left\{ 1, 2, \ldots, M \right\}$.

Proof. $(1).$ By the first three equations of (1.3), we obtain

for all $j \in \left\{ 1, 2, \ldots, M - 1 \right\}$,

for all $j \in \left\{ 1, 2, \ldots, M - 1 \right\}$ and

for all $j \in \left\{ 1, 2, \ldots, M - 1 \right\}$.

$(2).$ By our assumption of non-negative initial values $S_{h, 1}, I_{h, 1}, R_{h, 1}$ and by the three equations (3.3)–(3.5), we conclude that

hold for all $j \in \left\{ 1, 2, \ldots, M \right\}$.

$(3).$ By non-negativity and conservation of total population size (3.1) from Theorem 2, it follows that

hold for all $j \in \left\{ 1, 2, \ldots, M \right\}$.

$(4).$ Hence, all possible solutions of (1.3) remain non-negative and bounded.

3.3.   Unique solvability

It directly follows that our numerical solution scheme (1.3) is uniquely solvable for all $j \in \left\{ 1, 2, \ldots, M \right\}$.

Theorem 4. The discretization scheme (1.3), based on the explicit solutions (3.3)–(3.5), is well-defined and uniquely solvable for all $j \in \left\{ 1, 2, \ldots, M \right\}$.

Proof. (3.3)–(3.5) prove that our numerical solution scheme (1.3) is well-defined. Unique solvability is a direct consequence of these equations as well. Hence, our assertion is proven.

Remark 1. We briefly want to note why we call solutions (3.3)–(3.5) explicit. It obviously holds

for all $j \in \left\{ 1, 2, \ldots, M \right\}$ explicitly. Plugging (3.6) into (3.4) yields

for all $j \in \left\{ 1, 2, \ldots, M \right\}$ as an explicit formula. Finally, by using (3.7), we get

for all $j \in \left\{ 1, 2, \ldots, M \right\}$. For these reasons, we talk about an explicit numerical solution scheme.

3.4.   Monotonicity properties and long-time behavior

We show that monotonicity properties and long-time behavior of the time-continuous problem formulation (1.1) translates to our time discretization (1.3).

Theorem 5. 1). The sequence $\left\{ S_{h, j} \right\}_{j = 1}^{M}$ is monotonically decreasing and there exists a non-negative constant $S_{h}^{\star}$ such that $\lim\limits_{j \to \infty} S_{h, j} = S_{h}^{\star}$ holds.

2). The sequence $\left\{ R_{h, j} \right\}_{j = 1}^{M}$ is monotonically increasing and there exists a non-negative constant $R_{h}^{\star}$ such that $\lim\limits_{j \to \infty} R_{h, j} = R_{h}^{\star}$ holds.

3). It holds $\lim\limits_{j \to \infty} I_{h, j} = 0$ if also $T \to \infty$ and $h_{j + 1} \not = 0$ are assumed.

Proof. $(1).$ Since $\left\{ S_{h, j} \right\}_{j = 1}^{M}$ is non-negative and bounded above by the total population size $N$, the inequality

from (3.6) implies that this sequence is also monotonically decreasing which proves our first claim.

$(2).$ Since $\left\{ R_{h, j} \right\}_{j = 1}^{M}$ is non-negative and bounded above by the total population size $N$, the inequality

from (3.5) yields that this sequence is monotonically increasing which shows our second assertion.

$(3).$ We conclude

from (3.5) and this shows that our third claim is also true.

Remark 2. Both Theorems 1 and 5 imply that our time-continuous problem formulation (1.1) and our time-discrete problem formulation (1.3) converge towards the disease-free equilibrium state in the long run.

3.5.   Error analysis

We start this subsection by listing all our assumptions for our main result regarding linear convergence of the time-discrete solution towards the time-continuous one:

($A_{1}1$) We consider the time interval $\left[0, T \right]$ where

holds;

($A_{1}2$) The initial conditions of the time-continuous and of the time-discrete problems coincide;

($A_{1}3$) The time-continuous solution functions $S, I, R \, \colon \, \left[0, T \right] \longrightarrow \mathbb{R}$ should be twice continuously differentiable;

($A_{1}4$) The time-varying transmission rate $\alpha \, \colon \, \left[0, T \right] \longrightarrow \left[0, \infty \right)$ and the time-varying recovery rate $\beta \, \colon \, \left[0, T \right] \longrightarrow \left[0, \infty \right)$ should be once continuously differentiable;

($A_{1}5$) The time-varying transmission rate and the time-varying recovery rate should be bounded, i.e., there are non-negative constants $\alpha_{ \text{min}}, \alpha_{ \text{max}}, \beta_{ \text{min}}, \beta_{ \text{max}}$ such that

and

hold for all $t \in \left[0, T \right]$;

($A_{1}6$) Set $h : = \max\limits_{p \in \mathbb{N}} h_{p}$.

Theorem 6. Let all assumptions ($A_{1}1$)–($A_{1}6$) be fulfilled. Then the time-discrete solution convergences linearly towards the time-continuous solution if $h \to 0$.

Proof. We shortly explain our proof's strategy because it is relatively technical. If we first assume that our time-discrete solution coincides with the time-continuous solution at a certain time point $t_{p} \in \left[0, T \right]$ for an arbitrary $p \in \left\{ 1, 2, \ldots, M - 1 \right\}$, then we investigate the local error propagation to the next time point $t_{p + 1}$. Secondly, we consider error propagation in time and finally, we examine accumulation of these errors. We denote time-discrete solutions, for example, by $S_{h, p}$ and time-continuous solutions by $S \left(t_{p} \right)$ at the same time point $t_{p}$.

1) For investigation of local errors, we assume that

hold for an arbitrary $p \in \left\{ 1, 2, \ldots, M - 1 \right\}$ on the time-interval $\left[t_{p}, t_{p + 1} \right]$. Since we only take one time step into account, we denote corresponding time-discrete solutions by $\widetilde{S_{h, p + 1}}$, $\widetilde{I_{h, p + 1}}$ and $\widetilde{R_{h, p + 1}}$.

1.1) We have by (3.6) from Remark 1 that

holds. This implies

By the mean value Theorem, there exists $\xi_{S, 1} \in \left(t_{p}, t_{p + 1} \right)$ such that we obtain

This yields

For $II_{S, 1}$, we notice that

and hence,

holds due to boundedness of time-continuous solutions. Using again the mean value theorem for $III_{S, 1}$, we see that there is $\xi_{\alpha, 1} \in \left(t_{p}, t_{p + 1} \right)$ such that we obtain

Hence, we conclude

for $III_{S, 1}$. Plugging (3.15) into (3.13) yields

Thus, we finally obtain

1.2) By (3.7) from Remark 1, we get

Hence, it follows

By the mean value theorem, there exists $\xi_{I, 1} \in \left(t_{p}, t_{p + 1} \right)$ such that

is valid. This implies

For $II_{I, 1}$, we observe

and can conclude

By boundedness and the mean value theorem, we receive the following estimates

for $VI_{I, 1}$, $V_{I, 1}$ and $IV_{I, 1}$. By boundedness, the triangle inequality and the mean value theorem, we can deduce

and consequently, this implies

Hence, we conclude that

holds. By defining

we finally deduce

from (3.27).

1.3) By (3.8) from Remark 1, we get

Hence, it follows

By using similar argument as in previous steps with the help of the mean value theorem and (3.29), we can deduce

1.4) By setting

we conclude

2) Normally, the points $\left(t_{p}, S_{p} \right)$, $\left(t_{p}, I_{p} \right)$ and $\left(t_{p}, R_{p} \right)$ do not exactly lie on the graph of the time-continuous solution. We then must investigate how procedural errors such as $\left| S_{h, p} - S \left(t_{p} \right) \right|$, $\left| I_{h, p} - I \left(t_{p} \right) \right|$ or $\left| R_{h, p} - R \left(t_{p} \right) \right|$ propagate to the $\left(p + 1 \right)$-th time step.

2.1) We observe that

holds where $I_{S, 2}$ is given and estimated by

We have

and

Plugging (3.37) and (3.38) into (3.36) yields

Using (3.35), we obtain

2.2) We see that

holds and get

Plugging (3.40) into (3.42) yields

2.3) We obtain

2.4) For our summary, we define

for all $p \in \left\{ 0, 1, \ldots, M - 1 \right\}$ and

for all $p \in \left\{ 1, \ldots, M \right\}$. By (3.40), (3.43) and (3.44), we define the constants

and

Hence, it holds

3) We obtain

and inductively

which proves linear convergence of the time-discrete solution towards the time-continuous solution.

3.6.   Algorithmic summary

Our algorithm for our proposed non-standard finite-difference-method for (1.1) is portrayed in Algorithm 2 below.

4.   A simple parameter identification algorithm

4.1.   Problem formulation and solvability

We consider the parameter identification problem (1.3) for all $j \in \left\{ 1, 2, \ldots, M - 1 \right\}$ where we take the following assumptions into account:

($A_{2}1$) The sequences $\left\{ S_{h, p} \right\}_{p = 1}^{M}$, $\left\{ I_{h, p} \right\}_{p = 1}^{M}$ and $\left\{ R_{h, p} \right\}_{p = 1}^{M}$ are given;

($A_{2}2$) It holds $0 \leq S_{h, p} \leq N$ for all $p \in \left\{ 1, \ldots, M \right\}$ and the sequence $\left\{ S_{h, p} \right\}_{p = 1}^{M}$ is monotonically decreasing;

($A_{2}3$) It holds $0 \leq I_{h, p} \leq N$ for all $p \in \left\{ 1, \ldots, M \right\}$;

($A_{2}4$) It holds $0 \leq R_{h, p} \leq N$ for all $p \in \left\{ 1, \ldots, M \right\}$ and the sequence $\left\{ R_{h, p} \right\}_{p = 1}^{M}$ is monotonically increasing;

($A_{2}5$) We assume $S_{h, p} \not = 0$, $I_{h, p} \not = 0$ and $R_{h, p} \not = 0$ for all $p \in \left\{ 1, \ldots, M \right\}$.

Our idea now is to reduce (1.3) to

for all $p \in \left\{ 1, \ldots, M - 1 \right\}$. Hence, we obtain the sought time-varying transmission rate and the sought time-varying recovery rate by

for all $p \in \left\{ 1, \ldots, M - 1 \right\}$. Finally, we have explicit formulas for the time-varying transmission and recovery rates at every time point which we can evaluate by the known data from the numbers of susceptible, infected and recovered people. We have the following theorem.

Theorem 7. Let all assumptions ($A_{2}1$)–($A_{2}5$) be fulfilled. (4.1) for parameter identification is uniquely solvable for all $p \in \left\{ 1, \ldots, M - 1 \right\}$ and it holds $\alpha_{p + 1} \geq 0$ and $\beta_{p + 1} \geq 0$ for all $p \in \left\{ 1, \ldots, M - 1 \right\}$.

Proof. 1) Unique solvability is a direct consequence of (4.2).

2) $\left\{ S_{h, p} \right\}_{p = 1}^{M}$ is monotonically decreasing. Hence, $S_{h, p} - S_{h, p + 1} \geq 0$ for all $p \in \left\{ 1, \ldots, M - 1 \right\}$ and we have

3) $\left\{ R_{h, p} \right\}_{p = 1}^{M}$ is monotonically increasing. Hence, $R_{h, p + 1} - R_{h, p} \geq 0$ for all $p \in \left\{ 1, \ldots, M - 1 \right\}$ and we have

This finishes our claim's proof.

4.2.   Algorithmic summary

Our algorithmic summary for our proposed parameter identification problem (4.2) reads:

5.   Numerical examples

In this section, we want to strengthen our theoretical findings by two numerical examples. We apply GNU OCTAVE for our implementation ^[43].

5.1.   Example 1: Theoretical results

The initial conditions for our first example read $S \left(0 \right) = 5000$, $I \left(0 \right) = 5000$, $R \left(0 \right) = 0$, $N = 10,000$, $\alpha = 0.5$, $\beta = 0.1$ and our simulation interval is given by $\left[0, T \right]$ with $T = 70$. We compare three numerical algorithms, namely our proposed non-standard finite-difference-method (NSFDM), the explicit Eulerian scheme (EE) and an explicit Runge-Kutta scheme with two stages (RK2) and use equidistant time step sizes $h$ for our computations.

As it can be clearly seen in Figures 1 and 2, only our novel non-standard finite-difference-method is unconditionally non-negativity-preserving and stable. Both explicit standard time-discretizations suffer from time-step restrictions to obtain non-negativity preservation. In Figure 3, we observe that all three numerical algorithms possess same properties for smaller time step sizes $h$. On the positive side for all algorithms, they conserve total population size $N$ as depicted in Figure 4.

Finally, underlining our convergence result from Theorem 6, we compare results of our non-standard finite-difference-method (NSFDM) for different time step sizes $h$ with a fine scale Runge-Kutta scheme of order 4 (RK4) at $T = 70$ where a fine time step size of $0.0001$ is applied for (RK4). The errors

● Error for $S \left(t \right)$: $\left| S_{ \text{NSFDM}, M} - S_{ \text{RK4}, M} \right|$;

● Error for $I \left(t \right)$: $\left| I_{ \text{NSFDM}, M} - I_{ \text{RK4}, M} \right|$;

● Error for $R \left(t \right)$: $\left| R_{ \text{NSFDM}, M} - R_{ \text{RK4}, M} \right|$

are used for this purpose. It can be clearly seen from Table 1 that the theoretical convergence order is achieved by our non-standard finite-difference-method.

5.2.   Example 2: COVID-19 in Spain (01 March 2020 to 15 May 2020)

We consider real-world data from Spain ^[44]. $N = 47,370,000$ is the total population size of Spain. The cumulative cases of infected people are modified such

hold from the given data. We take $\alpha \left(t \right) = 0.52 \cdot \text{e}^{-0.03 \cdot t}$ and $\beta \left(t \right) = 0.045$. For our simulation, we additionally assume $h = 0.75$ as an equidistant time step size.

Our estimated transmission rate $\alpha \left(t \right)$ and constant recovery rate $\beta \left(t \right)$ are user-chosen and adapted to the data in Figures 6 and 7, while the modified data for infected and recovered people are depicted in Figure 5. We solely consider the beginning of the pandemic because, later, vaccines have been available and reinfections have occurred after virus modifications. In Figure 8, we see results of our numerical simulations with $h = 0.75$. Results are relatively sensitive with respect to time step sizes and estimates rates, e.g., time-varying transmission rates. Hence, we can only conclude that good data acquisition and different predictions in different settings are crucial for future pandemics to make predictions for courses of such future epidemics.

6.   Concluding remarks

In this article, we first summarized some analytical properties of a time-continuous epidemiological model (1.1) as given in ^[15]. Afterwards, we designed a new non-standard finite-difference-method (1.3) as a time-discrete version of the continuous model. In Section 3, we contributed our main results on our proposed non-standard finite-difference-method (1.3). Namely, we investigated unconditionally non-negativity preservation, conservation of total population size, convergence towards a disease-free equilibrium point and linear convergence by letting the time size $h \to 0$. Conclusively for our theoretical findings, we developped a simple parameter identification algorithm in Section (4). In Section 5, we gave two numerical examples to underline our theoretical findings. Comparing Algorithms 1 and 2, we conclude that the latter one is definitely easier to implement and it can serve as a starting point for constructing higher order time discretization schemes for (1.1). Regarding practical applications, our proposed time-discretization by a non-standard finite-difference-method preserves all desired properties from the time-continuous model which are expected from our real-word intuition. Additionally, the simple parameter identification algorithm produce fast graphs from given data where trends in time-varying transmission rates such as contact restrictions might be seen from.

There might be different directions to extend our work:

● It seems straightforward to integrate further compartments into our model such that we can additionally investigate influences such as, for example, vaccines or the simple extension to the SEIR model. For example, we readily can construct

as one variant of a non-standard finite-difference-method for the SEIR model and it follows directly that this scheme conserves total population size $N$ and can be written in an explicit fashion as our proposed non-standard finite-difference-method for our considered SIR model by

to see that it also preserves non-negativity assuming positive initial conditions. One disadvantage is that there are different possibilities to construct such non-standard finite-difference-methods by non-local approximations and a suitable choice depends on the time-continuous properties which should be transferred to the time-discrete case. However, as seen by this example, it is flexible to construct different possibilities of algorithms and they are easy to implement;

● It might be of interest to construct higher order non-standard finite-difference-methods from our first-order non-standard finite-difference-method in order to get higher accuracy although, especially, our second example demonstrated that predictions of such models need to be considered with care. One straightforward way to improve accuracy, and hence, the convergence order of numerical solution schemes is application of Richardson's extrapolation ^[45];

● One can use our discretization model for the discretization of optimal control problems for minimizing certain epidemic aspect such as eradication time ^[46].

Code and data availability

Our code and used data from ^[44] can be found under https://github.com/bewa87/2023-SIR-FirstOrder-NSFDM.

Acknowledgments

We thank the editor and both anonymous reviewers for their valuable comments which definitely helped us to improve quality and readability of our manuscript.

Conflict of interest

The authors declare there is no conflict of interest.