Construction and numerical analysis of a fuzzy non-standard computational method for the solution of an SEIQR model of COVID-19 dynamics

1.   Introduction

The COVID-19 belongs to a large class of deadly viruses that have infected millions of people around the world and seriously challenged not only their lifestyles but also the economies and GDP of countries themselves. It was first identified as the cause of many pneumonia cases in December 2019, when the first case of respiratory infection was reported in Wuhan, China ^[1]. Transmission of the coronavirus typically occurs when an exposed person catches infected droplets released by infected people while sneezing, exhaling, or coughing. In general, 80 to 85% of infected people recover from the COVID-19 pandemic without specific treatment ^[2]. Confirmed cases of COVID-19 infection worldwide have risen to nearly 110 million, with more than 2.4 million deaths due to the disease. According to Worldometers ^[3], the total number of active cases in the world on March 19, 2020, was 32,129 and the total number of deaths worldwide was recorded as 4.6 million. The maximum number of deaths due to COVID-19 in the United States, Brazil, India, Russia, Mexico, and Pero were recorded as 0.6, 0.58, 0.44, 0.26, and 0.19 million, respectively. Worldometers data shows that on Jan 17, 2022, the total number of cases in the world was 328,677,751, a total number of deaths were recorded as 5.5 million and deaths due to COVID-19 in the United States, Brazil, India, Russia, Mexico, and Pero were recorded as 0.87, 0.62, 0.48, 0.30 and 0.20 million, respectively.

Like many other research questions regarding COVID-19 disease, reliable estimation of transmission dynamics is an important part of the research. The novel COVID-19 is still causing great panic among people around the world. The world is facing the fifth wave of this new epidemic. Various approaches are being considered to combat this deadly disease. These approaches have been offered in the form of modeling, data analysis, control of disease spread and clinical insights. Saif Ullah and Khan developed a mathematical model to study the transmission dynamics and potential control of the COVID-19 pandemic in Pakistan by formulating a mathematical model without optimal control and time-dependent control variables ^[4]. It was concluded that the implementation of strict social distancing and contact tracing for those exposed to quarantine is the most effective strategy to minimize the disease burden. Naveed et al. studied the mathematical analysis of the coronavirus model with time delay effect by analyzing the reproduction number of the model. Social distancing, travel restrictions, quarantine, isolation, and hospitalization were seen as delaying factors and it was concluded that more delaying tactics ultimately led to control of the pandemic ^[5]. Shatanawi et al. studied the stochastic coronavirus model by developing Euler Maruyama, stochastic Euler and stochastic RK methods respectively ^[6]. A non-standard finite difference approach has been developed which preserves dynamic consistency, boundedness and positivity. Moreover, the results were compared with the deterministic model. N. Shahid et al. have developed a Spatio-temporal epidemic model of the COVID-19 with advection and diffusion process and discussed the consistency, stability, and positivity of the developed model ^[7]. Naveed et al. have developed a mathematical model for the transmission of COVID-19 ^[8]. Routh Hurwitz criteria and the Lasalle invariance principle were used for showing the stability of the model. The effect of the delay factor on the reproduction number was also discussed. Naveed et al. have proposed a fractional-order mathematical model for COVID-19 and explained that the fractional epidemic model provides a better understanding and biologically more information about disease dynamics ^[9]. Gao et al. studied the 2019-nCoV infection system with a Caputo-defined non-local operator and obtained results using the fractional natural decomposition method ^[10]. Khan et al. studied the interaction between bats and unknown hosts, between humans and the reservoir of infection in detail ^[11]. The mathematical results are presented and a fractional model is formulated. Rafiq et al. proposed a susceptible infected treatment and recovered (SITR) dynamical model of COVID-19 ^[12]. The model is studied numerically and convergence and dynamical consistency of the model is discussed. The effect of treatment on the basic reproduction number is also examined.

