Next Article in Journal
Effects of Relative Density and Grading on the Particle Breakage and Fractal Dimension of Granular Materials
Next Article in Special Issue
A New Incommensurate Fractional-Order Discrete COVID-19 Model with Vaccinated Individuals Compartment
Previous Article in Journal
Particle Swarm Optimization Fractional Slope Entropy: A New Time Series Complexity Indicator for Bearing Fault Diagnosis
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

On the Modeling of COVID-19 Transmission Dynamics with Two Strains: Insight through Caputo Fractional Derivative

1
Department of Mathematics, Faculty of Science and Technology, Universitas Airlangga, Surabaya 60115, Indonesia
2
Department of Mathematics and Applied Mathematics, University of Johannesburg, Auckland Park 2006, South Africa
3
Division of Infectious Diseases and Global Public Health, University of California, San Diego, CA 92093, USA
*
Author to whom correspondence should be addressed.
Fractal Fract. 2022, 6(7), 346; https://doi.org/10.3390/fractalfract6070346
Submission received: 30 May 2022 / Revised: 14 June 2022 / Accepted: 15 June 2022 / Published: 21 June 2022

Abstract

:
The infection dynamics of COVID-19 is difficult to contain due to the mutation nature of the SARS-CoV-2 virus. This has been a public health concern globally with the impact of the pandemic on the world’s economy and mode of living. In the present work, we formulate and examine a fractional model of COVID-19 considering the two variants of concern on the disease transmission pathways, namely SARS-CoV-2 and D614G on our model formulation. The existence and uniqueness of our model solutions were analyzed using the fixed point theory. Mathematical analyses were presented, and the model’s basic reproduction numbers R 01 and R 02 were determined. The model has three equilibria: the disease-free equilibrium, that endemic for strain 1, and that endemic for strain 2. The locally asymptotic stability of the equilibria was established based on the R 01 and R 02 values. Caputo fractional operator was used to simulate the model to study the dynamics of the model solution. Results from numerical simulations envisaged that an increase in the transmission parameters of strain 1 leads to an increase in the number of infected individuals. On the other hand, an increase in the strain 2 transmission rate gives rise to more infection. Furthermore, it was established that there is an increased number of infections with a negative impact of strain 1 on strain 2 dynamics and vice versa.

1. Introduction

Coronavirus disease (COVID-19) is an infectious illness caused by SARS-CoV-2. The first case was announced in Wuhan, China on 31 December 2019 [1]. COVID-19 transmission can occur from human to human through splashes from the nose or mouth that come out through coughing, sneezing, talking, singing, or breathing heavily. In addition, contagion can also occur when a healthy person touches the surface of droplet-contaminated objects around an infected person [1]. The COVID-19 disease has a negative effect on people’s lives, including psychological, economic, and social effects. The COVID-19 pandemic also increases the risk of suicide due to lockdown, social distancing, and financial crisis. Measures that can minimize the spread of COVID-19 include educating the public about the efforts to prevent the spread of COVID-19 [2]. Other preventive measures include the use of masks when around others, washing hands as often as possible, keeping at least 6 feet away from others, staying home and isolating yourself from others when sick, avoiding crowds, avoiding poorly ventilated rooms, covering your mouth with tissues when coughing and sneezing, and routinely cleaning surfaces that are frequently touched [3]. According to the World Health Organisation (WHO), there were 85,091,012 cases of COVID-19 with 1,861,005 deaths worldwide as of 6 January 2021 [4].
Before the COVID-19 outbreak, there were six different coronaviruses species that could infect humans, namely HCoV-229E, HCoV-NL63, HCoV-OC43, HCoV-HKU1, SARS-CoV, and MERS-CoV [5]. Common symptoms of COVID-19 infection incorporate acute respiratory disorders such as fever, cough, and shortness of breath [6]. SARS-CoV-2 is a virus that can mutate and produce new strains. One mutation of SARS-CoV-2 is D614G. The D614G strains is an amino acid mutation in position 614 that results in a change from D (aspartic acid) to G (glycine) [7]. The D614G strains is the most dominant variant in the world because it spreads faster than other types and tends to be more contagious [8]. The SARS-CoV-2 variants have been classified by CDC and the U.S government SARS-CoV-2 Interagency Group (SIG) into four categories [9]: (i) Variant being monitored, which consists of Alpha, Beta, Gamma, Delta, Epsilon, Eta, Iota, Kappa, Mu, and Zeta (P.2), (ii) Variant of interest, (iii) Variant of concern: Omicron and (iv) Variant of high consequence. The Omicron variant (B.1.1.529) was declared a variant of concern by the WHO on 26 November 2021 [10]. Since the detection of the Omicron variant in South Africa via genome sequencing [11], it has also spread rapidly across the globe with an increase in the new number of daily COVID-19 cases, which is linked to the Omicron variant [4]. In particular, the Omicron variant has about 60 mutations from the original wild-type variant that emanated originally from Wuhan, China, in 2020, with more than 30 in the spike protein [12]. These spike mutations of Omicron have prominent significant concerns about current vaccines’ potential effectiveness and the level of immunity conferred from previous infections of COVID-19. In addition, it also has increased transmissibility, hospitalization rate, mortality, disease severity, and decreased vaccine-induced protection from severe disease as well as neutralization by antibodies generated during previous infection or after vaccination [9].
Since COVID-19 was declared by the WHO as a global pandemic, a number of mathematical models have been elaborated to conceive the spread of COVID-19. Some of the research work done in the field of mathematical modeling related to the COVID-19 problem applied both integer [13,14,15,16,17,18,19,20] and fractional order in their modeling frameworks [21,22,23,24,25,26,27]. Fractional calculus, which uses fractional derivatives (FD), is a new rapidly growing area in mathematical science that picks up a memory effect and hereditary characteristics of various physical and natural phenomena occurring in engineering, technology, physical and natural sciences [28,29,30]. The most widely used FD includes the Caputo, Caputo–Fabrizio (CF) and Atangana–Baleanu–Caputo fractional-order operators. Other researchers have also analyzed COVID-19 by exploring possible control measures to mitigate the spread of the virus through the application of fractional derivative theory. The authors in [24] proposed a fractional-order model which considers the effect of the pathogens in the environment and vaccination on the transmission dynamics of COVID-19 in Indonesia. Their finding suggests that a decrease in the pathogen level will reduce the number of infected individuals in the human population, which helps reduce the spread of the disease.
Recent work on COVID-19 modeling which incorporates two strains includes models with vaccination [31], optimal control [32], vaccination scale-up [33], with pre-symptomatic individuals [34] and non-pharmaceutical interventions [35] in the form integer model. It is essential to explore the COVID-19 transmission caused by two different types of variants using a mathematical framework. Limited evidence indicates that those who recover from the disease can be reinfected again when there is a new strain but at a rate less severe than the previous infection. Therefore, we investigate the potential impact of two variants on the COVID-19 transmission pathway by considering a case scenario where it was caused by two different strains of the SARS-CoV-2 variants. We are now motivated in the present study to modify the mathematical model of the spread of COVID-19 that was developed by Rezapour et al. [21] by including a transmission path with two strains. The model initially did not consider strains of the COVID virus in the modeling approach. Therefore, to further refine the work in [21], we considered two strains of the virus, namely SARS-CoV-2 and D614G, in our modeling approach by applying the Caputo operator. The fractional model leads to more interesting and unique findings than the integer model due to providing a better understanding of disease dynamics. For example, the authors in [24] indicated that the Caputo model gives a preferable approximation to the reported infected cases.
Furthermore, the model modifications made include dividing the population of exposed individuals into two, based on the strain of the virus, that is, the exposed population to strain 1 and the exposed population to strain 2. In addition, the population of infected individuals is divided into those infected with strain 1 and those infected with strain 2. The CF basic mathematical preliminaries is discussed in the next section.

2. Basic Caputo Fractional-Order Preliminaries

