Prediction of COVID-19 spread with models in different patterns: A case study of Russia

Mehmet Akif Cetin; Seda Igret Araz

doi:10.1515/phys-2024-0009

Article Open Access

Prediction of COVID-19 spread with models in different patterns: A case study of Russia

Mehmet Akif Cetin and Seda Igret Araz

Published/Copyright: April 29, 2024

Published by

Become an author with De Gruyter Brill

Submit Manuscript Author Information Explore this Subject

From the journal Open Physics Volume 22 Issue 1

Abstract

This study deals with a mathematical model that examines the spread of Coronavirus disease (COVID-19). This model has been handled with different processes such as deterministic, stochastic, and deterministic–stochastic. First of all, a detailed analysis is presented for the deterministic model, which includes the positivity of the solution, the basic reproduction number, the disease, and endemic equilibrium points. Then, for the stochastic model, we investigate under which conditions, the solution exists and is unique. Later, model is reconsidered with the help of the piecewise derivative, which can combine deterministic and stochastic processes. Numerical simulations are presented for all these processes. Finally, the model has been modified with the rate indicator function. The model presenting these four different situations is compared with the real data in Russia. According to the results obtained from these situations, the model that is obtained by adding the rate indicator function predicts the COVID-19 outbreak in Russia more accurately. Thus, it is concluded that the model with the rate indicator function presents more realistic approach than the previous ones.

Keywords: COVID-19 spread; fractional differentiation and integration; rate indicator function

1 Introduction

The novel Coronavirus disease (COVID-19) was first encountered in a group of patients who developed respiratory symptoms (fever, cough, shortness of breath) in late December in Wuhan Province, China [1]. As a result of research, it was a virus identified on January 13, 2020. The epidemic was initially detected in seafood and animal markets in Wuhan and then spread from person to person, spreading to other provinces in China, especially Wuhan, and then to the whole world. Although the countries sometimes take measures, the increases in cases have not been prevented. As a result, all countries have faced the waves, even most countries have witnessed the fifth wave. During this long time, numerous mathematical models are studied and epidemics in many countries are tried to be predicted. Sometimes deterministic models and sometimes stochastic models are used to predict and interpret epidemics [2–6]. In addition to classical order models, fractional order models are also discussed and the dynamics of the epidemic have been examined [7–12].

Some studies in the literature are as follows: Masandawa et al. have considered a mathematical model that combines both public control measures and health workers [13]. Ariffin et al. have determined the severity of the disease and the approximate number of days needed for the disease trend to decline in Malaysia, which has been severely affected by the COVID-19 outbreak [14]. Premarathna et al. have made estimations of the future number of COVID-19 cases in Sri Lanka [15]. Lotfi et al. developed a regression-based robust optimization (RO) approach with a mean absolute deviation objective function to estimate the number of patients affected by the COVID-19 outbreak in Iran [16]. Asempapa et al., using the COVID-19 data from Brazil and South African countries, created a mathematical model between the two subgroups that they classified as low risk and high risk and showed that it is possible to prevent the disease with drug-free interventions [17]. Samui et al., using India’s COVID-19 epidemic data within a certain period, have found a compartmental mathematical model to find the causes of disease transmission and control the epidemic [18]. Alqarni et al. have constructed a new mathematical model for the transmission dynamics of COVID-19 using disease cases reported in the Kingdom of Saudi Arabia over a given time period [19]. Grassly et al. developed a mathematical model that aims to investigate the potential impact of different polymerase chain reaction tests and isolation strategies on the transmission of the severe acute respiratory syndrome SARS-CoV-2 [20].

The objective of this research is to design a mathematical model that can effectively forecast the number of COVID-19 cases in Russia within a given time frame. There are three situations to consider: a deterministic process for the epidemic, a stochastic process for the epidemic, and a combination of both. In addition to these situations, the last one is the situation where the rate indicator function, which is very successful in estimating the numbers of infected, is added to a model that will be presented in the next section. In this study, we compare all of these patterns that will estimate the number of cases in Russia [1] which has faced a new wave recently, and determine the best prediction for the infected cases in this country.

Some definitions of fractional derivatives will be now presented. We start with the definition of the fractional derivative with power-law kernel [7] is given as

(1) D t α 0 C f ( t ) = 1 Γ ( 1 − α ) ∫ 0 t f ′ ( τ ) ( t − τ ) − α d τ ,

where 0 < α ≤ 1 . The definition of the fractional derivative with exponential decay kernel [8] is presented as

(2) D t α 0 C F f ( t ) = M ( α ) 1 − α ∫ 0 t f ′ ( τ ) exp − α 1 − α ( t − τ ) d τ ,

where 0 < α ≤ 1 . Here, M ( α ) is the normalization function such that M ( 0 ) = M ( 1 ) = 1 . The definition of the fractional derivative with Mittag–Leffler kernel [9] is defined by

(3) D t α 0 A B C f ( t ) = A B ( α ) 1 − α ∫ 0 t f ′ ( τ ) E α − α 1 − α ( t − τ ) α d τ ,

where 0 < α ≤ 1 and

(4) E α ( z ) = ∑ k = 0 ∞ z k Γ ( α k + 1 ) ,

is the Mittag–Leffler function. Here,

(5) A B ( α ) = 1 − α + α Γ ( α )

is the normalization function, noting that A B ( 0 ) = A B ( 1 ) = 1 . The definitions for fractional integrals, incorporating power law, exponential decay, and Mittag–Leffler kernel are as follows:

(6) I t α 0 C f ( t ) = 1 Γ ( α ) ∫ 0 t f ( τ ) ( t − τ ) α − 1 d τ , I t α 0 C F f ( t ) = 1 − α M ( α ) f ( t ) + α M ( α ) ∫ 0 t f ( τ ) d τ , I t α 0 A B C f ( t ) = 1 − α A B ( α ) f ( t ) + α A B ( α ) Γ ( α ) ∫ 0 t f ( τ ) ( t − τ ) α − 1 d τ .

2 Model formulation

In this section, we consider a mathematical model to understand the dynamics of the COVID-19 pandemic in Russia. Employing this model [21], we will be able to have an idea about the rate of spread of the disease and predict the number of future cases. The mathematical equations under investigation are represented by

(7) S ˙ = ρ − β s ( 1 − ψ ξ ) S I s − β a ( 1 − ψ ξ ) S I a − μ S , E ˙ = β s ( 1 − ψ ξ ) S I s + β a ( 1 − ψ ξ ) S I a − ( φ + μ ) E , I ˙ s = ( 1 − θ ) φ E − ( α s + δ s + λ s + μ ) I s , I ˙ a = θ φ E − ( α a + λ a + μ ) I a , Q ˙ = α s I s + α a I a − ( λ q + δ q + μ ) Q , R ˙ = λ s I s + λ a I a + λ q Q − μ R , F ˙ = δ s I s + δ q Q ,

where initial conditions are as follows:

(8) S ( 0 ) = S 0 , E ( 0 ) = E 0 , I s ( 0 ) = I s 0 , I a ( 0 ) = I a 0 , Q ( 0 ) = Q 0 , R ( 0 ) = R 0 , and F ( 0 ) = F 0 .

The model discussed in this article will be divided into seven compartments. Let N ( t ) be the total human population to be considered. It is important to note that the nonlinear terms S*I_s and S*I_s are assumed to be divided by N for simplicity in notation. At any given time, the number of susceptible is called S , the number of exposed E , the number of symptomatically infected I s , the number of asymptomatically infected I a , the number of quarantined Q , the number of recovered R , and the number of death F . Now, let us talk about the parameters that we will use in the model and what they mean. The recruitment rate into the susceptible compartment ρ , the effective contact rate β s , the effectiveness of social distancing β a , proportion of people wearing face masks ψ , the effectiveness of face masks ξ , proportion of exposed individuals showing clinical symptoms after the incubation period 1 − θ , progression rate from exposed compartment to infectious compartment φ , isolation rate for individuals in symptomatically infected compartment α s , recovery rate of individuals in symptomatically infected compartment λ s , COVID-19 disease mortality rate for individuals in the infectious compartment δ s , isolation rate of asymptomatically contagious individuals α a , recovery rate of asymptomatically infectious individuals λ a , recovery rate of individuals in the quarantined compartment λ q , COVID-19 disease mortality rate for individuals in quarantined compartment δ q and natural death rate of all individuals is denoted by μ .

The COVID-19 model is considered in a feasible region Ω such that

(9) Ω = ( S ( t ) , E ( t ) , I s ( t ) , I a ( t ) , Q ( t ) , R ( t ) , F ( t ) ) ∈ R t 7 : N ( t ) ≤ ρ μ .

2.1 Positivity and boundness

In this subsection, we shall show that the solutions of the model are nonnegative and bounded if the initial conditions are positive for all t > 0 . Assuming that the classes S ( t ) and I s ( t ) have same signs and ψ ξ < 1 , for E ( t ) class, we can have the following inequality:

(10) E ( t ) ≥ E 0 e − ( φ + μ ) t , ∀ t > 0 .

Since E ( t ) is positive, we write

(11) I ˙ s = ( 1 − θ ) φ E − ( α s + δ s + λ s + μ ) I s ≥ − ( α s + δ s + λ s + μ ) I s .

Then, we have

(12) I s ( t ) ≥ I s 0 e − ( α s + δ s + λ s + μ ) t , ∀ t > 0 .

Using the same routine, we have

(13) I a ( t ) ≥ I a 0 e − ( α a + λ a + μ ) t , ∀ t > 0 , Q ( t ) ≥ Q 0 e − ( λ q + δ q + μ ) t , ∀ t > 0 , R ( t ) ≥ R 0 e − μ t , ∀ t > 0 .

Assuming that the classes I s ( t ) and Q ( t ) are integrable, then we have

