Transmission dynamics and stability of fractional order derivative model for COVID-19 epidemic with optimal control analysis

Nomenclature: SEIR: susceptible–exposed–infectious–recovered; WHO: World Health Organization; DEs: differential equations; FODEs: Fractional order differential equations; ML: Mittag–Leffler; LT: Laplace transform; DFE: Disease free equilibrium; EE: endemic equilibrium; BRN: basic reproduction number; ABM-PCS: Adam–Bashforth–Moulton predictor–corrector scheme; FBSM: forward–backward sweep method; FOCP: fractional optimal control problem.

1.   Introduction

Wuhan, China, first reported COVID-19 in December 2019. It is an infectious respiratory disease caused by a novel Coronavirus. The disease swiftly disseminated worldwide, resulting in a multitude of fatalities. Numerous factors, such as close contact, coughing, and respiratory emissions from individuals infected with COVID-19, played a significant role in its swift spread. Given its unpredictable nature, containing the disease proved to be challenging ^[1]. Researchers worldwide are dedicated to devising practical strategies aimed at combatting the COVID-19 pandemic. To unravel the intricate dynamics of the current pandemic, diverse methodologies have been employed. Vaccination is a powerful preventative method that guards against infectious diseases. In order to reduce the global incidence of individuals infected with COVID-19, vaccination is crucial, as is the stringent application of non-pharmaceutical therapies. Currently, a surge of research papers has emerged in the literature, presenting mathematical models across diverse scientific domains ^[2,3]. Fractional calculus (FC) and fractional order differential operators are essential tools when modeling a wide range of systems developed in various fields of science and technology ^[4,5,6]. The successful formulation and analysis of these models underscore their significance in advancing understanding and innovation across multiple scientific and engineering domains. In recent decades, mathematical modeling has found applications in predicting the future dynamics of diverse emerging and re-emerging infectious diseases. The results obtained from these model-driven evaluations have been extremely helpful in directing policymakers and public health administrations.

The power of fractional calculus to depict and capture phenomena that traditional integer order calculus is unable to sufficiently capture is one of its main advantages. FC enables more precise modeling and analysis of systems using fractional derivatives and integrals. Presently, an abundance of mathematical frameworks have emerged within the realm of FC, encompassing diverse models for Chikungunya ^[7], malaria ^[8], dengue ^[9], ebola ^[10,11], lassa ^[12], HIV/AIDS ^[13], tuberculosis (TB) ^[14], Hepatitis B and E ^[15,16], population dynamics ^[17], cancer epidemiology ^[18], and predator–prey population model ^[19]. At present, fractional optimal control problems (FOCPs) have surfaced as a formidable asset in the domain of epidemiological control. These strategies can ultimately contribute to disease eradication by controlling the spread within individuals. FOCPs aim to optimize a cost functional within a dynamic system governed by fractional differential equations. FOCPs extend classical optimal control problems, offering a versatile framework for systems with memory effects or non-local interactions ^[20]. Notably, researchers have employed FOCPs to model and design control strategies for diseases such as HIV/AIDS ^[21] and TB ^[22]. FOCPs have also been successfully applied to model and control a variety of other diseases (see ^{[23,24,25,26,27,28]}) and the references cited therein.

In this context, numerous researchers have recently crafted a plethora of models to explore the transmission dynamics of COVID-19. Many scholars have delved into COVID-19 modeling, utilizing data from various countries to refine their analyses. Among these, mathematical models stand out as formidable tools, effectively elucidating various facets of the COVID-19 phenomenon. As the disease emerged, nations worldwide implemented a spectrum of measures to curb its transmission, resulting in a notable reduction in its spread.

1.1.   Motivation

Motivated by the successful application of fractional order differential equations in modeling real-world problems, particularly infectious diseases' transmission dynamics, this study extends the existing epidemic model ^[29] by incorporating the well-established Caputo derivative operator to capture the memory effects present in real-world phenomena. Notably, this work builds upon our prior research conducted in ^[30,31,32].

1.2.   The novelty and significance of this study