Here, we provide the basic definitions for fractional calculus that will be employed to consider the COVID-19 model with two strains.
Definition 1.
According to [36], let χ C n be a function; then, the Caputo derivative with a fractional order σ ( n 1 , n ) for n N is defined to be
C D t σ ( χ ( t ) ) = 1 Γ ( n σ ) 0 t d n d t n χ ( φ ) ( t φ ) n σ 1 d φ .
Thus, it is obviously that C D t σ ( χ ( t ) ) χ ( t ) if σ 1 .
Definition 2.
According to Definition 1, the integral operator for the function χ : R B with order σ 0 is given by
C I t σ ( χ ( t ) ) = 1 Γ ( σ ) 0 t χ ( φ ) ( t φ ) ( σ 1 ) d φ .
for 0 < σ < 1 , t > 0 .
Lastly, we give the definition which denotes the equilibrium for Caputo fractional model. Hence, it is established as follows
Definition 3.
According to [37], if χ * is the equilibrium point, then
C D t σ ( χ ( t ) ) = f ( t , χ ( t ) )
for σ ( 0 , 1 ) if and only if f ( t , χ * ) = 0
The remaining of the paper is arranged as follows. Section 3 explains details of the model formulation followed by model analyses in Section 4. Meanwhile, Section 5 presents the numerical simulations and Section 6 summarizes the main finding of the paper.

3. The Model

In this section, we describe a mathematical model of the spread of COVID-19 with two strains. The model consists of six subpopulations, namely, the susceptible class ( S ) , the class of individuals exposed to strain 1 ( E 1 ), the class of individuals exposed to strains 2 ( E 2 ), infected with strains 1 ( I 1 ), infected with strain 2 ( I 2 ), and recovered class from either strain 1 or 2 ( R ) . Thus, the total population at any time t is given by N ( t ) = S ( t ) + E 1 ( t ) + E 2 ( t ) + I 1 ( t ) + I 2 ( t ) + R ( t ) . The explicit notation and description of each variable and parameter used in this mathematical model are given in Table 1.
Our model is built based on the following assumptions in terms of COVID-19 with two strains; (1) both the exposed and infected with strains 1 and 2 populations can transmit the disease; (2) the transmission rate of the infected individual with strain 2 is greater than that of the infected individuals with strain 1; (3) the transmission rate of individuals exposed to strain 2 is greater than that of individuals exposed to strain 1; (4) only the infected individuals die from the disease; (5) the recovered individuals cannot be reinfected.
Based on the notation and assumptions, a transmission diagram of the mathematical model of the spread of COVID-19 with two strains can be constructed as shown in Figure 1.
Combining the transmission diagram shown in Figure 1 and the model assumptions, we have the following ordinary differential equation (ODE) governing our mathematical model for the spread of COVID-19 with two strains, which we shall explain in detail. Equation (1) represents the rate of change of the susceptible individual population to the unit of time. Susceptible populations can increase due to the recruitment rate and decrease due to interactions between susceptible and exposed individuals to strain 1 and individuals infected with strain 1. In addition, susceptible populations can also decrease due to interactions between the susceptible and exposed individuals to strain 2 and individuals infected with strain 2. The population can also be reduced due to the natural death rate.
d S d t = Λ α 1 E 1 + β 1 I 1 S α 2 E 2 + β 2 I 2 S μ S .
Equation (2) represents the rate of change of the population of individuals exposed to strain 1. This population increases because of the interaction between susceptible individuals and individuals exposed to strain 1 and infected with strain 1. This population can decrease due to the rate of development from individuals exposed to strain 1 to individuals infected with strain 1. The population can also decrease due to the natural death rate.
d E 1 d t = α 1 E 1 + β 1 I 1 S p E 1
where p = ( μ + δ ) .
Equation (3) represents the rate of change of the population of individuals infected with strain 1. This population increases because of the rate of transition of individuals exposed to strain 1 to individuals infected with strain 1. The population can decrease due to the natural death rate and death rate due to strain 1. In addition, the population can also decrease due to the recovery rate of individuals infected with strain 1.
d I 1 d t = δ E 1 q I 1
for q = ( μ + σ 1 + θ 1 ) .
Equation (4) represents the rate of change of the population of individuals exposed to strain 2. This population increases because of the interaction between susceptible individuals and individuals exposed to strain 2 and individuals infected with strain 2. This population can decrease due to the rate of development from individuals exposed to strain 2 to individuals infected with strain 2. It can also decrease due to the natural death rate.
d E 2 d t = α 2 E 2 + β 2 I 2 S r E 2
in which r = ( μ + E ) .
Equation (5) represents the rate of change of the population of individuals infected with strain 2 with respect to time. This population increases because of the rate of development of individuals exposed to strain 2 to individuals infected with strain 2. The population can decrease due to the natural death rate and death rate due to strain 2. In addition, the population can also decrease due to the recovery rate of individuals infected with strain 2.
d I 2 d t = E E 2 w I 2
for w = ( μ + σ 2 + θ 2 ) .
Equation (6) represents the rate of change of population of individuals recovered from both strains. This population increases because of the recovery rate of individuals infected with strain 1 and the rate of recovery of individuals infected with strain 2. The population can decrease due to the natural death rate.
d R d t = θ 1 I 1 + θ 2 I 2 μ R .
Therefore, combining Equations (1)–(6), we have the general system of ODE describing COVID-19 transmission pathways with two trains:
d S d t = Λ α 1 E 1 + β 1 I 1 S α 2 E 2 + β 2 I 2 S μ S , d E 1 d t = α 1 E 1 + β 1 I 1 S p E 1 , d I 1 d t = δ E 1 q I 1 , d E 2 d t = α 2 E 2 + β 2 I 2 S r E 2 , d I 2 d t = E E 2 w I 2 , d R d t = θ 1 I 1 + θ 2 I 2 μ R ,
subject to the initial conditions
S ( 0 ) > 0 , E 1 ( 0 ) 0 , I 1 ( 0 ) 0 , E 2 ( 0 ) 0 , I 2 ( 0 ) 0 , R ( 0 ) 0 .
We transform the model Equation (7) into a fractional differential model, which yields
C D t σ S ( t ) = 1 Γ ( 1 σ ) 0 t ( t φ ) σ S ( φ ) d φ , C D t σ E 1 ( t ) = 1 Γ ( 1 σ ) 0 t ( t φ ) σ E 1 ( φ ) d φ , C D t σ I 1 ( t ) = 1 Γ ( 1 σ ) 0 t ( t φ ) σ I 1 ( φ ) d φ , C D t σ E 2 ( t ) = 1 Γ ( 1 σ ) 0 t ( t φ ) σ E 2 ( φ ) d φ , C D t σ I 2 ( t ) = 1 Γ ( 1 σ ) 0 t ( t φ ) σ I 2 ( φ ) d φ , C D t σ R ( t ) = 1 Γ ( 1 σ ) 0 t ( t φ ) σ R ( φ ) d φ .
Using the Caputo operator, we covert system (7) to a fractional ODE system:
0 C D t σ S ( t ) = Λ α 1 E 1 ( t ) + β 1 I 1 ( t ) S ( t ) α 2 E 2 + β 2 I 2 ( t ) S ( t ) μ S ( t ) , 0 C D t σ E 1 ( t ) = α 1 E 1 ( t ) + β 1 I 1 ( t ) S ( t ) p E 1 ( t ) , 0 C D t σ I 1 ( t ) = δ E 1 ( t ) q I 1 ( t ) , 0 C D t σ E 2 ( t ) = α 2 E 2 ( t ) + β 2 I 2 ( t ) S ( t ) r E 2 ( t ) , 0 C D t σ I 2 ( t ) = E E 2 ( t ) w I 2 ( t ) , 0 C D t σ R ( t ) = θ 1 I 1 ( t ) + θ 2 I 2 ( t ) μ R ( t ) ,
subjective to initial conditions
0 C D t σ S ( 0 ) > 0 , 0 C D t σ E 1 ( 0 ) 0 , 0 C D t σ I 1 ( 0 ) 0 , 0 C D t σ E 2 ( 0 ) 0 , 0 C D t σ I 2 ( 0 ) 0 , 0 C D t σ R ( 0 ) 0 .
The parameters of the model and state variables are assumed to be positive over the modeling time horizon.