The concepts of susceptible, exposed, infectious, quarantined and recovered are uncertain due to varying degrees of susceptibility, exposure, infectivity, quarantine and recovery among individuals of the population. Differences may arise if considered population groups with different habits, customs and age brackets have different resistances, etc. More realistic models are needed to account for these different levels of individuals. According to Mishra et al. epidemic systems, especially those dealing with infectious diseases, need to be treated differently due to high uncertainties. These uncertainties are due to the fact that the epidemic strength of an infectious agent depends, among other things, on the proportion of susceptible nodes and the proportion of infectious nodes. Sensitivity and infectivity are inherently fuzzy terms and therefore ideal arguments for discussing fuzzy logic ^[13]. Since the parameters used in epidemic models are uncertain. This uncertainty can be described by introducing the fuzzy theory. The use of fuzzy logic in biological systems has great potential but is used less frequently. Zadeh introduced the fuzzy theory in 1965 ^[14]. Barros et al. proposed an SI model using fuzzy theory and considering the transmission coefficient as a fuzzy set ^[15]. The average number of people infected is compared to the average change in virus load and the basic reproduction value is analyzed. Mondal et al. have concentrated on the plague models having the transmission coefficient as a fuzzy set and formulated an SIS model ^[16]. The fuzzy basic reproduction number is examined and a pathogenic threshold state is achieved when the system undergoes transcritical branching. Verma et al. studied a fuzzy SIR model by introducing new parameters ^[17]. In addition, a comparative study of equilibrium points for classical and fuzzy models is performed. The fuzzy reproduction number was also examined for groups of people with different viral loads. Ortega et al. ^[18] employed the fuzzy logic for the prediction in the epidemiology related problems in the infectious disease. A model of rabies among the partially vaccinated dogs was discussed. A comparison between the fuzzy linguistic rules and the classical differential equation was also done. Renu Verma et al. ^[19] have studied the dynamics of Ebola virus disease by employing fuzziness in all biological parameters. SEIR and SEIRHD models were prosed for the transmission trajectories of the Ebola outbreak. The existence of the equilibria and their stability was studied by employing triangular fuzzy numbers. The stability of the equilibria was related to the basic reproduction number which was calculated with the help of the next-generation matrix and the numerical methods were used to support theoretical work. Das and Pal developed a SIR model with imprecise biological parameters ^[20]. The existence of equilibrium points was discussed and the numerical simulation was done to support the analytical results. Sadhukhan et al. ^[21] studied about food chain model with optimal harvesting in a fuzzy environment. The stability of the system is studied using the Lyapunov function and the existence of a bioeconomic equilibrium is discussed. Jafelice et al. introduced a model for the development of the HIV positive population and the manifestation of AIDS focusing on the mode of transmission from HIV to AIDS ^[22]. Due to its uncertain nature, the transmission rate is considered a fuzzy set of viral load. Allaoui et al. ^[23] presented a fuzzy model to giving up smoking. Moreover, the results obtained with the fuzzy theory are compared with those of the classical case. Jessica and Filipe suggested a method of dermoscopy images using fuzzy numbers and concluded that fuzzy numbers could be applied to slide images of skin lesions. The use of pre-processing steps to minimize potential artwork or noise effects in the image was also important in increasing the result rate ^[24]. Boventura and Gonzaga proposed edge detection in grayscale images based on fuzzy number theory to remove uncertainties about the shades of gray that make up the image and then calculate the suitability of pixels to an area homogeneous around the image. Pixels that do not belong to the region are classified as edge pixels ^[25]. Al-Amin et al. discussed techniques for solving non-homogeneous linear fuzzy difference equations ^[26]. Moreover, the equilibrium and stability analysis has been performed. Mishra and Prajapati developed a fuzzy SEIQRS model for the transmission of malicious codes in computer networks by simulating the results for different parameters and analyzed the model's stability ^[27]. The effectiveness of antivirus software and node blocking due to malicious code attacks has also been analyzed. Lefevre et al. introduced a fuzzy application of epidemiological models related to the prevalence of HIV in a sample of individuals with injectable drugs and different fuzzy scenarios were analyzed for a different number of users and a different number of HIV test samples per year for the samples used in trials varied from case to case ^[28].

Since each community changes with the evolution of the environment, even the biological parameters used in mathematical models are not always fixed ^[29]. Global warming is the main cause of the many problems attributed to the increase in the earth's average temperature. The change in temperature also affects the rate of transmission of the virus in the population. David et al. studied the relationship between average temperatures and confirmed cases of COVID-19 in Brazil ^[30]. The results showed that temperatures were in a negative linear relationship with the number of confirmed cases. Irfan et al. investigated the relationship between temperature and COVID-19 transmission in different provinces of Pakistan ^[31]. Low-temperature provinces were found to have strong associations between temperature and COVID-19 transmissibility. In this sense, fuzzy mathematical models are more meaningful than crisp models.