In this study, we extend the conventional integer order framework by integrating a fractional order model for COVID-19, utilizing the Caputo fractional derivative. We utilize authentic data from Trivandrum city, the capital of Kerala, from January 1, 2022, to February 24, 2022. The aim was to construct an adequate model that enhances comprehension of the virus's dynamics. This article focuses on assessing a mathematical framework designed to control a COVID-19 epidemic. To better reflect the real-world complexity of COVID-19 transmission, our study adopts a model based on the Caputo fractional derivative, which provides enhanced accuracy and a more nuanced understanding of the effects of different vaccination approaches. The model integrates vaccination strategies and incorporates the concept of time expansion to evaluate their efficacy within a public health context. Our study aims to bridge the gap between classical epidemiological modeling and real-world complexity, offering valuable insights into disease control measures. The originality of the proposed endeavor include the following aspects:

● We investigate the interplay between partial and full vaccination and COVID-19 disease transmission, utilizing a Caputo fractional SEIR framework.

● We establish the existence and uniqueness results of the model, utilizing the context of a fixed-point approach.

● We delve into the stability of the DFE and EE points using the BRN and perform a detailed sensitivity analysis of the proposed model.

● Using a fractional ABM-PCS along with a FBSM, we solve the FOCP.

1.3.   Structure of the paper

The following is an overview of the paper's setup: The basic concepts and preliminary studies on fractional derivatives in the Caputo sense are presented in Section 2. In Section 3, we describe the Caputo fractional SEIR model for COVID-19. The existence and uniqueness results, as well as the positivity and boundedness of the model, are established in Section 4. The equilibrium points of the model and their basic reproduction number are given in Section 5. Section 6 delves into an exploration of the stability analysis and the stability criteria of the model. We conduct a sensitivity analysis for the suggested SEIR model in Section 7. We analyze the optimality of the system in Section 8 using Pontryagin's maximum principle. Section 9 utilizes MATLAB to perform numerical simulations and analyze the results. The conclusion of our work is finally included in Section 10.

2.   Preliminaries

This section presents some definitions necessary for the subsequent discussion.

Definition 1. ^[33] The Riemann–Liouville fractional derivative of order $\eta$ for a continuous function $f(t)$ is defined by

Definition 2. ^[33] The Riemann–Liouville fractional integral of order $\eta$ for the function $f:\mathbb{J}\to\mathbb{R}$ is defined as

where $n-1 < \eta\leq n$, $n\in \mathbb{N}$ and $\Gamma(.)$ is the Euler gamma function.

Definition 3. ^[33] The Caputo fractional derivative for function $f(t)\in C^{n}$ of order $\eta$ is defined as

where $n\in \mathbb{N}$ is such that $n-1 < \eta\leq n$.

Definition 4. ^[33] The Laplace transform of the Caputo fractional differential operator of order $\eta \in (n-1, n]$ is given by,

Definition 5. ^[33,34] The Mittag–Leffler function $\mathit{\boldsymbol{E}}_{\eta, \beta}$ is defined by the power series

and satisfies the equality given by ^[34]

where $\Gamma(.)$ is the Euler gamma function.

3.   Formulation of the Caputo fractional SEIR model

In this article, we extend the study in ^[29] with the Caputo fractional model. Here the total population is denoted by $N(t)$ is partitioned into four distinct groups, namely: susceptible individuals $S(t)$, exposed individuals $E(t)$, infectious individuals $I(t)$ and recovered individuals $R(t)$ at time $t$. The SEIR model of the Caputo derivative of fractional order $0 < \eta\leq 1$ is as follows:

Here, $S(0) = S_0\geq0, E(0) = E_0\geq0, I(0) = I_0\geq 0$, and $R(0) = R_0\geq0.$ In the abovementioned model (3.1), B denotes the recruitment rate and $\mu$ is the death rate in all compartments. The parameter $\beta$ represents the effective contact rate, $\delta$ is the rate of transmission from the exposed to the infected population, $P_{v}$ denotes the partially vaccinated individuals, $F_{v}$ is the completely vaccinated individuals, $\mu_{1}$ is the death rate due to infection, and $\gamma$ is rate of recovery from infection. Figure 1 represents the schematic diagram for the SEIR model.

