*3.2. Local Stability of the Equilibrium Points*

**Theorem 3.** *The disease-free equilibrium point P*<sup>0</sup> *of the system (6) is locally asymptotically stable if* R<sup>0</sup> < 1.

**Proof.** The eigenvalues of the Jacobian matrix of system (6) evaluated at point *P*0, are given as the roots of polynomial

$$(-\delta\_I - \lambda)(-\mu\_0 - \lambda)\left[\lambda^2 + \left(\delta\_I + \mathcal{C} + \frac{\beta\Lambda}{\mu\_0}\right)\lambda + \left(\delta\_I\mathcal{C} + \frac{\delta\_I\beta\Lambda}{\mu\_0} - \frac{\beta e^{-\delta\_{I\_E}\Lambda}N\delta\_I\Lambda}{\mu\_0}\right)\right] = 0.$$

Therefore, the first two eigenvalues of the Jacobian matrix evaluated at *P*<sup>0</sup> are *λ*<sup>1</sup> = −*δ<sup>I</sup>* < 0 and *λ*<sup>2</sup> = −*μ*<sup>0</sup> < 0.

Next, since R<sup>0</sup> < 1 then the coefficients of equation

$$
\lambda^2 + \left(\delta\_I + \mathbb{C} + \frac{\beta \Lambda}{\mu\_0}\right) \lambda + \frac{\delta\_I}{\mu\_0} (\mathbb{C}\mu\_0 + \beta \Lambda)(1 - \mathcal{R}\_0) = 0,\tag{23}
$$

are positives. Thus, since there is no sign change between its terms, and by Descartes' sign rule it is concluded that Equation (23) does not have positive roots. Now, if *λ* is replaced by −*λ* in (23) then

$$
\lambda^2 - \left(\delta\_I + \mathcal{C} + \frac{\beta \Lambda}{\mu\_0}\right)\lambda + \frac{\delta\_I}{\mu\_0}(\mathcal{C}\mu\_0 + \beta \Lambda)(1 - \mathcal{R}\_0) = 0. \tag{24}
$$

Thus, if R<sup>0</sup> < 1 Equation (24) has two sign changes in its terms, and by Descartes' sign rule it is concluded that there are exactly two negative roots in Equation (23). Therefore, *P*<sup>0</sup> = 0, 0, 0, <sup>Λ</sup> *μ*0 is locally asymptotically stable if R<sup>0</sup> < 1.

**Theorem 4.** *The P*<sup>∗</sup> *endemic point of the system (6) is locally asymptotically stable if* R<sup>0</sup> > 1*.*

**Proof.** We note that *<sup>R</sup>*<sup>1</sup> <sup>=</sup> *Ne*−*δIE*<sup>Δ</sup> <sup>−</sup> <sup>1</sup> <sup>&</sup>gt; 0. Thus, the characteristic equation is given by

$$\begin{pmatrix} -\delta\_I - \lambda & \frac{\mathsf{C}\mathscr{E}^{-\delta\_{\mathsf{I}}\mathsf{A}}}{R\_1} & \frac{\beta \Lambda R\_1 e^{-\delta\_{\mathsf{I}}\mathsf{A}}}{\mathsf{C}} - \mu\_0 e^{-\delta\_{\mathsf{I}}\mathsf{A}} \\\\ N\delta\_I & -\mathsf{C} - \frac{\mathsf{C}}{R\_1} - \lambda & \mu\_0 - \frac{\beta \Lambda R\_1}{\mathsf{C}} \\\\ 0 & -\frac{\mathsf{C}}{R\_1} & -\frac{\beta \Lambda R\_1}{\mathsf{C}} - \lambda \end{pmatrix} = 0$$

Therefore, the first eigenvalue of the Jacobian matrix will be *λ*<sup>1</sup> = −*δIE* < 0, and the other, are the roots of polynomial

$$
\lambda^3 + \left(\delta\_l + \mathbb{C} + \frac{\mathbb{C}}{R\_1} + \frac{\beta \Lambda R\_1}{\mathbb{C}}\right) \lambda^2 + \left(\frac{\delta\_l \beta \Lambda R\_1}{\mathbb{C}} + \beta \Lambda R\_1 + \frac{\mathbb{C} \mu\_0}{R\_1}\right) \lambda + \delta\_l (\mathbb{C} \mu\_0 + \beta \Lambda)(\mathbb{R}\_0 - 1) = 0. \tag{25}
$$