(14) F ( t ) ≥ F 0 + ∫ ( δ s I s + δ q Q ) d t , ∀ t > 0 .

Finally, we will discuss the positivity for S ( t ) class. For this, we need to define the following norm:

(15) ∥ f ∥ ∞ = sup t ∈ D f ∣ f ∣ .

Hence, for the class S , we have

(16) S ˙ ( t ) = ρ − β s ( 1 − ψ ξ ) S I s − β a ( 1 − ψ ξ ) S I a − μ S ≥ − ( β s ( 1 − ψ ξ ) I s + β a ( 1 − ψ ξ ) I a + μ ) S ≥ − ( β s ( 1 − ψ ξ ) ∣ I s ∣ + β a ( 1 − ψ ξ ) ∣ I a ∣ + μ ) S ≥ − ( β s ( 1 − ψ ξ ) sup t ∈ D I s ∣ I s ∣ + β a ( 1 − ψ ξ ) sup t ∈ D I a ∣ I a ∣ + μ ) S ≥ − ( β s ( 1 − ψ ξ ) ∥ I s ∥ ∞ + β a ( 1 − ψ ξ ) ∥ I a ∥ ∞ + μ ) S ≥ − ϰ S ,

where

(17) ϰ = β s ( 1 − ψ ξ ) ∥ I s ∥ ∞ + β a ( 1 − ψ ξ ) ∥ I a ∥ ∞ + μ .

Thus, we write

(18) S ( t ) ≥ S 0 e − ϰ t , ∀ t > 0 .

However, we know that the addition of all equations of model leads to

(19) lim t → ∞ N ( t ) = ρ μ .

Hence, the positive solutions of the model are bounded.

2.2 Equilibrium points of the model

In this subsection, we present the disease-free and endemic equilibrium points of Model (7). To do this, we set the right-hand side of each equation in Model (7) equal to zero. Thus, the disease-free equilibrium point is found as

(20) E 0 = ρ μ , 0 , 0 , 0 , 0 , 0 , 0 ,

and the endemic equilibrium point is calculated as

(21) S ∗ = ρ β s ( 1 − ψ ξ ) I s + β a ( 1 − ψ ξ ) I a + μ , E ∗ = β s ( 1 − ψ ξ ) S I s + β a ( 1 − ψ ξ ) S I a φ + μ , I s ∗ = ( 1 − θ ) φ E α s + δ s + λ s + μ , I a ∗ = θ φ E α a + λ a + μ , Q ∗ = α s I s + α a I a λ q + δ q + μ , R ∗ = λ s I s + λ a I a + λ q Q μ .

2.3 Reproduction number

In this subsection, the next-generation method is used to obtain reproduction number, which serves to understand the dynamics of the disease [22]. To achieve this, we consider the classes E , I s , and I a . Employing the next-generation method, the following Jacobian matrices are evaluated:

(22) F ˜ = J F = 0 β s ( 1 − ψ ξ ) ρ μ β a ( 1 − ψ ξ ) ρ μ 0 0 0 0 0 0

and

(23) V ˜ − 1 = J V − 1 = 1 φ + μ 0 0 ( 1 − θ ) φ ( φ + μ ) k 1 1 k 1 0 θ φ ( φ + μ ) k 2 0 1 k 2 ,

where

(24) k 1 = α s + δ s + λ s + μ , k 2 = α a + λ a + μ .

Then, the spectral radius of the ( F ˜ V ˜ − 1 ) , in other words, the reproduction number, is derived as

(25) R 0 = ρ μ β s ( 1 − ψ ξ ) ( 1 − θ ) φ ( φ + μ ) k 1 + β a ( 1 − ψ ξ ) θ φ ( φ + μ ) k 2 .

3 Existence and uniqueness of solution of the model

In this section, the conditions under which the solution of the COVID-19 model with the Atangana–Baleanu fractional derivative exists and is unique will be investigated. We start with examining the existence of the solution of the associated model.

3.1 Existence of the solution

In this section, we show that the solution of the considered model exists and is unique. To do this, we consider the following model with Mittag–Leffler kernel:

(26) D t α 0 A B C S = ρ − β s ( 1 − ψ ξ ) S I s − β a ( 1 − ψ ξ ) S I a − μ S , D t α 0 A B C E = β s ( 1 − ψ ξ ) S I s + β a ( 1 − ψ ξ ) S I a − ( φ + μ ) E , D t α 0 A B C I s = ( 1 − θ ) φ E − ( α s + δ s + λ s + μ ) I s , D t α 0 A B C I a = θ φ E − ( α a + λ a + μ ) I a , D t α 0 A B C Q = α s I s + α a I a − ( λ q + δ q + μ ) Q , D t α 0 A B C R = λ s I s + λ a I a + λ q Q − μ R , D t α 0 A B C F = δ s I s + δ q Q ,

with the initial conditions

(27) S ( 0 ) = S 0 , E ( 0 ) = E 0 , I s ( 0 ) = I s 0 , I a ( 0 ) = I a 0 , Q ( 0 ) = Q 0 , R ( 0 ) = R 0 , and F ( 0 ) = F 0 .

To avoid complexity, we will use some notations

(28) Y = S E I s I a Q R F , F ˜ = F 1 F 2 F 3 F 4 F 5 F 6 F 7 = ρ − β s ( 1 − ψ ξ ) S I s − β a ( 1 − ψ ξ ) S I a − μ S β s ( 1 − ψ ξ ) S I s + β a ( 1 − ψ ξ ) S I a − ( φ + μ ) E ( 1 − θ ) φ E − ( α s + δ s + λ s + μ ) I s θ φ E − ( α a + λ a + μ ) I a α s I s + α a I a − ( λ q + δ q + μ ) Q λ s I s + λ a I a + λ q Q − μ R δ s I s + δ q Q .

Theorem 1

The kernels F 1 , F 2 , … , F 7 in the aforementioned system hold the Lipschitz condition and contraction if the following inequalities satisfy

(29) 0 ≤ β s ( 1 − ψ ξ ) K 1 + β a ( 1 − ψ ξ ) K 2 + μ < 1 , 0 ≤ ( φ + μ ) < 1 , 0 ≤ ( α s + δ s + λ s + μ ) < 1 , 0 ≤ ( α a + λ a + μ ) < 1 , 0 ≤ ( λ q + δ q + μ ) < 1 , 0 ≤ μ < 1 , 0 ≤ K 7 ′ < 1 .

Proof

For proof, we start with the function F 1 . Let S 1 and S 2 be two functions. Then, we write

(30) ∥ F 1 ( t , S 1 ) − F 1 ( t , S 2 ) ∥ = ∥ − ( β s ( 1 − ψ ξ ) I s + β a ( 1 − ψ ξ ) I a + μ ) ( S 1 − S 2 ) ∥ ≤ ( β s ( 1 − ψ ξ ) ∥ I s ∥ + β a ( 1 − ψ ξ ) ∥ I a ∥ + μ ) ∥ S 1 − S 2 ∥ ≤ ( β s ( 1 − ψ ξ ) K 1 + β a ( 1 − ψ ξ ) K 2 + μ ) ∥ S 1 − S 2 ∥ ≤ K 1 ′ ∥ S 1 − S 2 ∥ ,

where

(31) K 1 = sup t ∈ D I ∥ I s ∥ , K 2 = sup t ∈ D I ∥ I a ∥ K 1 ′ = ( β s ( 1 − ψ ξ ) K 1 + β a ( 1 − ψ ξ ) K 2 + μ ) .

Similarly,

(32) ∥ F 2 ( t , E 1 ) − F 2 ( t , E 2 ) ∥ ≤ ( φ + μ ) ∥ E 1 − E 2 ∥ ≤ K 2 ′ ∥ E 1 − E 2 ∥ , ∥ F 3 ( t , I s 1 ) − F 3 ( t , I s 2 ) ∥ ≤ ( α s + δ s + λ s + μ ) ∥ I s 1 − I s 2 ∥ ≤ K 3 ′ ∥ I s 1 − I s 2 ∥ , ∥ F 4 ( t , I a 1 ) − F 4 ( t , I a 2 ) ∥ ≤ ( α a + λ a + μ ) ∥ I a 1 − I a 2 ∥ ≤ K 4 ′ ∥ I a 1 − I a 2 ∥ , ∥ F 5 ( t , Q 1 ) − F 5 ( t , Q 2 ) ∥ ≤ ( λ q + δ q + μ ) ∥ Q 1 − Q 2 ∥ ≤ K 5 ′ ∥ Q 1 − Q 2 ∥ , ∥ F 6 ( t , R 1 ) − F 6 ( t , R 2 ) ∥ ≤ μ ∥ R 1 − R 2 ∥ ≤ K 6 ′ ∥ R 1 − R 2 ∥ , ∥ F 7 ( t , F 1 ) − F 7 ( t , F 2 ) ∥ ≤ K 7 ′ ∥ F 1 − F 2 ∥ .

In the general form, System (26) can be converted to

(33) Y ( t ) = Y 0 + 1 − α A B ( α ) F i ( t , Y ) + α A B ( α ) Γ ( α ) ∫ 0 t ( t − τ ) α − 1 F 1 ( τ , Y ) d τ .

Using recursive formula yields

(34) Y n ( t ) = Y 0 + 1 − α A B ( α ) F i ( t , Y n − 1 ) + α A B ( α ) Γ ( α ) ∫ 0 t ( t − τ ) α − 1 F i ( τ , Y n − 1 ) d τ .

Taking the difference between two successive iterations leads to

(35) Y ∗ ( t ) = Y n ( t ) − Y n − 1 ( t ) = 1 − α A B ( α ) [ F i ( t , Y n − 1 ) − F 1 ( t , Y n − 2 ) ] + α A B ( α ) Γ ( α ) ∫ 0 t ( t − τ ) α − 1 [ F 1 ( τ , Y n − 1 ) − F 1 ( τ , Y n − 2 ) ] d τ ,

