Numerical Simulation and Parameter Estimation of Fractional-Order Dynamic Epidemic Model for COVID-19

Rong Kang, Tianzeng Li^*

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

DOI: 10.4236/jamp.2024.1210206                

Abstract

The outbreak of COVID-19 in 2019 resulted in numerous infections and deaths. In order to better study the transmission of COVID-19, this article adopts an improved fractional-order SIR model. 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 and combined with the improved MH-NMSS-PSO parameter estimation method to fit the real data of Delhi, India from April 1, 2020 to June 30, 2020. The results show that the fitting effect is quite ideal. Finally, long-term predictions were made on the number of infections. We accurately estimate that the peak number of infections in Delhi, India, can reach around 2.1 million. This paper also compares the fitting performance of the integer-order COVID-19 model and the fractional-order COVID-19 model using the real data from Delhi. The results indicate that the fractional-order model with different orders, as we proposed, performs the best.

Keywords

Parameter Estimation, COVID-19, Infectious Disease Model, Fractional-Order Derivative

1. Introduction

On January 30, 2020, the first case of COVID-19 was detected in Kerala state. The first case was reported in Delhi on March 2 of that year, and the number of cases in Delhi has since grown exponentially. The government has also taken steps to strengthen healthcare services [1]. 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 [2]. 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. Assumptions about the ability of the virus to spread and strategies for preventing epidemics are embodied in the mathematical rules for the transfer of a population from one group to another, resulting in the calculation of the infection rate and peak intervals of the virus [3]. Taimoor et al. applied a time-series-based SIR model to simulate a COVID-19 pandemic in Pakistan [4]. In their analysis of COVID-19, Uçar et al. utilized the SIR model and employed the Taylor matrix method. However, there is a need for further improvement in the accuracy of this method in generating the final simulated data [5]. Meanwhile, Cooper et al. applied the SIR model to predict the transmission of COVID-19 across different communities [6]. On the other hand, Yicheng Chen and Bingen Lu proposed a time-dependent SIR model to observe the spread and recovery rates at specific time points [7]. Additionally, Vinitsky et al. developed a mathematical model that effectively strikes a balance between accurately characterizing the pandemic using the SIR model and maintaining simplicity in practical estimation [8].

Fractional order calculus, a generalized version of integer order calculus, has been widely explored in various fields [9]. 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 [10]. Majee et al. investigated the fractional-order complex dynamics of 2019 coronavirus disease in the context of SIR systems and proposed and analyzed a fractional-order SIR-type epidemiological model incorporating a saturated processing function [11]. Ndairou and Torres [12] used fractional-order modeling to mathematically analyze the neocoronavirus pneumonia in Wuhan, Spain, and Portugal. Lu et al. [13] investigated an epidemiological model based on the CTRW and distributional delay for abnormal diffusion epidemic modeling. Additionally, Joshi et al. [14] analyzed the stability of a non-singular fractional order COVID-19 model with nonlinear treatment and incidence rates. Finally, Li T.Z. et al. [15] summarized the effectiveness of the fractional order viral model in solving practical problems.

In this research, our primary objective is to investigate an advanced model for studying the SIR system of COVID-19. Our specific focus is on utilizing Caputo’s fractional order derivatives. Initially, we conducted an examination of COVID-19 systems displaying fractional orders with identical magnitudes. Subsequently, we propose alternative COVID-19 systems with varying fractional orders. To address this system, we put forward an improved format known as GMMP, resulting in an implicit difference equation as the solution. By implementing Newton’s method, we can effectively solve this nonlinear equation and acquire the numerical solution. Furthermore, we tackle the inverse problem of parameter estimation by employing the MH-NMSS-PSO algorithm in order to estimate the relevant parameters. This innovative algorithm integrates authentic data from the outbreak of New Crown Pneumonia in Delhi, India, spanning from April 1, 2020 to June 30, 2020 [16]. Our findings indicate that by utilizing the new fractional order and parameter values, we are able to obtain a numerical solution for the SIR fractional order COVID-19 system that closely aligns with the present data and exceeds the capabilities of existing models in the literature. Furthermore, our long-term predictions display a significant resemblance to the precise data.

The structure of this article is as follows: Section 3 introduces fractional derivatives and the fractional-order COVID-19 system. In Section 4, some properties of the fractional-order COVID-19 system are presented, including the feasible domain, bounded solutions, and equilibrium points of the system. Section 5 describes the method of obtaining numerical solutions for the COVID-19 system using the GMMP format. Section 6 introduces the improved parameter estimation method, MH-NMSS-PSO. In Section 7, real data from Delhi state in India is used for fitting and long-term predictions. Finally, Section 8 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) [17].