If R<sup>0</sup> > 1 all the coefficients of Equation (25) are positive, i.e., there is no sign change between their terms, and by Descartes's sign rule it is concluded that there are no positive roots. Now if *λ* is replaced by −*λ* in (25) it gives us that

$$-\lambda^3 + \left(\delta\_I + \mathbb{C} + \frac{\mathbb{C}}{R\_1} + \frac{\beta \Lambda R\_1}{\mathbb{C}}\right)\lambda^2 - \left(\frac{\delta\_I \beta \Lambda R\_1}{\mathbb{C}} + \beta \Lambda R\_1 + \frac{\mathbb{C}\mu\_0}{R\_1}\right)\lambda + \delta\_I (\mathbb{C}\mu\_0 + \beta \Lambda)(\mathbb{R}\_0 - 1) = 0\tag{26}$$

Then, if R<sup>0</sup> > 1 the polynomial (26) has three sign changes between its terms, and by Descartes's sign rule it is concluded that there are exactly three negative roots of Equation (25). Thus, *P*<sup>∗</sup> is locally asymptotically stable if R<sup>0</sup> > 1.

## *3.3. Global Stability Analysis of the Mathematical Model*

Since the variable *IE* does not appear in the other three equations, without loss of generality we will only consider the following three-dimensional system,

$$\begin{aligned} \frac{dI(t)}{dt} &= \beta T(t-\Delta)V(t-\Delta)e^{-\delta\_{\parallel\_{E}}\Delta} - \delta\_{I}I(t) \\ \frac{dV(t)}{dt} &= N\delta\_{I}I(t) - (\mathbb{C} + \beta T(t))V(t) \\ \frac{dT(t)}{dt} &= \Lambda - \beta T(t)V(t) - \mu\_{0}T(t) \end{aligned} \tag{27}$$

To analyze the global stability of the equilibrium points of the system (27), we use the method of the Lyapunov's functions, and we will use the Volterra function

$$G(\mathbf{x}) = \mathbf{x} - \mathbf{1} - \ln \mathbf{x} \tag{28}$$

for *x* > 0, which is no negative for any *x* > 0 and *G*(*x*) = 0 if and only if *x* = 1.

**Theorem 5.** *If* R<sup>0</sup> < 1 *then the disease-free equilibrium point Pf* = 0, 0, <sup>Λ</sup> *μ*0 *of system (27) is globally asymptotically stable.*

**Proof.** We define the Lyapunov functional

$$\mathcal{V}(t) = e^{-\delta\_{l\_{\mathbb{E}}}\Lambda} T^0 G\left(\frac{T(t)}{T^0}\right) + \frac{1}{N} V(t) + I(t) + e^{-\delta\_{l\_{\mathbb{E}}}\Lambda} \int\_0^\Lambda \beta T(t-\theta) V(t-\theta) \,d\theta \,d\theta$$

Now, calculating the time derivative of V(*t*) along the solution of model (27), one gets that

$$\begin{split} \frac{d\mathcal{V}(t)}{dt} &= e^{-\delta\_{\mathbb{E}}\Delta} \frac{T(t) - T^{0}}{T(t)} \frac{dT(t)}{dt} + \frac{1}{N} \frac{dV(t)}{dt} + \frac{dI(t)}{dt} + e^{-\delta\_{\mathbb{E}}\Delta} \frac{d}{dt} \int\_{0}^{\Delta} \beta T(t - \theta) V(t - \theta) \,d\theta \\ &= -\mu\_{0} e^{-\delta\_{\mathbb{E}}\Delta} \frac{\left(T(t) - T^{0}\right)^{2}}{T(t)} - e^{-\delta\_{\mathbb{E}}\Delta} \beta T(t) V(t) + e^{-\delta\_{\mathbb{E}}\Delta} \beta V(t) T^{0} + \delta\_{I} I(t) - \frac{\mathbb{C}}{N} V(t) - \frac{\beta T(t) V(t)}{N} \,d\theta \\ &+ \beta T(t - \Delta) V(t - \Delta) e^{-\delta\_{\mathbb{E}}\Delta} - \delta\_{I} I(t) - \beta T(t - \Delta) V(t - \Delta) e^{-\delta\_{\mathbb{E}}\Delta} - \delta\_{I} I(t) + e^{-\delta\_{\mathbb{E}}\Delta} \beta T(t) V(t) \\ &\leq -\mu\_{0} e^{-\delta\_{\mathbb{E}}\Delta} \frac{(T(t) - T^{0})^{2}}{T(t)} + \frac{\mathbb{C}}{N} (\Re\_{0} - 1) . \end{split}$$

