A nonstandard finite difference scheme for a time-fractional model of Zika virus transmission

1.   Introduction

The Zika virus (ZIKV) is an emerging arbovirus that is transmitted by several so-called vectors, the most important being the Aedes aegypti mosquito. Vectors are living organisms that can transmit infectious pathogens between humans, or from animals to humans. ZIKV was first isolated from a macaca monkey in the Zika forest in Uganda in 1947, giving the virus its name, cf. ^[1,2].

The first major ZIKV epidemic began 2007 on the Yap archipelago in the Federated States of Micronesia, where a high number of cases were recorded in about 75% of the population within a few months ^[3,4]. Later, a worldwide epidemic occurred in French Polynesia (2013–2014) with approximately 28,000 cases (about 11%) of the total population ^[5]. In 2015, ZIKV was reported in Brazil via viremic travelers or infected mosquitoes ^[6], it also began to spread in Mexico ^[7]. Messina ^[8] showed that up to 2.17 billion people live in "risk areas" (tropical and subtropical regions).

The ZIKV infection is associated with mild symptoms: Fever, headache, rash, myalgia, and conjunctivitis, similar to other arboviruses (dengue or chikungunya) ^[9] and no deaths have been reported to date. Nevertheless, ZIKV has emerged as a major cause of the development of the Guillain-Barré syndrome ^[10]. Also, there is uncertainty about the outcome of co-infections with other arboviruses such as Dengue fever. Furthermore, there is no available treatment for ZIKV infection. Patient care is based on symptomatic treatment with a combination of acetaminophen and antihistamine medications ^[4].

Several mathematical models have been developed to address different categories in epidemiology, such as prediction of disease outbreaks and evaluation of control strategies ^{[11,12,13,14]}. The first mathematical epidemic model dates back to Kermack and McKendrick (1927), who were concerned with mass events in the susceptible, infected, and recovered (SIR) disease transmission cycle ^[15]. Manore and Hyman ^[16] proposed a mathematical model for ZIKV representing disease transition and population dynamics Gao ^[17] developed a model of ZIKV transmission through bites of Aedes mosquitoes and also through sexual contact. Lee and Pietz ^[18] developed a mathematical model for Zika virus using logistic growth in human populations. Nishiura et al. ^[19] proposed a mathematical Zika model that exhibits the same dynamics as Dengue fever.

Fractional order approaches were used in COVID-19 transmission models by using fractional order Caputo derivative ^[20], the analysis of semi-analytical solutions of a hepatitis B epidemic model using the Caputo-Fabrizio operator ^[21], the study of stability and Lyapunov functions for HIV/AIDS epidemic models with the Atangana-Baleanu-Caputo derivative ^[22], the mathematical modeling of the measles epidemic with optimized fractional order under the classical Caputo differential operator ^[23].

Let us briefly mention most recent research directions. Wang et al. ^[24] developed and analyzed a partial differential equation (PDE) model with periodic delay to understand the transmission dynamics of the Zika virus. This reaction-diffusion model takes into account spatial structure, seasonality and the temperature sensitivity of the incubation period, factors that play a crucial role in vector-borne diseases like Zika. The authors gave insights into the estimation of the basic reproduction number ($R_0$), the impact of spatial averaging on $R_0$ in periodic systems, the potential underestimation of $R_0$ in the absence of seasonality, and the significance of shortening the incubation period in reducing the risk of disease transmission.

Ibrahim and Dénes ^[25] created a compartmental model for Zika virus transmission, emphasizing microcephaly. It determines the basic reproduction number ($R_0$), assesses global stability conditions, and fits the model to Colombian data (2015–2017). The findings underscore the significance of preventing mosquito bites, controlling mosquito populations, and using protection during sexual contact to minimize Zika-related microcephaly cases.

Murugappan et al. ^[26] developed and analyzed a non-linear mathematical model to understand the dynamics of Zika virus transmission. The study focused on the comparative analysis of transmission dynamics in males, females and children. The calculation of the reproductive ratio provides insights into the spread of the Zika virus. They analytically determine the equilibrium points and their stability, supported by numerical simulations. Their research aims to identify the population most affected by Zika transmission, contributing valuable information for targeted intervention strategies.

Next, we want to mention two recent related works (on other diseases) with useful approaches for analysis of stability analysis of the transmission dynamics. Oguntolu et al. ^[27] developed a deterministic model to study the dynamics of soil-transmitted helminth diseases. The model includes two equilibrium points: the disease-free equilibrium (DFE) and the endemic equilibrium (EEP). Stability analysis of these points provides insights into disease persistence and potential eradication. The authors apply optimal control theory to identify effective intervention strategies related to hygiene awareness rates in susceptible and infectious classes.

Peter et al. ^[28] formulated and analyzed an epidemic model for COVID-19, considering both the first and second vaccination doses. The control reproduction number is derived and the stability of the system is assessed, with the COVID-free equilibrium being stable when the reproduction number $R_0$is less than one. The model is calibrated using real data from Malaysia, and a global sensitivity analysis identifies key parameters influencing disease dynamics. The results highlight the importance of adherence to preventive measures and increased vaccination rates in reducing the spread of the disease.

