key: cord-0955909-sxawafk2
authors: Shatanawi, Wasfi; Raza, Ali; Arif, Muhammad Shoaib; Abodayeh, Kamaledin; Rafiq, Muhammad; Bibi, Mairaj
title: Design of nonstandard computational method for stochastic susceptible–infected–treated–recovered dynamics of coronavirus model
date: 2020-09-18
journal: Adv Differ Equ
DOI: 10.1186/s13662-020-02960-y
sha: 98275223bab8aa38893cc2ee7fe939d4a41c663d
doc_id: 955909
cord_uid: sxawafk2

The current effort is devoted to investigating and exploring the stochastic nonlinear mathematical pandemic model to describe the dynamics of the novel coronavirus. The model adopts the form of a nonlinear stochastic susceptible-infected-treated-recovered system, and we investigate the stochastic reproduction dynamics, both analytically and numerically. We applied different standard and nonstandard computational numerical methods for the solution of the stochastic system. The design of a nonstandard computation method for the stochastic system is innovative. Unfortunately, standard computation numerical methods are time-dependent and violate the structure properties of models, such as positivity, boundedness, and dynamical consistency of the stochastic system. To that end, convergence analysis of nonstandard computational methods and simulation with a comparison of standard computational methods are presented.

Humanity is enduring many diseases of variable lethality since its birth. Ebola, HIV, and Lassa fever are just a few of them. Ebola is a devastating disease that is transmitted to humans from carrier nonhuman primates, named fruit bats, and destroyed races of humankind. HIV is transferred from cross-species of chimpanzee to humans. In the earlier 19th century it was unknown. It spread rapidly to five continents of the world, killing 300,000 people, including children and women, because its signs and symptoms did not accompany any transmission. Lassa fever has its severity and history of destroying humankind. It is transmitted to humans via rats. Not only this fever, but many other diseases are also present that prove themselves devastating for the living being. Therefore, scientists tried very hard to build instruments to encounter the adverse effects of ailments and to produce possible treatment via vaccine or medicine. Among several other vicious diseases, COVID-19 has uprooted the humanity by killing many of people and is still consuming many lives to date. It was first discovered in the city Wuhan, Province Hubei in China [1, 2] . It is a pandemic that causes respiratory disorder and is transmitted through sneezing droplets of infected individuals. These droplets can fall on the objects around the effected and enter a healthy individual through contact. The number of cases is surging dramatically, raping developed and undeveloped countries together. It is frequently spreading so that it becomes impossible for the world's aristocrats to overcome [3] . Currently, every continent is a sufferer, and among them, China, Iran, the UK, the USA, Spain, Italy are considered the most effected countries. Major symptoms of this disease include lethargy, dry cough, followed by fever [4] . Few patients may develop pains and aches, running nose with nasal congestion, diarrhea, and sore throat. Some individuals have developed these as mild symptoms, while others may show their severe forms. Those affected with mild symptoms recover through special treatment. Every individual out of six showed recovery. This virus shows severe symptoms in older individuals and in those who are already indulged with some sort of disease, like diabetes, cardiovascular disorders, cancer, chronic respiratory disorders, etc. The outbreak of COVID-19 raised many questions like how much time is required to formulate the vaccine or medicine to treat this disease effectively [5] . When would the world be free of all the deadly viruses? How much recession and life loss will occur, and will the world overcome this loss? No answer is available right now. Even though the recovery ratio is much larger compared to the death ratio, the number of fatalities is surging [6, 7] . Khan and Atangana [8] have suggested a fractional-order COVID-19 model with the Atangana-Baleanu-Caputo operator and have introduced this model to evaluate the infection in Wuhan. In the absence of efficient vaccines and treatment to mitigate the COVID-19 pandemic, the function of lock-down is studied [9] . In formulating the proposed mathematical model, the author used the new fractional operator's method. In [10] Khan analyzes the dynamics of the current chaotic system, i.e., using Caputo-Fabrizio and Atangana-Baleanu derivatives. Some other articles on Atangana-Baleanu derivative applications can be found in [11, 12] . A coronavirus model has recently been mathematically considered in [13] . The authors used Pakistan's actual data and discussed the potential control and infection removal from Pakistan. In [14] , where potential removal of the virus was discussed, the data from Ghana and its study using a mathematical model were taken into consideration.