Thus, *<sup>d</sup>*V(*t*) *dt* <sup>&</sup>lt; 0 when <sup>R</sup><sup>0</sup> <sup>&</sup>lt; 1. Therefore, by Lyapunov–LaSalle Invariance Principle, the infection-free equilibrium *Ef* is globally asymptotically stable.

## **4. Design of a** *NSFD* **Scheme for the Mathematical Model**

The use of differential equations in the modeling of the transmission of infectious diseases has represented a versatile tool to understand better the dynamics of a variety of infectious diseases [7,59,60,74–76]. Mathematical models based on differential equations have been useful to study how to reduce the burden of infectious diseases. The models allow the determination of optimal controls and estimate the impact of a variety of virus on the disease dynamics [67,75,77]. One advantage of mathematical models is that different simulations can be performed, and this allows us to analyze different main driven factors of epidemics under a variety of complex scenarios [8,11,59]. However, there are no general formulas that allow the obtaining of precise analytical solutions for many differential equation systems. These solutions exist only occasionally and are often difficult to find, so good approximations are necessary that preserve the qualitative properties of said solution, for which numerical methods have been used, see [24,25,31,48–50,78–83].

Discrete epidemic models generated by numerical methods contain additional parameters to those that already exist in differential equations, such as the time and space steps. Variations in these additional parameters can generate solutions to the discrete equations that do not correspond to any solution of the original differential equations, producing fictitious bifurcations, artificial chaos, spurious solutions, and false stable states [24–26,45,83,84]. Therefore, we must choose numerical discrete schemes that guarantee the qualitative properties of the mathematical models. There are several methods that can be used to obtain accurate solutions. For instance, the Richardson extrapolation on uniform and nonuniform grids or *NSFD* schemes have been used for that end [45,85–90] .

Another, important aspect where a robust numerical method plays an important role is when solving inverse problems to estimate the parameters of the model [48,51,56,57,91,92]. Thus, for mathematical models of a variety of virus is of paramount importance to have a robust and efficient numerical method for solving the differential equations [25,49,51]. Usually, when a differential equation is solved numerically a certain tolerance is prescribed and this has an impact in the success in estimating the parameters [48,49,93,94]. In this paper, we deal with a mathematical model that is based on system of differential equations with discrete time delay [17,32,60,72,95–97]. There are different numerical methods to deal with this type of equations, and some are analogous to the ones used for ordinary differential equations but with additional issues [50,98–102]. One particular numerical method that we are interested in is by using *NSFD* schemes to guarantee some properties of the continuous mathematical model. Some previous works using this methodology have been developed for linear and nonlinear delay differential equations [103–107].

For the construction of a discrete numerical scheme that allows us to efficiently approximate the solutions of the system (6), we use the methodology proposed by Ronald Mickens, see [24–27]. In that order of ideas, for the discrete approximation of the time derivative of a function *<sup>X</sup>*(*t*) <sup>∈</sup> *<sup>C</sup>*1(R), we define the non-standard derivative as

$$\frac{d\_N X(t)}{dt} = \frac{X(t+h) - X(t)}{\varrho(h)} + \mathcal{O}(h),\tag{29}$$

where *<sup>ϕ</sup>*(*h*) is a real positive valued function that satisfies *<sup>ϕ</sup>*(*h*) = *<sup>h</sup>* <sup>+</sup> <sup>O</sup>(*h*2), and *<sup>N</sup>* is to denote the non-standard derivative.

Although there is no general algorithm for constructing an *NSFD* schema that approximates the solutions of a given system of differential equations, the following general rules are often useful to correctly design these schemes.