4. Mathematical Analyses

In this section, we establish some basic properties of our model (9).

4.1. Basic Properties of the Fractional Model

We verify that the fractional model (9) is positive, bounded, and biologically meaningful over the modeling time for t 0 .
Proof. 
The total dynamic of the fractional model (9) is given by the following equation
0 C D t τ N ( t ) = 0 C D t τ S ( t ) + 0 C D t τ E 1 ( t ) + 0 C D t τ I 1 ( t ) + 0 C D t τ E 2 ( t ) + 0 C D t τ I 2 ( t ) + 0 C D t τ R ( t ) , = Λ μ N ( q + w ) , Λ μ N .
The solutions of (10) yield
N ( t ) N ( 0 ) E σ , 1 ( μ t σ ) + Λ t σ E σ , σ + 1 ( μ t σ ) ,
where E σ , σ + 1 is the Mittag–Leffler function defined by E σ , σ + 1 ( t ) = Σ k = 0 t k Γ ( σ k + ( σ + 1 ) ) .
We obtain N ( t ) Λ μ when t .
The biological feasible region Γ ¯ of the COVID-19 model with two strains is expressed by
Γ ¯ ¯ = S ( t ) , E 1 ( t ) , I 1 ( t ) , E 2 ( t ) , I 2 ( t ) , R ( t ) R + 6 Λ μ .
Hence, the fractional model (9) which incorporates two strains is biologically feasible and bounded for all time t 0 .  □

4.2. Existence and Uniqueness of Solutions

Now, we examine and prove the existence and uniqueness of the solutions of the model (9) with the initial condition using the fixed point theory.
Let G ( Q ) be the Banach space of the real valued continuous function over Q = [ 0 , c ] with the norm given by
| | S | | = sup t Q | S ( t ) | , | | E 1 | | = sup t Q | E 1 ( t ) | , | | I 1 | | = sup t Q | I 1 ( t ) | , | | E 2 | | = sup t Q | E 2 ( t ) | ,
| | I 2 | | = sup t Q | I 2 ( t ) | and | | R | | = sup t Q | R ( t ) | .
Further, we also have the norm
| | S ( t ) , E 1 ( t ) , I 1 ( t ) , E 2 ( t ) , I 2 ( t ) , R ( t ) | | = | | S ( t ) | | + | | E 1 ( t ) | | + | | I 1 ( t ) | | + | | E 2 ( t ) | | +
| | I 2 ( t ) | | + | | R ( t ) | | .
Proof. 
We transform the model (9) using the Caputo integral
S ( t ) = S ( 0 ) + C D t σ Λ α 1 E 1 + β 1 I 1 S α 2 E 2 + β 2 I 2 S μ S , E 1 ( t ) = E 1 ( 0 ) + C D t σ α 1 E 1 + β 1 I 1 S p E 1 , I 1 ( t ) = I 1 ( 0 ) + C D t σ δ E 1 q I 1 , E 2 ( t ) = E 2 ( 0 ) + C D t σ α 2 E 2 + β 2 I 2 S r E 2 , I 2 ( t ) = I 2 ( 0 ) + C D t σ E E 2 w I 2 R ( t ) = R ( 0 ) + C D t σ θ 1 I 1 + θ 2 I 2 μ R
Utilizing Definition (2), we obtain
S ( t ) = S ( 0 ) + 1 Γ ( σ ) 0 t ( t φ ) σ 1 F 1 ( φ , S ( φ ) ) d φ , E 1 ( t ) = E 1 ( 0 ) + 1 Γ ( σ ) 0 t ( t φ ) σ 1 F 2 ( φ , E 1 ( φ ) ) d φ , I 1 ( t ) = I 1 ( 0 ) + 1 Γ ( σ ) 0 t ( t φ ) σ 1 F 3 ( φ , I 1 ( φ ) ) d φ , E 2 ( t ) = E 2 ( 0 ) + 1 Γ ( σ ) 0 t ( t φ ) σ 1 F 4 ( φ , E 2 ( φ ) ) d φ , I 2 ( t ) = I 2 ( 0 ) + 1 Γ ( σ ) 0 t ( t φ ) σ 1 F 5 ( φ , I 2 ( φ ) ) d φ , R ( t ) = R ( 0 ) + 1 Γ ( σ ) 0 t ( t φ ) σ 1 F 6 ( φ , R ( φ ) ) d φ ,
with the following kernels:
F 1 ( t , S ( t ) ) = Λ α 1 E 1 + β 1 I 1 S α 2 E 2 + β 2 I 2 S μ S , F 2 ( t , E 1 ( t ) ) = α 1 E 1 + β 1 I 1 S p E 1 , F 3 ( t , I 1 ( t ) ) = δ E 1 q I 1 , F 4 ( t , E 2 ( t ) ) = α 2 E 2 + β 2 I 2 S r E 2 , F 5 ( t , I 2 ( t ) ) = E E 2 w I 2 F 5 ( t , R ( t ) ) = θ 1 I 1 + θ 2 I 2 μ R .
The expression F j in Equation (12) for j = 1 , , 6 fulfills the Lipschitz for the model variables ( S ( t ) , E 1 ( t ) , I 1 ( t ) , E 2 ( t ) , I 2 ( t ) , R ( t ) ) that hold an upper bound. Firstly, we provide the proof for the S variable. Let S and S ˜ be two different functions. We examine the proof for S ( t ) .
| | F 1 ( t , S ( t ) ) F 1 ( t , S ˜ ( t ) ) | | = | | α 1 E 1 + β 1 I 1 α 2 E 2 + β 2 I 2 μ ( S ( t ) S ˜ ( t ) ) | | .
Setting ζ 1 = ( | | α 1 E 1 + β 1 I 1 + α 2 E 2 + β 2 I 2 + μ | | ) , we get
| | F 1 ( t , S ( t ) ) F 1 ( t , S ¯ ( t ) ) | | ζ 1 | | S ( t ) S ¯ ( t ) | | .
Likewise, for the continuing of state variables, we carry that
| | F 2 ( t , E 1 ( t ) ) F 2 ( t , E 1 ˜ ( t ) ) | | ζ 2 | | E 1 ( t ) E 1 ˜ ( t ) | | , | | F 3 ( t , I 1 ( t ) ) F 3 ( t , I 1 ˜ ( t ) ) | | ζ 3 | | I 1 ( t ) I 1 ˜ ( t ) | | , | | F 4 ( t , E 2 ( t ) ) F 4 ( t , E 2 ˜ ( t ) ) | | ζ 4 | | E 2 ( t ) E 2 ˜ ( t ) | | , | | F 5 ( t , I 2 ( t ) ) F 5 ( t , I 2 ˜ ( t ) ) | | ζ 5 | | I 2 ( t ) I 2 ˜ ( t ) | | , | | F 6 ( t , R ( t ) ) F 6 ( t , R ˜ ( t ) ) | | ζ 6 | | R ( t ) R ˜ ( t ) | | ,
where ζ 1 , ζ 2 , ζ 3 , ζ 4 , ζ 5 and ζ 6 assert the Lipschitz constants for each kernel F j for j = 1 , , 6 respectively. Hence, the Lipschitz condition is fulfilled. To finish the uniqueness proof by employing the Banach theory, we also represent that system (11) is recursive in nature. Therefore
S n ( t ) = S ( 0 ) + 1 Γ ( σ ) 0 t ( t φ ) σ 1 F 1 ( φ , S n 1 ( φ ) ) d φ , E 1 n ( t ) = E 18 ( 0 ) + 1 Γ ( σ ) 0 t ( t φ ) σ 1 F 2 ( φ , E 1 ( n 1 ) ( φ ) ) d φ , I 1 n ( t ) = I 1 ( 0 ) + 1 Γ ( σ ) 0 t ( t φ ) σ 1 F 3 ( φ , I 1 ( n 1 ) ( φ ) ) d φ , E 2 n ( t ) = E 2 ( 0 ) + 1 Γ ( σ ) 0 t ( t φ ) σ 1 F 4 ( φ , E 2 ( n 1 ) ( φ ) ) d φ , I 2 n ( t ) = I 2 ( 0 ) + 1 Γ ( σ ) 0 t ( t φ ) σ 1 F 5 ( φ , I 2 ( n 1 ) ( φ ) ) d φ , R n ( t ) = R ( 0 ) + 1 Γ ( σ ) 0 t ( t φ ) σ 1 F 6 ( φ , R n 1 ( φ ) ) d φ .
After some algebraic calculations of the difference between the successive expression together with the initial conditions, we obtain
Φ S , n ( t ) = S n ( t ) S n 1 ( t ) = 1 Γ ( σ ) 0 t ( t φ ) σ 1 F 1 ( φ , S n 1 ( φ ) ) F 1 ( φ , S n 2 ( φ ) ) d φ , Φ E 1 , n ( t ) = E 1 n ( t ) E 1 n 1 ( t ) = 1 Γ ( σ ) 0 t ( t φ ) σ 1 F 2 ( φ , V n 1 ( φ ) ) F 2 ( φ , V n 2 ( φ ) ) d φ , Φ I 1 , n ( t ) = I 1 n ( t ) I 1 n 1 ( t ) = 1 Γ ( σ ) 0 t ( t φ ) σ 1 F 3 ( φ , E n 1 ( φ ) ) F 3 ( φ , E n 2 ( φ ) ) d φ , Φ E 2 , n ( t ) = E 2 n ( t ) E 2 n 1 ( t ) = 1 Γ ( σ ) 0 t ( t φ ) σ 1 F 4 ( φ , I n 1 ( φ ) ) F 4 ( φ , E 2 n 2 ( φ ) ) d φ , Φ I 2 , n ( t ) = I 2 n ( t ) I 2 n 1 ( t ) = 1 Γ ( σ ) 0 t ( t φ ) σ 1 F 5 ( φ , I 2 n 1 ( φ ) ) F 5 ( φ , I 2 n 2 ( φ ) ) d φ , Φ R , n ( t ) = R n ( t ) R n 1 ( t ) = 1 Γ ( σ ) 0 t ( t φ ) σ 1 F 6 ( φ , R n 1 ( φ ) ) F 6 ( φ , R n 2 ( φ ) ) d φ ,
where in (14)
S n ( t ) = j = 0 n Φ S , j ( t ) , E 1 n ( t ) = j = 0 n Φ E 1 , j ( t ) , I 1 n ( t ) = j = 0 n Φ I 1 , j ( t ) E 2 n ( t ) = j = 0 n Φ E 2 , j ( t ) , I 2 n ( t ) = j = 0 n Φ I 2 , j ( t ) , R n ( t ) = j = 0 n Φ R , j ( t )
Furthermore, let
Φ ( S , ( n j ) ) ( t ) = S n 1 ( t ) S n 2 ( t ) , Φ ( E 1 , ( n j ) ) ( t ) = E 1 n 1 ( t ) E 1 n 2 ( t ) , Φ ( I 1 , ( n j ) ) ( t ) = I 1 n 1 ( t ) I 1 n 2 ( t ) , Φ ( E 2 , ( n j ) ) ( t ) = E 2 n 1 ( t ) E 2 n 2 ( t ) , Φ ( I 2 , ( n j ) ) ( t ) = I 2 n 1 ( t ) I 2 n 2 ( t ) , Φ ( R , ( n j ) ) ( t ) = R n 1 ( t ) R n 2 ( t ) ,
then we have that
| | Φ S , n ( t ) | | ζ 1 Γ ( σ ) 0 t ( t φ ) σ 1 | | Φ S , n 1 ( φ ) | | d φ , | | Φ E 1 , n ( t ) | | ζ 2 Γ ( σ ) 0 t ( t φ ) σ 1 | | Φ E 1 , n 1 ( φ ) | | d φ , | | Φ I 1 , n ( t ) | | ζ 3 Γ ( σ ) 0 t ( t φ ) σ 1 | | Φ E 2 , n 1 ( φ ) | | d φ , | | Φ E 2 , n ( t ) | | ζ 4 Γ ( σ ) 0 t ( t φ ) σ 1 | | Φ I 2 , n 1 ( φ ) | | d φ , | | Φ I 2 , n ( t ) | | ζ 5 Γ ( σ ) 0 t ( t φ ) σ 1 | | Φ I 2 , n 1 ( φ ) | | d φ , | | Φ R , n ( t ) | | ζ 6 Γ ( σ ) 0 t ( t φ ) σ 1 | | Φ R , n 1 ( φ ) | | d φ .
Since the conditions for a Lipschitz function are all satisfied by S ( t ) , E 1 ( t ) , I 1 ( t ) , E 2 ( t ) , I 2 ( t ) , R ( t ) and the functions F j for j = 1 , , 6 , we deduce that it is bounded. □