4.   Qualitative study of the proposed model

4.1.   Uniqueness results

The unique solution to the system (3.1) is proven in this part. The model (3.1) is expressed as follows:

where

Now (3.1) is given by

Only if

By the integral representation of (4.3), which is equivalent to the model (3.1), is expressed as follows:

Consider the Banach space $\chi$ = $C\left([0, b], \mathbb{R}^4\right)$, which implies that $\left\|\upsilon\right\| = \sup\limits_{t\in J}|\upsilon(t)|$, where $|\upsilon(t)| = |S(t)|+|E(t)|+|I(t)|+|R(t)|$, and $S, E, I, R\in C([0, b]).$ Similar to ^[28], we use Krasnoselskii's fixed-point theorem to prove the existence of the solution for the model (3.1). To apply the assumptions $(\mathbb{A}1)$ and $(\mathbb{A}2)$ in ^[28], to the nonlinear functions $\kappa\in C([t\in J, \mathbb{R}^{4}])$ and $\kappa:t\in J\times\mathbb{R}^4 \rightarrow \mathbb{R}^{4}$, which is continuous and bounded.

$({\mathbb{A}1})$ There exists a constants $\psi\in C\left([0, b], \mathbb{R}^{4}_{+}\right) > 0$, such that $|\kappa(t, \upsilon)| \leq \psi(t),$ for all $(t, \upsilon)\in J\times\mathbb{R}^4$.

$({\mathbb{A}2})$ A constant $\mathcal{L}_{\kappa} > 0$ exists such that for $\forall\; t\in J$ and every ${\upsilon_{1}(t)}, {\upsilon_{2}(t)}\in C$, the following holds:

Lemma 1. ^[35] Let ${\bf T}$ be a closed, convex non empty subset of a Banach space $\chi$, and let ${\bf B}_1$ and ${\bf B}_2$ be two operators satisfying the following conditions:

(i) ${\bf B}_{1}\upsilon+{\bf B}_{2}\upsilon\in {\bf T}$, for all $\upsilon\in{\bf T},$

(ii) ${\bf B}_{1}$ is compact and continuous,

(iii) ${\bf B}_{2}$ is a contraction.

Then there $s\in {\bf T}$ exists such that $s = {\bf B}_{1}s + {\bf B}_{2}s.$

Theorem 1. Under the assumption $({\mathbb{A}1})$, the abovementioned model (3.1) has a solution provided, that $\mathcal{L}_{\kappa}\left\|\upsilon_{1}\left(t_{0}\right)-\upsilon_{2}\left(t_{0}\right)\right\| < 1.$

Proof. Consider $\sup\limits_{t \in J}|\psi(t)| = \|\psi\|$, with $\zeta\geq\|\upsilon_{0}\|+\Omega\|\psi\|,$ where $\Omega = {b^\eta}[{\Gamma({\eta+1}})^{-1}]$. We define ${\bf C}_{\zeta} = \{\upsilon \in \chi:\|\upsilon\| \leq \zeta\}$. The two operators ${\bf B}_1, {\bf B}_2$ on ${\bf C}_{\zeta}$ are given by

and

For any $\upsilon_{1}, \upsilon_{2} \in {\bf C}_{\zeta}$, it follows that

Therefore, ${\bf B}_1\upsilon_{1}+{\bf B}_2\upsilon_{2}\in {\bf C}_{\zeta}.$

We now prove that ${\bf B}_2$ is a contraction operator. For any $\upsilon_{1}, \upsilon_{2}\in{\bf C}_{\zeta}$, we obtain

which implies that $\left\|\left({\bf B}_{2}\upsilon_{1}\right)(t)- \left({\bf B}_{2}\upsilon_{2}\right)(t)\right\|\leq\mathcal{L}_{\kappa}\left\|\upsilon_{1}(t_0)-\upsilon_{2}(t_0)\right\|.$

Here, $\mathcal{L}_{\kappa}$ is a Lipschitz constant and $\mathcal{L}_{\kappa} < 1$. Hence, ${\bf B}_2$ is contraction operator.