Let us denote by *I<sup>n</sup> <sup>E</sup>*, *<sup>I</sup>n*, *<sup>V</sup><sup>n</sup>* and *<sup>T</sup><sup>n</sup>* the approximations of *IE*(*nh*), *<sup>I</sup>*(*nh*), *<sup>V</sup>*(*nh*) and *T*(*nh*), respectively, for *n* = 0, 1, 2..., and for *h* size step in time of the scheme. The value of *I <sup>n</sup>*+<sup>1</sup> *<sup>E</sup>* for *<sup>n</sup>* <sup>=</sup> 0, 1, ··· , is calculated using Equation (8) and a quadrature formula. For this case, we use

$$I\_{\mathbb{E}}^{n+1} = e^{-\delta\_{l\underline{k}}\Lambda(n+1)\hbar} \left[ I\_{\mathbb{E}}^{n} + \frac{m\_{1}\hbar}{2} \left( \beta T^{n+1} V^{n+1} e^{\delta\_{l\underline{k}}\Lambda(n+1)\hbar} + \beta T^{n+1-m\_{1}} V^{n+1-m\_{1}} e^{\delta\_{l\underline{k}}\Lambda(n+1-m\_{1})\hbar} \right) \right],\tag{30}$$

with Δ = *m*1*h*.

Next, we make the following non-local approximations of the terms on the right side of the system (27)

$$\begin{cases} \begin{array}{ccccc} \beta T(t)V(t) & \rightarrow & \beta T^n V^n \\ -\beta T(t-\Delta)V(t-\Delta)e^{-\delta l\_{\mathbb{E}}\Delta} & \rightarrow & -\beta T^{n-m\_1+1}V^{n-m\_1}e^{-\delta l\_{\mathbb{E}}\Delta} \\ & -\delta\_{\mathbb{E}}I\_E(t) & \rightarrow & -\delta\_{\mathbb{E}}I\_E^{n+1} \\ & -\delta\_I I(t) & \rightarrow & -\delta\_I I^{n+1} \\ & N\delta\_I I(t) & \rightarrow & N\delta\_I I^{n+1} \\ & -(\mathbb{C} + \beta T(t))V(t) & \rightarrow & -(\mathbb{C} + \beta T^n)V^{n+1} \\ & -\beta T(t)V(t) & \rightarrow & -\beta T^{n+1}V^n \\ & -\mu\_0 T(t) & \rightarrow & -\mu\_0 T^{n+1} \end{array} \tag{31}$$

Then, we approximate the derivatives on the left side of the system (27) as follows

$$\begin{cases} \quad \frac{d\_N I(t)}{dt} \quad \rightarrow \quad \frac{I^{n+1} - I^n}{\varrho(h)}\\ \quad \frac{d\_N V(t)}{dt} \quad \rightarrow \quad \frac{V^{n+1} - V^n}{\varrho(h)}\\ \quad \frac{d\_N T(t)}{dt} \quad \rightarrow \quad \frac{T^{n+1} - T^n}{\varrho(h)} \end{cases} \tag{32}$$

Consequently, the system (27) can be discretized as an implicit *NSFD* scheme given by

$$\begin{aligned} \frac{T^{n+1} - T^n}{\varphi(h)} &= \Lambda - \beta T^{n+1} V^n - \mu\_0 T^{n+1}, \\ \frac{I^{n+1} - I^n}{\varphi(h)} &= \beta T^{n-m\_1+1} V^{n-m\_1} e^{-\delta\_{I\_E}\Lambda} - \delta\_I I^{n+1}, \\ \frac{V^{n+1} - V^n}{\varphi(h)} &= N\delta\_I I^{n+1} - \mathbb{C}V^{n+1} - \beta T^n V^{n+1}. \end{aligned} \tag{33}$$

And the explicit form is given by

$$\begin{split} T^{n+1} &= \frac{\varrho(h)\Lambda + T^n}{1 + \varrho(h)\left(\beta V^n + \mu\_0\right)'}, \\ I^{n+1} &= \frac{\varrho(h)\beta T^{n-m\_1+1} V^{n-m\_1} \varepsilon^{-\delta\_{I\_c}\Delta} + I^n}{1 + \varrho(h)\delta\_I}, \\ V^{n+1} &= \frac{\varrho(h)N\delta\_I I^{n+1} + V^n}{1 + \varrho(h)\left(C + \beta T^n\right)'} \end{split} \tag{34}$$