4.3. Model Equilibrium and the Basic Reproduction Number

By utilizing Definition 3, the fractional COVID-19 model with two strains has three equilibria: namely, the disease-free equilibrium ( E 0 ), the endemic for strain 1 ( E 1 ), and the endemic for strain 2 ( E 2 ). The disease-free equilibrium (DFE) of the model is
E 0 = ( S , E 1 , I 1 , E 2 , I 2 , R ) = Λ μ , 0 , 0 , 0 , 0 , 0 .
Furthermore, using the Next-Generation Matrix method as in [38], the COVID-19 model with two strains has the basic reproduction number denoted by ( R 0 ) as
R 0 = max { R 01 , R 02 } ,
where
R 01 = Λ α 1 q + β 1 Λ δ μ p q and R 02 = Λ α 2 w + β 2 Λ E μ r w .
Importantly, in the above expression for our model R 0 , we have that R 01 is the contribution from strain 1, while R 02 is the contribution from stain 2 of the virus. Hence, R 0 is responsible as a key transmissibility parameter which is a function of transmission rates and other model parameters and is required to cause an average infectious in a susceptible population.
Furthermore, two endemic equilibria of the COVID-19 model are calculated and found to be E 1 = S * , E 1 * , I 1 * , 0 , 0 , R * and E 2 = S * * , 0 , 0 , E 2 * * , I 2 * * , R * * where
S * = p q q α 1 + β 1 δ , E 1 * = μ q R 01 1 ( α 1 q + β 1 δ ) , I 1 * = δ μ R 01 1 ( α 1 q + β 1 δ ) , R * = θ 1 δ R 01 1 ( α 1 q + β 1 δ )
and
S * * = r w w α 2 + β 2 E , E 2 * * = μ w R 02 1 ( α 2 w + β 2 E ) , I 2 * * = E μ R 02 1 ( α 2 w + β 2 E ) , R * * = θ 2 E R 02 1 ( α 2 w + β 2 E ) .
Thus, the existence conditions for the endemic equilibrium point E 1 and E 2 meet the conditions R 01 > 1 and R 02 > 1 , respectively.

4.3.1. Stability Analysis of the Equilibrium Point

The stability of the COVID-19 model with two strains can be analyzed as follows. The Jacobian matrix of the model is obtained as
J = Q 0 α 1 S * β 1 S * α 2 S * β 2 S * 0 Q 1 α 1 S * p β 1 S * 0 0 0 0 δ q 0 0 0 Q 2 0 0 α 2 S * r β 2 S * 0 0 0 0 E w 0 0 0 θ 1 0 θ 2 μ
where Q 0 = α 1 E 1 * + β 1 I 1 * α 2 E 2 * + β 2 I 2 * μ , Q 1 = α 1 E 1 * + β 1 I 1 * , Q 2 = α 2 E 2 * + β 2 I 2 * .

4.3.2. Local Stability of DFE

