Establishment of a Fractional Order COVID-19 Model and Its Feasibility Analysis

Rong Kang, Tianzeng Li^*, Yu Zhao

School of Mathematics and Statistics, Sichuan University of Science and Engineering, Zigong, China.

DOI: 10.4236/jcc.2024.1210006                

Abstract

This paper investigates an improved SIR model for COVID-19 based on the Caputo fractional derivative. Firstly, the properties of the model are studied, including the feasible domain and bounded solutions of the system. Secondly, the stability of the system is discussed, among other things. Then, the GMMP method is introduced to obtain numerical solutions for the COVID-19 system. Numerical simulations were conducted using MATLAB, and the results indicate that our model is valuable for studying virus transmission.

Keywords

Grid Approximation Method, COVID-19, Infectious Disease Model, Fractional-Order Derivative

1. Introduction

The emergence and continued spread of infectious diseases pose a threat to the health and lives of humans and other organisms, affecting high-quality development in various aspects such as the economy. Moreover, infectious diseases are not only diverse in classification but also may have their transmission characteristics altered due to mutations in their pathogenic agents. Studying the dynamics of infectious disease models helps predict the development trends of these diseases, providing strong theoretical support for formulating effective and reasonable control strategies. Therefore, to address the risks and challenges posed by various infectious diseases, it is particularly important and practical to use infectious disease dynamics to study the dynamic behavior of different infectious disease models [1].

In 2019, the World Health Organization (WHO) reported the first case of unknown pneumonia. This pneumonia was caused by a coronavirus, a virus that causes infections in both animals and humans, and was later named novel coronavirus (COVID-19) by WHO on February 11, 2020 [2]. COVID-19 can be transmitted through various routes such as coughing, respiratory droplets, nasal secretions from sneezing, aerosol transmission, salivary droplets, and surface physical contact. Signs and symptoms of COVID-19 through the contact route of transmission may appear 2 to 14 days after exposure to the virus. Common signs and symptoms may include fatigue, cough, fever, and early on there may be loss of taste and smell; other symptoms may include headache, nausea, weakness, diarrhea, and runny nose [3]. In severe cases of infection, breathing difficulties and inability to stay awake may occur. Hence, research into the spread of COVID-19 is necessary, and several mathematical models have been developed for this purpose. The SIR model is still the most prevalent, which presents a simple and effective approach that categorizes the population into susceptible individuals denoted as S, infected individuals denoted as I, and recovered individuals denoted as R [4].

Fractional order calculus, a generalized version of integer order calculus, has been widely explored in various fields [5]-[7]. While initially focusing on theoretical knowledge, recent years have witnessed the increasing use of fractional order calculus due to advancements in society. This is because fractional order calculus is better suited to describe certain physical phenomena. One advantage of fractional order calculus is the flexibility it provides in the choice of real derivatives or integral operators. Additionally, the fractional order derivative, being a nonlocal operator, incorporates a memory function [8]. This makes the functionality more versatile. Several researchers have already applied fractional order calculus to SIR modeling. Bialic et al. proposed a predictive model for the spread of COVID-19 in Italy and Spain, using a fractional order generalized susceptible infection recovery (SIR) epidemiological model. The accuracy of this model was confirmed by the results of the COVID-19 spread model prediction [9]. Multiple simulations using the SIR model were conducted to study the spread of COVID-19 in various countries. These simulations found that the predicted number of confirmed cases can be further improved for the most part [10]. This is mainly due to the fact that during the course of the outbreak, people who tested positive were quarantined and some asymptomatic infected people did not know whether they were carrying the virus or not, resulting in a shortage of active cases. Thus, we have introduced a novel parameter p that signifies the individuals who have undergone the test and tested positive, as well as those who were infected and quarantined.

In this article, the main contributions are as follows:

1) The improved fractional-order SITPR model is adopted, which is concise and accurate.

2) Analyzing the feasibility of the model, bounded solutions, computing the equilibrium points of the model, and examining the stability of these equilibrium points.

3) The GMMP method is used to obtain the numerical solution of the fractional-order COVID-19 system.

The structure of this article is as follows: Section 2 introduces fractional derivatives and the fractional-order COVID-19 system. In Section 3, some properties of the fractional-order COVID-19 system are presented, including the feasible domain, bounded solutions, and equilibrium points of the system. Section 4 introduced the GMMP numerical method and analyzed the dynamic epidemiology of the COVID-19 model. Section 5 provides a summary of the work presented in this article.

2. Basic Knowledge and the SIR Model

2.1. Factional Derivatives

The utilization of fractional order calculus finds extensive application in various disciplines and, thus, significantly influences the swift progress of scientific investigations. The three most common definitions of fractional order derivatives are: Grünwald-Letnikov (G-L), Caputo, and Riemann-Liouville (R-L) [11].

Definition 1. [12] Define the α-order Riemann-Liouville derivative of the function $f\left(t\right)$ as follows:

${}_a^{RL}{D}_t^{\alpha }f\left(t\right)=\frac{{\text{d}}^n}{\text{d}{t}^n}{}_{a}D{}_t^{-\left(n-\alpha \right)}f\left(t\right)=\frac{1}{\text{Γ}\left(n-\alpha \right)}\frac{{\text{d}}^n}{\text{d}{t}^n}{\displaystyle {\int }_a^t{\left(t-\tau \right)}^{n-\alpha -1}f\left(\tau \right)\text{d}\tau },$ (1)