In this paper, we aim to suggest and present mathematical analysis revealing the spread of such a deathly disease, and develop some prediction with real-world data [15, 16] . Also, the proposed structure-preserving nonstandard computational method, named the stochastic nonstandard finite difference (SNSFD), is applied for the given pandemic model [17, 18] . This is the critical point of this paper. The flow of the paper is based on the following sections: In Sect. 2, the formulation of the deterministic susceptible-infectedtreated-recovered model is given. In Sect. 3, we formulate the stochastic susceptibleinfected-treated-recovered model and study its threshold dynamics. In Sect. 4, we present different numerical methods for the stochastic model and develop convergence analysis. Finally, in Sect. 5, the conclusion is presented.

In this section, we consider the dynamics of a human population, which is divided into four components, or subpopulations. The population will be represented by N : [0, ∞) → R, as a function of the time t ≥ 0. Furthermore, each of the components of the population will be denoted by a nonnegative differentiable function S, I, T, R : [0, ∞) → R. The components of the human population are described as follows: S(t) denotes those who are not infected yet with the coronavirus, but have some serious diseases, I(t) denotes humans infected with the coronavirus, T(t) denotes those who are availed via the vaccination or precautionary measures or quarantine, and R(t) denotes those who have recovered from the coronavirus. The nonnegative constants of the model are presented as follows: B denotes the natural rate of new-borns or rate of individuals who have travelled from other countries, β is the rate of interaction between individuals with serious diseases and individuals infected with coronavirus, δ is the rate of individuals with serious diseases who have recovered directly from the coronavirus due to their immune system or quarantine or extensive use of precautions, α is the rate of individuals who may die due to coronavirus or naturally, μ is the rate of individuals infected with coronavirus who are treated, and ρ is the rate of treated individuals who are completely recovered due to quarantine or vaccination. The dynamics of the coronavirus model is given by the nonlinear system of coupled ordinary differential equations as follows:

Obviously, the identity N(t) = S(t) + I(t) + T(t) + R(t) is satisfied at all time instances t ≥ 0. Let S 0 , I 0 , T 0 , and R 0 be nonnegative real numbers that represent the initial sizes of each component of the population, respectively. More clearly, the following conditions are satisfied:

Lemma 1 For any given nonnegative initial conditions, there exists a unique solution S, I, T, R, respectively, for all t ≥ 0. Moreover, it satisfies the following inequality of boundedness:

Proof The total dynamics of the model (1)-(4) is obtained by adding the four equations as follows:

It follows that

Thus a solution of the model (1)-(4) exists for given initial conditions and is eventually bounded on every finite time interval.

Proof Using Eqs. (5a) and (5b), it follows that as the time approaches infinity, t → ∞, the population is bounded by a positive number so on the set ,

therefore, the set is positively invariant.

It is easy to see that there are three steady states of Eqs. (1) to (4) as follows:

• Trivial equilibrium (TE) = (S, I, T, R) = (0, 0, 0, 0);

Notice that the reproduction number R O is the spectral radius of G 1 G -1 2 [4] , where

More precisely, notice that

Let us consider the vector T , the possible changes in the given model are as presented in Table 1 . Table 1 Transition probabilities

Notice that the expectation and variance of the given model are as follows:

. So, the stochastic differential equation of the given model is as follows:

and where the Brownian motion is denoted by W (t).

This method can be applied to Eq. (6) as follows [19] :

where t is the time step size and W m = W t m+1 -W t m is a random variable having the standard normal distribution. It is normally distributed between stochastic drift and stochastic diffusion, i.e., W m ∼ N(0, 1).

In this section, we shall choose parameters from Eq. (1) to (4) and change them into random parameters with small noise as β dt = β dt + σ dW (t) as follows [20] :

The Wiener process is denoted by W k (t), while σ is the randomness of Eqs. (8) to (11) .

The system of Eqs. (8) to (11) is nonintegrable because of the Wiener process.

Let us introduce R S o = R d o -σ 2 2(α+μ) , where R S o denotes the stochastic reproduction number.

If the initial data satisfy (S(0), I(0), T(0), R(0)) ∈ R 4 + , then the system (8)-(11) has a unique solution, and the solution (S(t), I(t), T(t), R(t)) belongs to .

Definition If lim t→∞ I(t) = 0, then in the system (8)- (11) , infected individuals will experience extinction. 

More precisely, by integrating from [0, t], ∀t ≥ 0,

Notice that a local continuous martingale is defined as 