The Jacobian matrix of DFE E 0 = ( S , E 1 , I 1 , E 2 , I 2 , R ) = Λ μ , 0 , 0 , 0 , 0 , 0 is given by
J ( E 0 ) = μ α 1 Λ μ β 1 Λ μ α 2 Λ μ β 2 Λ μ 0 0 α 1 Λ μ p β 1 Λ μ 0 0 0 0 δ q 0 0 0 0 0 0 α 2 Λ μ r β 2 Λ μ 0 0 0 0 E w 0 0 0 θ 1 0 θ 2 μ
It follows that from the Jacobian matrix J ( E 0 ) , we have the characteristic equation
λ + μ 2 λ 2 + λ a 1 + a 2 λ 2 + λ b 1 + b 2 = 0
where the coefficients are
a 1 = ( q + p ) ( 1 R 1 )
a 2 = p q 1 R 01
b 1 = ( w + r ) ( 1 R 2 )
b 2 = r w 1 R 02 .
In addition, we have
R 1 = α 1 Λ μ q + p ,
and
R 2 = α 2 Λ μ w + r .
Based on Equation (15), the eigenvalues are λ 1 = λ 2 = μ containing a negative real part, and other eigenvalues are obtained from the solution of the following equations:
λ 2 + λ a 1 + a 2 = 0
and
λ 2 + λ b 1 + b 2 = 0 .
Using the Routh–Hurwitz criteria, the characteristic Equations (16) and (17) have negative real roots if the coefficients of a 1 , a 2 > 0 and b 1 , b 2 > 0 . The coefficients a 1 and a 2 are positive if R 1 < 1 and R 01 < 1 , respectively.
We look at the relationship between R 1 and R 01 .
R 01 R 1 = α 1 Λ q + β 1 Λ δ μ p q α 1 Λ μ q + p = α 1 Λ μ q 2 + β 1 Λ δ μ q + β 1 Λ δ μ p μ 2 p q q + p .
Since all the parameters are assumed to be positive, then R 01 R 1 > 0 R 01 > R 1 . Hence, we obtain R 1 < R 01 < 1 . Next, the conditions for the coefficients b 1 and b 2 will be determined. The coefficients b 1 and b 2 are positive if R 2 < 1 and R 02 < 1 , respectively. We now check the relationship between R 2 and R 02 .
R 02 R 2 = α 2 Λ w + β 2 Λ E μ r w α 2 Λ μ w + r = α 2 Λ μ w 2 + β 2 Λ E μ w + β 2 Λ E μ r μ 2 r w w + r .
Hence, we have the condition R 02 R 2 > 0 R 02 > R 2 . So, R 2 < R 02 < 1 . Using the approach in [39], the argument of the roots of Equation (15) are all larger than σ ϕ 2 if R 01 < 1 and R 02 < 1 . Thus, the DFE COVID-19 model with two strains will be locally asymptotically stable if R 01 < 1 and R 02 < 1 . This shows that if R 0 < 1 , then the COVID-19 disease with two strains can be reduced in the community.

4.3.3. Local Stability of the Endemic Equilibrium for Strain 1

The Jacobian matrix at E 1 is given as follows:
J ( E 1 ) = k 1 α 1 S * β 1 S * α 2 S * β 2 S * 0 k 2 k 3 β 1 S * 0 0 0 0 δ q 0 0 0 0 0 0 k 4 β 2 S * 0 0 0 0 E w 0 0 0 θ 1 0 θ 2 μ
with
k 1 = α 1 E 1 * β 1 I 1 * μ
k 2 = α 1 E 1 * + β 1 I 1 *
k 3 = α 1 S * p
k 4 = α 2 S * r
The characteristic equation of Jacobian matrix J ( E 1 ) is
λ + μ λ 2 + λ c 1 + c 2 λ 3 + λ 2 d 1 + λ d 2 + d 3 = 0
with
c 1 = α 2 Λ μ R 01 R 01 R 2 1
c 2 = r w R 02 R 01 R 01 R 02 1
d 1 = α 1 q 2 + β 1 δ q + β 1 δ p + α 1 q μ + β 1 δ μ + R 01 1 ( β 1 δ μ + α 1 μ q ) ( α 1 q + β 1 δ )
d 2 = R 01 1 μ q α 1 + δ μ β 1 p + q + μ α 1 q 2 + δ μ β 1 p + δ μ β 1 q α 1 q + β 1 δ
d 3 = R 01 1 ( α 1 p μ q 2 + β 1 p q δ μ ) α 1 q + β 1 δ .
The eigenvalue of matrix E 1 is λ 1 = μ and the other eigenvalues are the roots of the following characteristic equation:
λ 2 + λ c 1 + c 2 = 0
and
λ 3 + λ 2 d 1 + λ d 2 + d 3 = 0 .
By using the Routh–Hurwitz criterion, the solution of Equation (19) has a negative real part if the coefficients c 1 , c 2 > 0 . The coefficient c 1 and c 2 are positive, respectively, if they satisfy the following criteria
(i).
R 01 R 2 1 > 0 R 01 > R 2 .
(ii).
R 01 R 02 1 > 0 R 01 > R 02 .
Using the relation R 2 and R 02 , i.e., R 02 > R 2 , the two conditions stated in (i) and (ii) imply that R 01 > R 02 > R 2 .
Furthermore, we will examine the solution of (20) that has negative real roots if the coefficients d 1 , d 3 > 0 , and d 1 d 2 d 3 > 0 . It is clear that d 1 , d 3 are positive whenever R 01 > 1 . Next, the condition d 1 d 2 d 3 > 0 can be satisfied if R 01 1 > 0 R 01 > 1 . Thus, it follows from the Routh–Hurwitz criteria [39] that the COVID-19 model at E 1 is asymptotically stable if and only if R 01 > 1 and R 01 > R 02 > R 2 .

4.3.4. Local Stability of the Endemic Equilibrium for Strain 2

The Jacobian matrix at the endemic equilibrium E 2 is given by
J ( E 2 ) = m 1 α 1 S * * β 1 S * * α 2 S * * β 2 S * * 0 0 m 3 β 1 S * * 0 0 0 0 δ q 0 0 0 m 2 0 0 m 4 β 2 S * * 0 0 0 0 E w 0 0 0 θ 1 0 θ 2 μ
with
m 1 = α 2 E 2 * * β 2 I 2 * * μ
m 2 = α 2 E 2 * * + β 2 I 2 * *
m 3 = α 1 S * * p
m 4 = α 2 S * * r .
We obtain one eigenvalue λ = μ , and the remaining eigenvalues containing negative real parts can be determined through the equations given by
λ 2 + λ e 1 + e 2 λ 3 + λ 2 f 1 + λ f 2 + f 3 = 0 ,
where
e 1 = α 1 Λ μ R 02 R 02 R 1 1
e 2 = p q R 01 R 02 R 02 R 01 1
f 1 = α 2 ω 2 + α 2 μ w + β 2 E w + β 2 E r + β 2 E μ + R 02 1 ( α 2 μ w + β 2 E μ ) ( w α 2 + β 2 E )
f 2 = α 2 μ w 2 + β 2 E μ w + β 2 E μ r + R 02 1 w + r ( μ w α 2 + E μ β 2 ) ( α 2 w + β 2 E )
f 3 = R 02 1 ( α 2 r μ w 2 + β 2 r w E μ ) α 2 w + β 2 E
Based on the Equation (21), we will determine the roots of the following characteristic equation:
λ 2 + λ e 1 + e 2 = 0 ,
and
λ 3 + λ 2 f 1 + λ f 2 + f 3 = 0 .
Similarly, from the Routh–Hurwitz criteria, Equation (22) has negative real roots if the coefficients of e 1 , e 2 > 0 . It can be seen that e 1 > 0 and e 2 > 0 if R 02 > R 1 and R 02 > R 01 , respectively.
The characteristic Equation (23) has negative real roots if the coefficients f 1 , f 3 > 0 , and f 1 d 2 f 3 > 0 . It can be observed that f 1 , f 3 are positive whenever R 02 > 1 . Next, the condition f 1 f 2 f 3 > 0 can be satisfied if R 02 1 > 0 R 02 > 1 . Thus, it follows from the Routh–Hurwitz criteria [39] that the COVID-19 model at E 2 is asymptotically stable if and only if R 02 > 1 and R 02 > R 01 > R 1 .