Definition 1. 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. 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. 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 SIR Fractional Order COVID-19 System

The most extensively used infectious disease model is the SIR Model, recognized for its precision and brevity. Specifically, in the traditional SIR Model, the entire population is divided into three distinct groups. The susceptible individuals, denoted by S, have not yet been infected. The infected individuals, denoted by I, have the ability to transmit the disease. Lastly, the recovered individuals, denoted by R, have successfully overcome the infection and acquired immunity, thus being excluded from the infected population [18]. The overall population size, represented by N, remains constant throughout. As a result, the differential equation can be expressed as follows, as shown in Figure 1.

$\frac{\text{d}S\left(t\right)}{\text{d}t}=-\frac{\beta S\left(t\right)I\left(t\right)}{N},$

$\begin{array}{l}\frac{\text{d}I\left(t\right)}{\text{d}t}=\frac{\beta S\left(t\right)I\left(t\right)}{N}-\gamma I\left(t\right),\\ \frac{\text{d}R\left(t\right)}{\text{d}t}=\gamma I\left(t\right),\end{array}$ (6)

where

$N=S\left(0\right)+I\left(0\right)+R\left(0\right).$

Figure 1. SIR model flowchart.

The different parameters in this SIR model denote different meanings and can be expressed as follows:

(i) The infection rate, denoted as $\beta$ , is an important parameter in studying the spread of a disease.

(ii) On the other hand, the recovery rate, denoted as $\gamma$ , is also critical in understanding the dynamics of an epidemic.

Set the parameters to $\beta =0.583$ , $\gamma =0.51$ within the allowable range [19].

Table 1. Selected data on confirmed diagnoses reported from the state of Delhi in India [16].

Date	Cured	Deaths	Confirmed
2020/4/01	6	2	152
2020/4/03	8	4	219
2020/4/06	19	7	523
2020/4/13	27	24	1154
2020/4/23	724	48	2248
2020/4/30	1092	56	3439
2020/5/04	1362	64	4549
2020/5/12	2129	73	7233
2020/5/19	4485	168	10,054
2020/5/22	5567	194	11,659
2020/5/30	7846	398	17,386
2020/6/04	9542	606	23,645
2020/6/07	10,664	761	27,654
2020/6/13	13,398	1214	36,824
2020/6/20	23,569	2035	53,116
2020/6/21	31,294	2112	56,746
2020/6/24	39,313	2301	66,602
2020/6/27	47,091	2492	77,240
2020/6/29	52,607	2623	83,077

Utilizing data provided by the state of Delhi, India, up until April 1, 2020, as shown in Table 1, it is estimated that the total population of Delhi is approximately 190 million [20]. Utilizing this information, we can set the initial conditions as follows:

$S\left(0\right)=189999842,I\left(0\right)=152,R\left(0\right)=6.$

These initial conditions indicate that the number of infections is initially low and will increase if no intervention measures are implemented. To solve the non-linear differential Equation (6), the ODE45 function in Matlab was utilized. Examining Table 1, we observe the number of infections in Delhi, India from April 1 to June 29, 2020. The final results obtained from the simulation are displayed in Figure 2, with the root-mean-square relative error calculated as $g\left(U\right)=0.3937$ . However, it is important to note that further improvements can be made to achieve a more accurate representation of the spread of COVID-19 in Delhi, India.

Figure 2. In the state of Delhi, India, a comparison was made between the number of confirmed cases of neo-coronavirus, denoted as $I\left(t\right)$ , and the numerical results derived from the model (6). The root-mean-square relative error was computed, resulting in a value of $g\left(U\right)=0.3937$ .

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 [21]. 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 [22]. The flow chart of the model is shown in Figure 3 and is expressed by the differential equation as

Figure 3. Flowchart of the SIR model with the introduction of new parameters.

$\{\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),\hfill \\ 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),\hfill \\ T{P}^{\prime }\left(t\right)=pI\left(t\right)-\gamma TP\left(t\right)-\left(\delta +d\right)TP\left(t\right),\hfill \\ {R^{\prime }}_1\left(t\right)=\gamma I\left(t\right)-\delta R_1\left(t\right),\hfill \\ {R^{\prime }}_2\left(t\right)=\gamma TP\left(t\right)-\delta R_2\left(t\right),\hfill \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).$

The meaning of the parameters in the system is shown in Table 2. In the system (7), $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 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 [23]. Therefore, a fractional order COVID-19 system of the same order is described as follows:

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

$\{\begin{array}{l}{}_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),\\ {}_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),\\ {}_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),\\ {}_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)$ .