where $n-1<\alpha <n$ , $n\in Z^{+}$ , and $\text{Γ}\left(z\right)={\displaystyle {\int }_0^{\infty }{t}^{z-1}{\text{e}}^{-t}\text{d}t}$ is a gamma function.

While the R-L definition has a very high status in theoretical analysis, the Caputo derivative with initial value conditions is more applicable in practical applications.

Definition 2. [12] Define the α-order Caputo derivative of the function $f\left(t\right)$ as follows:

${}_a^C{D}_t^{\alpha }f\left(t\right)={}_{0}D{}_t^{-\left(n-\alpha \right)}\frac{{\text{d}}^n}{\text{d}{t}^n}f\left(t\right)=\frac{1}{\text{Γ}\left(n-\alpha \right)}{\displaystyle {\int }_a^t{\left(t-\tau \right)}^{n-\alpha -1}{f}^{\left(n\right)}\left(\tau \right)\text{d}\tau },$ (2)

where $n-1<\alpha <n$ , $n\in Z^{+}$ .

Definition 3. [12] Define the α-order Grünwald-Letnikov derivative of the function $f\left(t\right)$ as follows:

$\begin{array}{c}{}_a^{GL}{D}_t^{\alpha }f\left(t\right)=\underset{\begin{array}{c}h\to 0\\ \varrho h=t\end{array}}{\lim }{h}^{-\alpha }{\displaystyle \sum _{r=0}^{\varrho }{\left(-1\right)}^r}\left(\begin{array}{l}\alpha \hfill \\ r\hfill \end{array}\right)f\left(t-rh\right)\\ ={\displaystyle \sum _{k=0}^n\frac{f^{\left(k\right)}\left(\alpha \right){\left(t-\alpha \right)}^{k-\alpha }}{\Gamma \left(k+1-\alpha \right)}}+\frac{1}{\Gamma \left(\varrho +1-\alpha \right)}\cdot {\displaystyle {\int }_a^t{\left(t-v\right)}^{\varrho -\alpha }}{f}^{\left(\varrho +1\right)}\left(v\right)\text{d}\nu ,\end{array}$ (3)

the function $f\left(t\right)$ possesses continuous derivatives up to order n within the interval $\left[a,t\right]$ .

The R-L derivative is equivalent to the G-L derivative based on the definition of the fractional order derivative. However, it is important to note that the Caputo derivative differs from the R-L derivative and can be represented in the following manner:

${}_a^{RL}{D}_t^{\alpha }f\left(t\right)={}_0{}^{C}D{}_t^{\alpha }f\left(t\right)+{\displaystyle \sum }_{k=0}^{n-1}{r}_k^{\alpha }\left(t\right)f^{\left(k\right)}\left(a\right),$ (4)

here, $n-1<\alpha <n$ , $n\in Z^{+}$ , and $f^{\left(k\right)}\left(a\right)$ with $k=0,1,\cdots ,n-1$ represent the initial conditions, and $r_k^{\alpha }=\frac{t^{k-\alpha }}{\text{Γ}\left(k+1-\alpha \right)}$ .

In this article, our preference lies in employing the Caputo operator and examining and discussing the scenario for the case of $n=1$ , where $\alpha \in \left(0,1\right)$ . The aforementioned Equation (4) can be written as follows:

${}_a^{RL}{D}_t^{\alpha }f\left(t\right)={}_0{}^{C}D{}_t^{\alpha }f\left(t\right)+r_0^{\alpha }f\left(a\right).$ (5)

Here, we consider the case where $r_0^{\alpha }=\frac{t^{-\alpha }}{\text{Γ}\left(1-\alpha \right)}$ .

2.2. The New SITPR Fractional Order COVID-19 System

Multiple simulations using the SIR model were conducted to study the spread of COVID-19 in various countries. These simulations found that the predicted number of confirmed cases can be further improved for the most part [10]. This is mainly due to the fact that during the course of the outbreak, people who tested positive were quarantined and some asymptomatic infected people did not know whether they were carrying the virus or not, resulting in a shortage of active cases. Thus, we have introduced a novel parameter p that signifies the individuals who have undergone the test and tested positive, as well as those who were infected and quarantined, and thus an improved SIR model is extracted in this paper [13]. The expressed by the differential equation as