Since $\kappa$ is continuous, this implies that the operator ${\bf B}_{1}$ must also be continuous. Moreover, for any $t\in J$ and $\upsilon\in{\bf C}_{\zeta}$, we have

which implies that ${\bf B}_{1}$ is uniformly bounded.

We now establish the compactness of the operator ${\bf B}_{1}$. Consider $\sup\limits _{(t, \upsilon)\in J\times{\bf C}_{\zeta}}|\kappa(t, \upsilon(t))| = \kappa^{*}.$

Given any $t_1, t_2\in J$ such that $t_2\geq t_1$, we have

Thus, ${\bf B}_{1}$ is equicontinuous and relatively compact on ${\bf C}_{\zeta}.$ As a result, by the Arzelá–Ascoli theorem, ${\bf B}_{1}$ is compact on ${\bf C}_{\zeta}$, since the operator has already been shown to be uniformly bounded and continuous. Hence, the model (3.1) has at least one solution on $t\in J$ according to the fixed point theorem of Krasnoselskii's.

Theorem 2. Under the assumption $({\mathbb{A}2})$, the model (3.1) has a unique solution whenever $\Omega\mathcal{L}_{\kappa} < 1,$ where $\Omega = {b^\eta}[{\Gamma({\eta+1}})^{-1}]$.

Proof. Define the operator ${\bf B}:\chi\rightarrow \chi$ defined by

Here, ${\bf B}$ is well defined, with the unique solution of model the (3.1) being its fixed point. Consider $\sup\limits_{t\in J}\left\|\kappa(t, 0)\right\| = {\bf Q_1}$ and $\kappa\geq\left\|\upsilon_0\right\|+\Omega{\bf Q_1},$ where $\Omega = {b^\eta}[{\Gamma({\eta+1}})^{-1}].$ Here, $\mathbb{P}_{\kappa} = \left\{{\upsilon \in \chi:\left\|\upsilon\right\|\leq\kappa}\right\}$ is closed and convex.

For any $\upsilon\in\mathbb{P}_{\kappa}$, we obtain,

For any ${\upsilon}_{1}, {\upsilon}_{2}\in\chi$, we get

which implies that $\left\|({\bf B}\upsilon_{1})-({\bf B}\upsilon_{2})\right\|\leq\Omega\mathcal{L}_{\kappa}\left\|\upsilon_{1}-\upsilon_{2}\right\|.$

Hence, the model (3.1) has a unique solution by Banach's contraction principle.

4.2.   Positivity and boundedness

For the fractional order model (3.1) to maintain biological validity, it is anticipated that its solution will remain both positive and bounded throughout all time intervals.

Theorem 3. For all t $\geq 0$, the state variables are non-negative, and the region $\Theta = \{(S, E, I, R)\in\mathbb{R}^{4}_{+}:0\leq S+E+I+R \leq \frac{B}{\mu}\}$ is positive invariant for the model (3.1).

Proof. From the model (3.1), we have

We also have ${ }^C D^\eta N(t) = B-\mu N(t)-\mu_1 I(t)$ and $I(t) \geq 0$ for all $t \in[0, T]$. Knowing that $\mu_1 \geq 0$, it follows that $\mu_1 I(t) \geq 0$. Therefore, $-\mu_1 I(t) \leq 0$. Thus, $B-\mu N(t)-\mu_1 I(t) \leq$ $B-\mu N(t)$. Consequently, ${ }^C D^\eta N(t) \leq B-\mu N(t)$. The inequality ${ }^C D^\eta N(t) \leq B-\mu N(t)$ is valid if and only if $I(t) \geq 0$ for all $t \in[0, T]$. Since the term $\left(-\mu_1 I(t)\right)$ is non positive, it follows that the inequality can be written as:

Using Laplace and inverse Laplace transforms (see ^[33]), we have

If we apply the equality given in ^[34], similar to (^[28], Theorem 4.1), the described model (3.1) is bounded by $\frac{B}{\mu}.$ In conclusion, all state variables are non-negative, and hence $\Theta$ is positively invariant.