The above models are insufficient in connection to the development of fuzzy numerical and mathematical techniques. Keeping this in mind, we have studied an SEIQR model with fuzzy parameters. The use of fuzzy theory helps us to deal with issues of uncertainty quantification in mathematical disease modeling. Consequently, the use of fuzzy parameters helps us to explain the transmission of the COVID-19 virus more accurately. We calculated the fuzzy reproduction number and fuzzy equilibrium points for the studied model. In addition, a fuzzy NSFD method for the model is developed which is a fuzzy extension of Micken's NSFD scheme ^[32]. The fuzzy stability of the proposed method is also discussed.

The rest of this study is structured as the basic definitions are presented in Section 2 which will be used in this study and the fuzzy model formulation is discussed. Fuzzy reproduction number and fuzzy equilibrium points have also been discussed in this section. The numerical solution and the resulting simulation results are presented in Sections 3 and 4, respectively. Concluding remarks and future directions relevant to this research are summarized in Section 5.

2.   SEIQR model with Fuzzy Parameters

In this section, we introduce the extended SEIQR model with fuzzy parameters. First, let's mention some basic definitions which will be useful for this study.

2.1.   Fuzzy subset

A fuzzy subset $\mathrm{A}$ of the universal set $\mathrm{X}$ is represented symbolically by the membership function $\mathrm{A}:\; \mathrm{X}\;\to \;[0, \;1]$, where $\mathrm{A}\left(\mathrm{x}\right)$ indicates the degree of membership of $\mathrm{x}$ in the fuzzy set $\mathrm{A}$ ^[33].

2.2.   Fuzzy number

A fuzzy subset $A$ in $\mathbb{R}$ is called fuzzy number when ^[33]:

● All $\delta$-levels of $A$ are non-empty, with $0\le \delta \le 1$, that is, $A$ must be normal.

● All $\delta$-levels of $A$ are closed intervals of $\mathbb{R}$.

● The support of $A$, that is, $support = \left\{x\in \mathbb{R}:A\left(x\right) > 0\right\}$ is bounded.

2.3.   Triangular fuzzy number (TFN)

A fuzzy number $A = \left(a, b, c\right)$ is said to be triangular ^[33] if its membership function is given by

where $a\le b\le c.$

2.4.   Expected value of a fuzzy number

The expected value of a fuzzy number was introduced by B. Liu and Y. Liu ^[34]. It is denoted by $E\left[\zeta \right]$ and is defined by

where $Cr$ is the credibility measure and for any real number $r$, it can be defined as

2.5.   Expected value of a TFN

The expected value of a TFN $A = \left(a, b, c\right)$ is given by ^[35]

2.6.   Fuzzy basic reproductive number ${{R}_{c}}^{f}$

The fuzzy basic reproductive number ${{R}_{c}}^{f}$ is defined as ^[35]

where ${R}_{c}\left(\upsilon \right)$ is the reproduction number and $E\left[{R}_{c}\left(\upsilon \right)\right]$ is the expected value of a TFN defined in formula 2.

The SEIQR model consists of 5 segments, i.e., susceptible, exposed, infected, quarantined and recovered subpopulations respectively. Consider the following system of 5 differential equations describing the SEIQR model proposed by Hussain et al. ^[36].

The corresponding fuzzy SEIQR model can be written as

The detail of the model parameters is given in Table 1.

We suppose that the infection rate $\beta \left(\upsilon \right)$, the recovery rate $\gamma \left(\upsilon \right)$ and the mortality rate $\xi \left(\upsilon \right)$ due to COVID-19 are fuzzy numbers that depend on the individual virus load. Let $\beta = \beta \left(\upsilon \right)$ be the possibility that virus transmission occurs during an encounter between a susceptible person and an infected person. It was introduced by Barros et al. ^[15] and can be defined as

The $\beta \left(\upsilon \right)$ will be maximum if the $\upsilon$ is maximum and it will be negligible when the $\upsilon$ is minimum. ${\upsilon }_{min}$ is the minimum virus load needed for disease transmission and the disease transmission is maximum at ${\upsilon }_{M}$, where it is equal to 1. The membership function of $\beta \left(\upsilon \right)$ is shown in Figure 1.

The recovery rate $\gamma = \gamma \left(\upsilon \right)$ introduced by Verma et al ^[17] is also a fuzzy number and can be defined as

The minimum recovery rate is ${\gamma }_{0} > 0$.

The death rate $\xi = \xi \left(\upsilon \right)$ introduced by Verma et al ^[17] can also be considered as a fuzzy number as it increases due to the increase of the infection of the disease and can be defined as