$\{\begin{array}{l}{S}^{\prime }\left(t\right)=A-\frac{\beta S\left(t\right)I\left(t\right)}{N\left(t\right)}-\delta S\left(t\right),\\ I^{\prime }\left(t\right)=\frac{\beta S\left(t\right)I\left(t\right)}{N\left(t\right)}-pI\left(t\right)-\gamma I\left(t\right)-\left(\delta +d\right)I\left(t\right),\\ T{P}^{\prime }\left(t\right)=pI\left(t\right)-\gamma TP\left(t\right)-\left(\delta +d\right)TP\left(t\right),\\ {R^{\prime }}_1\left(t\right)=\gamma I\left(t\right)-\delta R_1\left(t\right),\\ {R^{\prime }}_2\left(t\right)=\gamma TP\left(t\right)-\delta R_2\left(t\right),\end{array}$ (6)

where

$N\left(0\right)=S\left(0\right)+I\left(0\right)+TP\left(0\right)+R_1\left(0\right)+R_2\left(0\right).$

In the system (6), $TP\left(t\right)$ represents infected individuals who test positive and are thus quarantined, $R_1\left(t\right)$ represents the number of individuals who have been removed (recovered + dead) from the group I, and $R_2\left(t\right)$ represents the number of individuals removed (recovered + dead) from TP. In recent times, the use of fractional differential equations has become increasingly widespread. Many researchers have discovered that models incorporating fractional derivatives can produce more accurate measurement data than traditional classical models, demonstrating a higher level of consistency [14]. Therefore, a fractional order COVID-19 system of the same order is described as follows. The meaning of the parameters is shown in Table 1.

Table 1. Explanation of model parameters and variables.

Variable	Explanation
S	The number of susceptible individuals
I	The number of un-isolated infected
TP	Infected individuals who test positive and are thus quarantined
$R_1$	The number of individuals who have been removed (recovered + dead) from the group I
$R_2$	Number of individuals removed (recovered + dead) from TP
p	Quarantine rate or proportion of infected individuals tested
$\gamma$	Recovery rate
$\beta$	Infection rate
A	Recruitment rate of the susceptible population
d	Mortality of infected individuals
$\delta$	Natural death rate
N	Total population

$\{\begin{array}{l}{\lambda }_{\alpha }{}_0{}^{C}D{}_t^{\alpha }S\left(t\right)=A-\frac{\beta S\left(t\right)I\left(t\right)}{N\left(t\right)}-\delta S\left(t\right),\\ {\lambda }_{\alpha }{}_0{}^{C}D{}_t^{\alpha }I\left(t\right)=\frac{\beta S\left(t\right)I\left(t\right)}{N\left(t\right)}-pI\left(t\right)-\gamma I\left(t\right)-\left(\delta +d\right)I\left(t\right),\\ {\lambda }_{\alpha }{}_0{}^{C}D{}_t^{\alpha }TP\left(t\right)=pI\left(t\right)-\gamma TP\left(t\right)-\left(\delta +d\right)TP\left(t\right),\\ {\lambda }_{\alpha }{}_0{}^{C}D{}_t^{\alpha }{R}_1\left(t\right)=\gamma I\left(t\right)-\delta R_1\left(t\right),\\ {\lambda }_{\alpha }{}_0{}^{C}D{}_t^{\alpha }{R}_2\left(t\right)=\gamma TP\left(t\right)-\delta R_2\left(t\right),\end{array}$ (7)

where $N\left(0\right)=S\left(0\right)+I\left(0\right)+TP\left(0\right)+R_1\left(0\right)+R_2\left(0\right)$ , and $\alpha \in \left(0,1\right)$ .

To keep the system concise, we define $R_i=R_1+R_2$ . Consequently, we obtain

$\{\begin{array}{l}{\lambda }_{\alpha }{}_0{}^{C}D{}_t^{\alpha }S\left(t\right)=A-\frac{\beta S\left(t\right)I\left(t\right)}{N\left(t\right)}-\delta S\left(t\right),\\ {\lambda }_{\alpha }{}_0{}^{C}D{}_t^{\alpha }I\left(t\right)=\frac{\beta S\left(t\right)I\left(t\right)}{N}-pI\left(t\right)-\gamma I\left(t\right)-\left(\delta +d\right)I\left(t\right),\\ {\lambda }_{\alpha }{}_0{}^{C}D{}_t^{\alpha }TP\left(t\right)=pI\left(t\right)-\gamma TP\left(t\right)-\left(\delta +d\right)TP\left(t\right),\\ {\lambda }_{\alpha }{}_0{}^{C}D{}_t^{\alpha }{R}_i\left(t\right)=\gamma I\left(t\right)+\gamma TP\left(t\right)-\delta R_i\left(t\right),\end{array}$ (8)

where $N\left(0\right)=S\left(0\right)+I\left(0\right)+TP\left(0\right)+R_i\left(0\right)$ , and $\alpha \in \left(0,1\right)$ .

It can be observed that the Equations (7) and (8) exhibit dissimilar units on the left and right sides. More precisely, their respective units are expressed as (days)^−α and (days)⁻¹. By dividing the left-hand side of the equations in systems (7) and (8) by ${\lambda }_{\alpha }$ , the dimensions become (days)^α−1, thereby enabling the attainment of identical units on both sides of the system equations, namely (days)⁻¹. Generally, in a fractional order system, it becomes feasible to acquire ${\lambda }_{\alpha }=1$ . The purpose of this paper is to utilize the MH-NMSS-PSO approach in order to identify an $\alpha$ appropriate collection of fractional orders as well as parameters. Subsequent subsections will establish the essential prerequisites for ensuring the boundedness of the system and determining the feasible domain.

3. The Properties of the Model

3.1. Feasible and Bounded Solutions

Denote ${ℝ}_{+}^4:=\left\{\vartheta \in {ℝ}^4:\vartheta \ge 0\right\}$ and $\vartheta \left(t\right)={\left(S\left(t\right),I\left(t\right),TP\left(t\right),R\left(t\right)\right)}^{\text{T}}$ .

To establish the primary outcome, it is necessary to demonstrate that the positivity and invariance of ${ℝ}_{+}^4$ are imperative, which is very helpful for us to prove the boundedness of the solution.

Theorem 3.1. System (8) has a unique solution, which remains within the set of positive real numbers (${ℝ}_{+}^4$ ).

Proof. The uniqueness of the existence of a solution for system (8) is obtained from theorem 3.1 [15]. Subsequently, we will demonstrate that ${ℝ}_{+}^4$ maintains positivity throughout, since

$\begin{array}{l}{{}^C{D}_{0,t}^{\alpha }S\left(t\right)|}_{S=0}=A\ge 0,\\ {{}^C{D}_{0,t}^{\alpha }I\left(t\right)|}_{I=0}=0\ge 0,\\ {{}^C{D}_{0,t}^{\alpha }TP\left(t\right)|}_{TP=0}=pI\ge 0,\\ {{}^C{D}_{0,t}^{\alpha }R\left(t\right)|}_{R=0}=\gamma I+\gamma Tp\ge 0.\end{array}$ (9)

By using the equation above (9) and theorem 1 [16], we get that ${ℝ}_{+}^4$ is positively invariant.

Next, we will obtain the principal outcome of this specific subsection concerning the feasibility range and boundedness of the solution.

Theorem 3.2. The solution for the system (8) that is both feasible and bounded is as follows:

$\mathcal{D}:=\left\{\left(S,I,TP,R\right)\in {ℝ}_{+}^4:0<N\left(t\right)=S\left(t\right)+I\left(t\right)+TP\left(t\right)+R\left(t\right)\le \frac{A}{\delta }\right\}.$ (10)

Proof. For the derivation of the conclusion, we use the linear property of the Caputo derivative and then add up the four equations of the system (8) to get

$\begin{array}{l}{}^{C}D{}_{0,t}^{\alpha }S\left(t\right)+{}^{C}D{}_{0,t}^{\alpha }I\left(t\right)+{}^{C}D{}_{0,t}^{\alpha }TP\left(t\right)+{}^{C}D{}_{0,t}^{\alpha }R\left(t\right)={}^{C}D{}_{0,t}^{\alpha }N\left(t\right)\\ =A-\delta N\left(t\right)-d\left(I\left(t\right)+TP\left(t\right)\right)\le A-\delta N\left(t\right).\end{array}$ (11)

Hence

${}^c{D}_{0,t}^{\alpha }N\left(t\right)\le A-\delta N\left(t\right).$ (12)

Using the lemma 3 [17], we get

$N\left(t\right)\le \left(N\left(0\right)-\frac{A}{\delta }\right)E_{\alpha }\left(-\delta t^{\alpha }\right)+\frac{A}{\delta },$ (13)

where $\alpha \in \left(0,1\right)$ and $E_{\alpha }\left(-\delta t^{\alpha }\right)$ is the Mittag-Leffler function. Since, $E_{\alpha }\left(-\delta t^{\alpha }\right)\ge 0$ and ${\lim }_{t\to \infty }{E}_{\alpha }\left(-\delta t^{\alpha }\right)=0$ , then the following can be obtained from the above equation:

$\underset{t\to \infty }{\lim \sup }N\left(t\right)\le \frac{A}{\delta }.$ (14)

Based on the derivation and proof presented above, we can conclude that a feasible domain and bounded solutions of the system can be determined by the set $\mathcal{D}$ .

3.2. Equilibrium Points and Their Stability Analysis

Given that the initial three equations in the system (8) do not rely on the last variable R, we can assess the system’s equilibrium and stability using the subsequent methodology.

$\{\begin{array}{l}{}^{C}D{}_{0,t}^{\alpha }S\left(t\right)=A-\frac{\beta S\left(t\right)I\left(t\right)}{N\left(t\right)}-\delta S\left(t\right),\\ {}^{C}D{}_{0,t}^{\alpha }I\left(t\right)=\frac{\beta S\left(t\right)I\left(t\right)}{N\left(t\right)}-pI\left(t\right)-\gamma I\left(t\right)-\left(\delta +d\right)I\left(t\right),\\ {}^{C}D{}_{0,t}^{\alpha }TP\left(t\right)=pI\left(t\right)-\gamma TP\left(t\right)-\left(\delta +d\right)TP\left(t\right),\\ {}^{C}D{}_{0,t}^{\alpha }N\left(t\right)=A-\delta N\left(t\right)-d\left(I\left(t\right)+TP\left(t\right)\right).\end{array}$ (15)

Now, to find the stable point of the system (15), let

$\{\begin{array}{l}{}^{C}D{}_{0,t}^{\alpha }S\left(t\right)=0,\\ {}^{C}D{}_{0,t}^{\alpha }I\left(t\right)=0,\\ {}^{C}D{}_{0,t}^{\alpha }TP\left(t\right)=0,\\ {}^{C}D{}_{0,t}^{\alpha }N\left(t\right)=0.\end{array}$ (16)

Next, the system’s equilibrium points can be represented as $E_0=\left(\frac{A}{\delta },0,0,\frac{A}{\delta }\right)$ and $E_1=\left(\overline{S},\overline{I},\overline{TP},\overline{N}\right)$ , where

$\begin{array}{l}\overline{S}=\frac{-A\left(\gamma +\delta \right)}{\delta \left(d-\beta +\frac{p\beta }{p+\gamma +\delta +d}\right)},\\ \overline{I}=\frac{A}{p+\gamma +\delta +d}+\frac{A\left(\gamma +\delta \right)}{d\left(p+\gamma +\delta +d\right)-\beta \left(\gamma +\delta +d\right)},\\ \overline{TP}=\frac{p}{\gamma +\delta +d}\left[\frac{A}{p+\gamma +\delta +d}+\frac{A\left(\gamma +\delta \right)}{d\left(p+\gamma +\delta +d\right)-\beta \left(\gamma +\delta +d\right)}\right],\\ \overline{N}=\frac{-A\beta \left(\gamma +\delta \right)}{\delta \left[d\left(p+\gamma +\delta +d\right)-\beta \left(\gamma +\delta +d\right)\right]}.\end{array}$ (17)

Jacobian matrix is given by

$J=\left(\begin{array}{cccc}-\frac{\beta I}{N}-\delta & -\frac{\beta S}{N}& 0& \frac{\beta SI}{N^2}\\ \frac{\beta I}{N}& \frac{\beta S}{N}-\left(p+\gamma +\delta +d\right)& 0& -\frac{\beta SI}{N^2}\\ 0& p& -\left(\gamma +\delta +d\right)& 0\\ 0& -d& -d& -\delta \end{array}\right).$ (18)

Theorem 3.3. If $\frac{\beta }{p+\gamma +\delta +d}<1$ , the disease-free equilibrium $E_0=\left(\frac{A}{\delta },0,0,\frac{A}{\delta }\right)$ of the system is locally asymptotically stable; otherwise, if $\frac{\beta }{p+\gamma +\delta +d}>1$ it is unstable.

Proof. The Jacobian matrix is obtained by following these steps at $E_0=\left(\frac{A}{\delta },0,0,\frac{A}{\delta }\right)$

$J\left(E_0\right)=\left(\begin{array}{cccc}-\delta & -\beta & 0& 0\\ 0& \beta -\left(p+\gamma +\delta +d\right)& 0& 0\\ 0& p& -\left(\gamma +\delta +d\right)& 0\\ 0& -d& -d& -\delta \end{array}\right).$ (19)

By solving the characteristic equation of the Jacobian matrix $J\left(E_0\right)$ , the eigenvalues of $J\left(E_0\right)$ can be determined. This leads to the equation.

${\left(\lambda +\delta \right)}^2\left(\lambda +\gamma +\delta +d\right)\left[\lambda -\left(\beta -p-\gamma -\delta -d\right)\right]=0.$ (20)

Therefore, the eigenvalues of the Jacobian matrix $J\left(E_0\right)$ are

$\begin{array}{l}{\lambda }_1=-\delta ,\\ {\lambda }_2=-\delta ,\\ {\lambda }_3=\beta -\left(p+\gamma +\delta +d\right),\\ {\lambda }_4=-\left(\gamma +\delta +d\right).\end{array}$ (21)

Based on the analysis presented above, we can observe that if $\frac{\beta }{p+\gamma +\delta +d}<1$ ,

all the eigenvalues will have negative real parts. Consequently, the disease-free equilibrium point $E_0$ will be locally asymptotically stable. However, if

$\frac{\beta }{p+\gamma +\delta +d}>1$ , the disease-free equilibrium point $E_0$ will be unstable. In this context, the ratio $\frac{\beta }{p+\gamma +\delta +d}$ is referred to as the basic reproduction number,

represented by $R_0$ . From a biological perspective, $R_0$ can be interpreted as follows: if $R_0$ is less than one, the infection will eventually disappear; whereas if it is greater than one, the infection will persist. Therefore, we can now proceed to discuss the asymptotic stability of the endemic equilibrium in the given system.

Now, let $M=p+\gamma +\delta +d$ . The obtained Jacobian matrix is denoted as $J\left(E_1\right)$ ,

$J\left(E_1\right)=\left(\begin{array}{cccc}\frac{\delta \left[dM-\beta \left(M-p\right)\right]}{M\left(\gamma +\delta \right)}& -M& 0& \frac{\delta \left[dM-\beta \left(M-p\right)\right]}{-\beta \left(\gamma +\delta \right)}-\frac{M\delta }{\beta }\\ \frac{\delta \left[dM-\beta \left(M-p\right)\right]}{-M\left(\gamma +\delta \right)}-\delta & 0& 0& \frac{\delta \left[dM-\beta \left(M-p\right)\right]}{\beta \left(\gamma +\delta \right)}+\frac{M\delta }{\beta }\\ 0& p& p-M& 0\\ 0& -d& -d& -\delta \end{array}\right).$ (22)

Let $Q=\frac{dM-\beta \left(M-p\right)}{M\left(\gamma +\delta \right)}$ . Then, the characteristic equation for the Jacobian matrix $J\left(E_1\right)$ is to be simplified,

$T\left(\lambda \right)=\left(\delta +\lambda \right)X\left(\lambda \right)\text{and}X\left(\lambda \right)={\lambda }^3+D_1{\lambda }^2+D_2\lambda +D_3,$ (23)

where

$\begin{array}{l}{D}_1=M-p-\delta Q,\\ D_2=\delta Q\left(p-M\right)+\delta \left(Q+M\right)\frac{d-\beta }{\beta },\\ D_3=\delta \left(Q+M\right)\left[\left(p-M\right)+\frac{dM}{\beta }\right].\end{array}$ (24)

Regarding the polynomial $X\left(\lambda \right)$ , we can have the following definition:

$D\left(X\right)=18D_1{D}_2{D}_3+D_1^2{D}_2^2-4D_3{D}_1^3-4D_2^3-27D_3^3.$ (25)

Using Proposition 1 given in the article [18], if we observe that the stable point $E_1$ is asymptotically steady, then $D\left(X\right)>0$ , $D_1>0$ , $D_3>0$ , and $D_1{D}_2>D_3$ .

4. Numerical Solutions for Fractional-Order Differential Equations

4.1. Further Precision of the COVID-19 Model

By combining the equations for $TP\left(t\right)$ and $R_2\left(t\right)$ outlined in system (8), the equation for $C\left(t\right)$ , representing the total number of confirmed cases at time t, can be derived. Now classifying everyone as susceptible, and therefore with a negligible recruitment rate A for the susceptible population, the exact system is as follows:

$\{\begin{array}{l}{\lambda }_{\alpha }{}_0{}^{C}D{}_t^{\alpha }S\left(t\right)=-\frac{\beta S\left(t\right)I\left(t\right)}{N}-\delta S\left(t\right),\\ {\lambda }_{\alpha }{}_0{}^{C}D{}_t^{\alpha }I\left(t\right)=\frac{\beta S\left(t\right)I\left(t\right)}{N}-pI\left(t\right)-\gamma I\left(t\right)-\left(\delta +d\right)I\left(t\right),\\ {\lambda }_{\alpha }{}_0{}^{C}D{}_t^{\alpha }C\left(t\right)=pI\left(t\right),\\ {\lambda }_{\alpha }{}_0{}^{C}D{}_t^{\alpha }{R}_1\left(t\right)=\gamma I\left(t\right)-\delta R_1\left(t\right),\end{array}$ (26)

where $N\left(0\right)=S\left(0\right)+I\left(0\right)+C\left(0\right)+R_1\left(0\right)$ , and $\alpha \in \left(0,1\right)$ .

4.2. Description of the GMMP Method

A number of numerical methods are currently used to find fractional order systems, including the prediction correction method, the Millin transformation method, etc. In this research paper, we employ the GMMP method to compute the numerical solution of the system (26). The system (26) can be expressed as: $\lambda \odot {}_a{}^{C}D{}_t^{\alpha }z\left(t\right)=f\left(t,z\left(t\right)\right)$ , here we can observe that $z\left(t\right)={\left(S\left(t\right),I\left(t\right),C\left(t\right),R_1\left(t\right)\right)}^{\text{T}}$ . Assuming the adoption of a uniform grid over the interval $\left[a,t\right]$ , as follows: $a=t_0<t_1<t_2<\cdots <t_N=t$ , and $\text{Δ}t=t_{i+1}-t_i=h$ . Also assume that $z\left(t\right)$ is sequential in each limited interval $\left(a,t\right)$ .

First, we rely on the classical finite-difference notation [19] to proceed with our analysis and $\alpha \in \left(0,1\right)$ .

$\frac{1}{h^{\alpha }}{\text{Δ}}_h^{\alpha }z\left(t\right)=\frac{1}{h^{\alpha }}\left(z\left(x_n\right)-{\displaystyle \sum }_{k=0}^n{c}_k^{\alpha }z\left(x_{n-k}\right)\right).$ (27)

In the Equation (27), $c_k^{\alpha }={\left(-1\right)}^k\left(\begin{array}{l}\alpha \hfill \\ k\hfill \end{array}\right)$ represents the most commonly used binomial coefficients.

Both R-L and G-L fractional order derivatives can be expressed in the following format according to Equation (27).

${}_a{}^{RL}D{}_t^{\alpha }z\left(t\right)={}_a{}^{GL}D{}_t^{\alpha }z\left(t\right)=\underset{h\to 0}{\lim }\frac{1}{h^{\alpha }}{\Delta }_h^{\alpha }z\left(t\right)\approx \frac{1}{h^{\alpha }}{\Delta }_h^{\alpha }z\left(t\right).$ (28)

The error term $r_0^{\alpha }$ in Equation (5) can also be represented using a uniform grid, leading to the following expression for the Caputo fractional order derivative [20]

${}_a^C{D}_t^{\alpha }z\left(t\right)\approx \frac{1}{h^{\alpha }}\left(z\left(x_n\right)-{\displaystyle \sum }_{k=0}^n{c}_k^{\alpha }z\left(x_{n-k}\right)\right)-r_n^{\alpha }z\left(a\right),$ (29)

where $z\left(a\right)$ is the initial condition, and $r_n^{\alpha }=r_0^{\alpha }\left(t_n\right)={\omega }_{0,-1}^{\alpha }{(n)}^{-\alpha }$ . The function $\omega$ is given by ${\omega }_{\mu ,\nu }^{\alpha }=\frac{\text{Γ}\left(\mu \alpha +1\right)}{\text{Γ}\left(v\alpha +1\right)}$ , $\mu ,\nu \in {\mathbb{N}}_0\cup \left\{-1\right\}$ . And $r_n^{\alpha }z\left(a\right)$ tends to zero when $n\to \infty$ [21].

Table 2 below shows the fundamental procedure of GMMP (Gorenflo-Mainardi-Moretti-Paradisi Programme).

Table 2. The fundamental procedure of GMMP.

Step	Specifics
Step 1	The initial stage of this process involves obtaining the difference grid fromthe solution region, following which a finite number of grids are utilized to represent the infinite solution domain of the differential equation.
Step 2	The calculation is employed to compute the difference quotient,we know that the differential terms in a differential equation can beexpressed by finite differences.
Step 3	The finite number of unknown variables can be found in the differenceequation for the discrete case.

In a general sense, the fractional order nonlinear equation can be expressed in terms of the Caputo operator as follows:

$\begin{array}{l}\lambda \odot {}_a{}^{C}D{}_t^{\alpha }z\left(t\right)=f\left(t,z\left(t\right)\right),0\le t\le T,\\ z^{\left(k\right)}\left(a\right)=z_0^{\left(k\right)},k=0,1,\cdots ,n-1.\end{array}$ (30)

Finally, combining (29) and (30), The equation can be obtained through deduction as follows:

$z\left(x_n\right)=h^{\alpha }\oslash \lambda \odot f\left(t_n,z\left(x_n\right)\right)+{\displaystyle \sum }_{k=0}^n{c}_k^{\alpha }z\left(x_{n-k}\right)+h^{\alpha }{r}_n^{\alpha }{z}_0.$ (31)

The reasoning presented above leads to the ultimate formulation (31) of the fractional order nonlinear equation, which takes the form of an implicit difference scheme. One can interpret (31) as an equation concerning the unfamiliar variable $z\left(x_n\right)$ , which appears on both ends of the equation. Hence, Equation (31) can be utilized to numerically solve any initial value problem associated with a Caputo fractional order differential equation. To effectively tackle such problems, it is necessary to introduce the relevant equations and employ Matlab software (R2020a) for their resolution. Since the Equation (31) possesses an unfamiliar variable $z\left(x_n\right)$ on both ends, the Newton’s algorithm is selected to obtain the value of $z\left(x_n\right)$ through Equation (31).

5. Discussion

The purpose of the numerical simulation is to examine the influence of order and parameter change on the dynamic action of the system. We set the initial values and parameters within an appropriate range as follows:

$N\left(0\right)=20000,S\left(0\right)=19400,I\left(0\right)=500,C\left(0\right)=100,R\left(0\right)=0.$

$\beta =0.8465,\gamma =0.15,d=0.0002,\delta =0.06,p=\mathrm{0.388.}$

Figure 1. The dynamic changes of the system when the order α = 0.7.

Figure 2. The dynamic changes of the system when the order α = 0.75.

Figure 3. The dynamic changes of the system when the order α = 0.8.

Figure 4. The impact of changing the order on (S).

Figure 5. The impact of changing the order on (I).

Figure 6. The impact of changing the order on (C).

Figure 7. The impact of changing the order on (R).

Using Matlab software to generate the plots, we set the order α to 0.70, 0.75, and 0.80, with the results shown in Figures 1-3. It can be observed that each compartment in the model tends to stabilize over time, indicating that under reasonable control, the virus may eventually be eradicated. A smaller value of α leads to a faster trend of change. Finally, we also discussed the specific impacts of different values of α on each compartment. The results consistently stabilize, as shown in Figures 4-7, demonstrating that our model and numerical method are highly feasible, which is also significant for studying the spread of other viruses.

6. Conclusion

In this paper, we study the same order of fractional order COVID-19 mathematical model. Furthermore, an investigation into the SIR model is conducted, which incorporates a novel parameter p into the fundamental SIR model. Our research encompasses an analysis of the positive invariance of ${ℝ}_{+}^4$ , the derivation of the system’s solution’s boundedness and uniqueness, and an examination of the system’s basic reproductive number, denoted as $R_0$ . Additionally, we investigate the equilibrium point and stability of the system. This paper uses the GMMP method to obtain numerical solutions for the system and generates results using Matlab software. However, the paper does not discuss the impact of changes in each parameter on each compartment, which would be a valuable research topic. In the future, more complex models such as SEIR or SEAIRQ can be used to study the spread of the epidemic more comprehensively, taking into account factors like isolation and recovery. Additionally, combining methods such as Particle Swarm Optimization (PSO) for parameter optimization can enhance the accuracy and predictive capability of the model. This combination will contribute to a deeper understanding of epidemic dynamics and the development of more effective response strategies.

Acknowledgements

This work was supported by College Student Innovation and Entrepreneurship Training Program of Sichuan University of Science and Engineering (grant numbers Y2023335), the Fund of the Center for Mathematics and Finance Research (Grant No. SCMF202301).

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Brauer, F. (2017) Mathematical Epidemiology: Past, Present, and Future. Infectious Disease Modelling, 2, 113-127. https://doi.org/10.1016/j.idm.2017.02.001
[2]	Li, Z., Jones, C., Ejigu, G.S., George, N., Geller, A.L., Chang, G.C., et al. (2021) Countries with Delayed COVID-19 Introduction—Characteristics, Drivers, Gaps, and Opportunities. Globalization and Health, 17, Article No. 28. https://doi.org/10.1186/s12992-021-00678-4
[3]	Rezaei, N., Ashkevarian, S., Fathi, M.K., Hanaei, S., Kolahchi, Z., Ladi Seyedian, S., et al. (2021) Introduction on Coronavirus Disease (COVID-19) Pandemic: The Global Challenge. In: Rezaei, N., Ed., Coronavirus Disease—COVID-19, Springer International Publishing, 1-22. https://doi.org/10.1007/978-3-030-63761-3_1
[4]	Uçar, D. and Çelik, E. (2022) Analysis of Covid 19 Disease with SIR Model and Taylor Matrix Method. AIMS Mathematics, 7, 11188-11200. https://doi.org/10.3934/math.2022626
[5]	Jiang, Y., Zhang, B., Shu, X. and Wei, Z. (2020) Fractional-order Autonomous Circuits with Order Larger than One. Journal of Advanced Research, 25, 217-225. https://doi.org/10.1016/j.jare.2020.05.005
[6]	Miao, Y., Gao, Z. and Yang, C. (2022) Adaptive Fractional-Order Unscented Kalman Filters for Nonlinear Fractional-Order Systems. International Journal of Control, Automation and Systems, 20, 1283-1293. https://doi.org/10.1007/s12555-021-0163-4
[7]	Elwy, O., AbdelAty, A.M., Said, L.A., Madian, A.H. and Radwan, A.G. (2020) Two Implementations of Fractional-Order Relaxation Oscillators. Analog Integrated Circuits and Signal Processing, 106, 421-432. https://doi.org/10.1007/s10470-020-01640-x
[8]	Barros, L.C.d., Lopes, M.M., Pedro, F.S., Esmi, E., Santos, J.P.C.D. and Sánchez, D.E. (2021) The Memory Effect on Fractional Calculus: An Application in the Spread of Covid-19. Computational and Applied Mathematics, 40, Article No. 72. https://doi.org/10.1007/s40314-021-01456-z
[9]	Kozioł, K., Stanisławski, R. and Bialic, G. (2020) Fractional-Order SIR Epidemic Model for Transmission Prediction of COVID-19 Disease. Applied Sciences, 10, Article 8316. https://doi.org/10.3390/app10238316
[10]	Malavika, B., Marimuthu, S., Joy, M., Nadaraj, A., Asirvatham, E.S. and Jeyaseelan, L. (2021) Forecasting COVID-19 Epidemic in India and High Incidence States Using SIR and Logistic Growth Models. Clinical Epidemiology and Global Health, 9, 26-33. https://doi.org/10.1016/j.cegh.2020.06.006
[11]	Scherer, R., Kalla, S.L., Tang, Y. and Huang, J. (2011) The Grünwald-Letnikov Method for Fractional Differential Equations. Computers & Mathematics with Applications, 62, 902-917. https://doi.org/10.1016/j.camwa.2011.03.054
[12]	Pan, W., Li, T. and Ali, S. (2021) A Fractional Order Epidemic Model for the Simulation of Outbreaks of Ebola. Advances in Difference Equations, 2021, Article No. 161. https://doi.org/10.1186/s13662-021-03272-5
[13]	Singh, A.K., Mehra, M. and Gulyani, S. (2021) A Modified Variable-Order Fractional SIR Model to Predict the Spread of COVID-19 in India. Mathematical Methods in the Applied Sciences, 46, 8208-8222. https://doi.org/10.1002/mma.7655
[14]	Diethelm, K. (2012) A Fractional Calculus Based Model for the Simulation of an Outbreak of Dengue Fever. Nonlinear Dynamics, 71, 613-619. https://doi.org/10.1007/s11071-012-0475-2
[15]	Lin, W. (2007) Global Existence Theory and Chaos Control of Fractional Differential Equations. Journal of Mathematical Analysis and Applications, 332, 709-726. https://doi.org/10.1016/j.jmaa.2006.10.040
[16]	Odibat, Z.M. and Shawagfeh, N.T. (2007) Generalized Taylor’s Formula. Applied Mathematics and Computation, 186, 286-293. https://doi.org/10.1016/j.amc.2006.07.102
[17]	Li, H., Zhang, L., Hu, C., Jiang, Y. and Teng, Z. (2016) Dynamical Analysis of a Fractional-Order Predator-Prey Model Incorporating a Prey Refuge. Journal of Applied Mathematics and Computing, 54, 435-449. https://doi.org/10.1007/s12190-016-1017-8
[18]	Ahmed, E., El-Sayed, A.M.A. and El-Saka, H.A.A. (2006) On Some Routh-Hurwitz Conditions for Fractional Order Differential Equations and Their Applications in Lorenz, Rössler, Chua and Chen Systems. Physics Letters A, 358, 1-4. https://doi.org/10.1016/j.physleta.2006.04.087
[19]	Gorenflo, R., Mainardi, F., Moretti, D. and Paradisi, P. (2002) Time Fractional Diffusion: A Discrete Random Walk Approach. Nonlinear Dynamics, 29, 129-143. https://doi.org/10.1023/a:1016547232119
[20]	Quintana Murillo, J. and Bravo Yuste, S. (2009) On Three Explicit Difference Schemes for Fractional Diffusion and Diffusion-Wave Equations. Physica Scripta, 136, Article ID: 014025. https://doi.org/10.1088/0031-8949/2009/t136/014025
[21]	Zhang, Y., Yang, X. and Li, X. (2018) Error Estimates of Finite Element Methods for Nonlinear Fractional Stochastic Differential Equations. Advances in Difference Equations, 2018, Article No. 215. https://doi.org/10.1186/s13662-018-1665-0