5.   Equilibrium points and BRN

In this part, we evaluate the equilibrium points of the model (3.1). Here, the DFE is $\mathbb{E}^{0}$ = $(S^{0}, E^{0}, I^{0}, R^{0}) = \Big[\frac{B}{\mu+P_{v}+F_{v}}, 0, 0, \frac{F_{v}B}{\mu(\mu+P_{v}+F_{v})}\Big]$.

We apply the next-generation matrix method from ^[36] to evaluate the BRN for (3.1). The BRN, represented by $\mathcal{R}_{0}$, is determinedd from the principal eigenvalue of the matrix $\mathcal{F}\mathcal{V}^{-1}$.

Here

Hence, the BRN of the model (3.1) is

The EE of the model (3.1) is $\mathbb{E}^{*} = (S^{*}, E^{*}, I^{*}, R^{*})$, where

where $I^{*}$ is the positive solution of the equation

where $X = (\mu+\delta)(\mu+\mu_{1}+\gamma)\beta, Y = (\mu+\delta)(\mu+\mu_{1}+\gamma)(\mu+F_{v}+P_{v})-B\beta\delta$, and $Z = -P_{v}B\delta.$ Here, $Y$ can be written in terms of $\mathcal{R}_{0}$

$I^{*}$ is obtained by solving the quadratic Eq (5.1). Substituting the value of $I^{*}$ in the endemic equation $S^{*}, E^{*}, R^{*}$, the EE is determined. Clearly, the value of EE is based on the sign of B, which can found by the value of $\mathcal{R}_{0}$ in B, whether $\mathcal{R}_{0} > 1$, $\mathcal{R}_{0} = 1$ or $\mathcal{R}_{0} < 1$. Hence, by Descartes' rule of signs, only one unique positive EE point exists.

6.   Local stability analysis

Theorem 4. The DFE point $\mathbb{E}^{0}$ is locally asymptotically stable (LAS), if $\mathcal{R}_{0} < 1$ and is unstable if $\mathcal{R}_{0} > 1$.

Proof. The Jacobian matrix $\mathbb{J}$ of the system (3.1) is obtained as follows:

Here, $\Big|\mathbb{J}({\mathbb{E}^{0}})-\lambda I\Big|$ = 0 is given by

Here, $\lambda_{1} = -\mu$ which is negative, whereas the other eigenvalues can be derived from the cubic equation as follows:

where

By the Routh–Hurwitz condition, the DFE point, $\mathbb{E}^{0}$, is LAS if $s_{i}\; (i = 1, 2, 3)$ are non-negative and $s_{1}s_{2} > s_{3}$ when $\mathcal{R}_{0} < 1$ if and only if $(\mu+\delta)(\mu+\mu_{1}+\gamma)-\frac{B\beta\delta}{\mu+P_{v}+F_{v}} > 0$ and $(\mu+P_{v}+F_{v})(\mu+\delta)(\mu+\mu_{1}+\gamma)-B\beta\delta > 0$. Therefore, $\mathbb{E}^{0}$ is LAS if $\mathcal{R}_{0} < 1$ and is unstable if $\mathcal{R}_{0} > 1$.

Theorem 5. The EE point $\mathbb{E}^{*}$ of the COVID-19 model is LAS if $\mathcal{R}_{0} > 1$.

Proof. The Jacobian matrix of the system (3.1) at $\mathbb{E}^{*}$ is obtained as follows:

Here, $\Big|\mathbb{J}({\mathbb{E}^{*}})-\lambda I\Big|$ = 0 is given by

Here, obviously, one of the eigenvalue is $\lambda_{1} = -\mu$ which is negative, whereas the other eigenvalues can be derived from the cubic equation.

where

By the Routh–Hurwitz condition, the endemic equilibrium $\mathbb{E}^{*}$ is LAS if $r_{i}\; (i = 1, 2, 3)$ are non-negative and $r_{1}r_{2}-r_{3} > 0$ when $\mathcal{R}_{0} > 1$ if and only if $(\mu+P_{v}+F_{v})(\mu+\delta)-\beta\delta S^{*} < 0$ and $(\mu+P_{v}+F_{v})(\mu+\delta)(\mu+\mu_{1}+\gamma)-\beta\delta(\mu+P_{v}+F_{v})S^{*} < 0$. Therefore, $\mathbb{E}^{*}$ is LAS if $\mathcal{R}_{0} > 1.$