In this work we derive a new nonstandard finite difference scheme (NSFD) for a recent SEIR (susceptible-exposed-infectious-recovered) model ^[13] that describes the spread of the Zika virus using a human-mosquito compartmental model and a human-mosquito-monkey compartmental model. Despite the fact that this NSFD scheme has a nonlinear denominator function, this schemes has a couple of favourable properties: It is explicit and due to its construction it reproduces important properties of the solution, like the number and location of fixed-points, the positivity and certain conservation laws. The goal of this work is to briefly demonstrate, in detail, how the NSFD methodology is to be applied to a system of coupled ordinary differential equations (ODEs), where the discretizations are dynamical consistent with the basic properties of the continuous differential equations, e.g., positivity, asymptotic behavior, memory effects, etc.

This paper is organized as follows. In Section 2, we formulate the ZIKV transmission models. Section 3 includes the analysis, especially the boundedness of the solution and the stability analysis of the two considered models. In Section 4, we design the nonstandard finite difference method for the two proposed models and show how it can be extended to time-fractional variants of the models using the L1 method. In Section 4, we propose NSFD schemes for the conventional and the time-fractional version of our models. The numerical results of our novel schemes are shown in Section 5. Finally, Section 6 presents the conclusions and some outlook.

2.   The ZIKV transmission models

This section outlines an extended version of two mathematical compartmental models that describe the transmission dynamics of ZIKV, as presented in the work by Hamou Maamar et al. ^[13]. Our models incorporate varying population sizes for humans, vectors (mosquitoes) and nonhuman primates, expanding upon previous models. Additionally, we introduce fractional order operators, particularly the Caputo derivative, to generalize these models.

In areas without nonhuman primates, such as Yap State and French Polynesia, ZIKV is likely maintained in a human-mosquito-human cycle, suggesting that the virus has adapted to humans as reservoir hosts ^[29]. This setting will lead us the first model, formulated in a SEIR-SEI framework.

Boorman and Porterfield ^[30] showed in a laboratory setting that Monkeys can become infected with ZIKV. However, there is no evidence that ZIKV is transmitted to humans through contact with animals. On the other hand, the presence of specific antiviral antibodies in various nonhuman primates, suggesting that other reservoirs may play a role in the ZIKV transmission cycle, cf. ^[31]. For this reason, we also consider a second extended model.

2.1.   The model parameters and assumptions

The human population is divided into four classes (so-called 'compartments'): Susceptible, exposed (latently infected), infected and recovered (individuals who have acquired immunity). We denote the number in each compartment by $S_h$, $E_h$, $I_h$ and $R_h$. Accordingly, we divide the vector population (adult female mosquitoes) into three compartments: Susceptible, exposed and infected, with the analogous notation $S_v$, $E_v$ and $I_v$. Next, we define the total number of populations as

The parameters and assumptions in our model are based on biological evidence and epidemiological insights, cf. ^[13]:

● $B$ is the average number of bites per mosquito per day.

● $\beta_{vh}$ is the probability rate that a bite from an infectious vector will infect a human, the product $B\beta_{vh}$ is the number of disease-transmitting bites per infectious mosquito per day, and the product $B\beta_{vh}I_v(t)$ is the number of disease-transmitting bites per day in the entire mosquito population at time $t$ (measured in days). However, multiplying $B\beta_{vh}I_v(t)$ by the proportion of susceptible people at time $t$ represents the number of disease-transmitting bites per day by infectious mosquitoes on susceptible people at time $t$ (the daily rate at which susceptible people are exposed).

● The parameter $\mu_h$ is the proportion of the human population that dies each day ('human mortality rate').

● $\nu_h$ is the daily rate at which exposed people become infected ('human infection rate').

● $\eta_h$ denotes the daily rate at which infected people become immune. ('human immunity rate').

● The parameter $\beta_{hv}$ is the probability rate that the bite of an infectious human will infect a mosquito; $B\beta_{hv}$ is the number of disease-transmitting bites per mosquito per day. Thus, the product $B\beta_{hv}S_v(t)$ is the number of bites per day that result in disease being transmitted by susceptible mosquitoes at time $t$. Multiplying $B\beta_{hv}S_v(t)$ by the proportion of infectious people at time $t$ the complete rate of disease-transmitting bites at time $t$ (the daily rate at which susceptible mosquitoes become infected).

● The parameter $\mu_v$ is the proportion of the mosquito population that dies each day ('mosquito mortality rate').

● $\nu_v$ denotes the daily rate at which exposed mosquitoes become infected ('mosquito infection rate').

We include a constant system inflow, the per-capita birth rates $\Lambda_h$, $\Lambda_v$ (e.g., birth of new individuals that can become infected).

2.2.   The human-mosquito model

Now, we are ready to formulate the first model. The system of ODEs has the following form