Table 2. 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
R1	The number of individuals who have been removed (recovered + dead) from the group I
R2	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

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

Theorem 3.1. System (8) has a unique solution, which remains within the set of positive real numbers (${ℝ}_{+}^4$ ). 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}}$ .

Proof. The uniqueness of the existence of a solution for system (8) is obtained from theorem 3.1 [24]. 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 [25], 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 [26], 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}{llll}-\delta \hfill & -\beta \hfill & 0\hfill & 0\hfill \\ 0\hfill & \beta -\left(p+\gamma +\delta +d\right)\hfill & 0\hfill & 0\hfill \\ 0\hfill & p\hfill & -\left(\gamma +\delta +d\right)\hfill & 0\hfill \\ 0\hfill & -d\hfill & -d\hfill & -\delta \hfill \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 [27], 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

To obtain precise numerical results, we have chosen to utilize data exclusively from the Indian state of Delhi in order to accurately analyze the spread of COVID-19. Consequently, 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. The data is in Table 1 and is publicly available. The total population of Delhi is 190 million [24]. 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}{}_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),\\ {}_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),\\ {}_0{}^{C}D{}_t^{\alpha }C\left(t\right)=pI\left(t\right),\\ {}_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 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}{\text{lim}}\frac{1}{h^{\alpha }}{\text{Δ}}_h^{\alpha }z\left(t\right)\approx \frac{1}{h^{\alpha }}{\text{Δ}}_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 [28]