where

(36) Y ( t ) = ∑ n = 0 ∞ Y n ∗ ( t ) .

Applying the norm on both sides of the aforementioned equation, we can have

(37) ∥ Y ∗ ( t ) ∥ = ∥ Y n ( t ) − Y n − 1 ( t ) ∥ ≤ 1 − α A B ( α ) ∥ F i ( t , Y n − 1 ) − F i ( t , Y n − 2 ) ∥ + α A B ( α ) Γ ( α ) ∫ 0 t ( t − τ ) α − 1 [ F i ( τ , Y n − 1 ) − F i ( τ , Y n − 2 ) ] d τ .

Since the functions F 1 , F 2 , … , F 7 hold the Lipschitz condition, we have the following:

(38) ∥ Y n ∗ ( t ) ∥ ≤ 1 − α A B ( α ) K i ′ ∥ Y n − 1 ∗ ( t ) ∥ + α A B ( α ) Γ ( α ) K i ′ ∫ 0 t ( t − τ ) α − 1 ∥ Y n − 1 ∗ ( τ ) ∥ d τ ,

which completes the proof.□

Theorem 2

System (26) has a solution under the conditions that we have from Eq. (38):

(39) K i ′ 1 − α A B ( α ) + T α A B ( α ) Γ ( α ) < 1 , i = 1 , 2 , … , 7 .

Proof

From (38), we write

(40) ∥ Y n ∗ ( t ) ∥ ≤ ∥ S ( 0 ) ∥ K 1 ′ 1 − α A B ( α ) + T α A B ( α ) Γ ( α ) n .

We now show that the functions in the aforementioned equation are solutions of the model. To do this, we assume

(41) Y ( t ) − Y ( 0 ) = Y n ( t ) − Y ^ n ( t ) ,

where Y ^ n ( t ) are the reminder terms of series solution. Then, we need to show

(42) ∥ Y ^ n ( t ) ∥ → 0 .

Then, we write

(43) ∥ Y ^ n ( t ) ∥ ≤ 1 − α A B ( α ) ∥ F i ( t , Y ) − F i ( t , Y n − 1 ) ∥ + α A B ( α ) Γ ( α ) ∫ 0 t ( t − τ ) α − 1 ∥ F i ( τ , Y ) − F i ( τ , Y n − 1 ) ∥ d τ ≤ K i ′ 1 − α A B ( α ) + t α A B ( α ) Γ ( α ) ∥ Y − Y n − 1 ∥ .

Applying this recursively yields

(44) ∥ Y ^ n ( t ) ∥ ≤ ( K i ′ ) n 1 − α A B ( α ) + t α A B ( α ) Γ ( α ) n + 1 ∥ Y − Y n − 1 ∥ ≤ ( K i ′ ) n 1 − α A B ( α ) + T α A B ( α ) Γ ( α ) n + 1 ∥ Y − Y n − 1 ∥ .

Applying lim n → ∞ on both sides of the aforementioned equations, we can have

(45) ∥ Y ^ n ( t ) ∥ → 0 ,

which completes the proof. We can say that the system has a solution.□

3.2 Uniqueness of the solution

To show that solution of the system is unique, we assume that the system has two solutions such that Y 1 = Y i 1 ( t ) and Y 2 = Y i 2 ( t ) , i = 1 , 2 , … , 7 . To complete our proof, we need to show Y i 1 ( t ) = Y i 2 ( t ) . For this aim, we consider the following integral equations:

(46) Y i 1 ( t ) = Y i , 0 1 ( t ) + 1 − α A B ( α ) F i ( t , Y i 1 ) + α A B ( α ) Γ ( α ) ∫ 0 t ( t − τ ) α − 1 F i ( τ , Y i 1 ) d τ

and

(47) Y i 2 ( t ) = Y i , 0 2 ( t ) + 1 − α A B ( α ) F i ( t , Y i 2 ) + α A B ( α ) Γ ( α ) ∫ 0 t ( t − τ ) α − 1 F i ( τ , Y i 2 ) d τ .

Taking difference of the aforementioned equalities

(48) Y i 1 ( t ) − Y i 2 ( t ) = 1 − α A B ( α ) [ F i ( t , Y i 1 ) − F i ( t , Y i 2 ) ] + α A B ( α ) Γ ( α ) ∫ 0 t ( t − τ ) α − 1 [ F i ( τ , Y i 1 ) − F i ( τ , Y i 2 ) ] d τ .

Taking absolute value on both sides of Eq. (48) yields

(49) ∣ Y i 1 − Y i 2 ∣ ≤ 1 − α A B ( α ) ∣ F i ( t , Y i 1 ) − F i ( t , Y i 2 ) ∣ + α A B ( α ) Γ ( α ) ∫ 0 t ( t − τ ) α − 1 ∣ F i ( τ , Y i 1 ) − F i ( τ , Y i 2 ) ∣ d τ ,

and consider the Lipschitz condition

(50) ∣ Y i 1 ( t ) − Y i 2 ( t ) ∣ ≤ 1 − α A B ( α ) L i ∣ Y i 1 − Y i 2 ∣ + t α A B ( α ) Γ ( α ) L i ∣ Y i 1 − Y i 2 ∣ .

Arranging the aforementioned inequality leads to

(51) ∣ Y i 1 ( t ) − Y i 2 ( t ) ∣ 1 − 1 − α A B ( α ) + T α A B ( α ) Γ ( α ) L i ≤ 0 .

Hence, we have

(52) ∣ Y i 1 ( t ) − Y i 2 ( t ) ∣ = 0 ⇒ Y i 1 ( t ) = Y i 2 ( t ) , i = 1 , 2 , … , 7 .

Then, we can conclude that the model has a unique solution.

4 Numerical Solution of COVID-19 model with Mittag–Leffler kernel

Throughout this section, we construct a numerical scheme with Newton polynomial for the model of COVID-19 spread where the kernel is Mittag–Leffler

(53) D t α 0 A B C S ( t ) = ρ − β s ( 1 − ψ ξ ) S I s − β a ( 1 − ψ ξ ) S I a − μ S , D t α 0 A B C E ( t ) = β s ( 1 − ψ ξ ) S I s + β a ( 1 − ψ ξ ) S I a − ( φ + μ ) E , D t α 0 A B C I s ( t ) = ( 1 − θ ) φ E − ( α s + δ s + λ s + μ ) I s , D t α 0 A B C I a ( t ) = θ φ E − ( α a + λ a + μ ) I a , D t α 0 A B C Q ( t ) = α s I s + α a I a − ( λ q + δ q + μ ) Q , D t α 0 A B C R ( t ) = λ s I s + λ a I a + λ q Q − μ R , D t α 0 A B C F ( t ) = δ s I s + δ q Q .

We can facilitate the aforementioned system as

(54) D t α 0 A B C S ( t ) = f 1 ( S , E , I s , I a , Q , R , F , t ) , D t α 0 A B C E ( t ) = f 2 ( S , E , I s , I a , Q , R , F , t ) , D t α 0 A B C I s ( t ) = f 3 ( S , E , I s , I a , Q , R , F , t ) , D t α 0 A B C I a ( t ) = f 4 ( S , E , I s , I a , Q , R , F , t ) , D t α 0 A B C Q ( t ) = f 5 ( S , E , I s , I a , Q , R , F , t ) , D t α 0 A B C R ( t ) = f 6 ( S , E , I s , I a , Q , R , F , t ) , D t α 0 A B C F ( t ) = f 7 ( S , E , I s , I a , Q , R , F , t ) .

Applying the associated integral, we obtain the following:

(55) Y ( t ) = 1 − α A B ( α ) F i ( Y , t ) + α A B ( α ) Γ ( α ) × ∫ 0 t ( t − τ ) α − 1 F i ( Y , τ ) d τ , i = 1 , … , 7 .

If the functions f 1 , f 2 , f 3 , f 4 , f 5 , f 6 , and f 7 can be approximated within the interval [ t k , t k + 1 ] using Newton polynomial, the following scheme can be achieved

(56) Y ( t n + 1 ) = Y ( 0 ) + 1 − α A B ( α ) F i ( Y n , t n ) + α ( Δ t ) α A B ( α ) Γ ( α + 1 ) ∑ k = 2 n F i ( Y k − 2 , t k − 2 ) A 1 + α ( Δ t ) α A B ( α ) Γ ( α + 2 ) ∑ k = 2 n [ F i ( Y k − 1 , t k − 1 ) − F i ( Y k − 2 , t k − 2 ) ] A 2 + α ( Δ t ) α 2 A B ( α ) Γ ( α + 3 ) ∑ k = 2 n [ F i ( Y k , t k ) − 2 F i ( Y k − 1 , t k − 1 ) + F i ( Y k − 2 , t k − 2 ) ] A 3 ,

where

(57) A 1 = [ ( n − k + 1 ) α − ( n − k ) α ] , A 2 = ( n − k + 1 ) α ( n − k + 3 + 2 α ) − ( n − k ) α ( n − k + 3 + 3 α ) , A 3 = ( n − k + 1 ) α ( 2 ( n − k ) 2 + ( 3 α + 10 ) ( n − k ) + 2 α 2 + 9 α + 12 ) − ( n − k ) α ( 2 ( n − k ) 2 + ( 5 α + 10 ) ( n − k ) + 6 α 2 + 18 α + 12 ) .

5 Numerical simulation

In this section, we present the numerical simulations for our model with the following Atangana–Baleanu fractional derivative:

(58) D t α 0 A B C S ( t ) = ρ − β s ( 1 − ψ ξ ) S I s − β a ( 1 − ψ ξ ) S I a − μ S , D t α 0 A B C E ( t ) = β s ( 1 − ψ ξ ) S I s + β a ( 1 − ψ ξ ) S I a − ( φ + μ ) E , D t α 0 A B C I s ( t ) = ( 1 − θ ) φ E − ( α s + δ s + λ s + μ ) I s , D t α 0 A B C I a ( t ) = θ φ E − ( α a + λ a + μ ) I a , D t α 0 A B C Q ( t ) = α s I s + α a I a − ( λ q + δ q + μ ) Q , D t α 0 A B C R ( t ) = λ s I s + λ a I a + λ q Q − μ R , D t α 0 A B C F ( t ) = δ s I s + δ q Q .

Here, the initial conditions are considered as

(59) S ( 0 ) = 10,000 , E ( 0 ) = 1,750 , I s ( 0 ) = 1,950 , I a ( 0 ) = 1,965 , Q ( 0 ) = 90 , R ( 0 ) = 0 , F ( 0 ) = 10 ,

and the parameters are

(60) ρ = 5 , β s = 0.1 , ψ = 0.1 , ξ = 0.5 , β a = 0.1 , μ = 3.6529 × 1 0 − 5 , φ = 1 6 , θ = 0.5 , α s = 0.2 , δ s = 0.015 , λ s = 1 10 , α a = 0.2 , δ a = 0.015 , λ a = 1 10 , λ q = 0.05 , δ q = 0.015 .

For the deterministic system of COVID-19 spread, we present the numerical simulations in Figures 1, 2, 3, 4.

Figure 1

Numerical visualization for the class of susceptible and exposed individuals.

Figure 2

Numerical visualization for the class of symptomatically and asymptomatically infected individuals.

Figure 3

Numerical visualization for the class of quarantined and recovered individuals.

Figure 4

Numerical visualization for the class of died individuals.

6 Stochastic model of COVID-19 spread

In this section, we consider the following stochastic COVID-19 model:

(61) d S = [ ρ − β s ( 1 − ψ ξ ) S I s − β a ( 1 − ψ ξ ) S I a − μ S ] d t + σ 1 S ( t ) d B 1 ( t ) , d E = [ β s ( 1 − ψ ξ ) S I s + β a ( 1 − ψ ξ ) S I a − ( φ + μ ) E ] d t + σ 2 E ( t ) d B 2 ( t ) , d I s = [ ( 1 − θ ) φ E − ( α s + δ s + λ s + μ ) I s ] d t + σ 3 I s ( t ) d B 3 ( t ) , d I a = [ θ φ E − ( α a + λ a + μ ) I a ] d t + σ 4 I a ( t ) d B 4 ( t ) , d Q = [ α s I s + α a I a − ( λ q + δ q + μ ) Q ] d t + σ 5 Q ( t ) d B 5 ( t ) , d R = [ λ s I s + λ a I a + λ q Q − μ R ] d t + σ 6 R ( t ) d B 6 ( t ) , d F = [ δ s I s + δ q Q ] d t + σ 7 F ( t ) d B 7 ( t ) ,

where σ i and B i ( t ) , { i = 1 , … , 7 } are the stochastic component and standard Brownian function, respectively.

6.1 Existence and uniqueness of the model

In this subsection, we prove that the solution of Model (26) exists and is unique. For simplicity, let us rewrite Model (26) as follows:

(62) d S = u 1 ( t , S ) d t + C 1 ( t , S ) d B 1 ( t ) , d E = u 2 ( t , E ) d t + C 2 ( t , E ) d B 2 ( t ) , d I s = u 3 ( t , I s ) d t + C 3 ( t , I s ) d B 3 ( t ) , d I a = u 4 ( t , I a ) d t + C 4 ( t , I a ) d B 4 ( t ) , d Q = u 5 ( t , Q ) d t + C 5 ( t , Q ) d B 5 ( t ) , d R = u 6 ( t , R ) d t + C 6 ( t , R ) d B 6 ( t ) , d F = u 7 ( t , F ) d t + C 7 ( t , F ) d B 7 ( t ) .

For simplicity, we shall take as C i ( t , x i ) = σ i x i , i = 1 , … , 7 . To achieve our aim, we need to show that linear growth condition and the Lipschitz condition are satisfied, i.e.,

(63) ∣ u i ( t , x i ) ∣ 2 , ∣ C i ( t , x i ) ∣ 2 < l i ( 1 + ∣ x i ∣ 2 ) , ∀ t ∈ [ 0 , T ] ,

and

(64) ∣ u i ( t , x i 1 ) − u i ( t , x i 2 ) ∣ , ∣ C i ( t , x i 1 ) − C i ( t , x i 2 ) ∣ < l ¯ i ∣ x i 1 − x i 2 ∣ , ∀ t ∈ [ 0 , T ] .

First, we shall prove the linear growth condition. For the function u 1 ( t , S ) , we have

(65) ∣ u 1 ( t , S ) ∣ 2 = ∣ ρ − β s ( 1 − ψ ξ ) S I s − β a ( 1 − ψ ξ ) S I a − μ S ∣ 2 ≤ 4 ( ρ 2 + β s 2 ( 1 − ψ ξ ) 2 ∣ S ∣ 2 ∣ I s ∣ 2 + β a 2 ( 1 − ψ ξ ) 2 ∣ S ∣ 2 ∣ I a ∣ 2 + μ 2 ∣ S ∣ 2 ) = 4 ( ρ 2 + ( β s 2 ( 1 − ψ ξ ) 2 ∣ I s ∣ 2 + β a 2 ( 1 − ψ ξ ) 2 ∣ I a ∣ 2 + μ 2 ) ∣ S ∣ 2 ) ≤ 4 ( ρ 2 + ( β s 2 ( 1 − ψ ξ ) 2 sup 0 ≤ t ≤ T ∣ I s 2 ∣ + β a 2 ( 1 − ψ ξ ) 2 sup 0 ≤ t ≤ T ∣ I a 2 ∣ + μ 2 ) ∣ S ∣ 2 ) ≤ 4 ( ρ 2 + ( β s 2 ( 1 − ψ ξ ) 2 ∥ I s ∥ ∞ 2 + β a 2 ( 1 − ψ ξ ) 2 ∥ I a ∥ ∞ 2 + μ 2 ) ∣ S ∣ 2 ) = 4 ρ 2 1 + ( β s 2 ( 1 − ψ ξ ) 2 ∥ I s ∥ ∞ 2 + β a 2 ( 1 − ψ ξ ) 2 ∥ I a ∥ ∞ 2 + μ 2 ) ρ 2 ∣ S ∣ 2 ,

under the condition

(66) m 1 = ( β s 2 ( 1 − ψ ξ ) 2 ∥ I s ∥ ∞ 2 + β a 2 ( 1 − ψ ξ ) 2 ∥ I a ∥ ∞ 2 + μ 2 ) ρ 2 < 1 .

Then, the following inequality holds:

(67) ∣ u 1 ( t , S ) ∣ 2 ≤ l 1 ( 1 + ∣ S ∣ 2 ) .

For the function u 2 ( t , E ) , we obtain

(68) ∣ u 2 ( t , E ) ∣ 2 = ∣ β s ( 1 − ψ ξ ) S I s + β a ( 1 − ψ ξ ) S I a − ( φ + μ ) E ∣ 2 ≤ 3 ( β s 2 ( 1 − ψ ξ ) 2 ∣ S ∣ 2 ∣ I s ∣ 2 + β a 2 ( 1 − ψ ξ ) 2 ∣ S ∣ 2 ∣ I a ∣ 2 + ( φ + μ ) 2 ∣ E ∣ 2 ) ≤ 3 β s 2 ( 1 − ψ ξ ) 2 sup 0 ≤ t ≤ T ∣ S 2 ∣ sup 0 ≤ t ≤ T ∣ I s 2 ∣ + β a 2 ( 1 − ψ ξ ) 2 sup 0 ≤ t ≤ T ∣ S 2 ∣ sup 0 ≤ t ≤ T ∣ I a 2 ∣ + ( φ + μ ) 2 ∣ E ∣ 2 ≤ 3 ( β s 2 ( 1 − ψ ξ ) 2 ∥ S ∥ ∞ 2 ∥ I s ∥ ∞ 2 + β a 2 ( 1 − ψ ξ ) 2 ∥ S ∥ ∞ 2 ∥ I a ∥ ∞ 2 + ( φ + μ ) 2 ∣ E ∣ 2 ) ≤ 3 ( 1 − ψ ξ ) 2 ∥ S ∥ ∞ 2 [ β s 2 ∥ I s ∥ ∞ 2 + β a 2 ∥ I a ∥ ∞ 2 ] × 1 + ( φ + μ ) 2 ( 1 − ψ ξ ) 2 ∥ S ∥ ∞ 2 [ β s 2 ∥ I s ∥ ∞ 2 + β a 2 ∥ I a ∥ ∞ 2 ] ∣ E ∣ 2 ,

such that

(69) m 2 = ( φ + μ ) 2 ( 1 − ψ ξ ) 2 ∥ S ∥ ∞ 2 [ β s 2 ∥ I s ∥ ∞ 2 + β a 2 ∥ I a ∥ ∞ 2 ] < 1 ,

such that ψ ξ ≠ 1 . Thus, we have

(70) ∣ u 2 ( t , E ) ∣ 2 ≤ l 2 ( 1 + ∣ E ∣ 2 ) .

Similarly,