5. Numerical Simulations

This section is devoted to carrying out the simulation numerically to quantify some of the analytical findings and give predictions using our model. We now give the parameter values and method of estimation in the following subsection.

5.1. Parameter Estimation

The parameter values for our simulation are obtained from the published literature, while some were estimated when the model reached its endemic steady state. All the parameterized are now given in Table 2 below.

5.2. Local Sensitivity Analysis

Sensitivity analysis is needed to determine the most influential parameters in the spread of the disease. The parameter values used during the sensitivity analysis are presented in Table 2. Sensitivity analysis aims to determine the parameters that have a major influence on a model [40]. This can be known through the sensitivity index of each parameter. The parameter sensitivity index ( i n d m ) is given by
i n d m = R 0 m m R 0
where m is the parameter to be analyzed. In this work, R 0 is used as follows:
R 01 = α 1 Λ q + δ β 1 Λ μ p q , R 02 = α 2 Λ w + E β 2 Λ μ r w
The results of the parameter sensitivity index calculation can be seen in Table 3 below.
Based on Table 3, it can be seen that the parameters Λ , α 1 , β 1 , ( Λ , α 2 , β 2 ) have positive indices showing that R 01 ( R 02 ) increases with the parameters. The rest of the parameters have negative sensitivity indices, which demonstrate that R 01 and R 02 decrease with the increase in these parameters. For example, the sensitivity index of the parameter β 1 is 0.9969, which mean that increasing (decreasing) β 1 10% will increase (decrease) the value R 01 by 9.969%. In contrast, increasing (decreasing) parameter θ 1 10% will decrease (increase) R 01 by 9%. Similarly, the sensitivity index of parameter β 2 is 0.9959, which implies increasing (decreasing) β 2 10% will increase (decrease) the value R 02 by 9.959%. Thus, the increasing (decreasing) parameter θ 2 10% will decrease (increase) R 02 by 8.71%.
Figure 2 depicts the impact of the parameters β 1 , θ 1 , β 2 , and θ 2 on R 0 . Based on the contour plot in Figure 2, it can be seen that R 0 is increasing with the increase of β 1 and β 2 , while it is decreased with the increase in θ 1 and θ 2 .

5.3. Simulation Results Using Caputo Operator

Here, we devote the numerical simulations for the fractional model (9) in the Caputo operator framework. The fractional model is solved numerically by taking the method as mentioned in [41] and the parameter values set out in Table 2.
Figure 3a characterizes susceptible individuals ( S ) as the fractional-order derivative moves from 0.6 to the integer-order 1. From this figure, we can see that the number of susceptible individuals decreases as time increases. Naturally, this situation occurs in many epidemiological models as more individuals become infected a few weeks after the outbreak and leave the ( S ) class. Figure 3b depicts individuals exposed to the strain 1 virus. The population of E 1 increases steadily as the fractional-order derivative increases toward the classical value. As it approaches the value 1, the increase becomes faster. This result is expected because more susceptible individuals contact the disease and move into the exposed class after the first few weeks of the outbreak. Figure 3c illustrates changes in the number of strain-1 infected individuals with time. The number of strain 1-infected individuals increases as the fractional-order moves toward the classical value 1. Here, the fractional order is very sensitive and drives the infection. The increase becomes faster as the fractional order becomes very close to the value 1. This is motivated by the fact that most of the exposed individuals move to the symptomatic and infectious stage a few weeks after infection. Figure 3d depicts the dynamics of individuals exposed to the strain of two viruses. As the fractional-order derivative moves toward the integer part, the population of individuals exposed to the strain two virus increases. This result is similar to the result obtained in Figure 3b, showing that the fractional order has the same effect on the exposed population regardless of the strain that the population is exposed to during the time of the infection. Similar explanations can be used to relate the similarity between the result of Figure 3c and that of Figure 3e. Figure 3f shows the changes in the number of recovered individuals from both strain one and strain two infections as the fractional order changes. The number of recovered individuals increases as the fractional-order derivative moves up toward the classical value. This steady increase originates from the fact that individuals are recovering from both strains simultaneously. The recovery of infected individuals is usually encouraged to facilitate disease containment. Increasing the fractional order is thus suggested here as we observe on the graph that the increase in the population of the R class is much faster when the fractional order is high.
Figure 4a,b depict the changes that occur in the class I 1 when the infection rates are varied. We observe a steady increase in the number of individuals in the I 1 class as the strain one infection rate increases. This increase is expected because more susceptible individuals come into contact with strain one infectious individuals and contract the virus. On the other hand, increasing the strain two infection rate gives the reverse results. The steady fall in the number of individuals in class I 1 as the strain two infection rate increases is only an indication that increasing the strain two infection rate and leaving the strain one infection rate constant only leads to an increase in the strain two infectious population, giving time for the strain one infectious population to decrease. Figure 4c,d, on the other hand, show the effects of changes in the recovery rates on the population of class I 1 . In Figure 4c, we observe that when the rate of recovery from the strain one infection decreases from 0.16978 to 0.15978, there is a corresponding steady increase in the population of class I 1 . However, a decrease in the recovery rate from the strain one infection from 0.15978 to 0.14978 skyrockets the population of class I 1 . This can be explained by the fact that when the rate of recovery from the disease is very low, the prevalence of infectious individuals and hence the chances of disease spread will increase. On the contrary, an increase in the number of individuals in class I 1 is attained upon increasing the recovery rate from the strain two infections. The results obtained in Figure 5 can be explained in a similar fashion.

6. Conclusions

In this paper, we proposed and analyzed a fractional model of COVID-19 to track the impact of multiple stains of the SARS-CoV-2 variants during the epidemic. Mathematical analyses of the model, which envisage the positivity of solutions, basic reproduction number, and the local stability of the model, were all presented. Our results suggest that an increase in the number of infections has two variants driving the disease than only one variant. This result agrees that most countries have also experienced higher mortality and infection rates during the second or third wave than the first wave of the COVID-19 pandemic.
We acknowledge that our model is not without shortcomings. Some of the limitations are as follows: first, the model was not fitted to historic epidemiological COVID-19 data for a specific setting. Second, few of the parameter values were available in the published literature, and some were estimated where possible. Lastly, our model did not take into account the asymptomatic individuals; this may lead to not accounting for all the possible routes of infectivity in the modeling or the transmissions of the disease during the epidemic. Notwithstanding these limitations, we suggest that our findings remain valid and applicable in the control of the COVID-19 infection generated by the combination of two different strains of the SARS-CoV-2 virus if implemented by policymakers. The model can be extended by considering specific mutations of the variants at different stages and developing a multi-scale model to analyze the effects of these variants on the human population.

Author Contributions

Conceptualization, F. and E.Y.; methodology, F. and E.Y.; software, F. and C.W.C.; validation, F., C.A. and C.W.C.; formal analysis, F., E.Y. and C.W.C.; investigation, F. and C.W.C.; data curation, F. and C.W.C.; writing—original draft preparation, F., E.Y., M.L.J. and C.W.C.; writing—review and editing, F. and C.W.C.; visualization, F. and C.W.C.; supervision, F.; project administration, F.; funding acquisition, F. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

No data were used for the production of this manuscript.

Acknowledgments