The dynamical system described by (2.2) is depicted in Figure 1. We note that by a convention in epidemiology models all parameters in (2.2) are assumed to be positive.

Summing up the equations in (2.2) gives immediately the ODE system for the time evolution of the total populations of humans and mosquitos

that can be solved easily, cf. Section 4.4. Since the Zika virus transmission has a faster dynamic than the human birthrate and the human natural mortality, $N_h(t)$ can be regarded as a conserved quantity of the above ODE system, if we set for simplicity $\Lambda_h = \mu_h = 0$. Note that this is not the case for the vector (mosquito) which has a comparable dynamic and the asymptotic behavior

This well-known limiting behavior can be exploited for a further simplification of the model (2.2) (so-called 'limiting model') by removing $I_v$, and thus the remaining vector components $S_v$, $E_v$ can be plotted in a 2D phase diagram, cf. ^[32]. We return later to this property when designing the numerical scheme.

2.3.   The human-mosquito-monkey model

Accordingly, we define the total monkey population as

Next, we introduce similar parameters for the monkey population, cf. ^[13]:

● $\beta_{vm}$ is the probability rate that a bite from an infectious mosquito will infect a monkey.

● The parameter $\mu_m$ is the proportion of the monkey population that dies each day.

● $\nu_m$ is the daily rate at which exposed monkeys become infected.

● $\eta_m$ the daily rate at which infected monkeys become immune.

The corresponding system of ODEs for the temporal evolution of the human, vector and monkey population has the following form

The dynamical system described by Eq (2.6) is depicted in Figure 2.

Again, summing up the equations in (2.6) yields for the total populations

with simple exact solutions, see Section 4.4. Analogously, $N_h(t)$ and $N_m(t)$ can be regarded as a conserved quantity of the above ODE system, if we set $\Lambda_h = \mu_h = 0$ and $\Lambda_m = \mu_m = 0$.

Using standard arguments (see e.g., ^[33]) it can be easily shown that both ODE systems (2.2), (2.6) preserve the positivity of the solution. This basic property should be respected by any reasonable numerical method and yields as a byproduct the stability of the scheme.

2.4.   A Fractional-order human-vector model

The fractional-order dynamics of the transmission of the Zika virus to human and vector populations is given by the following system

with the initial conditions

where $^{C}D^\alpha X(t)$ is the Caputo derivative and it is defined as:

Adding the equations of the system (2.8) yields the fractional ODEs

In the model given above, we modified the right-hand sides parameters $\mu_h^\alpha$, $B^\alpha$, $\nu_h^\alpha$, $\eta_h^\alpha$, $\mu_v^\alpha$ and $\nu_v^\alpha$ using the procedure described in Diethelm ^[34] in order to adjust the dimensions because the dimension of the left-hand sides of the equations is (time)$^{-\alpha }$. Note that in the limit case $\alpha\to1$, the system (2.8) reduces to the classical system given in (2.2).

3.   Analysis of the models

3.1.   Non-negativity and boundedness of solutions

The positivity and boundedness of the solutions of an epidemiological system are essential properties. Therefore, it is important to prove that all subpopulations in the systems (2.2), (2.6) and (2.8) are non-negative and bounded for all times $t\ge 0$. The following results show how to confirm these two properties.

We now focus on the human-mosquito system (2.2) and prove the following theorem, which confirms the positivity and boundedness of the system.

Theorem 3.1. The closed region

is a positively invariant set for the system (2.2).

Proof. Let $S_h(0) > 0$, then

By using the Comparison Lemma ^[35], we have

Similarly, it can be shown that

From Eqs (4.4) and (4.5) the quantities $N_h(t)$ and $N_v(t)$ are non-negative for all $t\ge0$, and

Thus, $S_h(t)$, $E_h(t)$, $I_h(t)$, $R_h(t)$, $S_v(t)$, $E_v(t)$ and $I_v(t)$ are bounded.

The corresponding proof for the human-mosquito-monkey system (2.6) is analogous.

Remark 3.2. The properties of positivity and boundedness in (fractional) epidemiological systems are essential for maintaining biological relevance, ensuring model consistency, facilitating stability analysis, supporting accurate numerical simulations and providing reliable insights for public health decision making. These properties contribute to the overall robustness and applicability of epidemiological models in the understanding and management of infectious diseases. In more detail they read:

● Biological interpretation:

– Positivity: In epidemiological models, the state variables typically represent populations, such as the number of susceptible, infected, or recovered individuals. Positivity ensures that these quantities remain non-negative, consistent with the biological reality that populations cannot have negative values.

– Boundedness: Bounded solutions imply that population levels do not grow indefinitely. This is consistent with the finite nature of populations in the real world, where resources and environmental factors impose limits on the growth of individuals within a given ecosystem.

● Model consistency:

– Positivity: The requirement of positivity maintains meaningful interpretation of model variables. Negative values would lack biological significance and could lead to unrealistic scenarios that are inconsistent with the principles of epidemiology.