(71) ∣ u 3 ( t , I s ) ∣ 2 = ∣ ( 1 − θ ) φ E − ( α s + δ s + λ s + μ ) I s ∣ 2 ≤ 2 ( ( 1 − θ ) 2 φ 2 ∣ E ∣ 2 + ( α s + δ s + λ s + μ ) 2 ∣ I s ∣ 2 ) ≤ 2 ( ( 1 − θ ) 2 φ 2 sup 0 ≤ t ≤ T ∣ E 2 ∣ + ( α s + δ s + λ s + μ ) 2 ∣ I s ∣ 2 ) ≤ 2 ( 1 − θ ) 2 φ 2 ∥ E ∥ ∞ 2 × 1 + ( α s + δ s + λ s + μ ) 2 ( 1 − θ ) 2 φ 2 ∥ E ∥ ∞ 2 ∣ I s ∣ 2 ≤ l 3 ( 1 + ∣ I s ∣ 2 ) ,

under the condition

(72) m 3 = ( α s + δ s + λ s + μ ) 2 ( 1 − θ ) 2 φ 2 ∥ E ∥ ∞ 2 < 1 and θ ≠ 1 .

Next,

(73) ∣ u 4 ( t , I a ) ∣ 2 = ∣ θ φ E − ( α a + λ a + μ ) I a ∣ 2 ≤ 2 ( θ 2 φ 2 sup 0 ≤ t ≤ T ∣ E 2 ∣ + ( α a + λ a + μ ) 2 ∣ I a ∣ 2 ) ≤ 2 θ 2 φ 2 ∥ E 2 ∥ ∞ 1 + ( α a + λ a + μ ) 2 θ 2 φ 2 ∥ E 2 ∥ ∞ ∣ I a ∣ 2 ≤ l 4 ( 1 + ∣ I a ∣ 2 ) ,

where

(74) m 4 = ( α a + λ a + μ ) 2 θ 2 φ 2 ∥ E 2 ∥ ∞ < 1 .

Next,

(75) ∣ u 5 ( t , Q ) ∣ 2 = ∣ α s I s + α a I a − ( λ q + δ q + μ ) Q ∣ 2 ≤ 3 ( α s 2 sup 0 ≤ t ≤ T ∣ I s 2 ∣ + α a 2 sup 0 ≤ t ≤ T ∣ I a 2 ∣ + ( λ q + δ q + μ ) 2 ∣ Q ∣ 2 ) ≤ 3 ( α s 2 ∥ I s 2 ∥ ∞ + α a 2 ∥ I a 2 ∥ ∞ ) × 1 + ( λ q + δ q + μ ) 2 α s 2 ∥ I s 2 ∥ ∞ + α a 2 ∥ I a 2 ∥ ∞ ≤ l 5 ( 1 + ∣ Q ∣ 2 ) ,

under the condition

(76) m 5 = ( λ q + δ q + μ ) 2 α s 2 ∥ I s 2 ∥ ∞ + α a 2 ∥ I a 2 ∥ ∞ < 1 .

For the function u 6 ( t , R ) , we write

(77) ∣ u 6 ( t , R ) ∣ 2 = ∣ λ s I s + λ a I a + λ q Q − μ R ∣ 2 ≤ 4 ( λ s 2 sup 0 ≤ t ≤ T ∣ I s 2 ∣ + λ a 2 sup 0 ≤ t ≤ T ∣ I a 2 ∣ + λ q 2 sup 0 ≤ t ≤ T ∣ Q 2 ∣ + μ 2 ∣ R ∣ 2 ) ≤ 4 ( λ s 2 ∥ I s 2 ∥ ∞ + λ a 2 ∥ I a 2 ∥ ∞ + λ q 2 ∥ Q 2 ∥ ∞ ) × 1 + μ 2 λ s 2 ∥ I s 2 ∥ ∞ + λ a 2 ∥ I a 2 ∥ ∞ + λ q 2 ∥ Q 2 ∥ ∞ ∣ R ∣ 2 ≤ l 6 ( 1 + ∣ R ∣ 2 ) ,

where

(78) m 6 = μ 2 λ s 2 ∥ I s 2 ∥ ∞ + λ a 2 ∥ I a 2 ∥ ∞ + λ q 2 ∥ Q 2 ∥ ∞ < 1 .

Finally, we have

(79) ∣ u 7 ( t , F ) ∣ 2 = ∣ δ s I s + δ q Q ∣ 2 ≤ 2 ( δ s 2 ∥ I s 2 ∥ ∞ + δ q 2 ∥ Q 2 ∥ ∞ ) ( 1 + ∣ F ∣ 2 ) ≤ l 7 ( 1 + ∣ F ∣ 2 ) .

For the stochastic part,

(80) ∣ C 1 ( t , S ) ∣ 2 = σ 1 2 ∣ S ∣ 2 ≤ σ 1 2 ( 1 + ∣ S ∣ 2 ) , ∣ C 2 ( t , E ) ∣ 2 ≤ σ 2 2 ( 1 + ∣ E ∣ 2 ) , ∣ C 3 ( t , I s ) ∣ 2 ≤ σ 3 2 ( 1 + ∣ I s ∣ 2 ) , ∣ C 4 ( t , I a ) ∣ 2 ≤ σ 4 2 ( 1 + ∣ I a ∣ 2 ) , ∣ C 5 ( t , Q ) ∣ 2 ≤ σ 5 2 ( 1 + ∣ Q ∣ 2 ) , ∣ C 6 ( t , R ) ∣ 2 ≤ σ 6 2 ( 1 + ∣ R ∣ 2 ) , ∣ C 7 ( t , F ) ∣ 2 ≤ σ 7 2 ( 1 + ∣ F ∣ 2 ) .

Thus, we verified the first condition if the following is satisfied:

(81) max { m 1 , m 2 , m 3 , m 4 , m 5 , m 6 } < 1 .

We now prove the Lipschitz condition

(82) ∣ u 1 ( t , S 1 ) − u 1 ( t , S 2 ) ∣ 2 = ∣ ( − β s ( 1 − ψ ξ ) I s − β a ( 1 − ψ ξ ) I a − μ ) ( S 1 − S 2 ) ∣ 2 ≤ 3 ∣ ( 1 − ψ ξ ) 2 ( β s 2 ∣ I s 2 ∣ + β a 2 ∣ I a 2 ∣ ) + μ 2 ∣ ∣ S 1 − S 2 ∣ 2 ≤ 3 ∣ ( 1 − ψ ξ ) 2 ( β s 2 sup 0 ≤ t ≤ T ∣ I s 2 ∣ + β a 2 sup 0 ≤ t ≤ T ∣ I a 2 ∣ + μ 2 ) ∣ ∣ S 1 − S 2 ∣ 2 ≤ 3 ∣ ( 1 − ψ ξ ) 2 ( β s 2 ∥ I s 2 ∥ ∞ + β a 2 ∥ I a 2 ∥ ∞ + μ 2 ) ∣ ∣ S 1 − S 2 ∣ 2 ≤ l ¯ 1 ∣ S 1 − S 2 ∣ 2 ,

(83) ∣ u 2 ( t , E 1 ) − u 2 ( t , E 2 ) ∣ 2 = ∣ − ( φ + μ ) ( E 1 − E 2 ) ∣ 2 ≤ ∣ ε 1 2 + ( φ + μ ) 2 ∣ ∣ E 1 − E 2 ∣ 2 ≤ l ¯ 2 ∣ E 1 − E 2 ∣ 2 ,

(84) ∣ u 3 ( t , I s 1 ) − u 3 ( t , I s 2 ) ∣ 2 = ∣ − ( α s + δ s + λ s + μ ) ( I s 1 − I s 2 ) ∣ 2 ≤ ∣ ε 2 2 + ( α s + δ s + λ s + μ ) 2 ∣ ∣ I s 1 − I s 2 ∣ 2 ≤ l ¯ 3 ∣ I s 1 − I s 2 ∣ 2 ,

(85) ∣ u 4 ( t , I a 1 ) − u 4 ( t , I a 2 ) ∣ 2 = ∣ − ( α a + λ a + μ ) ( I a 1 − I a 2 ) ∣ 2 ≤ ∣ ε 3 2 + ( α a + λ a + μ ) 2 ∣ ∣ I a 1 − I a 2 ∣ 2 ≤ l ¯ 4 ∣ I a 1 − I a 2 ∣ 2 ,

(86) ∣ u 5 ( t , Q 1 ) − u 5 ( t , Q 2 ) ∣ 2 = ∣ − ( λ q + δ q + μ ) ( Q 1 − Q 2 ) ∣ 2 ≤ ∣ ε 4 2 + ( λ q + δ q + μ ) 2 ∣ ∣ Q 1 − Q 2 ∣ 2 ≤ l ¯ 5 ∣ Q 1 − Q 2 ∣ 2 ,

(87) ∣ u 6 ( t , R 1 ) − u 6 ( t , R 2 ) ∣ 2 = ∣ − μ ( R 1 − R 2 ) ∣ 2 ≤ ∣ ε 5 2 + μ 2 ∣ ∣ R 1 − R 2 ∣ 2 ≤ l ¯ 6 ∣ R 1 − R 2 ∣ 2 ,

(88) ∣ u 7 ( t , F 1 ) − u 7 ( t , F 2 ) ∣ 2 ≤ l ¯ 7 ∣ F 1 − F 2 ∣ 2 .

Also, we write

(89) ∣ C 1 ( t , S 1 ) − C 1 ( t , S 2 ) ∣ 2 ≤ 3 2 σ 1 2 ∣ S 1 − S 2 ∣ 2 , ∣ C 2 ( t , E 1 ) − C 2 ( t , E 2 ) ∣ 2 ≤ 3 2 σ 2 2 ∣ E 1 − E 2 ∣ 2 , ∣ C 3 ( t , I s 1 ) − C 3 ( t , I s 2 ) ∣ 2 ≤ 3 2 σ 3 2 ∣ I s 1 − I s 2 ∣ 2 , ∣ C 4 ( t , I a 1 ) − C 4 ( t , I a 2 ) ∣ 2 ≤ 3 2 σ 4 2 ∣ I a 1 − I a 2 ∣ 2 , ∣ C 5 ( t , Q 1 ) − C 5 ( t , Q 2 ) ∣ 2 ≤ 3 2 σ 5 2 ∣ Q 1 − Q 2 ∣ 2 , ∣ C 6 ( t , R 1 ) − C 6 ( t , R 2 ) ∣ 2 ≤ 3 2 σ 6 2 ∣ R 1 − R 2 ∣ 2 , ∣ C 7 ( t , F 1 ) − C 7 ( t , F 2 ) ∣ 2 ≤ 3 2 σ 7 2 ∣ F 1 − F 2 ∣ 2 .