The death rate $\xi \left(\upsilon \right)$ will be higher when the amount of virus $\upsilon$ is the highest i.e., ${\upsilon }_{0} < \upsilon$ and the maximum death is $1-\zeta$, $(\zeta \ge 0)$.

Because the amount of virus is different for each group of people. To make the model more realistic, we only consider human individuals from a given group $N$ with a classification (weak, medium and strong) given by an expert, which can be seen as a linguistic variable with a membership function ${\Gamma }\left(\upsilon \right)$ and is given by

The membership functions of $\beta \left(\upsilon \right)$, $\gamma \left(\upsilon \right)$, $\xi \left(\upsilon \right)$ and ${\Gamma }\left(\upsilon \right)$ are shown in Figures 1–4 respectively.

2.7.   The fuzzy basic reproductive number ${{R}_{\mathit{\boldsymbol{c}}}}^{\mathit{\boldsymbol{f}}}$

We find ${R}_{c}$ by incorporating the next generation matrix approach ^[37] and it is given by

Since ${R}_{c}$ is a function of the amount of virus $\upsilon$ and we analyze it for different amounts of virus.

Case 1: If $\upsilon < {\upsilon }_{min}$, then from Eq (6), we have $\beta \left(\upsilon \right) = 0$ and consequently from Eq (10) we obtain,

Case 2: If ${\upsilon }_{min} < \upsilon \le {\upsilon }_{M}$, then from Eq (6), we have $\beta \left(\upsilon \right) = \frac{\upsilon -{\upsilon }_{min}}{{\upsilon }_{M}-{\upsilon }_{min}}$ and consequently from Eq (10) we obtain,

Case 3: If ${\upsilon }_{M} < \upsilon < {\upsilon }_{max}$, then from Eq (6), we have $\beta \left(\upsilon \right) = 1$ and consequently from Eq (10) we obtain,

The basic reproduction number ${R}_{c}\left(\upsilon \right)$ is an increasing function of the parasitic virus load $\upsilon$ and it is well-defined as a fuzzy variable. Consequently, the expected value of ${R}_{c}\left(\upsilon \right)$ is well-defined and it can be expressed as a triangular fuzzy number as:

Now by using formulas (2) and (3), we find the fuzzy reproduction number as follows:

2.8.   Fuzzy equilibrium analysis

Case 1: If $\upsilon < {\upsilon }_{min}$, then from Eq (6), we have $\beta \left(\upsilon \right) = 0$, Substituting it the in system (5), we get $S = \frac{{\Lambda }}{\mu }$, $E = 0$, $I = 0$, $Q = 0$ and $R = 0$. Therefore, we obtain:

which is the DFE point. This is the situation where there is no virus in the population. Biologically, the disease is eradicated when the virus level is below the minimum level required in the population for disease transmission.

Case 2: If ${\upsilon }_{min} < \upsilon \le {\upsilon }_{M}$, then from Eq (6), we have $\beta \left(\upsilon \right) = \frac{\upsilon -{\upsilon }_{min}}{{\upsilon }_{M}-{\upsilon }_{min}}$. In this case, we obtain the equilibrium point

Case 3: If ${\upsilon }_{M} < \upsilon < {\upsilon }_{max}$, then from Eq (6), we have $\beta \left(\upsilon \right) = 1$. In this case, we obtain the equilibrium point

The equilibrium points ${E}^{*}({S}^{*}, {E}^{*}, {I}^{*}, {Q}^{*}, {R}^{*})$ and ${E}^{**}({S}^{**}, {E}^{**}, {I}^{**}, {Q}^{**}, {R}^{**})$ are called endemic equilibrium points. These equilibria occur when the virus is larger than the minimum required and persists in the population.

3.   Numerical modeling of fuzzy SEIQR model

In this segment, we will examine a fuzzy nonstandard finite difference method (FNSFD) based on Micken's theory ^[32] to solve the dynamical system 5.

3.1.   Fuzzy NSFD scheme

To develop a fuzzy NSFD scheme for the studied model, first-order time derivatives are described by using forward difference approximation i.e., $ḟ\left(t\right)$ can be approximated as

${s}^{n}$, ${e}^{n}$, ${i}^{n}$, ${q}^{n}$ and ${r}^{n}$ are the approximations of $S\left(nh\right)$, $E\left(nh\right)$, $I\left(nh\right)$, $Q\left(nh\right)$ and $R\left(nh\right)$ for $n = \mathrm{0, 1}, \mathrm{2, 3}, \dots .$, where $h$ is the time step size. Thus the FNSFD scheme for system 5 can be written as