– Boundedness: Bounded solutions ensure that model predictions remain within realistic and feasible limits. Unbounded solutions could lead to predictions of population sizes or infection levels that are physically implausible.

● Stability analysis:

– Positivity: Positivity is often a critical criterion for stability in dynamical systems. It ensures that solutions do not exhibit erratic behavior and that the system does not diverge into unrealistic states.

– Boundedness: Bounded solutions contribute to the stability of the system by preventing uncontrolled growth or decay. Stability is essential for making reliable predictions about the long-term behavior of the epidemiological system.

● Numerical simulations:

– Positivity: Preserving positivity is critical when using numerical methods to solve differential equations. Algorithms that preserve positivity help prevent numerical artifacts and ensure that simulated solutions remain physically meaningful.

– Boundedness: Numerical stability is often linked to the boundedness of solutions. Unbounded solutions can lead to numerical instability, making it difficult to obtain accurate and reliable results from simulation.

● Public health implications:

– Positivity and boundedness: From a public health perspective, positivity and boundedness ensure that model predictions are realistic and actionable. Decision makers rely on models to guide interventions, and unrealistic predictions could lead to misinformed public health strategies.

The following theorem highlights the positivity and boundedness of the fractional-order human-vector model (2.8):

Theorem 3.3. The region $\Omega ^\alpha = \Big\{ (S_h, E_h, I_h, R_h, S_v, E_v, I_v) \in \mathbb{R}_+^7 \colon 0\le N_{\alpha, h}$ $\le \frac{\Lambda_h^\alpha}{\mu_h^\alpha} \mathit{\text{ and }} 0\le N_{\alpha, v}\le \frac{\Lambda_v^\alpha}{\mu_v^\alpha}\Big\}$ is a non-negative invariant for the model (2.8) for $t\ge 0$.

Proof. We have

and using the Laplace transform, we obtain

then

and applying the inverse Laplace transform, we get

where $E_{\alpha, \alpha +1}$ denotes the Mittag-Leffler function

Using the well-known recurrence relation for the Mittag-Leffler function ^[36] for $\beta = 1$,

we may write the Eq (3.1) as

and thus

We proceed similarly to derive the equation of $N_{\alpha, v}(t)$,

and conclude that

As a result, the functions $S_h$, $E_h$, $I_h$, $R_h$, $S_v$, $E_v$ and $I_v$ are all non-negative.

3.2.   Stability analysis

System (2.8) always has a disease-free equilibrium (DFE) at:

where

The infection components considered in (2.8) model consist of $E_h$, $I_h$, $E_v$ and $I_v$. Using the next generation approach ^[34], the basic reproduction number of the model (2.8) is defined as $R_0^\alpha = \rho (F_{\alpha} V_{\alpha}^{-1})$, where $F_\alpha$ represents the new infection matrix and $V_\alpha$ represents the transition matrix. The values for $F_\alpha$ and $V_\alpha$ are provided below:

and

Thus,

Remark 3.4. The basic reproduction number $R_0$ is a critical epidemiological parameter that represents the average number of secondary infections produced by an infected individual in a fully susceptible population. In the context of Zika control, understanding the biological significance of $R_0$ is critical to developing effective strategies to manage and mitigate the spread of the virus. Here are some key aspects of its biological significance:

● Transmission dynamics: $R_0$ provides insight into the potential for Zika virus transmission within a population. When $R_0$ is greater than 1, each infected individual infects, on average, more than one other person, indicating the potential for sustained transmission. Conversely, if $R_0$ is less than 1, the virus is likely to die out in the population.

● Epidemic potential: A high $R_0$ indicates a higher risk of a Zika epidemic. The higher the $R_0$, the more difficult it will be to control the spread of the virus. Understanding epidemic potential helps public health officials strategically allocate resources and implement interventions.

● Effectiveness of control measures: Control interventions such as vector control, vaccination, and public health campaigns can affect $R_0$. If interventions are successful in reducing the transmission rate, they can reduce $R_0$ below the critical threshold of 1, leading to a decrease in the number of new infections and eventual control of the outbreak.

● Population susceptibility: $R_0$ accounts for the proportion of the population that is susceptible to the virus. As susceptibility decreases due to factors such as prior exposure or vaccination, $R_0$ is effectively reduced, helping to control Zika transmission.

● Targeted vaccination strategies: Understanding $R_0$ helps prioritize and target vaccination efforts. High-risk populations or areas with a high potential for transmission, as indicated by a high $R_0$, can be prioritized for vaccination campaigns to achieve maximum impact in controlling the spread of Zika.

● Vector control strategies: $R_0$ is influenced by factors related to the mosquito vector (e.g., Aedes mosquitoes). Implementing vector control measures, such as eliminating breeding sites or using insecticides, can reduce the transmission rate, thereby lowering $R_0$ and preventing Zika transmission.

● Monitoring and surveillance: Continuous monitoring of $R_0$ provides a real-time assessment of the dynamics of Zika transmission. This information is valuable for public health authorities to adapt control measures based on the evolving epidemiological situation.

3.3.   Local and global stability of DFE