For each N N, define the set I N = {0, 1, 2, . . . , N}. In this section, we will provide and analyze a discretization of the system (8)- (11) . To that end, we consider the temporal period of length T > 0. Fix a uniform partition of the temporal interval [0, T] consisting of N subintervals, and let k = T N . Defining t m = mk, for each m I N , we will employ the notation S m , I m , T m , and R m to represent the numerical approximation to the values of the functions S, I, T, and R, respectively, at time t m .

Also, we will apply different explicit and implicit methods to the given system. As expected, we will employ the discrete initial data (S 0 , I 0 , T 0 , R 0 ), where S 0 = S(0), I 0 = I(0), T 0 = T(0), R 0 = R(0).

This method can be applied to the system of Eqs. (8) to (11) as follows [21] :

where k is time step size and W m = W t m+1 -W t m .

This method can be applied to the system of Eqs. 

A 2 = k B -β S m + A 1 2 I m + B 1 2 -(α + δ) S m + A 1 2 -S m + A 1 2 I m + B 1 2 σ W m , B 2 = k β S m + A 1 2 I m + B 1 2 -(α + μ) I m + B 1 2 + S m + A 1 2 I m + B 1 2 σ W m , C 2 = k μ I m + B 1 2 -α T m + C 1 2 + δ S m + A 1 2 -ρ T m + C 1 2 , D 2 = k ρ T m + C 1 2 -α R m + D 1 2 .

Stage 4 

where k is the time step size and W m = W t m+1 -W t m .

This method can be applied to the system of Eqs. (8) to (11) as follows:

where k is the time step size and W m = W t m+1 -W t m .

In this section, we shall present the following theorems for positivity, boundedness, consistency, and stability.

Theorem 2 For any given initial value (S m (0), I m (0), T m (0), R m (0)) ∈ R 4 + , the system of Eqs. (17) to (20) has a unique positive solution (S m , I m , T m , R m ) ∈ R 4 + for any m ≥ 0.

Proof All the parameters and initial conditions must be nonnegative due to biological reasoning. The proof is straightforward.

α } for all m ≥ 0 is a positively invariant feasible region for Eqs. (17) to (20) .

Proof We rewrite the system (17) to (20) as follows:

Theorem 4 For any m ≥ 0, the system of discrete dynamical Eqs. (17) to (20) has the same steady states as that of the continuous dynamical Eqs. (8) to (11) .

Proof Solving the Eqs. (17) to (20) , we get three states as follows:

• 

Theorem 5 For any m ≥ 0, the proposed computational method is stable if the eigenvalues of Eqs. (17) to (20) lie inside the unit circle.

Proof Consider the right-hand sides of Eqs. (17)- (20) as functions F, G, H, and L as follows:

The Jacobian matrix is defined as 

,

,

Now we want to linearize the model about the steady-state of the model for virus-free state V 1 = ( B α+δ , 0, 0, 0) and R o < 1. The given Jacobian there is Table 1 were used, and the spectral radius was obtained for various temporal time step sizes k in [0, 1000] The eigenvalues of this Jacobian matrix are

So, the given model is stable around V 2 as shown in Fig. 1 . This shows that all the eigenvalues of the Jacobian matrix lie inside the unit circle. So, the given model is stable around V 1 . We plot the largest eigenvalues of the Jacobian matrix associated with V 2 . Notice that the spectral radius is always less than one, as desired.

In this section, we present a comparison analysis of existing numerical methods with the nonstandard computational method by using different values of the parameters as presented in [22] , B = 0.3, β = 0.3, α = 0.3, δ = 0.001, ρ = 0.03, σ = 0.2, for R S 0 < 1. For R S 0 > 1, B = 0.3, β = 0.5, α = 0.3, δ = 0.001, ρ = 0.03, σ = 0.2. Moreover, we used S 0 = 0.7, I 0 = 0.05, T 0 = 0.2, R 0 = 0.05.

Example 1 (Simulation for the free of coronavirus state) For given real data, the reproductive value is R S 0 = 0.6975 < 1. So, the system converges to V 1 = ( B α+δ , 0, 0, 0). However, in Fig. 2 , we just take a susceptible component; the numerical solution of the system using the nonstandard and standard computational methods is presented. In Fig. 2, (a) , (c), and (e), the results converge to V 1 for certain small temporal step size. On the other hand, in Fig. 2, (b), (d) , and (f), the results violated the structural properties of the system. The comparison graph exhibits the fact that the proposed nonstandard computational method is capable of preserving the structural properties of the solutions in the sense of biological reasoning, as desired.