On simplification the Eqs (17)–(21), we obtain

3.2.   Fuzzy stability and convergence analysis of the NSFD scheme

In this part, the stability and convergence of the developed FNSFD scheme at DFE point ${E}^{0}({S}^{0}, {E}^{0}, {I}^{0}, {Q}^{0}, {R}^{0})$ and endemic equilibrium points ${E}^{\mathrm{*}}({S}^{\mathrm{*}}, {E}^{\mathrm{*}}, {I}^{\mathrm{*}}, {Q}^{\mathrm{*}}, {R}^{\mathrm{*}})$ and ${E}^{\mathrm{*}\mathrm{*}}({S}^{\mathrm{*}\mathrm{*}}, {E}^{\mathrm{*}\mathrm{*}}, {I}^{\mathrm{*}\mathrm{*}}, {Q}^{\mathrm{*}\mathrm{*}}, {R}^{\mathrm{*}\mathrm{*}})$ respectively are discussed here.

The systems (22)–(26) can be written as:

The Jacobean matrix corresponding to the systems (27)–(31) is

The Jacobean at the DFE ${E}^{0}({S}^{0}, {E}^{0}, {I}^{0}, {Q}^{0}, {R}^{0}) = \left(\frac{{\Lambda }}{\mu }, 0, 0, 0, 0\right)$ is

The FNSFD scheme will be unconditionally convergent iff absolute eigenvalues of the above Jacobean matrices at the DFE point ${E}^{0}\left({S}^{0}, {E}^{0}, {I}^{0}, {Q}^{0}, {R}^{0}\right)$ are less than unity, i.e., $\left|{\mathrm{\lambda }}_{\mathrm{i}}\right| < 1, i = 1, 2, 3, 4, 5.$

Eigenvalues of the above Jacobean matrix are ${\lambda }_{1} = {\lambda }_{2} = \frac{1}{1+h\mu } < 1$, ${\lambda }_{3} = \frac{1}{1+h\alpha +h\mu } < 1$, ${\lambda }_{4} = \frac{1}{1+h\varphi +h\mu } < 1$ and ${\lambda }_{5} = \frac{1}{1+h(\mu +\gamma (\upsilon )+\xi (\upsilon \left)\right)} < 1$. Since all eigenvalues are less than one which proves the desired result.

The Jacobean at the first endemic equilibrium ${E}^{*}({S}^{*}, {E}^{*}, {I}^{*}, {Q}^{*}, {R}^{*})$ is given by

The Jacobean at the second endemic equilibrium point ${E}^{**}({S}^{**}, {E}^{**}, {I}^{**}, {Q}^{**}, {R}^{**})$ is given by

Again the FNSFD scheme will be unconditionally convergent iff absolute eigenvalues of the above Jacobean matrices at the EE points ${E}^{\mathrm{*}}({S}^{\mathrm{*}}, {E}^{\mathrm{*}}, {I}^{\mathrm{*}}, {Q}^{\mathrm{*}}, {R}^{\mathrm{*}})$ and ${E}^{\mathrm{*}\mathrm{*}}({S}^{\mathrm{*}\mathrm{*}}, {E}^{\mathrm{*}\mathrm{*}}, {I}^{\mathrm{*}\mathrm{*}}, {Q}^{\mathrm{*}\mathrm{*}}, {R}^{\mathrm{*}\mathrm{*}})$ are less than unity, i.e., $\left|{\mathrm{\lambda }}_{\mathrm{i}}\right| < 1, i = 1, 2, 3, 4, 5.$

The largest eigenvalues of the Jacobeans $J\left({E}^{\mathrm{*}}\right)$ and $J\left({E}^{\mathrm{*}\mathrm{*}}\right)$ have been plotted by using the MATLAB database and shown in Figure 5a and Figure 5b respectively.

These figures show that the spectral radius in both cases remains less than 1 for each time step. This implies that all the eigenvalues of the Jacobin matrices at the endemic equilibrium points lie in a unit circle, which guarantees that the FNSFD scheme is convergent for each time step used.

4.   Results and discussions

In this section, we presented the graphical results of the proposed FNSFD method.