The following theorem discuss the local stability of DFE.

Theorem 3.5. The disease-free equilibrium of the proposed fractional-order model is locally asymptotically stable if $R_0^\alpha < 1$ and is unstable if $R_0^\alpha > 1$.

Proof. The Jacobian matrix of system (2.8) at DFE is given by,

the characteristic polynomial is then given by,

then $X_1 = X_2 = -\mu_v^\alpha$ and $X_3 = -\mu_v^\alpha$ are three eigenvalues. The remaining eigenvalues correspond to the roots of the following polynomial

where

and

The polynomial $q(X)$ has the following Hurwitz matrix :

So, by the Routh-Hurwitz criterion the roots of $q(X)$ have negative real parts if and only if all principal minors are strictly positive, that is,

and

In order to prove the global stability of the equilibrium points, we need to recall the following result:

Lemma 3.6 (See ^[37]). Let $X(t) \in \mathbb{R}$ be a continuous and differentiable function. Then, for any time instant $t\ge0$

where $g(x) = x-1-\ln{x}$.

Note that for $\alpha = 1,$ the inequality in (3.4) becomes equality. Now, taking into account the Lyapunov direct method, we provide the global stability of the equilibria in the following theorem.

Theorem 3.7. If $(R_0^\alpha)^2 < \frac{N_{\alpha, v}(0) N_{\alpha, h}(0) }{N_{\alpha, h}^* N_{\alpha, v}^*} < 1$, then the DFE is globally asymptotically stable.

Proof. We consider the following Lyapunov function

where

with

Now using Lemma 3.6, the derivative of $V$ in the Caputo sense with respect to $t$ is given by:

and thus

which implies

and have

thus

This implies that if $(R_0^\alpha)^2 < \frac{N_{\alpha, v}(0) N_{\alpha, h}(0) }{N_{\alpha, h}^* N_{\alpha, v}^*}$ then $^{C}D^\alpha V(t) < 0$ for all $(S_h, E_h, I_h, R_h, S_v, \; E_v, I_v)\neq E_{DF}$ and $^{C}D^\alpha V(t) = 0$ for $(S_h, E_h, I_h, R_h, S_v, E_v, I_v) = E_{DF}$. Therefore, by LaSalle's invariance principle, the DFE is globally asymptotically stable.

4.   The nonstandard finite difference method

In this section, we explain the technique of nonstandard finite difference schemes (NSFDs). A NSFD scheme is constructed to satisfy the positivity condition and the conservation laws. Consequently, the solutions are bounded, i.e., stable. Also, only the fixed-points of the ODE systems (2.2) and (2.6) appear in the NSFD scheme. The specific full details are not given; we refer to the book of Mickens ^[38] for the discretization strategy.

4.1.   Nonstandard finite difference schemes

NSFD methods for the numerical integration of differential equations had their origin in a paper by Mickens published in 1989 ^[38]. In this section, an NSFD scheme is constructed to satisfy the essential positivity condition and the conservation law for $\Lambda_h = \mu_h = 0$, $\Lambda_v = \mu_v = 0$ and $\Lambda_m = \mu_m = 0$ which leads as a byproduct to the stability of the scheme. We will also check that the equilibrium points of the ODE model also appear in the proposed NSFD scheme.

Let us recall that schemes such as those based on Runge-Kutta methods can yield wrong negative solutions (see ^[39,40]) can produce 'false' or 'spurious' fixed-points, which are not fixed points of the original ODE system, cf. ^[41].

Finally, we will determine in Section 4.4 the so-called denominator function $\phi(h)$, such that we obtain the correct long-time behaviour. We refer to ^[42,43], where we established an NSFD scheme for a similar compartment model as here.

We remind the reader that a numerical scheme for a system of first-order differential equations is called NSFD scheme if at least one of the following conditions ^[38] is satisfied:

● The orders of the discrete derivatives should be equal to the orders of the corresponding derivatives appearing in the differential equations.

● Discrete representations for derivatives must, in general, have nontrivial denominator functions.

– The first-order derivatives in the system are approximated by the generalized forward difference method (forward Euler method)

with the numerical approximation $u^n\approx u(t_n)$, $n = 0, 1, 2\dots$ on the uniform grid $t_n = n\, h$ with the step size $h = \Delta t$.

– Here, $\phi\equiv\phi(h) > 0$ is the so-called denominator function such that $\phi(h) = h+\mathcal{O}(h^2)$. This function $\phi(h)$ is chosen so that the discrete solution has the same asymptotic behaviour as the analytical solution, see Section 4.4.

● The nonlinear terms are approximated by non-local discrete representations, for instance by a suitable function of several points of a mesh, like $u^2(t_n) \approx u^n\, u^{n+1}$ or $u^3(t_n) \approx (u^n)^2 \, u^{n+1}$.

● Special conditions that hold for either the ODE and/or its solutions should also apply to the difference equation model and/or its solution, e.g., positivity of the solution, convexity of the solution (in finance) and equilibrium points of the ODE system, including their local asymptotic stability properties.