Then, we can conclude that the system has a unique solution.

6.2 Numerical solution of stochastic COVID-19 model with Mittag–Leffler kernel

In this subsection, we present a numerical algorithm based on the Newton interpolation polynomial to solve stochastic model with Mittag–Leffler kernel:

(90) D t α 0 A B C S ( t ) = F 1 ( S , E , I s , I a , Q , R , F , t ) + σ 1 S ( t ) d B 1 ( t ) , D t α 0 A B C E ( t ) = F 2 ( S , E , I s , I a , Q , R , F , t ) + σ 2 E ( t ) d B 2 ( t ) , D t α 0 A B C I s ( t ) = F 3 ( S , E , I s , I a , Q , R , F , t ) + σ 3 I s ( t ) d B 3 ( t ) , D t α 0 A B C I a ( t ) = F 4 ( S , E , I s , I a , Q , R , F , t ) + σ 4 I a ( t ) d B 4 ( t ) , D t α 0 A B C Q ( t ) = F 5 ( S , E , I s , I a , Q , R , F , t ) + σ 5 Q ( t ) d B 5 ( t ) , D t α 0 A B C R ( t ) = F 6 ( S , E , I s , I a , Q , R , F , t ) + σ 6 R ( t ) d B 6 ( t ) , D t α 0 A B C F ( t ) = F 7 ( S , E , I s , I a , Q , R , F , t ) + σ 7 F ( t ) d B 7 ( t ) ,

subject to

(91) σ i = ( σ 1 , σ 2 , σ 3 , σ 4 , σ 5 , σ 6 , σ 7 ) , i = 1 , … , 7 , B i ( t ) = ( B 1 ( t ) , B 2 ( t ) , B 3 ( t ) , B 4 ( t ) , B 5 ( t ) , B 6 ( t ) , B 7 ( t ) ) .

The aforementioned model can be solved using the Newton interpolation polynomial as follows:

(92) Y ( t n + 1 ) = Y ( 0 ) + 1 − α A B ( α ) F i ( Y n , t n ) + 1 − α A B ( α ) σ 1 Y i ( t n ) ( B i ( t n + 1 ) − B i ( t n ) ) + α ( Δ t ) α A B ( α ) Γ ( α + 1 ) ∑ k = 2 n F i ( Y k − 2 , t k − 2 ) A 1 + α ( Δ t ) α A B ( α ) Γ ( α + 2 ) ∑ k = 2 n [ F i ( Y k − 1 , t k − 1 ) − F i ( Y k − 2 , t k − 2 ) ] A 2 + α ( Δ t ) α 2 A B ( α ) Γ ( α + 3 ) ∑ k = 2 n × F i ( Y k , t k ) − 2 F i ( Y k − 1 , t k − 1 ) + F i ( Y k − 2 , t k − 2 ) A 3 + α ( Δ t ) α − 1 A B ( α ) Γ ( α + 1 ) ∑ k = 0 n σ i Y i ( t k ) ( B i ( t k + 1 ) − B i ( t k ) ) A 1 .

The numerical simulations will be depicted using the same parameters in Eq. (60). Also, we consider the following initial data and stochastic constants:

(93) S ( 0 ) = 2,150 , E ( 0 ) = 1,750 , I s ( 0 ) = 3,930 , I a ( 0 ) = 1,965 , Q ( 0 ) = 94 , R ( 0 ) = 3,766 , F ( 0 ) = 10 , σ 1 = 0.0003 , σ 2 = 0.0011 , σ 3 = 0.0012 , σ 4 = 0.00021 , σ 5 = 0.0008 , σ 6 = 0.0003 , σ 7 = 0.0007 .

The numerical simulations for stochastic system of COVID-19 spread are shown in Figures 5, 6, 7, 8.

Figure 5

Numerical visualization for the class of susceptible and infected individuals.

Figure 6

Numerical visualization for the class of symptomatically and asymptomatically infected individuals.

Figure 7

Numerical visualization for the class of quarantined and recovered individuals.

Figure 8

Numerical visualization for the class of died individuals.

6.3 Numerical solution of COVID-19 model with piecewise derivative

In this section, we combine two concepts called deterministic and stochastic for the considered model using piecewise derivative. To do this, we consider a piecewise process where the spread displays three processes, namely, classical behaviors, stochastic behaviors, and finally nonlocal behaviors. In this case, if we consider T as the last time of the spread, this is to say, the last day where a new infection occurs, then, for the first period of time ranging from 0 to T 1 , the mathematical model will be constructed with classical differential operator, the second-phrase stochastic approach will be used, the model will be with the Atangana–Baleanu differential operator and finally at the last phase. The mathematical model explaining this dynamic is then presented as follows:

(94) d Y i d t = F ( t , Y ) , if 0 ≤ t ≤ T 1 , Y i ( 0 ) = Y i , 0 , i = 1 , 2 , … , 7 , D t α T 1 A B C Y i = F ( t , Y ) , if T 1 ≤ t ≤ T 2 , Y i ( T 1 ) = Y i , 1 , d Y ( t ) = F ( t , Y ) d t + σ i Y i d B i ( t ) , if T 2 ≤ t ≤ T , Y i ( T 2 ) = Y i , 2 .

Using the same routine, the numerical solution can be obtained as

(95) Y i n 1 = Y i ( 0 ) + ∑ j 1 = 2 n 1 23 12 F ( t j 1 , Y j 1 ) − 4 3 F ( t j 1 − 1 , Y j 1 − 1 ) + 5 12 F ( t j 1 − 2 , Y j 1 − 2 ) Δ t , 0 ≤ t ≤ T 1 , Y i n 2 = Y i ( T 1 ) + 1 − α A B ( α ) F ( t n 2 , Y n 2 ) + α Δ t A B ( α ) Γ ( α + 1 ) ∑ j 2 = n 1 + 3 n 2 F ( t j 2 − 2 , Y j 2 − 2 ) { ( n 2 − j 2 + 1 ) α − ( n 2 − j 2 ) α } + α Δ t A B ( α ) Γ ( α + 2 ) ∑ j 2 = n 1 + 3 n 2 [ F ( t j 2 − 1 , Y j 2 − 1 ) − F ( t j 2 − 2 , Y j 2 − 2 ) ] ( n 2 − j 2 + 1 ) α ( n 2 − j 2 + 3 + 2 α ) − ( n 2 − j 2 ) α ( n 2 − j 2 + 3 + 3 α ) + α Δ t 2 A B ( α ) Γ ( α + 3 ) ∑ j 2 = n 1 + 3 n 2 [ F ( t j 2 , Y j 2 ) − 2 F ( t j 2 − 1 , Y j 2 − 1 ) + F ( t j 2 − 2 , Y j 2 − 2 ) ] × ( n 2 − j 2 + 1 ) α 2 ( n 2 − j 2 ) 2 + ( 3 α + 10 ) ( n 2 − j 2 ) + 2 α 2 + 9 α + 12 − ( n 2 − j 2 ) α 2 ( n 2 − j 2 ) 2 + ( 5 α + 10 ) ( n 2 − j 2 ) + 6 α 2 + 18 α + 12 , T 1 ≤ t ≤ T 2 , Y i n 3 = Y i ( T 2 ) + ∑ j 3 = n 2 + 3 n 3 23 12 F ( t j 3 , Y j 3 ) − 4 3 F ( t j 3 − 1 , Y j 3 − 1 ) + 5 12 F ( t j 3 − 2 , Y j 3 − 2 ) Δ t + σ i ∑ j = i n 3 Y i j 3 ( B i j 3 + 1 − B i j 3 ) , T 2 ≤ t ≤ T .

The numerical simulations of model with piecewise differential and integral operators are performed for α = 0.9 in Figures 9, 10, 11, 12.

Figure 9

Numerical visualization for the class of susceptible and infected individuals.

Figure 10

Numerical visualization for the class of symptomatically and asymptomatically infected individuals.

Figure 11

Numerical visualization for the class of quarantined and recovered individuals.

Figure 12

Numerical visualization for the class of died individuals.

For another case, we can consider the following system:

(96) d Y ( t ) = F ( t , Y ) d t + σ i Y i d B i ( t ) , if 0 ≤ t ≤ T 1 , Y i ( 0 ) = Y i , 0 , i = 1 , 2 , … , 7 , D t α T 1 A B C Y i = F ( t , Y ) , if T 1 ≤ t ≤ T 2 , Y i ( T 1 ) = Y i , 1 , d Y ( t ) = F ( t , Y ) d t + σ i Y i d B i ( t ) , if T 2 ≤ t ≤ T , Y i ( T 2 ) = Y i , 2 .

The numerical simulations for the aforementioned model are given in Figures 13, 14, 15, 16.

Figure 13

Numerical visualization for the class of susceptible and infected individuals.

Figure 14

Numerical visualization for the class of symptomatically and asymptomatically infected individuals.

Figure 15

Numerical visualization for the class of quarantined and recovered individuals.

Figure 16

Numerical visualization for the class of died individuals.

7 Modification of COVID-19 model with the rate indicator function