The authors would like to thank their respective universities and the anonymous potential reviewers for the production of this paper.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. World Health Organization (WHO). Coronavirus Disease (COVID-19) Pandemic. Available online: https://www.who.int/emergencies/diseases/novel-coronavirus-2019 (accessed on 30 April 2022).
  2. Thakur, V.; dan Jain, A. COVID-19-Suicides: A Global Psychological Pandemic. Brain Behav. Immun. 2020, 88, 952–953. [Google Scholar] [CrossRef] [PubMed]
  3. CDC. Your Health. 2020. Available online: https://www.cdc.gov/coronavirus/2019-ncov/your-health/index.html (accessed on 28 April 2022).
  4. World Health Organisation. COVID-19 Weekly Epidemiological Update Edition 70. Available online: https://www.who.int/docs/defaultsource/coronaviruse/situationreports/20211214_weekly_epi_update_70.pdf?sfvrsn=ad19bf83_3 (accessed on 14 May 2022).
  5. Guruprasad, L. Human Coronavirus Spike Protein-Host Receptor Recognition. Prog. Biophys. Mol. Biol. 2021, 161, 39–53. [Google Scholar] [CrossRef]
  6. WHO. Media Statement: Knowing the Risks for COVID-19. Available online: https://www.who.int/indonesia/news/detail/08-03-2020-knowing-the-risk-for-covid-19 (accessed on 29 May 2022).
  7. Korber, B.; Fischer, W.M.; Gnanakaran, S.; Yoon, H.; Theiler, J.; Abfalterer, W.; Hengartner, N.; Giorgi, E.E.; Bhattacharya, T.; Foley, B.; et al. Evaluating the Effects of SARS-CoV-2 Spike Mutation D614G on Transmissibility and Pathogenicity. Cell 2021, 184, 64–75. [Google Scholar]
  8. Korber, B.; Fischer, W.M.; Gnanakaran, S.; Yoon, H.; Theiler, J.; Abfalterer, W.; Hengartner, N.; Giorgi, E.E.; Bhattacharya, T.; Foley, B.; et al. Tracking changes in SARS-CoV-2 spike: Evidence that D614G increases infectivity of the COVID-19 virus. Cell 2020, 182, 812–827. [Google Scholar] [CrossRef] [PubMed]
  9. CDC. SARS-CoV-2 Variant Classifications and Definitions. Available online: https://www.cdc.gov/coronavirus/2019-ncov/variants/variant-classifications.html (accessed on 25 May 2022).
  10. World Health Organisation. Classification of Omicron (B.1.1.529): SARS-CoV-2 Variant of Concern. Available online: https://www.who.int/news/item/26-11-2021-classification-of-omicron-(b.1.1.529)-sars-cov-2-variant-of-concern (accessed on 25 April 2022).
  11. World Health Organisation. Network for Genomic Surveillance South Africa (NGS-SA). SARS-CoV-2 Sequencing. Available online: https://www.krisp.org.za/manuscripts/25Nov2021_B.1.1.529_Media.pdf (accessed on 25 September 2021).
  12. UK Health Security Agency. SARS-CoV-2 Variants of Concern and Variants under Investigation in England: Technical Briefing 31. Available online: https://assets.publishing.service.gov.uk/government/uploads/system/uploads/attachment_data/file/1042367/technical_briefing-31-10-december-2021.pdf (accessed on 10 March 2022).
  13. Wang, X.; Tang, T.; Cao, L.; Aihara, K.; Guo, Q. Inferring key epidemiological parameters and transmission dynamics of COVID-19 based on a modified SEIR model. Math. Model. Nat. Phenom. 2020, 15, 74. [Google Scholar] [CrossRef]
  14. Zeb, A.; Alzahrani, E.; Erturk, V.S.; Zaman, G. Mathematical model for coronavirus disease 2019 (COVID-19) containing isolation class. BioMed Res. Int. 2020, 2020, 3452402. [Google Scholar] [CrossRef]
  15. Aldila, D. Analyzing the impact of the media campaign and rapid testing for COVID-19 as an optimal control problem in East Java, Indonesia. Chaos Solitons Fractals 2020, 141, 109953. [Google Scholar] [CrossRef]
  16. Khan, M.A.; Atangana, A.; Alzahrani, E.; Fatmawati. The dynamics of COVID-19 with quarantined and isolation. Adv. Differ. Equations 2020, 1, 1–22. [Google Scholar] [CrossRef]
  17. Mishra, A.M.; Purohit, S.D.; Owolabi, K.M.; dan Sharma, Y.D. A Nonlinear Epidemiological Model Considering Asymtotic and Quarantine Classes for SARS-CoV-2 Virus. Chaos Solitons Fractals 2020, 138, 109953. [Google Scholar] [CrossRef]
  18. Aldila, D.; Ndii, M.Z.; Samiadji, B.M. Optimal control on COVID-19 eradication program in Indonesia under the effect of community awareness. Math. Biosci. Eng. 2020, 17, 6355–6389. [Google Scholar] [CrossRef]
  19. Alqahtani, R.T.; Ajbar, A. Study of Dynamics of a COVID-19 Model for Saudi Arabia with Vaccination Rate, Saturated Treatment Function and Saturated Incidence Rate. Mathematics 2021, 9, 3134. [Google Scholar] [CrossRef]
  20. Mushanyu, J.; Chukwu, W.; Nyabadza, F.; Muchatibaya, G. Modelling the potential role of super spreaders on COVID-19 transmission dynamics. Int. J. Math. Model. Numer. Optim. 2022, 12, 191–209. [Google Scholar]
  21. Rezapour, S.; Mohammadi, H.; Samei, M.E. SEIR epidemic model for COVID-19 transmission by Caputo derivative of fractional order. Adv. Differ. Equ. 2020, 1, 1–19. [Google Scholar] [CrossRef] [PubMed]
  22. Naik, P.A.; Yavuz, M.; Qureshi, S.; Zu, J.; Townley, S. Modeling and analysis of COVID-19 epidemics with treatment in fractional derivatives using real data from Pakistan. Eur. Phys. J. Plus 2020, 135, 795. [Google Scholar] [CrossRef] [PubMed]
  23. Boudaoui, A.; El hadj Moussa, Y.; Hammouch, Z.; Ullah, S. A fractional-order model describing the dynamics of the novel coronavirus (COVID-19) with nonsingular kernel. Chaos Solitons Fractals 2021, 146, 110859. [Google Scholar] [CrossRef]
  24. Chukwu, C.W.; Fatmawati. Modelling fractional-order dynamics of COVID-19 with environmental transmission and vaccination: A case study of Indonesia. AIMS Math. 2022, 7, 4416–4438. [Google Scholar] [CrossRef]
  25. Li, X.P.; DarAssi, M.H.; Khan, M.A.; Chukwu, C.W.; Alshahrani, M.Y.; Al Shahrani, M.; Riaz, M.B. Assessing the potential impact of COVID-19 Omicron variant: Insight through a fractional piecewise model. Results Phys. 2022, 38, 105652. [Google Scholar] [CrossRef]
  26. Bonyah, E.; Juga, M.; Fatmawati. Fractional dynamics of coronavirus with comorbidity via Caputo-Fabrizio derivative. Commun. Math. Biol. Neurosci. 2022, 2022, 12. [Google Scholar]
  27. Gao, W.; Veeresha, P.; Cattani, C.; Baishya, C.; Baskonus, H.M. Modified Predictor–Corrector Method for the Numerical Solution of a Fractional-Order SIR Model with 2019-nCoV. Fractal Fract. 2022, 6, 92. [Google Scholar] [CrossRef]
  28. Qureshi, S.; Bonyah, E.; Shaikh, A.A. Classical and contemporary fractional operators for modeling diarrhea transmission dynamics under real statistical data. Phys. Stat. Mech. Its Appl. 2022, 535, 122496. [Google Scholar] [CrossRef]
  29. Bonyah, E.; Juga, M.L.; Chukwu, C.W.; Fatmawati. A fractional order dengue fever model in the context of protected travelers. Alex. Eng. J. 2022, 61, 927–936. [Google Scholar] [CrossRef]
  30. Bonyah, E.; Chukwu, C.W.; Juga, M.; Fatmawati. Modeling fractional order dynamics of Syphilis via Mittag-Leffler law. AIMS Math. 2021, 6, 8367–8389. [Google Scholar] [CrossRef]
  31. Ade León, U.A.P.; Avila-Vales, E.; Huang, K.L. Modeling COVID-19 dynamic using a two-strain model with vaccination. Chaos Solitons Fractals 2022, 157, 111927. [Google Scholar] [CrossRef] [PubMed]
  32. Arruda, E.F.; Das, S.S.; Dias, C.M.; Pastore, D.H. Modelling and optimal control of multi strain epidemics, with application to COVID-19. PLoS ONE 2021, 16, e0257512. [Google Scholar] [CrossRef] [PubMed]
  33. Li, R.; Li, Y.; Zou, Z.; Liu, Y.; Li, X.; Zhuang, G.; Shen, M.; Zhang, L. Evaluating the impact of SARS-CoV-2 variants on the COVID-19 epidemic and social restoration in the United States: A mathematical modelling study. Front. Public Health 2021, 9, 2067. [Google Scholar] [CrossRef]
  34. González-Parra, G.; Arenas, A.J. Qualitative analysis of a mathematical model with presymptomatic individuals and two SARS-CoV-2 variants. Comput. Appl. Math. 2021, 40, 1–25. [Google Scholar] [CrossRef]
  35. Layton, A.T.; Sadria, M. Understanding the dynamics of SARS-CoV-2 variants of concern in Ontario, Canada: A modeling study. Sci. Rep. 2022, 12, 1–16. [Google Scholar] [CrossRef]
  36. Podlubny, I. Fractional Differential Equations: An Introduction to Fractional Derivatives, Fractional Differential Equations, to Methods of Their Solution and Some of Their Applications; Academic Press: San Diego, CA, USA, 1999. [Google Scholar]
  37. Vargas-De-Le´on, C. Volterra-type Lyapunov functions for fractional-order epidemic systems. Commun. Nonlinear Sci. Numer. Simul. 2015, 24, 75–85. [Google Scholar] [CrossRef]
  38. Van den Driessche, P.; Watmough, J. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. Bull. Math. Biol. 2008, 70, 29–48. [Google Scholar] [CrossRef]
  39. Ahmed, E.; El-Sayed, A.M.A.; El-Saka, H.A.A. On some Routh-Hurwitz conditions for fractional order differential equations and their applications in Lorenz, Rossler, Chua and Chen systems. Phys. Lett. 2006, 1, 1–4. [Google Scholar] [CrossRef]
  40. Chitnis, N.; Hyman, J.M.; Cushing, J.M. Determining important parameters in the spread of malaria through the sensitivity analysis of a mathematical model. Bull. Math. Biol. 2008, 70, 1272–1296. [Google Scholar] [CrossRef] [PubMed]
  41. Diethelm, K.; Freed, A.D. The FracPECE subroutine for the numerical solution of differential equations of fractional order. Forsch. Und Wiss. Rechn. 1999, 57–71. [Google Scholar]