In NSFD schemes, derivatives must be modeled by discrete analogues that take the form, cf. ^[38]

where $t_n = n\, h$, $u^n$ is the approximation of $u(t_n)$, and $\psi(h) = 1 + \mathcal{O}(h)$. The purpose of this more general time discretization (4.1) in NSFD schemes, is to properly model the asymptotic long-time behavior of the solution.

4.2.   NSFD scheme for the human-mosquito model

Next, we propose the following NSFD discretization for solving the ODE system (2.2)

with the denominator functions for each subsystem

The exact solutions of $N_h$ and $N_v$ are given by

and

Thus

Let us briefly comment on the discretizations of the nonlinear (here: quadratic) terms. For example, in the first line (4.2) we have discretized the nonlinear contact term $\beta_{vh}I_v(t)S_h(t)$ in (2.2) by $\beta_{vh}I_v^n S_h^{n+1}$ rather than, say, $I_v^n S_h^n$ or $I_v^{n+1} S_h^{n+1}$. The rule is that exactly one factor of the variable appearing in the time derivative (here $S_h$) must be taken at the new time level $n+1$. This is needed to obtain a positivity preserving scheme, see (4.6). In order not to destroy the explicit sequential evaluation, all other variables are taken from the previous time level, unless they are already known from a previous step, like $I_h^{n+1} S_v^{n+1}$ in the sixth line. If possible, discrete conservation properties (here: Total population of humans, vectors) must also be taken into account.

Observe that although the initial scheme (4.2) can be considered implicit, the variables at the ($n+1$)-th discrete-time level can be explicitly calculated in terms of the previously known variable values as given in the sequence of the equations above, i.e. we can rewrite it as an explicit form

The calculation must be done in exactly this order. All parameters appearing in these type of epidemic models are always non-negative. This is the convention used in fields related to the spread of diseases. From the explicit representation (4.6) it is easy to deduce that this scheme preserves the positivity, given some natural conditions on the parameters.

4.3.   NSFD scheme for the human-mosquito-monkey model

Correspondingly, the NSFD discretization for solving the ODE system (2.6) reads

Accordingly, we rewrite the scheme (4.7) in an explicit sequential formulation

4.4.   The denominator function

Finally, it only remains to correctly determine the denominator function $\phi(h)$. To do so, we reconsider the combined total population $N = N_h, N_v$ or $N_m$ of the ODE systems (2.2) and (2.6)), now without neglecting the birthrates and the natural mortality. Here, we introduce accordingly the combined values $\Lambda = \Lambda_h, \Lambda_v$ or $\Lambda_m$, $\mu = \mu_h, \mu_v$ or $\mu_m$ for the system (2.2) and the extended system (2.6). At a first glance, it looks inappropriate to add the populations of humans, mosquitos and monkeys, but this has purely mathematical reasons: it is used for the asymptotic behaviour that later leads to the denominator function $\phi(h)$, which must be the same for all components of the ODE system.

Adding the equations of (2.2) or (2.6), we easily obtain the following differential equation describing the dynamics of the combined total population $N$

It is solved by

with $N(0) = N_h(0)+N_v(0)+N_m(0)$. From (4.10) we immediately deduce that we have in the long term $\lim_{t\to\infty}N(t) = \Lambda/\mu$. Let us briefly note that this link between the transient dynamics and their 'natural' limiting systems can be used to reduce the dimension of this model, cf. ^[32].

Next, adding the equations in the discrete NSFD model (4.2) yields

i.e.,

The denominator function can be derived by comparing Eq (4.11) with the discrete version of Eq (4.10), that is

such that the (positive) denominator function is defined by

i.e.,

Note that the conservation property requires all the denominator functions $\phi(h)$ for the compartments to be the same. Otherwise, it would be impossible to obtain a discrete analogue like (4.11) which is also needed for stability reasons.

Remark 4.1. An even more accurate way to compute the denominator function would take into account the transition rate $\Upsilon_i$ at which the $i^{th}$ compartment is entered by individuals for all model compartments $\mathcal{K}_i$, $i = 1, 2, \dots$ (e.g., $\beta_{vh}$, $\nu_h$, $\eta_h$, $\nu_v , \dots$), cf. ^[44]. In this case the parameter $\mu$ occurring in the denominator function in Eq (4.15) would be replaced by a parameter $1/T^*$. $T^*$ could be determined as the minimum of the inverse transition parameters:

4.5.   A NSFD scheme for a time-fractional model

Again, let us consider a uniform temporal grid $t_0 = 0 < t_1 < \cdots < t_{N_T} = T$, $t_n = nT/ N_T$, where $N_T\in\mathbb{N}$. Next, we present a numerical approximation of the Caputo derivative using the NSFD method. We have

We discretize the term $\frac{\text{d}{X(\tau)}}{\text{d}{\tau}}$ on the interval $[t_j, t_{j+1}]$ as

where $X^j = X(t_j)$ and $\phi_{\alpha}(h)$ from (4.15).

