Design of nonstandard computational method for stochastic susceptible–infected–treated–recovered dynamics of coronavirus model

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.


Introduction
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.

Formulation of deterministic SITR model
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:

Basic properties
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.

Lemma 2 The closed set = {S(t), E(t), T(t), R(t) R
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.

Steady states of the model
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

Stochastic SITR model
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).

Euler-Maruyama method
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).

Parametric noise in SITR model
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.

Stochastic reproduction dynamics
Let us introduce R S o = R d o -σ 2 2(α+μ) , where R S o denotes the stochastic reproduction number.
More precisely, by integrating from [0, t], ∀t ≥ 0, Notice that a local continuous martingale is defined as

Numerical methodology
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 .

Stochastic Euler method
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 .

Stochastic Runge-Kutta method
This method can be applied to the system of Eqs.

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

Nonstandard computational method
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 .

Convergence analysis
In this section, we shall present the following theorems for positivity, boundedness, consistency, and stability.
Proof All the parameters and initial conditions must be nonnegative due to biological reasoning. The proof is straightforward.
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.
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.

Conclusion
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.