In this section, we modify model presented in this study by the rate indicator function introduced by Atangana and Araz [23]. First, we shall present the definition of the rate indicator function. The rate between two days can be calculated as

(97) r i = f d i + 1 − f d i d i + 1 − d i ,

where f d i + 1 is the observed data of day d i + 1 . The daily rate indicator function can be defined by

(98) r ( t ) = f d i + r i ( t − d i ) = f d i + f d i + 1 − f d i d i + 1 − d i ( t − d i ) .

The considered model with classical derivative can be modified with the rate indicator function as follows:

(99) S ′ ( t ) = ρ − β s ( 1 − ψ ξ ) S I s − β a ( 1 − ψ ξ ) S I a − μ S , E ′ ( t ) = β s ( 1 − ψ ξ ) S I s + β a ( 1 − ψ ξ ) S I a − ( φ + μ ) E , I s ′ ( t ) = ( 1 − θ ) φ E − r ( t ) ( α s + δ s + λ s + μ ) I s , I a ′ ( t ) = θ φ E − ( α a + λ a + μ ) I a , Q ′ ( t ) = α s I s + α a I a − ( λ q + δ q + μ ) Q , R ′ ( t ) = λ s I s + λ a I a + λ q Q − μ R , F ′ ( t ) = δ s I s + δ q Q ,

such that α s + δ s + λ s + μ = 1 .

7.1 Comparison between real data and model with rate function

In this subsection, we compare our model with experimental data in Russia since it is well known that mathematical models contribute to understanding and analyzing real-world problems. In order to compare infected cases in Russia with the model, the weekly data have been used from March 2, 2020, to August 29, 2022 [1].

Situation 1. Comparison between real data and the deterministic model

To achieve our aim, we consider the deterministic model

(100) D t α 0 A B C S ( t ) = ρ − β s ( 1 − ψ ξ ) S I s − β a ( 1 − ψ ξ ) S I a − μ S , D t α 0 A B C E ( t ) = β s ( 1 − ψ ξ ) S I s + β a ( 1 − ψ ξ ) S I a − ( φ + μ ) E , D t α 0 A B C I s ( t ) = ( 1 − θ ) φ E − r ( t ) ( α s + δ s + λ s + μ ) I s , D t α 0 A B C I a ( t ) = θ φ E − ( α a + λ a + μ ) I a , D t α 0 A B C Q ( t ) = α s I s + α a I a − ( λ q + δ q + μ ) Q , D t α 0 A B C R ( t ) = λ s I s + λ a I a + λ q Q − μ R , D t α 0 A B C F ( t ) = δ s I s + δ q Q ,

with the initial conditions

(101) S ( 0 ) = 144,000,000 , E ( 0 ) = 130,000,000 , I s ( 0 ) = 5 , I a ( 0 ) = 5 , Q ( 0 ) = 94 , R ( 0 ) = 6 ; F ( 0 ) = 0 .

While performing the simulations, the following parameters are taken as

(102) ρ = 12,500,000 , β s = 2 , ψ = 0.1 , ξ = 0.5 , β a = 0.5 , μ = 3.6529 × 1 0 − 6 , α s = 0.34 , δ s = 0.015 , λ s = 1 50 , α a = 0.3 , δ a = 0.015 , λ a = 1 20 , φ = 1 450 , θ = 0.06 , λ q = 0.0005 , δ q = 0.0015 , N = 144,000,000 .

The comparison between experimental data in Russia and the deterministic model is presented in Figure 17.

Figure 17

Comparison between the deterministic model and accumulative data in Russia.

When Figure 17 is examined, although the deterministic model is successful in predicting the future of the process, it has not been observed that it is successful enough in capturing the cumulative cases exactly.

Situation 2. Comparison between real data and stochastic model

We now consider the following model with stochastic constant:

(103) D t α 0 C S ( t ) = ρ − β s ( 1 − ψ ξ ) S I s − β a ( 1 − ψ ξ ) S I a − μ S + σ 1 S d B 1 ( t ) , D t α 0 C E ( t ) = β s ( 1 − ψ ξ ) S I s + β a ( 1 − ψ ξ ) S I a − ( φ + μ ) E + σ 2 E d B 2 ( t ) , D t α 0 C I s ( t ) = ( 1 − θ ) φ E − r ( t ) ( α s + δ s + λ s + μ ) I s + σ 3 I s d B 3 ( t ) , D t α 0 C I a ( t ) = θ φ E − ( α a + λ a + μ ) I a + σ 4 I a d B 4 ( t ) , D t α 0 C Q ( t ) = α s I s + α a I a − ( λ q + δ q + μ ) Q + σ 5 Q d B 5 ( t ) , D t α 0 C R ( t ) = λ s I s + λ a I a + λ q Q − μ R + σ 6 R d B 6 ( t ) , D t α 0 C F ( t ) = δ s I s + δ q Q + σ 7 F d B 7 ( t ) ,

with the initial conditions

(104) S ( 0 ) = 144,000,000 , E ( 0 ) = 130,000,000 , I s ( 0 ) = 5 , I a ( 0 ) = 5 , Q ( 0 ) = 94 , R ( 0 ) = 6 , F ( 0 ) = 0 .

While performing the simulations, the following parameters are taken as

(105) ρ = 12,500,000 , β s = 2 , ψ = 0.1 , ξ = 0.5 , β a = 0.5 , μ = 3.6529 × 1 0 − 6 , α s = 0.34 , δ s = 0.015 , λ s = 1 50 , α a = 0.3 , δ a = 0.015 , λ a = 1 20 , φ = 1 450 , θ = 0.06 , λ q = 0.0005 , δ q = 0.0015 , N = 144,000,000 .

The comparison between experimental data in Russia and stochastic model is presented in Figure 18.

Figure 18

Comparison between stochastic model and real data in Russia.

When Figure 18 is examined, it is observed that the model with stochastic constant captures more cumulative cases after the 50th week. Although the model with stochastic constant, which is successful in exhibiting non-stationary processes is not able to fully capture the cumulative cases, it is seen that it gives better results than the deterministic model.

Situation 3. Comparison between real data and piecewise model

In this case, we deal with the following model:

(106) Y ′ ( t ) = F ( t , Y ) , if 0 ≤ t ≤ T 1 Y ( 0 ) = Y i , 0 , i = 1 , 2 , … , 7 , D t α T 1 C F Y i = F ( t , Y ) , if T 1 ≤ t ≤ T 2 , Y ( T 1 ) = Y i , 1 , D t α T 2 A B C Y i = F ( t , Y ) , if T 2 ≤ t ≤ T 3 , Y ( T 2 ) = Y i , 2 , D t α T 3 C F Y i = F ( t , Y ) , if T 3 ≤ t ≤ T 4 , Y i ( T 3 ) = Y i , 3 , Y ′ ( t ) = F ( t , Y ) , if T 4 ≤ t ≤ T , Y i ( T 4 ) = Y i , 4 .

with the initial conditions

(107) S ( 0 ) = 144,000,000 , E ( 0 ) = 130,000,000 , I s ( 0 ) = 5 , I a ( 0 ) = 5 , Q ( 0 ) = 94 , R ( 0 ) = 6 , F ( 0 ) = 0 .

While performing the simulations, the following parameters are taken as

(108) ρ = 12,500,000 , β s = 2 , ψ = 0.1 , ξ = 0.5 , β a = 0.5 , μ = 3.6529 × 1 0 − 6 , α s = 0.34 , δ s = 0.015 , λ s = 1 50 , α a = 0.3 , δ a = 0.015 , λ a = 1 20 , φ = 1 450 , θ = 0.06 , λ q = 0.0005 , δ q = 0.0015 , N = 144,000,000 .

The comparison between experimental data in Russia and the deterministic model with piecewise derivative is given in Figure 19.

Figure 19

Comparison between deterministic model and accumulative data in Russia.

Now, let us compare the modified model by adding the stochastic constant to the piecewise model and the real data. Such model is represented by

(109) d Y ( t ) = F ( t , Y ) d t + σ i Y i d B i ( t ) , if 0 ≤ t ≤ T 1 , Y ( 0 ) = Y i , 0 , i = 1 , 2 , … , 7 , D t α T 1 C F Y i = F ( t , Y ) , if T 1 ≤ t ≤ T 2 , Y ( T 1 ) = Y i , 1 , D t α T 2 A B C Y i = F ( t , Y ) + σ i ˜ Y i d B i ( t ) , if T 2 ≤ t ≤ T 3 , Y ( T 2 ) = Y i , 2 , D t α T 3 C F Y i = F ( t , Y ) , if T 3 ≤ t ≤ T 4 , Y i ( T 3 ) = Y i , 3 , Y ′ ( t ) = F ( t , Y ) , if T 4 ≤ t ≤ T , Y i ( T 4 ) = Y i , 4 ,

where α = 0.98 , σ i = 0.1 , and σ i ˜ = 0.0013 . The comparison between experimental data in Russia and stochastic–deterministic model is shown in Figure 20.

Figure 20

Comparison between deterministic-stochastic model and real data in Russia.

In Figure 19, it can be clearly seen that the piecewise model performs quite well for cumulative cases compared to the previous two situations. In order to capture more data points, it is concluded from Figure 20 that the model where stochastic is added in the first and third intervals reach more reality.

Situation 4. Comparison between real data and model with the rate indicator function

Here, we modify the considered model by the rate indicator function, which is given by

(110) S ′ ( t ) = ρ − β s ( 1 − ψ ξ ) S I s − β a ( 1 − ψ ξ ) S I a − μ S , E ′ ( t ) = β s ( 1 − ψ ξ ) S I s + β a ( 1 − ψ ξ ) S I a − ( φ + μ ) E , I s ′ ( t ) = ( 1 − θ ) φ E − r ( t ) ( α s + δ s + λ s + μ ) I s , I a ′ ( t ) = θ φ E − ( α a + λ a + μ ) I a , Q ′ ( t ) = α s I s + α a I a − ( λ q + δ q + μ ) Q , R ′ ( t ) = λ s I s + λ a I a + λ q Q − μ R , F ′ ( t ) = δ s I s + δ q Q ,