where

Each equation in (2.8) can be written as

at the point $t = t_{n+1}$, we have

Now, we apply the scheme (4.16) to the system (2.8), we obtain

Setting $n = 0$, equations of $N_{\alpha, h}^{n+1}$ and $N_{\alpha, v}^{n+1}$ in (4.17) give

and

The exact solution of the Eqs (3.2) and (3.3) can be rewritten as

and

The denominator function $\phi_{\alpha, h}(h)$ ($\phi_{\alpha, v}(h)$ respectively) can be derived by comparing the exact version (4.20) ((4.21) respectively) with the discrete version (4.18) ((4.19) respectively), that is

It is not difficult to show that $\phi_{\alpha, h}(h)$ and $\phi_{\alpha, v}(h)$ reduce to the classical $\phi_h(h)$ and $\phi_v(h)$ in (4.3) when $\alpha = 1$.

5.   Numerical results

In this section, we present the numerical solution of the systems (2.2) and (2.6) using the NSFD schemes (4.6) and (4.8). Then, we compare it with the solution computed by the $\texttt{ode45}$ solver of Matlab.

5.1.   The human-mosquito model

We denote by $Y$ the matrix of order $N_T\times 7$ that contains the approximated solution determined by the $\texttt{ode45}$ solver which is given by

The parameters used to simulate the model are listed in the Table 1. The initial conditions are always set to

The following Figures 3–8 represent the trajectories in the three dimensional space of the human and the vector populations, respectively. They show that the NSFD method remains stable and approaches the disease-free equilibrium (DFE) or endemic equilibrium (EE) points.

The Figures 10 and 9 show that the approximate solutions obtained by the NSFD method and $\texttt{ode45}$ method are very closed to each other. However, the solution $Y$ obtained by the $\texttt{ode45}$ solver becomes negative for some values of $t$. This does not figure clearly in the curves because the smallest negative value of $Y$ is $-4.02\times10^{-7}$.

The Table 2 presents the percentage of negative values in the matrix $Y$ simulating the human-mosquito model (2.2) with the $\texttt{ode45}$ solver using the parameters for the disease-free point in the Table 1. The results given in Table 2 show that the NSFD preserves the positivity for all step sizes in $[0, T]$, which is a desirable modeling property. On the other side, the $\texttt{ode45}$ method yields solutions that becomes negative for some value of $t$.

5.2.   The human-mosquito-monkey model

Now we simulate the system for the data given in Tables 1 and 2. The initial conditions are always set to

Figures 11–14 show that the numerical solution approximates very well the solution of the continuous system by preserving positivity and converging towards the equilibrium points DFE or EE. Table 4 gives the percentage of negative values for the NSFD method and the $\texttt{ode45}$ solver. It can easily be seen that NSFD preserves the positivity of the continuous system where the $\texttt{ode45}$ solver failed in some cases.

5.3.   The time-fractional model

In this section, we provide some numerical simulations of the discrete model (4.17) with different values of fractional order $\alpha$. To proceed with the simulation, we use the parameter values in Table 1 and the initial conditions in (5.1). The numerical simulation results for the NSFD fractional order obtained for different values of $\alpha$ are displayed in Figures 15–20. These figures show two different scenarios:

$\textbf{Case 1 DFE.}$ The dynamical behavior of system for different values of $\alpha$ is shown in Figures 15–17 for $R_0^\alpha < 1$ which implies that it converges to the DFE. It is noticeable that due to the memory property of the Caputo fractional derivatives, the evolution of the system becomes slower each time the $\alpha$ decreases. Therefore, the system decays to the equilibrium like $t^{-\alpha}$, as previously established in ^[45].

$\textbf{Case 2 EE.}$ For $R_0^\alpha > 1$, Figures 18–20 show the impact of changing the Caputo fractional order $\alpha$ on Zika dynamics. The observed behavior from these figures demonstrates that the EE is shifted towards $EE, EE_{\alpha_1}, EE_{\alpha_2}$ and $EE_{\alpha_3}$ when $\alpha$ is decreasing.

The numerical results above show the memory effect for the fractional dynamical system, which does not occur in the ODE system as already proved by ^[46,47]. They also show that the new approach is very effective, preserves the positivity of the system, applies easier and can be used as an alternate method for solving fractional differential problems.

6.   Conclusions and outlook

In this work, we introduced a groundbreaking nonstandard finite difference (NSFD) method tailored to numerically solve a SEIR model elucidating the spread of the Zika virus. To evaluate its effectiveness, we compared the approximate solution with the $\texttt{ode45}$ solver solution in the absence of an exact solution.

Our numerical simulations robustly demonstrate that the discrete system not only converges to the same equilibrium points as the continuous system, but also establish the superiority of the NSFD scheme in preserving positivity, a critical aspect often compromised by standard ODE solvers.

A notable contribution of this research is the applicability of the NSFD methodology beyond the specific SEIR model studied. Our results indicate that this method can be seamlessly extended to various epidemic models, providing a versatile and powerful tool for understanding and predicting the spread of disease.