Figure 1. Transmission diagram of the COVID-19 model with two strains.
Figure 1. Transmission diagram of the COVID-19 model with two strains.
Fractalfract 06 00346 g001
Figure 2. Contour plot of (a) θ 1 versus β 1 and (b) θ 2 versus β 2 as a function of R 0 .
Figure 2. Contour plot of (a) θ 1 versus β 1 and (b) θ 2 versus β 2 as a function of R 0 .
Fractalfract 06 00346 g002
Figure 3. Numerical simulation for the fractional model (9) when σ = 1 , 0.9 , 0.8 , 0.7 and 0.6 , where (a) Susceptible individuals, (b) Exposed individuals to strain 1, (c) Infected individuals to strain 1, (d) Exposed individuals to strain 2, (e) Infected individuals to strain 2, and (f) Recovered individuals.
Figure 3. Numerical simulation for the fractional model (9) when σ = 1 , 0.9 , 0.8 , 0.7 and 0.6 , where (a) Susceptible individuals, (b) Exposed individuals to strain 1, (c) Infected individuals to strain 1, (d) Exposed individuals to strain 2, (e) Infected individuals to strain 2, and (f) Recovered individuals.
Fractalfract 06 00346 g003
Figure 4. Effect of the model parameters (a) β 1 , ( b ) β 2 , ( c ) θ 1 , and (d) θ 2 on I 1 individuals with fractional order σ = 0.95 . .
Figure 4. Effect of the model parameters (a) β 1 , ( b ) β 2 , ( c ) θ 1 , and (d) θ 2 on I 1 individuals with fractional order σ = 0.95 . .
Fractalfract 06 00346 g004
Figure 5. Effect of the model parameters (a) β 1 , ( b ) β 2 , ( c ) θ 1 , and (d) θ 2 on I 2 individuals with fractional order σ = 0.95 . .
Figure 5. Effect of the model parameters (a) β 1 , ( b ) β 2 , ( c ) θ 1 , and (d) θ 2 on I 2 individuals with fractional order σ = 0.95 . .
Fractalfract 06 00346 g005aFractalfract 06 00346 g005b
Table 1. Parameters of the COVID-19 model with two strains.
Table 1. Parameters of the COVID-19 model with two strains.
SymbolDescriptionUnits
Λ Recruitment rate I n d i v i d u a l t i m e
μ Natural death rate 1 t i m e
α 1 Transmission rate by E 1 1 t i m e
α 2 Transmission rate by E 2 1 t i m e
β 1 Transmission rate by I 1 1 t i m e
β 2 Transmission rate by I 2 1 t i m e
δ Progression rate from E 1 to I 1 1 t i m e
E Progression rate from E 2 to I 2 1 t i m e
σ 1 Death rate due to strain 1 1 t i m e
σ 2 Death rate due to strain 2 1 t i m e
θ 1 Recovery rate of infected with strain 1 1 t i m e
θ 2 Recovery rate of infected with strain 2 1 t i m e
Table 2. Parameter values of the COVID-19 model with two strains.
Table 2. Parameter values of the COVID-19 model with two strains.
ParametersValueSource
Λ 1930[15]
μ 1 70 × 365 [15]
α 1 2 × 10 11 [21]
β 1 1.2644 × 10 8 [15]
δ 0.0784 [15]
σ 1 0.015 [15]
θ 1 0.13978 [15]
α 2 2.65 × 10 11 Assumed
β 2 1.3 × 10 8 Assumed
E 0.0564 Assumed
σ 2 0.0143 Assumed
θ 2 0.1 Assumed
Table 3. The sensitivity index of the parameters.
Table 3. The sensitivity index of the parameters.
ParameterSensitivity Index ( R 01 )ParameterSensitivity Index ( R 02 )
Λ 0.9999 Λ 1
μ −1 μ −1
α 1 0.0031 α 2 0.0041
β 1 0.9969 β 2 0.9959
δ −0.0026 ε −0.0034
σ 1 −0.0966 σ 2 −0.1245
θ 1 −0.9 θ 2 −0.871
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Fatmawati; Yuliani, E.; Alfiniyah, C.; Juga, M.L.; Chukwu, C.W. On the Modeling of COVID-19 Transmission Dynamics with Two Strains: Insight through Caputo Fractional Derivative. Fractal Fract. 2022, 6, 346. https://doi.org/10.3390/fractalfract6070346

AMA Style

Fatmawati, Yuliani E, Alfiniyah C, Juga ML, Chukwu CW. On the Modeling of COVID-19 Transmission Dynamics with Two Strains: Insight through Caputo Fractional Derivative. Fractal and Fractional. 2022; 6(7):346. https://doi.org/10.3390/fractalfract6070346

Chicago/Turabian Style

Fatmawati, Endang Yuliani, Cicik Alfiniyah, Maureen L. Juga, and Chidozie W. Chukwu. 2022. "On the Modeling of COVID-19 Transmission Dynamics with Two Strains: Insight through Caputo Fractional Derivative" Fractal and Fractional 6, no. 7: 346. https://doi.org/10.3390/fractalfract6070346

Article Metrics

Back to TopTop