where α s + δ s + λ s + μ = 1 . The initial conditions are as follows:

(111) S ( 0 ) = 144,000,000 , E ( 0 ) = 13 , I s ( 0 ) = 2 , I a ( 0 ) = 4 , Q ( 0 ) = 2 , R ( 0 ) = 2 , F ( 0 ) = 1 .

While performing the simulations, the following parameters are considered as

(112) ρ = 0.0000195 , β s = 0.01 , ψ = 0.1 , ξ = 0.5 , β a = 0.5 , μ = 0.3 , α s = 0.3 , δ s = 0.2 , λ s = 0.2 , α a = 0.3 , δ a = 0.0015 , λ a = 1 20 , φ = 1 450 , θ = 0.7 , λ q = 0.0005 , δ q = 0.0015 , N = 144,000,000 .

The comparison between experimental data in Russia and model with the rate indicator function is performed in Figure 21.

Figure 21

Comparison between deterministic model with rate function and real data in Russia.

As shown in Figure 21, the model to which the rate indicator function is added has been quite successful in predicting infected cases in Russia.

8 Conclusion

In this study, a mathematical model depicting a COVID-19 outbreak is modified using different differential operators. The model is considered in four different situations, first deterministic, then stochastic, then piecewise derivative, and finally with the rate indicator function. For these models, various analyses are presented, including basic reproduction number, equilibrium points, existence, and uniqueness. Numerical simulations are depicted for each situation of the model. Finally, the best scenarios have been examined by comparing the weekly and cumulative data from March 2, 2020, to August 29, in Russia, with the deterministic model, stochastic model, piecewise model, and model with the rate indicator function. In this study, our focus was on the fact that the deterministic model is less reliable when dealing with cumulative data for Russia, and thus, the stochastic model is preferred for capturing a greater number of points. Still, when we evaluate real cumulative data, this study reveals that it has been noted that as the number of cases starts to level out, there is a subsequent rise in the overall data, which can be more accurately portrayed using piecewise model. Our analysis showed that the three cases used for predicting cumulative data were insufficient in predicting daily data. To overcome this, we considered the case incorporating the rate indicator function proposed by Atangana and Araz. As a result, we have analyzed the functionality of the rate indicator function, and it has been proven to have a high level of accuracy in predicting daily data. Briefly, according to the obtained results, the model with piecewise derivative gives better results than the deterministic and stochastic model for cumulative cases. However, none of the models were successfully enough for the weekly data. However, upon the examination of the presented simulations, it is evident that the model incorporating the rate indicator function has been the most successful in accurately representing both cumulative and weekly data. It is clearly seen that the rate indicator function will open new doors in mathematical modeling.

Funding information: The authors state no funding involved.
Author contributions: Mehmet Akif Çetin: data curation, writing – original draft, formal analysis, validation, writing – review and editing, investigation. Seda İğret Araz: conceptualization, methodology, writing – review and editing, supervision, formal analysis, validation, investigation, software. All authors have accepted responsibility for the entire content of this manuscript and approved its submission.
Conflict of interest: The authors state no conflict of interest.

References

[1] World Health Organization, WHO Coronavirus Disease (COVID-19) Dashboard. 2020. https://covid19.who.int/table. Search in Google Scholar

[2] Ozair M, Hussain T, Hussain M, Awan AU, Baleanu D, Abro KA. A mathematical and statistical estimation of potential transmission and severity of COVID-19: a combined study of Romania and Pakistan. BioMed Res Int. 2020;2020:5607236 10.1101/2020.05.02.20088989Search in Google Scholar

[3] Hussain T, Awan AU, Abro KA, Ozair M, Manzoor M. A mathematical and parametric study of epidemiological smoking model: a deterministic stability and optimality for solutions. Europ Phys J Plus. 2021;136:1–23. 10.1140/epjp/s13360-020-00979-4Search in Google Scholar

[4] Akbulut Arik I, Sari H, İIğret Araz S. Numerical simulation of COVID-19 model with integer and non-integer order: the effect of environment and social distancing. Results Phys. 2023;51. 10.1016/j.rinp.2023.106725Search in Google Scholar

[5] Atangana A, İğret Araz S. A novel COVID-19 model with fractional differential operators with singular and non-singular kernels: analysis and numerical scheme based on Newton polynomial. Alexandria Eng J. 2021;60(4):3781–806. 10.1016/j.aej.2021.02.016Search in Google Scholar

[6] Atangana A, Araz SI. Advanced analysis in epidemiological modeling: detection of waves. AIMS Math. 2022;7(10):18010–30. 10.3934/math.2022992Search in Google Scholar

[7] Caputo M, Linear model of dissipation whose Q is almost frequency independent. II. Geophys J Int. 1967;13(5):529–39. 10.1111/j.1365-246X.1967.tb02303.xSearch in Google Scholar

[8] Caputo M, Fabrizio M. On the notion of fractional derivative and applications to the hysteresis phenomena. Mechanica. 2017;52(13):3043–52. 10.1007/s11012-017-0652-ySearch in Google Scholar

[9] Atangana A, Baleanu D. New fractional derivatives with non-local and non-singular kernel. Theory Appl Heat Transfer Model Thermal Sci. 2016;20(2):763–9. 10.2298/TSCI160111018ASearch in Google Scholar

[10] Akbulut Arik I, Igret Araz S. Crossover behaviors via piecewise concept: A model of tumor growth and its response to radiotherapy. Results Phys. 2022;34. 10.1016/j.rinp.2022.105894Search in Google Scholar

[11] Atangana A, Igret Araz S. New concept in calculus: piecewise differential and integral operators. Chaos Solitons Fractals. 2021;145. 10.1016/j.chaos.2020.110638. Search in Google Scholar

[12] Alderremy AA, Gómez-Aguilar JF, Aly S, Saad KM. A fuzzy fractional model of coronavirus (COVID-19) and its study with Legendre spectral method. Results Phys. 2021;21:103773. 10.1016/j.rinp.2020.103773Search in Google Scholar PubMed PubMed Central

[13] Masandawa L, Mirau SS, Mbalawata IS. Mathematical modeling of COVID-19 transmission dynamics between healthcare workers and community. Results Phys. 2021;29:104731. 10.1016/j.rinp.2021.104731Search in Google Scholar PubMed PubMed Central

[14] Ariffin MRK, Gopal K, Krishnarajah I, Ilias ISC, Adam MB, Arasan J, et al. Mathematical epidemiologic and simulation modelling of first wave COVID-19 in Malaysia. Scientific Reports. 2021;11:20739. 10.1038/s41598-021-99541-0Search in Google Scholar PubMed PubMed Central

[15] Premarathna IHK, Srivastava HM, Juman ZAMS, AlArjani A, Uddin MS, Sana SS. Mathematical modeling approach to predict COVID-19 infected people in Sri Lanka. AIMS Mathematics. 2021;7(3):4672–99. 10.3934/math.2022260Search in Google Scholar

[16] Lotfi R, Kheiri K, Sadeghi A, Tirkolaee EB. An extended robust matematical model to project to course of COVID-19 epidemic in Iran. Ann Operat Res. 2022;2022. 10.1007/s10479-021-04490-6. Search in Google Scholar PubMed PubMed Central

[17] Asempapa R, Oduro B, Apenteng OO, Magagula VM. A COVID-19 mathematical model of at-risk populations with non-pharmaceutical preventive measures: the case of Brazil and South Africa. Infect Disease Model. 2022;7(1):45–61. 10.1016/j.idm.2021.11.005Search in Google Scholar PubMed PubMed Central

[18] Samui P, Mondal J, Khajanchi S. A mathematical model for COVID-19 transmission dynamics with a case study of India. Chaos Solitons Fractals. 2020;140:110173. 10.1016/j.chaos.2020.110173Search in Google Scholar PubMed PubMed Central

[19] Alqarni MS, Alghamdi M, Muhammad T, Alshomrani AS, Khan MA. Mathematical modeling for novel coronavirus (COVID-19) and control. Numer Methods Partial Differ Equ. 2020;38(4):760–76. 10.1002/num.22695Search in Google Scholar PubMed PubMed Central

[20] Grassly NC, Pons-Salort M, Parker EPK, White PJ, Ferguson NM. Comparison of molecular testing strategies for COVID-19 control: a mathematical modelling study. Lancet Infectious Diseases. 2020;20(12):1381–9. 10.1016/S1473-3099(20)30630-7Search in Google Scholar PubMed PubMed Central

[21] Riyapan P, Shuaib SE, Intarasit A. A mathematical model of COVID-19 pandemic: a case study of Bangkok, Thailand. Comput Math Methods Med. 2021;2021:6664483. 10.1155/2021/6664483Search in Google Scholar PubMed PubMed Central

[22] Driessche P, Watmough J. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. Math Biosci. 2002;180(1):29–48. 10.1016/S0025-5564(02)00108-6Search in Google Scholar

[23] Atangana A, Igret Araz S. Step forward in epidemiological modeling: Introducing the indicator function to capture waves. Results Phys. 2022;38:105638. 10.1016/j.rinp.2022.105638Search in Google Scholar

Received: 2023-10-31

Revised: 2024-02-16

Accepted: 2024-03-13

Published Online: 2024-04-29

This work is licensed under the Creative Commons Attribution 4.0 International License.

Articles in the same Issue

https://doi.org/10.1515/phys-2024-0009

Keywords for this article

COVID-19 spread; fractional differentiation and integration; rate indicator function

Creative Commons

BY 4.0