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–infected–treated–recovered model is given. In Sect. 3, we formulate the stochastic susceptible–infected–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 , as a function of the time . Furthermore, each of the components of the population will be denoted by a nonnegative differentiable function . The components of the human population are described as follows: denotes those who are not infected yet with the coronavirus, but have some serious diseases, denotes humans infected with the coronavirus, denotes those who are availed via the vaccination or precautionary measures or quarantine, and 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 is satisfied at all time instances . Let , , , and 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 . 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 is positively invariant.

Using Eqs. (5a) and (5b), it follows that as the time approaches infinity, , 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: where Notice that the reproduction number is the spectral radius of [4], where More precisely, notice that

;

;

,

Stochastic SITR model

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

Table 1

Transition probabilities

T_i=Transition	P_i=Probabilities
T_2=[−1100]^T	P_2=βSIΔt
T_3=[−100]^T	P_3=αSΔt
T_4=[−1010]^T	P_4=δSΔt
T_5=[0−100]^T	P_5=αIΔt
T_6=[0−110]^T	P_6=μIΔt
T_7=[00−10]^T	P_7=αTΔt
T_8=[00−11]^T	P_8=ρTΔt
T_9=[000−1]^T	P_9=αRΔt

Notice that the expectation and variance of the given model are as follows: Note that the stochastic drift is , while stochastic diffusion is .

So, the stochastic differential equation of the given model is as follows: with initial conditions , , and where the Brownian motion is denoted by .

Euler–Maruyama method

This method can be applied to Eq. (6) as follows [19]: where Δt is the time step size and is a random variable having the standard normal distribution. It is normally distributed between stochastic drift and stochastic diffusion, i.e., .

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 as follows [20]: The Wiener process is denoted by , 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 , where denotes the stochastic reproduction number.

Lemma 3

If the initial data satisfy , then the system (8)–(11) has a unique solution, and the solution belongs to Γ.

Definition

If , then in the system (8)–(11), infected individuals will experience extinction.

Theorem 1

If and , then the number of infected individuals in the system (8)–(11) exponentially tends to zero.

Assume that initial data satisfy and is a solution of system (8)–(11), with constant randomness σ and constant drift β, satisfying the stochastic differential equation for a Brownian motion W. Applying Itô’s lemma with gives More precisely, by integrating from , , Notice that a local continuous martingale is defined as with .

If , If , then so that, when we get as desired. Moreover,  □

Numerical methodology

For each , define the set . 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 . Fix a uniform partition of the temporal interval consisting of N subintervals, and let . Defining , for each , we will employ the notation , , , and to represent the numerical approximation to the values of the functions S, I, T, and R, respectively, at time .

Also, we will apply different explicit and implicit methods to the given system. As expected, we will employ the discrete initial data , where , , , .

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 .

Stochastic Runge–Kutta method

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

Stage 1

Stage 2

Stage 3

Stage 4

Final stage where k is the time step size and .

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 .

Convergence analysis

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

Theorem 2

For any given initial value , the system of Eqs. (17) to (20) has a unique positive solution for any .

Proof

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

Theorem 3

The region for all is a positively invariant feasible region for Eqs. (17) to (20).

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

Theorem 4

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

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

,

,

,

Theorem 5

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

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

The Jacobian matrix is defined as where Now we want to linearize the model about the steady-state of the model for virus-free state and .

The given Jacobian there is The eigenvalues of this Jacobian matrix are So, the given model is stable around 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 . We plot the largest eigenvalues of the Jacobian matrix associated with . Notice that the spectral radius is always less than one, as desired.  □

Figure 1

Graph of the spectral radius of the Jacobian matrix associated with the computational method (17)–(20) at the steady-state . The estimated parameter values of Table 1 were used, and the spectral radius was obtained for various temporal time step sizes k in

Computational results

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], , , , , , , for . For , , , , , , . Moreover, we used , , , .

Example 1

(Simulation for the free of coronavirus state)

For given real data, the reproductive value is . So, the system converges to . 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 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.

Figure 2

Comparison graphs of computational methods at the free of coronavirus state

Example 2

(Simulation for the existence of coronavirus state)

The reproductive value is . So, the system converges to . 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 for a certain small temporal step size. On the other hand, in Fig. 3, (b), (d), and (f), the results violated the structural properties of the system. Figure 3 exhibits the fact that the nonstandard computational method has kept the structural properties of the model such as positivity, boundedness, and dynamical consistency, as desired.

Figure 4

Graph of the effect of treatment on reproduction number

Figure 3

Comparison graph of computational methods at the existence of coronavirus state

Example 3

(Simulation for reproduction number with the effect of treatment)

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

Figure 5

Graph of the effect of quarantine on infected humans at the existence of a virus state

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.