${}_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 }{\left(n\right)}^{-\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$ [29].

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

Table 3. The fundamental procedure of GMMP.

Step	Specifics
Step 1	The initial stage of this process involves obtaining the difference grid from the 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 be expressed by finite differences.
Step 3	The finite number of unknown variables can be found in the difference equation 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 Equation (31) possesses an unfamiliar variable $z\left(x_n\right)$ on both ends, Newton’s algorithm is selected to obtain the value of $z\left(x_n\right)$ through Equation (31).

5. Estimating Parameters for Nonlinear Dynamic Systems with Fractional-Order

Estimation and correction of the parameters are necessary in order to bring the numerical method given by the fractional order dynamics model closer to the actual number of infections. The parameter estimation problem is an essential and difficult task. Estimating parameters becomes increasingly challenging when the parameters are constrained and the function f exhibits a high degree of non-linearity with respect to the unknown parameters. In the following, a parameter estimation technique is proposed where a fractional order dynamical system (26) with uncertain parameters can be represented as:

$\begin{array}{l}\lambda \odot {}_a{}^{C}D{}_t^{\alpha }z\left(t\right)=f\left(t,z\left(t\right),u_1,u_2,\cdots ,u_m\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}$ (32)

The variables $z={\left(z_1,z_2,\cdots ,z_n\right)}^{\text{T}}$ and $f={\left(f_1,f_2,\cdots ,f_n\right)}^{\text{T}}$ are n-dimensional vectors. Here, $f_i\left(i=1,2,\cdots ,n\right)$ representing the parameters that are subject to uncertainty. Additionally, m represents the number of parameters.

Next, the method used in this paper to identify the parameters, namely the MH-NMSS-PSO method, will be presented. As we know, this method is based on the optimisation of multimodal functions by combining the Nelder-Mead simplex search (NMSS) method and the PSO algorithm. The focus of the search parameters is different, with NMSS focusing on “development” and PSO on “exploration”. The first of these is the choice of initial points, which are predetermined for NMSS, whereas PSO is a random set of points. Finally, they also differ in the direction and condition of the steps, with the former moving away from the less effective function values and the latter moving towards the more functional points. Thus, they each have advantages, and for the best results, we need a combination of the two methods. This method has been widely applied in different fields [30]-[32]. In the following, we illustrate the steps of the MH-HMSS-PSO algorithm.

Now, suppose that $\left(u_1,u_2,\cdots ,u_m\right)\in H$ , H can be expressed as a Cartesian product of multiple closed intervals and is bounded.

$H=\left[u_1^{\left(\text{min}\right)},u_1^{\left(\text{max}\right)}\right]\times \left[u_2^{\left(\text{min}\right)},u_2^{\left(\text{max}\right)}\right]\times \cdots \times \left[u_m^{\left(\text{min}\right)},u_m^{\left(\text{max}\right)}\right].$ (33)

The RMSREF method calculates the estimate of the unknown parameter vector $\left(u_1^{\text{*}},u_2^{\text{*}},\cdots ,u_m^{\text{*}}\right)$ within the bounded domain H using the following procedure:

$g\left(U^{\text{*}}\right)=\underset{U\in H}{\text{min}}g\left(U\right)=\underset{p\in H}{\text{min}}\sqrt{\frac{{{\displaystyle \sum }}_{j=0}^N{\left(\left(z\left(t_j\right)-z_j\right)/z_j\right)}^2}{N+1}},$ (34)

where $z_j$ are real data, and $z\left(t_j\right)$ is one of the numerical solutions of Equation (32) for given parameters $U=\left(u_1,u_2,\cdots ,u_m\right)\in H$ (such as $C\left(t\right)$ ).

Step 1: To begin, create a population of $3m+1$ individuals. This population will be used to minimize the functions $g\left(U\right)$ , which involve m unknown parameters. To initialize the process, we will generate $m+1$ vertex points $U_i=\left(u_{1,i},u_{2,i},\cdots ,u_{m,i}\right)\in H$ , $\left(i=1,2,\cdots ,m+1\right)$ . These points will form the initial structure of a simplex with m dimensions,

$U_i=\left(u_{1,i},u_{2,i},\cdots ,u_{m,i}\right)\in H,i=1,2,\cdots ,m+1,$ (35)

where

$u_{j,i}=u_j^{\left(\text{min}\right)}+\left(i-1\right)\times \left(u_j^{\left(\text{max}\right)}-u_j^{\left(\text{min}\right)}\right)/\left(m+1\right),j=1,2,\cdots ,m;i=1,2,\cdots ,m+1.$ (36)

The PSO section involves generating two particles randomly in every dimension, as outlined below:

$U_i=\left(u_{1,i},u_{2,i},\cdots ,u_{m,i}\right)\in H,i=m+2,\cdots ,3m+1,$ (37)

where

$u_{j,i}=u_j^{\left(\text{min}\right)}+\text{Rand}\times \left(u_j^{\left(\text{max}\right)}-u_j^{\left(\text{min}\right)}\right),j=1,2,\cdots ,m;i=m+2,\cdots ,3m+1.$ (38)

The randomly generated number Rand lies within the interval $\left(0,1\right)$ . The initial velocities of the component in every dimension are determined using the following procedure:

$V_{j,i}=\left(V_j^{\left(\text{max}\right)}-V_i^{\left(\text{min}\right)}\right)/L_j,j=1,2,\cdots ,m;i=m+2,\cdots ,3m+1,$ (39)

where $L_j\left(j=1,2,\cdots ,m\right)$ are chosen integers.

Step 2: Assessment and evaluation: The goal function value $g\left(U\right)$ of every particle should be assessed and ranked accordingly. The ranking should be based on the objective function value in the following manner:

$g\left(U_1\right)\le g\left(U_2\right)\le \cdots \le g\left(U_{3m+1}\right).$ (40)

Step 3: The NMSS technique requires the application of an NMSS operator to the top $m+1$ st particle. This is done by calculating $U_o$ , which represents the center of gravity of all points except $U_{m+1}$ .

$U_o=\left(u_{1,o},u_{2,o},\cdots ,u_{m,o}\right)\in H,$ (41)

where $u_{j,o}=\frac{{{\displaystyle \sum }}_{i=1}^m{u}_{j,i}}{m},j=1,2,\cdots ,m$ . The following equations can be derived:

$U_r=\left(1+\kappa \right)U_o-\kappa U_{m+1},$ (42)

here, we need to point out a fact and $\kappa$ represents the reflection coefficient ($\kappa >0$ ). It has been recommended by Nelder and Mead to set $\kappa$ equal to 1 [33].

Then, we need to consider the three cases in Table 4:

Table 4. Three possible scenarios for step 3.

Case	Specifics
Case 1	If the following situations occur: $g\left(U_1\right)\le g\left(U_r\right)\le g\left(U_m\right)$ , then $U_r$ replaces $U_{m+1}$ .
Case 2	Secondly, if $g\left(U_r\right)\le g\left(U_1\right)$ , subsequently, calculate the expanded point using the following procedure: $U_e=\gamma U_r+\left(1-\gamma \right)U_o,$ if there is $g\left(U_e\right)\le g\left(U_1\right)$ , $U_e$ can be substituted for $U_{m+1}$ , otherwise $U_r$ replaces $U_{m+1}$ , on the contrary, $P_r$ replaces $P_{m+1}$ . In accordance with the recommendation of Nelder and Mead [34], it is advised to set $\gamma$ as 2.
Case 3	$g\left(U_r\right)\ge g\left(U_m\right)$ , if $g\left(U_r\right)\le g\left(U_{m+1}\right),U_r$ replaces $U_{m+1}$ . Calculate the contracted point in the following way: $U_c=\beta U_{m+1}+\left(1-\beta \right)U_c,$ if there is $g\left(U_c\right)\le g\left(U_{m+1}\right),U_r$ can be substituted for $U_{m+1}$ , otherwise, let: $U_i=\sigma U_i+\left(1-\sigma \right)U_1$ , $i=1,2,\cdots ,m+1$ .

Nelder and Mead propose the utilization of $\beta =0.5$ and $\sigma =0.5$ for optimization purposes [34].

Step 4: Utilize the PSO operator [35] to update 2m particles having the poorest value of the goal function.

Step 5: If the following situation arises: $S_c<\epsilon$ , then the interval is stopped. On the other hand, proceed to step 2. This criterion is determined using the following formula:

$S_c={\left[{\displaystyle \sum }_{i=1}^{m+1}\frac{{\left(\overline{g}-\sqrt{g_i}\right)}^2}{m+1}\right]}^{\frac{1}{2}},$ (43)

where $\overline{g}={\displaystyle \sum }_{i=1}^{m+1}\frac{g_i^{\text{*}}}{m+1}$, and $g_i^{\text{*}}=\sqrt{g_i}=\sqrt{g_i\left(u_1,u_2,\cdots ,u_m\right)}$ .

6. Numerical Simulation

In this section, we apply a modified algorithm combining the Nelder-Mead simplex search and particle swarm optimization (MH-NMSS-PSO) to determine the parameters of the SIR model for COVID-19, taking into account fractional order dynamics. Unlike conventional approaches for solving differential equations, addressing fractional order differential equations presents greater challenges. Employing traditional solution methods would yield only approximate solutions. Various techniques exist for solving fractional order differential formulas, including the straightforward power series approach and the fuzzy Mellin variation method [36]. By utilizing the parameters derived from the MH-NMSS-PSO algorithm, our results demonstrate that the newly obtained parameters and orders exhibit a closer fit to the actual data compared to the traditional SIR model. Furthermore, we examine the COVID-19 system of fractional orders, described by the Caputo fractional order derivatives, in Equation (44). Notably, the COVID-19 profiles with different orders exhibit significant consistency with the actual data from Delhi State, providing robust evidence for subsequent long-term infection predictions.

6.1. Fractional Order COVID-19 Epidemic Model of Identical Order

In this chapter, each parameter in the proposed new SIR model (26) has its specific biological significance. Based on the real data reported in India, the cumulative number of infections is more than 1.2 million cases and casualties are more than 17,000 as of May 3, 2021 [37]. In the appropriate range, we get the lethality of the virus $d=0.14$ . Similarly, we get $\delta =0.007$ based on the mortality rate of the population reported in India [38]. It can be seen in Figure 4. We still take p to be 0.388. From [16], we find that the parameters to be estimated are the fractional order $\alpha$ , the transmission rate of the novel coronavirus $\beta$ , and the recovery rate $\gamma$ . Write the vector of unknown parameters as $U=\left(\alpha ,\beta ,\gamma \right)$ . Then, it is equivalent to $\alpha =u_1$ , $\beta =u_2$ , $\gamma =u_3$ . And each parameter is meaningful only in a specific range of values. Thus, based on these ranges, the following intervals were chosen:

$0\le u_1\le 1,0\le u_2\le 1,0\le u_3\le 1,$

and

$h_1=h_2=h_3=\mathrm{0.01.}$

Next, employing the MH-NMSS-PSO technique for parameter estimation, we utilize the actual data with identical initial values and a period of 90 days. This approach leverages the numerical solution obtained from the modified GMMP method. The resulting value of $U^{*}$ is provided below:

$\alpha =u_1=0.4492,\beta =u_2=0.700,\gamma =u_3=\mathrm{0.0115.}$

The root-mean-square error in $g\left(U^{*}\right)$ outcome is measured as 0.3921. Analysis of the simulation results illustrated in Figure 5 implies that the values of the parameter $U^{*}$ , acquired via the MH-NMSS-PSO method, fittingly match the fractional-order COVID-19 system in comparison with Figure 2. Figure 6 was employed to compare Figure 2 and Figure 5, validating the efficacy of our approach. By examining the influence of each parameter on the variation in the number of infected individuals $C\left(t\right)$ , while preserving the remaining parameters constant, we can deduce that the parameters derived from the MH-NMSS-PSO parameter estimation method are optimal.

The above results indicate that the influence of different parameters, namely $\alpha$ , $\beta$ , and $\gamma$ , on $C\left(t\right)$ is also variable. The corresponding orders and the effects of the parameters are shown in Figures 7-9. This suggests that the parameters we estimated are indeed ideal.

Figure 4. Chart of mortality in India, 2000-2022.

Figure 5. Comparison of the MH-NMSS-PSO estimated parameter-based approach with the numerical results of the number of infected persons $C\left(t\right)$ for the 2020 outbreak in Delhi. The value of $g\left(U^{*}\right)$ is 0.3921, representing the root-mean-square relative error.

Figure 6. Comparison between system (6) and system (26), where $g\left(U\right)$ is 0.3937 for the former, and $g\left(U^{*}\right)$ is 0.3921 for the latter.

Figure 7. Based on the estimated parameter method MH-NMSS-PSO, the impact of $\alpha$ on the total number of infections $C\left(t\right)$ is investigated while keeping other variables constant.

Figure 8. Based on the estimated parameter method MH-NMSS-PSO, the impact of $\beta$ on the total number of infections $C\left(t\right)$ is investigated while keeping other variables constant.

Figure 9. Based on the estimated parameter method MH-NMSS-PSO, the impact of $\gamma$ on the total number of infections $C\left(t\right)$ is investigated while keeping other variables constant.

6.2. Fractional Order COVID-19 Epidemic Models with Different Orders

We can observe that several fractional order epidemic models share identical orders. In the present study, we apply the Caputo derivative to research the expression of the fractional order COVID-19 system with varying orders as follows:

$\{\begin{array}{l}{}_0{}^{C}D{}_t^{{\alpha }_1}S\left(t\right)=-\frac{\beta S\left(t\right)I\left(t\right)}{N}-\delta S\left(t\right),\\ {}_0{}^{C}D{}_t^{{\alpha }_2}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),\\ {}_0{}^{C}D{}_t^{{\alpha }_3}C\left(t\right)=pI\left(t\right),\\ {}_0{}^{C}D{}_t^{{\alpha }_4}{R}_1\left(t\right)=\gamma I\left(t\right)-\delta R_1\left(t\right).\end{array}$ (44)

In Equation (44), the parameters $\beta$ , $\delta$ , p, $\gamma$ , and d refer to the same values as before. N is the total number of people and $\alpha \in \left(0,1\right)$ . In the case of the fractional order COVID-19 system (44), we incorporate the parameter ${\lambda }_{{\alpha }_i}\left(i=1,2,3,4\right)$ to ensure that both sides of the equation possess the same dimension. To simplify the model, as elaborated in Section 5, we need to estimate the parameters ${\alpha }_1,{\alpha }_2,{\alpha }_3,{\alpha }_4,\beta$ , and $\gamma$ . We denote the vector of unknown parameters as $U=\left({\alpha }_1,{\alpha }_2,{\alpha }_3,{\alpha }_4,\beta ,\gamma \right)$ . Then, it is equivalent to ${\alpha }_1=u_1$ , ${\alpha }_2=u_2$ , ${\alpha }_3=u_3$ , ${\alpha }_4=u_4$ , $\beta =u_5$ , $\gamma =u_6$ . Each parameter holds significance only within a specific range of values. Therefore, it is crucial to choose appropriate parameter intervals that narrow down the target values. Considering these constraints, we select the following intervals:

$\begin{array}{l}0\le u_1\le 1,0\le u_2\le 1,0\le u_3\le 1,\\ 0\le u_4\le 1,0\le u_5\le 1,0\le u_6\le 1,\end{array}$

and

$h_1=h_2=h_3=h_4=h_5=h_6=\mathrm{0.01.}$

the obtained results $U^{*}$ are shown below:

${\alpha }_1=u_1=0.5475,{\alpha }_2=u_2=0.4297,{\alpha }_3=u_3=0.8672,$

${\alpha }_4=u_4=0.6312,\beta =u_5=0.6986,\gamma =u_6=\mathrm{0.1216.}$

The root-mean-square error $g\left(U^{*}\right)$ is 0.1640. The computational findings presented in Figure 10 reveal that the parameter values $U^{*}$ procured using the MH-NMSS-PSO method result in a better fit of the fractional-order COVID-19 system to the actual data, as compared to the results presented in Figure 2. As depicted in Figure 11, we have demonstrated that COVID-19 systems with different orders exhibit better performance than those with the same order. By investigating the impact of individual parameters on the fluctuation of the number of infected individuals $C\left(t\right)$ , while keeping other parameters constant, we can conclude that the parameters obtained using the MH-NMSS-PSO parameter estimation method are ideal, as the results indicate that the effect on $C\left(t\right)$ remains relatively unchanged when the value of ${\alpha }_4$ is varied. The effect of each order and parameter is illustrated in Figure 12.

Figure 10. Comparison of the MH-NMSS-PSO estimated parameter-based approach with the numerical results of the number of exposed persons $C\left(t\right)$ for the 2020 outbreak in Delhi. The value of $g\left(U^{*}\right)$ is 0.1640, representing the root-mean-square relative error.

Figure 11. Comparing the fitting performance of system (6), system (26), and system (44) with real data from Delhi state.

Figure 12. When the remaining parameters are constant, the impact of ${\alpha }_1$ , ${\alpha }_2$ , ${\alpha }_3$ , $\beta$ , and $\gamma$ on the count of COVID-19 cases $C\left(t\right)$ will be examined.

6.3. Long-Term Projections Based on the Number of Infections in Delhi

Figure 13 shows the long-term projections for Delhi state after April 1, 2020, assuming that the parameters take values in the COVID-19 system (44) of different orders. Thus, the cumulative count of verified instances in Delhi is likely to be around 2.1 million. The convergence of the curve is likely to start between June and July 2021. This is also in better agreement with the actual data and shows that our methodology is better.

Figure 13. Long-term projections of the number of infections in Delhi.

7. Conclusions

In this paper, we study the same order and different orders 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. In order to assess the order and unknown parameters, we utilized the MH-NMSS-PSO algorithm. The validity and accuracy of the recommended approach were examined using real data from the COVID-19 outbreak in Delhi state, India, spanning from April 1, 2020, to June 30, 2020. The numerical solution depicted in Figure 5 and Figure 10 exhibits a strong alignment with the actual data, substantiating the effectiveness of both the GMMP strategy and MH-NMSS-PSO scheme for parameter estimation. A comparison between the fitting effects illustrated in Figure 11 for the integer-order COVID-19 model (6) and the fractional-order COVID-19 models (26) and (44) with the real data reveals the superior fit achieved by the fractional-order model. Figures 7-9 and Figure 12 provide insights into the impact of individual parameters on the number of exposed individuals $C\left(t\right)$ while keeping other parameters constant. Our research surpasses the findings of the referenced literature [22], as indicated by our more accurate long-term predictions and the stabilization of COVID-19 infection cases in Delhi during July. Additionally, by analyzing the root-mean-square relative error $g\left(U^{*}\right)=0.3921$ for the fractional order model (26) with the same order and $g\left(U^{*}\right)=0.1640$ for the fractional order model (44) with different orders, it can be deduced that the fractional-order COVID-19 model with varying orders better matches the real data contrasted to the model with a uniform order.

Lastly, further research is necessary to explore the applicability of this model in the management of other infectious diseases.

Funding

This work was supported by the College Student Innovation and Entrepreneurship Training Program of Sichuan University of Science and Engineering [Grant Numbers: Y2023335].

Conflicts of Interest

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

References

[1]	Mangla, S., Pathak, A.K., Arshad, M., Ghosh, D., Sahoo, P.K., Garg, V.K., et al. (2021) Impact of Environmental Indicators on the COVID-19 Pandemic in Delhi, India. Pathogens, 10, Article No. 1003. https://doi.org/10.3390/pathogens10081003
[2]	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: Advances in Experimental Medicine and Biology, Springer International Publishing, 1-22.
[3]	Eroglu, E., Esenpnar, A.A., Çiċeek, M. and Tek, S. (2020) Sir Mathematical Model and Estimating of Covid-19 Epidemic Spreading. Fresenius Environmental Bulletin, 29, 9204-9210.
[4]	Taimoor, M., Ali, S., Shah, I. and Muwanika, F.R. (2022) COVID-19 Pandemic Data Modeling in Pakistan Using Time-Series SIR. Computational and Mathematical Methods in Medicine, 2022, Article ID: 6001876. https://doi.org/10.1155/2022/6001876
[5]	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
[6]	Cooper, I., Mondal, A. and Antonopoulos, C.G. (2020) A SIR Model Assumption for the Spread of COVID-19 in Different Communities. Chaos, Solitons & Fractals, 139, Article ID: 110057. https://doi.org/10.1016/j.chaos.2020.110057
[7]	Chen, Y., Lu, P., Chang, C. and Liu, T. (2020) A Time-Dependent SIR Model for COVID-19 with Undetectable Infected Persons. IEEE Transactions on Network Science and Engineering, 7, 3279-3294. https://doi.org/10.1109/tnse.2020.3024723
[8]	Vinitsky, S.I., Gusev, A.A., Derbov, V.L., Krassovitskiy, P.M., Pen’kov, F.M. and Chuluunbaatar, G. (2021) Reduced SIR Model of COVID-19 Pandemic. Computational Mathematics and Mathematical Physics, 61, 376-387. https://doi.org/10.1134/s0965542521030155
[9]	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
[10]	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 No. 8316. https://doi.org/10.3390/app10238316
[11]	Majee, S., Adak, S., Jana, S., Mandal, M. and Kar, T.K. (2022) Complex Dynamics of a Fractional-Order SIR System in the Context of Covid-19. Journal of Applied Mathematics and Computing, 68, 4051-4074. https://doi.org/10.1007/s12190-021-01681-z
[12]	Ndaïrou, F. and Torres, D.F.M. (2021) Mathematical Analysis of a Fractional COVID-19 Model Applied to Wuhan, Spain and Portugal. Axioms, 10, Article No. 135. https://doi.org/10.3390/axioms10030135
[13]	Lu, Z., Ren, G., Chen, Y., Meng, X. and Yu, Y. (2022) A Class of Anomalous Diffusion Epidemic Models Based on CTRW and Distributed Delay. International Journal of Biomathematics, 16, Article ID: 2250130. https://doi.org/10.1142/s1793524522501303
[14]	Joshi, H., Yavuz, M., Townley, S. and Jha, B.K. (2023) Stability Analysis of a Non-Singular Fractional-Order Covid-19 Model with Nonlinear Incidence and Treatment Rate. Physica Scripta, 98, Article ID: 045216. https://doi.org/10.1088/1402-4896/acbe7a
[15]	Chen, Y., Liu, F., Yu, Q. and Li, T. (2021) Review of Fractional Epidemic Models. Applied Mathematical Modelling, 97, 281-307. https://doi.org/10.1016/j.apm.2021.03.044
[16]	Kumar, S.R. (2020) Dataset on Novel Corona Virus Disease 2019 in India. https://www.kaggle.com/sudalairajkumar/covid19-in-india
[17]	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
[18]	Liao, Z., Lan, P., Liao, Z., Zhang, Y. and Liu, S. (2020) TW-SIR: Time-Window Based SIR for COVID-19 Forecasts. Scientific Reports, 10, Article No. 22454. https://doi.org/10.1038/s41598-020-80007-8
[19]	Sarkar, K., Khajanchi, S. and Nieto, J.J. (2020) Modeling and Forecasting the COVID-19 Pandemic in India. Chaos, Solitons & Fractals, 139, Article ID: 110049. https://doi.org/10.1016/j.chaos.2020.110049
[20]	Chandramouli, C. (2011) Census of India 2011, Provisional Population Totals. Government of India, 409-413.
[21]	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
[22]	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
[23]	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
[24]	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
[25]	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
[26]	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
[27]	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
[28]	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
[29]	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
[30]	Li, T., Wang, Y. and Pan, W. (2021) Parameter Estimation for the One-Term (multiterm) Fractional-Order SEIAR Models of Norovirus Outbreak. Advances in Mathematical Physics, 2021, Article ID: 5568897. https://doi.org/10.1155/2021/5568897
[31]	Qiao, Y., Wang, X., Qi, H. and Xu, H. (2021) Numerical Simulation and Parameters Estimation of the Time Fractional Dual-Phase-Lag Heat Conduction in Femtosecond Laser Heating. International Communications in Heat and Mass Transfer, 125, Article ID: 105355. https://doi.org/10.1016/j.icheatmasstransfer.2021.105355
[32]	Fakhouri, H.N., Hudaib, A. and Sleit, A. (2020) Hybrid Particle Swarm Optimization with Sine Cosine Algorithm and Nelder-Mead Simplex for Solving Engineering Design Problems. Arabian Journal for Science and Engineering, 45, 3091-3109. https://doi.org/10.1007/s13369-019-04285-9
[33]	Nelder, J.A. and Mead, R. (1965) A Simplex Method for Function Minimization. The Computer Journal, 7, 308-313. https://doi.org/10.1093/comjnl/7.4.308
[34]	Fan, S.S., Liang, Y. and Zahara, E. (2004) Hybrid Simplex Search and Particle Swarm Optimization for the Global Optimization of Multimodal Functions. Engineering Optimization, 36, 401-418. https://doi.org/10.1080/0305215041000168521
[35]	Qazza, A., Saadeh, R. and Salah, E. (2022) Solving Fractional Partial Differential Equations via a New Scheme. AIMS Mathematics, 8, 5318-5337. https://doi.org/10.3934/math.2023267
[36]	Azhar, N. and Iqbal, S. (2022) Solution of Fuzzy Fractional Order Differential Equations by Fractional Mellin Transform Method. Journal of Computational and Applied Mathematics, 400, Article ID: 113727. https://doi.org/10.1016/j.cam.2021.113727
[37]	COVID-19 India. https://www.covid19india.org/
[38]	World Population Prospects 2022. Department of Economic and Social Affairs Population Division. https://population.un.org/wpp/