7.   Sensitivity analysis

Here, we analyzed the sensitivity analysis on the basis of the BRN $\mathcal{R}_{0}$ by delving into the different aspects of its first derivative across diverse parameters. This strategy determines the extent to which each parameter's value contributes to the $\mathcal{R}_{0}$. Thus, finding the appropriate actions to stop the spread of COVID-19 can be facilitated by performing a parameter sensitivity analysis. The sensitivity index of $\mathcal{R}_{0}$ to the parameters is represented mathematically as follows:

Table 1 represents the sensitivity indices of $\mathcal{R}_{0}$ with respect to the various parameters.

The sensitivity indices presented in Table 1 show that as the values of B, $\beta$, and $\delta$ rise, and while the values of the other parameters remain constant, the value of $\mathcal{R}_{0}$ rises as well. Given that the indicators are showing signs of progress, this suggests that the disease is becoming more endemic. However, the value of $\mathcal{R}_{0}$ falls when the parameter values $\mu$, $P_v$, $F_v$, $\mu_1$, and $\gamma$ are reduced while the values of the remaining parameters remain unchanged. The fact that the indices are showing negative indicators indicates a decline in the disease's endemicity.

8.   Optimal control analysis

To eradicate an illness from the population, mathematical models with control are frequently employed. In the present section, a FOCP is formulated for the COVID-19 model (3.1) by introducing time-dependent control $u(t)$. Here, $u(t)$ is used to mitigate the transmission potential arising from interactions between susceptible individuals and those who are exposed or infected. The primary objective of the suggested section can be to increase the number of recovered individuals while perhaps reducing infections in exposed and infected compartments.

In order to do this, we determine the ideal control specification using Pontgrain's maximum principle ^[37,38]. The required findings for this procedure are provided. The objective functional for the model (3.1) is

Subject to (S.to)

with the initial conditions $S(0)\geq 0, E(0)\geq 0, I(0)\geq 0$, and $R(0)\geq 0$.

In (8.1), $C_{1}, C_{2}$ represents the positive weights and $C_{3}$ represent the relative cost measure of the control variable $u(t)$. We focus on determining the control parameter $u^{*}$ such that

Here, the control set is defined as

Remark 1. Here, we assume the control function

which means that the control function is Lebesgue-integrable and essentially bounded over the given time interval. The control function $u(t)$ is assumed to be in $L^\infty([0, T_f])$, with the $L^\infty$-norm defined as:

This ensures that the control function remains bounded and does not lead to unbounded interventions. The state variables $S, E, I$, and $R$ belong to the Banach space $L^2([0, T_f])$ with the standard $L^2$-norm given by

where $X(t)$ represents any of the state variables.

8.1.   Existence of an optimal control

For the existence of an optimal control, we take the optimal system (8.2) with the associated initial conditions at $t = 0.$ Our analysis reveals that the state system possesses positive bounded solutions whenever the initial conditions are non-negative and the controls are bounded and Lebesgue measurable. This finding on existence is formalized in the following theorem.

Theorem 6. An optimal control $u^{*}\in {\bf U}$ such that

S.to the systems (8.2).

Proof. As both the control and state variables are non-negative, the objective function (8.1) satisfies the convexity condition. Here, the control set is convex and closed, which shows that the optimal system satisfies the property of boundedness, thereby confirming the essential compactness necessary for the existence of optimal control. The integrand in the objective functional (8.1) is convex on ${\bf U}$, which shows that an optimal control $u^{*}\in {\bf U}$ exists for the system (8.2).

8.2.   Characterization of the optimal control functions

To establish the necessary conditions for optimality, we employ Pontryagin's maximum principle. Accordingly, consider the Hamiltonian $\mathbb{H}$ of the problem (8.1) as follows:

where $\lambda_{j}$, ($j = 1, 2\ldots 4$) represents the co-state variables, which are determined by solving the following equations:

8.3.   Optimality conditions

Theorem 7. ^[39] Let $u^{*}$ be an optimal control of the problems (8.1) and (8.2) with the corresponding states $S^{*}, E^{*}, I^{*}$, and $R^{*}$. In this case, $\lambda_{1}, \lambda_{2}, \lambda_{3}$, and $\lambda_{4}$ exists, satisfying the following:

With the transversality condition $\lambda_{1}(T_{f}) = \lambda_{4}(T_{f}) = 0$ and $\lambda_{2}(T_{f}) = -C_{1}, \lambda_{3}(T_{f}) = -C_{2}$. Furthermore, the optimal control $u^{*}$ is given by

Proof. The Hamiltonian function for our prescribed optimal control problem is as follows:

For $t\in [0, T_{f}]$, the adjoint equations and transversality conditions can be obtained by using the principle in ^[38] such that

The transversality conditions are $\lambda_{1}(T_{f}) = \lambda_{4}(T_{f}) = 0$ and $\lambda_{2}(T_{f}) = -C_{1}, \lambda_{3}(T_{f}) = -C_{2}$.

The condition for deriving the optimal control $u^*(t)$ is given by

The optimal control $u^{*}$ is derived by minimizing the cost functional ${\bf J}$, resulting in the following expression:

This completes the proof.

9.   Numerical results and discussion

In this section of the manuscript, we provide the numerical simulation results for the abovementioned COVID-19 model (3.1). The objective of Subsection 9.1 is to solve the system (3.1) by applying the ABM-PCS to the discretized form of the system. Furthermore, a description of FBSM employing the ABM-PCS to solve the optimality system (8.2) is provided in Subsection 9.2. These techniques provide precise numerical results over an extended period of time. For the simulations, we utilise MATLAB with the initial conditions and settings as listed below.

9.1.   ABM-PCS for the SEIR Model

For the initial values of the fractional order, the most widely used numerical method is the ABM-PCS method (see ^[40,41]), and it provides precise solutions over an extended period of time.

Consider the following general form of a differential equation:

where $A_{0}^b\in\mathbb{R}$, which is equal to the well-known Volterra integral equation

To compute the numerical solution of a nonlinear fractional order model (3.1), we apply the ABM-PCS. Now consider that the time interval $[0, T]$ can be divided into $N$ equal parts with a size of $h = \frac{T}{N}$, and each node can be represented as $t_n = nh, \; n = 0, 1, 2, \ldots, {N}$.

The corrector values are defined as:

The predictor values are given by

where

and

The initial values are $S(0) = S_{0} = 16.61989X10^5, E(0) = E_{0} = 0.16797X10^5$, $I(0) = I_{0} = 481,$ and $R(0) = R_{0} = 487$ along with the model parameters reported in ^[29] as follows.

Numerical simulations of the system (3.1) were performed for the integer and for the fractional orders $\eta = 1, 0.9, 0.8$, and $0.7$.

Figure 2 provides a visual representation of the dynamics in each compartment. This graphical representation emphasizes that the fractional derivative model provides complete information on disease dynamics. In Figure 2(a), we display the temporal dynamics of the susceptible population at various values of $\eta = 1, 0.9, 0.8$, and $0.7$. We can observe and confirm that the proportion of the susceptible population decreases with varying fractional orders of $\eta$. The decrease in the susceptible population is attributed to the vaccination rates $P_{v}$, $F_{v}$, and the natural infection rate, leading to a shift of individuals to the exposed class.