It is important to emphasize that our use of Caputo-type fractional derivatives to describe the temporal dynamics of epidemiological models addresses the crucial aspect of memory effects. This consideration becomes crucial in realistic systems, such as endemic models, where it captures the waning effects of vaccination or the biphasic decay behavior of infections and diseases. In essence, our study not only advances numerical methods for Zika virus modeling, but also contributes valuable insights applicable to a broader range of epidemiological scenarios.

By pursuing the following future research directions, we can further advance the field of epidemiological modeling, refine the NSFD methodology and contribute to more effective tools for understanding and managing infectious diseases:

● Enhancement of NSFD methodology: Investigate further refinements and optimizations to the NSFD method, exploring potential modifications or extensions that could enhance its computational efficiency and applicability to a broader range of epidemiological models. This could involve exploring different numerical schemes or considering adaptive strategies to handle varying model complexities.

● Comparative studies with other numerical methods: Conduct comparative studies with other state-of-the-art numerical methods beyond the $\texttt{ode45}$ solver. Assess the NSFD method against a spectrum of numerical solvers to provide a comprehensive understanding of its strengths and weaknesses in different contexts, ensuring a more robust validation of its superiority.

● Incorporation of real-world data: Extend the research by incorporating real-world epidemiological data related to the Zika virus or other infectious diseases. Validate the NSFD method against empirical data to assess its predictive accuracy and reliability in capturing complex dynamics observed in actual outbreaks.

● Generalization to other infectious agents: Explore the adaptability of the NSFD methodology to model the spread of other infectious agents. Investigate its effectiveness in diverse epidemiological scenarios involving different pathogens, transmission modes and population structures, thereby broadening its practical utility.

● Advanced memory effect modeling: Delve deeper into the modeling of memory effects in epidemiological systems. Investigate alternative fractional derivatives or memory kernels to capture more nuanced temporal dynamics, especially in scenarios with prolonged immunity, evolving vaccination strategies or changing patterns of human behavior.

● Multiscale modeling: Consider the incorporation of multiscale modeling approaches, examining how the NSFD method performs when integrating information across different spatial and temporal scales. This could lead to more comprehensive models that better reflect the intricate interactions within heterogeneous populations.

● Interdisciplinary collaboration: Foster interdisciplinary collaborations with experts in epidemiology, public health and mathematical modeling. Such collaborations can provide valuable insights, validate model assumptions and ensure the practical relevance of the NSFD methodology in addressing real-world challenges.

● Sensitivity analysis and uncertainty quantification: Conduct sensitivity analyses to identify key parameters influencing model outcomes and perform uncertainty quantification to assess the robustness of predictions. This can enhance the credibility of the NSFD method and contribute to its adoption in decision-making processes.

Finally, it will be an interesting aspect to use some published experimental data to illustrate the theoretical results. For this purpose we will develop special calibration routines that fit to NSFD schemes.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Conflict of interest

The authors declare there is no conflict of interest.

Appendix

A.1.   The human-mosquito model

The system (2.8) has a unique endemic equilibrium point that exists whenever $R_0^\alpha > 1$ and it is given by

A.2.   The human-mosquito-monkey model

The system (2.6) has two equilibrium points, the disease-free equilibrium $DFE = (\frac{\Lambda_h}{\mu_h}, 0, 0, 0, \frac{\Lambda_v}{\mu_v}, 0, 0, \frac{\Lambda_m}{\mu_m}, 0, 0, 0)^\top$ and the endemic equilibrium $EE = (S_h^{**}, E_h^{**}, I_h^{**}, R_h^{**}, S_v^{**}, E_v^{**}, I_v^{**}, S_m^{**}, E_m^{**}, I_m^{**}, \; R_m^{**})^\top$, where

$I_v^{**}$ is implicitly given as the zero of the following rational fraction expression

which is determined numerically. The basic reproduction number of (2.6) is

where

and

A.3.   The time-fractional model

The proof of the theorem requires the following lemma :

Lemma A.1. If $X^0, X^1, \dots, X^n\ge 0$ then

Proof. For $n\in\mathbb{N}^*$, we have

and

Thus

Theorem A.2 (Positivity of solution). Let the initial data $S_h^{0}, E_h^{0}, I_h^{0}, R_h^{0}, S_v^{0}, E_v^{0},$ and $I_v^{0} \ge0$, then all the components $S_h^{n+1}, E_h^{n+1}, I_h^{n+1}, R_h^{n+1}, S_v^{n+1}, E_v^{n+1},$ and $I_v^{n+1} \ge0$ in the system (4.17) are satisfied for all $n\in\mathbb{N}$.

Proof. We have for $n = 0$

We suppose that for $1, 2, \dots, n$, $S_h^n$, $E_h^n$, $I_h^n$, $R_h^n$, $S_v^n$, $E_v^n$ and $I_v^n\ge0$. The hypothesis of induction and Lemma A.1 allow for the statement for $n+1$, i.e.,