Example 2 (Simulation for the existence of coronavirus state) The reproductive value is R S 0 = 1.1958 > 1. So, the system converges to V 2 = (0.8000, 0.1480, 0.0473, 0.0047). Nevertheless, in Fig. 4 , we just take the infected component; the numerical solution of the system using the nonstandard and standard computational methods is presented. In Fig. 3, (a) , (c), and (e), the results converge to V 2 for a certain small temporal step size. On the other hand, Example 3 (Simulation for reproduction number with the effect of treatment) Let μ = 0.0488. Notice that the reproduction value decreases, moving the dynamics of the model from virus existence state to virus-free state. So, the virus-free state of the system is stable. However, Fig. 4 exhibits the fact that the increases in treatment strategy such as quarantine or vaccination can overcome the pandemic of coronavirus, as desired. Example 4 (Simulation for a component of the infected humans) Let us now for different values of μ (treatment rate) notice that the number of infected humans converges to zero. Eventually, the reported rate of infected humans has been controlled in certain scenarios. Consequently, Fig. 5 shows the fact that the treatment strategy has a vital role in the control of the pandemic of the coronavirus around the world.

This study focused on the stochastic susceptible-infected-treated-recovered model, which elaborates a simplified way to describe the dynamics of coronavirus in the human population. Also, we have applied different explicit and implicit computational methods to study the dynamics of the stochastic system. The proposed nonstandard computational method is unconditionally convergent as compared to other explicit numerical methods. This method preserves the structural properties of stochastic models, such as consistency, stability, positivity, and boundedness [12] . Moreover, this study provides a valuable proof that the confinement rules are vital to control the situation in a reasonable time. If the contact rate is dropped between people then stabilization of classes can be achieved significantly earlier. Therefore, there is no chance for governments to get rid of social distancing rules among individuals.

A mathematical model for simulation the phase-based transmissibility of a novel coronavirus

A conceptual model for the coronavirus disease 2019 (Covid-19) outbreak in Wuhan, China with individual reaction and government action

Covid-19 infection: origin, transmission and characteristics of human coronaviruses

Early dynamics of transmission and control of Covid-19: a mathematical modelling study

Modelling and optimal control of typhoid fever disease with cost effective strategies

Modeling the epidemic dynamics and control of Covid-19 outbreak in China

Dynamical behavior and sensitivity analysis of a delayed coronavirus epidemic model

Modeling the dynamics of novel coronavirus (2019-nCov) with fractional derivative

Modelling the spread of Covid-19 with new fractal-fractional operators: can the lock-down save mankind before vaccination?

The dynamics of a new chaotic system through the Caputo-Fabrizio and Atangana-Baleanu fractional operators

Modeling the impact of non-pharmaceutical interventions on the dynamics of novel coronavirus with optimal control analysis with a case study

Fractional model of HIV transmission with awareness effect

Modeling the impact of non-pharmaceutical interventions on the dynamics of novel coronavirus with optimal control analysis with a case study

Global stability and cost-effectiveness analysis of Covid-19 considering the impact of the environment: using data from Ghana

Stability behavior of mathematical model MERS corona virus spread in population

Transmission potential and severity of Covid-19 in South Korea

Competitive numerical analysis for stochastic HIV/AIDS epidemic model in a two-sex population

Construction of equivalent stochastic differential equation models

Numerical solution of stochastic SIR model via split-step forward Milstein method

A fundamental principle for constructing nonstandard finite difference schemes for differential equations

Convergence of the stochastic Euler scheme for locally Lipschitz coefficients

Design of a nonlinear SITR fractal model based on the dynamics of a novel coronavirus (Covid-19)

We always warmly thank the anonymous referees. We are also grateful to Vice-Chancellor, Air University, Islamabad, for providing an excellent research environment and facilities. The first and fourth author also thank Prince Sultan University for funding this work through the COVID19 Emergency Research Program.

This research project is funded by the Research and initiative centre RIC2020 # 65, Prince Sultan University.

All of the necessary data and the implementation details have been included in the manuscript.

No competing interests regarding this article.

The authors have equal contributions. All authors read and approved the final manuscript. 

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.Received: 11 June 2020 Accepted: 9 September 2020