We can examine the behavior of the fuzzy NSFD method for the SEIQR epidemic model for COVID-19 spread in the above graphs. The behavior of the graphs is investigated for various values of $h$. In Figure 6, the compartments of subpopulations susceptible, exposed, infected, quarantined and recovered have been shown at step size $h = 0.1$ at the DFE point. The results show the positive behavior and convergence of the proposed FNSFD method. From this, we conclude that the developed method can illustrate the actual behavior of the disease dynamics at the DFE point.

Figure 7 shows the convergence of the FNSFD method to the true equilibrium points of the continuous model at step size, $h = 0.1$ at first endemic equilibrium point. The figures show that the FNSFD method remains convergent and retains all the essential properties of the continuous dynamical system.

Figure 8 shows the convergence of the FNSFD method to the true equilibrium points of the continuous model for an increasing step size $h = 15$ at the first endemic equilibrium point. The figures show that the FNSFD method remains convergent for an increasing time step and preserves all the essential properties of the continuous dynamical system.

In Figure 9, the behavior of the FNSFD method is shown at the second endemic equilibrium point at a small step size $h = 0.1$. The results show that the method remains convergent and preserves all the essential properties of the continuous dynamical system.

In Figure 10, the behavior of the FNSFD method is shown at the second endemic equilibrium point for an increasing step size $h = 15$. The results show that the method is converging to the same equilibrium point and the step size does not affect it. Again we conclude that this method is a reliable tool for reflecting the actual behavior of the model.

Convergence, stability and positivity are the main features of epidemic models. It is clear from the graphs that the proposed FNSFD method converges to the same equilibrium points at different step sizes at all equilibrium points. The graphs also reveal that the method shows stable behavior and positivity at all equilibrium points. It is an explicit numerical scheme, easy to implement, shows stable behavior numerically and demonstrates a good agreement with analytic results possessed by continuous model. Moreover, it preserves all the essential features of a continuous dynamical system.

5.   Conclusions

In this study, we have considered a fuzzy SEIQR system that models the transmission of COVID-19. We assumed that the infected persons do not transmit the disease equally and each individual has a different degree of transfer of disease infectivity which depends on the quantity of virus possessed by an individual. Similarly, the recovery rate and disease induced death rate are also not same for each individual. We considered the disease transmission rate $\beta \left(\upsilon \right)$, the recovery rate $\gamma \left(\upsilon \right)$ and the additional death rate $\xi \left(\upsilon \right)$ as fuzzy variables depending on the amounts of virus. These parameters have fixed values in deterministic models and do not depend on the virus load directly. Therefore, the fuzzy model can be considered more flexible and balanced than the crisp system. The use of fuzzy theory helps to deal with issues of uncertainty quantification in mathematical disease modeling.

Since fuzzy variables are functions of the virus load which depend on the amounts of virus, we analyzed it for different amounts of virus. Keeping this in mind, we discussed fuzzy equilibrium points of the studied model by considering the amounts of virus in the population. We proved that if the amount of virus is less than the minimum amount required for disease transmission in the population, the disease-free equilibrium point is obtained. We reached the endemic equilibrium points, when the amounts of virus in the population were greater than the minimum amounts of virus required for disease transmission.

We calculated the basic reproduction number by utilizing the concept of the next-generation matrix method. Again, we analyzed it for different amounts of virus. We obtained fuzzy basic reproduction number by using the expected value of a fuzzy number.

Furthermore, we have developed a FNSFD scheme. The developed scheme has been analyzed for different amounts of virus. The positivity of the solutions of such dynamic population models must be preserved by the proposed numerical technique that can be visualized in this article. The proposed method showed positivity at small and large step sizes for different amounts of virus.

Fuzzy stability, another essential characteristic, is studied at equilibrium points and it is shown that all of the three equilibrium points have the same stability properties for the studied model. The proposed numerical scheme showed stable behavior numerically and demonstrated a good agreement with analytic results possessed by the continuous model.

We performed some numerical experiments to confirm our theoretical results. The simulation results also confirmed the numerical results. Thus, the developed method explored a structure-preserving technique for solving a mathematical system in epidemiology.

The present work mainly focuses on the construction, implementation and analysis of an explicit first-order numerical scheme in fuzzy non-standard finite difference environments. In the future, this work could be extended to fuzzy stochastic, fuzzy delayed, or fuzzy fractional dynamic senses. NSFD modeling theory could also be extended to age-structured fuzzy epidemic models.

Conflict of interest

The authors declare that they have no competing interests.