Figure 2(b) represent the exposed population with respect to time, which varies for various values of $\eta$. It can be observed that an exposed population is increased over time and slowly declines due to the values of the parameters $P_{v}$, $F_{v}$, and the fractional order $\eta$. The interplay between partial and full vaccination significantly affects the rate at which individuals transition from the exposed class to the infected or recovered compartments. Figure 2(c) shows the dynamic behavior of the infected population versus time. The results confirm that for smaller fractional orders $\eta$, the infection curve is prolonged, indicating a slower rate of decline. This observation aligns with theoretical expectations, as the fractional order introduces memory effects that delay recovery and prolong the infection's duration. The presence of the vaccination parameters $P_{v}$, $F_{v}$ effectively reduces the peak infection levels, demonstrating the critical role of vaccination strategies in mitigating disease spread. Figure 2(d) depicts the dynamic behavior of the recovered population versus time. We can confirm that the proportion of the recovered population increases significantly due to the parameters $P_{v}$, $F_{v}$, and the different values of the fractional order $\eta$. The impact of both partial and full vaccination is evident, as they contribute to a higher recovery rate and a reduction in the infected population. The fractional order further influences the rate of increase in recoveries, where lower values of $\eta$ lead to a more gradual rise in recovered cases, reflecting the long-term memory effects.

9.2.   FBSM based on ABM-PCS

This section presents a numerical simulation of the FOCP formulated in Section 8. When solving optimality systems, FBSM is a useful iterative technique. We extend FBSM to solve our FOCP using the following approach, which is based on the description of ABM-PCS mentioned above. The algorithm begins with a first estimation of the control variable (see ^[34]).

The numerical description of the FOCP for (8.2), as detailed in Subsection (9.2), was used to solve the optimality system, utilizing the initial state values $(S(0), E(0), I(0), R(0))$ = $(16.61989X10^5, 0.16797X10^5,481,487)$. The weight constant values were chosen as $C_{1} = 10, C_{2} = 200$, and $C_{3} = 0.001$, along with the parameter values in Table 2.

Control measures lead to an increase in the susceptible population compared with the lack of control (see Figure 3(a)). In Figure 3(b), (c) we observed that the impact of the control measure on the populations of exposed and infected people. The implementation of control parameters significantly reduces the number of exposed and infected individuals. As a result, this decline leads to an increase in the recovered population, as illustrated in Figure 3(d). The results indicate that increasing $u(t)$, particularly under optimal control strategies, effectively decreases the susceptible, exposed, and infected populations while enhancing the recovery rate. In conclusion, these results suggest that the combined implementation of the proposed control measures offers a superior and more impactful strategy. This approach effectively minimizes infection rates within the community, thereby safeguarding against future outbreaks.

Remark 2. This study extends the model in ^[29] by introducing control variables and evaluating their impact on reducing the spread of exposed and infected populations and increasing recovery. Although our study does not directly validate the model with real-world data, it lays the groundwork for designing optimal control strategies that can later be applied to refined and validated models. We agree that real-world applicability is essential, and we plan to combine the insights from ^[29] and ^[42] in future work to validate the control strategies with data from specific regions and address vaccine unavailability scenarios as well.

10.   Conclusions

In this study, a Caputo fractional derivative model for COVID-19 is presented in the SEIR context. As part of our study, we analyze the solution positivity as well as the boundedness of the system. Our stability analysis shows that when the BRN $\mathcal{R}_{0} < 1$, the system exhibits LAS around the DFE point ($\mathbb{E}^{0}$). On the other hand, the system stabilizes asymptotically around the EE point ($\mathbb{E}^{*}$) when $\mathcal{R}_{0} > 1$. A sensitivity analysis is conducted using $\mathcal{R}_{0}$ for this model. Pontryagin's maximum principle was used to derive the necessary optimality conditions for the FOCP. Using the ABM-PCS, we conducted numerical simulations, and developed the FBSM with ABM-PCS to solve the optimized system numerically. Our findings indicate that a larger proportion of completely vaccinated people correlates with a larger proportion of the population recovering. This suggests that comprehensive vaccination substantially mitigates the impact of infection, swiftly transitioning individuals to the recovered class. The COVID-19 virus can be effectively managed and inhibited by implementing control variables. In our study, we demonstrated that it is the most effective technique for eradicating diseases. In comparison with the classical model, our simulations indicate that the non-integer model better captures the dynamic behavior of the COVID-19 epidemiology model.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This work was supported by Chungbuk National University BK21's program (2024). This research was also supported by the Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (RS-2020-NR049604).

Conflict of interest

The authors declare there is no conflict of interest.