where *<sup>m</sup>*<sup>1</sup> <sup>=</sup> <sup>Δ</sup> *<sup>h</sup>* <sup>∈</sup> <sup>N</sup>. The initial conditions of scheme (34) are given by

$$T^j = \mathfrak{z}\_{1\prime}^j V^j = \mathfrak{z}\_{2\prime}^j j = -m\_{1\prime} - m\_1 + 1, \cdot, \cdot, 0, \cdot$$

The positivity of scheme (34) is trivially satisfied, since for *n* > 0 it holds that *Tn*, *In*, *V<sup>n</sup>* are positive.

**Theorem 6.** *Let* (*Tn*, *In*, *Vn*) *be a solution of system (34). Then is uniformly bounded for all n* > 0.

**Proof.** From the first equation of scheme (33), one gets that

$$\frac{T^{n+1} - T^n}{\varrho(h)} = \Lambda - \beta T^{n+1} V^n - \mu\_0 T^{n+1} \le \Lambda - \mu\_0 T^{n+1}.$$

When *<sup>n</sup>* <sup>→</sup> <sup>∞</sup> and since *<sup>ϕ</sup>*(*h*) = *<sup>h</sup>* <sup>+</sup> *<sup>O</sup>*(*h*2), then *<sup>ϕ</sup>*(*h*) coincide with 0 in the limit as *<sup>h</sup>* <sup>→</sup> 0. This implies that lim sup *<sup>n</sup>*→<sup>∞</sup> *<sup>T</sup><sup>n</sup>* <sup>≤</sup> Λ *μ*0 .

Next, let <sup>W</sup>*<sup>n</sup>* <sup>=</sup> *<sup>T</sup>n*−*m*<sup>1</sup> <sup>+</sup> *<sup>e</sup> <sup>δ</sup>IE*Δ*In*. From first and second equation of system (33), one obtains

$$\begin{aligned} \frac{\partial \mathcal{W}^{n+1} - \mathcal{W}^{n}}{\varrho} &= \frac{T^{n-m\_1+1} - T^{n-m\_1}}{\varrho(h)} + \varepsilon^{\delta\_{I\_E}\Delta} \frac{I^{n+1} - I^n}{\varrho(h)} = \Lambda - \mu\_0 T^{n-m\_1+1} - \delta\_I \varepsilon^{\delta\_{I\_E}\Delta} I^{n+1} \\ &\leq \Lambda - d\mathcal{W}^{n+1} .\end{aligned}$$

where *d* = min{*μ*0, *δI*}. Thus, lim sup *n*→∞ <sup>W</sup>*<sup>n</sup>* <sup>≤</sup> Λ *<sup>d</sup>* . Therefore, lim sup *n*→∞ *<sup>I</sup><sup>n</sup>* <sup>≤</sup> Λ*e* −*δIE*<sup>Δ</sup> *<sup>d</sup>* . Then, give <sup>&</sup>gt; 0 we can choose a *<sup>M</sup>* <sup>∈</sup> <sup>N</sup> such that *<sup>I</sup><sup>n</sup>* <sup>≤</sup> Λ*e* −*δIE*<sup>Δ</sup> *<sup>d</sup>* <sup>+</sup> for *<sup>n</sup>* <sup>≥</sup> *<sup>M</sup>*. Using the last equation of (33) it is concluded that

$$\frac{V^{n+1} - V^n}{\varphi(h)} \le N\delta\_I I^{n+1} - \mathcal{C}V^{n+1} \le N\delta\_I \left(\frac{\Lambda \varepsilon^{-\delta\_{I\_\mathbb{E}}\Lambda}}{d} + \epsilon\right) - \mathcal{C}V^{n+1}$$

for *n* ≥ *M* + 1. Then lim sup *n*→∞ *<sup>V</sup><sup>n</sup>* <sup>≤</sup> *Nδ<sup>I</sup> C* Λ*e* −*δIE*<sup>Δ</sup> *<sup>d</sup>* <sup>+</sup> , and as is for all > 0 it follows that lim sup *n*→∞ *<sup>V</sup><sup>n</sup>* <sup>≤</sup> *Nδ<sup>I</sup> C* Λ*e* −*δIE*<sup>Δ</sup> *d* . This completes the proof. Moreover, from Equation (27) it follows that *I<sup>n</sup> <sup>E</sup>* is bounded.
