Qualitative behavior of a fractional order disease model: Optimal control and cost-effectiveness analysis

Abstract

This study focuses on formulating and analyzing the classical fractional order S E I R (susceptible-exposed-infected-recovered) model by incorporating the vaccination compartment with Holling type-II incidence rate to describe the spread of infectious diseases, by using Caputo fractional derivative to capture memory-dependent processes. The mathematical foundation of the proposed S V E I R model has been laid by using some fundamental properties of the operator, along with essential lemmas and validated by existence and boundedness. To calculate the secondary infections, we calculate the reproduction number by using the next generation matrix method and disease persistence and eradication are also given by equilibrium analysis. The sensitivity analysis is performed to identify influential parameters for the reproduction number by using normalized sensitivity index method. The novelty of this work lies in integrating a fractional-order epidemic model with a detailed optimal control framework. Applying Pontryagin’s Maximum Principle to vaccination and treatment interventions, the study identifies optimal strategies that minimize both infection burden and implementation costs. Additionally, a cost-effectiveness assessment is carried out to compare the interventions based on their epidemiological impact and resource use.

Introduction

Infectious diseases have a profound effect on both society and the economy, often leading to widespread challenges. On a social level, they put immense pressure on healthcare systems, disrupt daily life and contribute to mental health issues due to fear, isolation, and loss. Education is also affected, as school closures hinder learning, while stigma against infected individuals creates social divisions. Economically, these diseases lower productivity, increase healthcare costs, and interrupt global trade and travel. Governments must allocate significant resources to manage outbreaks, while lower-income populations face greater financial hardships. Diseases such as COVID-19, tuberculosis, HIV/AIDS, malaria, influenza, and Ebola have historically caused major disruptions, influencing public health policies and economic stability worldwide. The study of infectious disease transmission patterns and the assessment of possible management strategies need the use of mathematical models. Mathematical models are essential for studying infectious disease transmission patterns and the assessment of possible management strategies. Traditional compartmental models, formulated using integer-order differential equations, have been extensively applied to describe disease spread. Nevertheless, these models frequently fail to take into consideration the long-term relationships and memory effects of biological systems. In recent years, fractional-order differential equations (FODEs) have gained prominence in disease modeling due to their ability to incorporate memory and hereditary properties. This approach offers a more flexible and accurate framework for capturing the complex nature of disease progression, especially when historical states influence current infection dynamics. The efficacy of this approach has been demonstrated in numerous recent studies, including models with fractal-fractional operators ¹ and those incorporating time delays and media coverage. ² This research focuses on formulating and analyzing a fractional-order compartmental model for infectious diseases, investigating its qualitative properties and numerical solutions to refine our grasp of pathogen dissemination and intervention techniques.

Mathematical modeling has been used by many researchers as a potent tool to comprehend the transmission of infectious diseases and to propose efficient vaccination and treatment plans.^3,4 Vaccination has consistently served as a fundamental approach in combating infectious diseases. Immunization programs are specifically designed to prevent and control viral illnesses, including poliomyelitis, influenza, pertussis, rabies, measles, encephalitis B, and hepatitis B. Many developing nations regularly provide these immunizations. According to the report issued by the WHO in 2005, the global eradication of smallpox—the last naturally occurring case was documented in 1977—was one of the most noteworthy achievements in public health achieved by vaccination. ⁵

The implementation of vaccination initiatives is critical in managing outbreaks and protecting community health. However, a key determinant of vaccine coverage is individual behavior. Since many immunization efforts are voluntary, personal decisions about vaccination often depend on factors such as public awareness, perceived severity of the disease and the community’s existing immunity level. Vardavas et al. ⁶ discussed two primary approaches for modeling human behavior and decision-making, underscoring the need for carefully structured and adequately funded Immunization campaigns to reach community-wide immunity. Although several studies^7–9 have underscored the significance of widespread vaccination in epidemic control, detailed evaluations of their viability are still scarce.

Vaccines stimulate the immune system by presenting a harmless form of the antigen, prompting B cells to generate antibodies that specifically target the pathogen. Importantly, the immune system retains this response, providing long-term protection, often lasting years, decades, or a lifetime against future infections. Many immunization protocols require several doses over time. As an illustration, in a hepatitis B vaccination research, subjects were given three doses of the recombinant vaccine at 0, 1, and 6 months 91.2% of those examined 11 years later still showed protective levels, and 99.8% generated anti-HBs antibodies after the last treatment. ¹⁰

Researchers around the globe utilize various incidence functions based on disease characteristics. Majee et al. ¹¹ formulated a fractional order SIR model employing a saturated incidence rate βSI/1 + aS + bI to capture both infection spread and treatment dynamics. Lu et al. ¹² examined the effect of media influence on epidemic progression using an SEI model that integrates media awareness. The infection function is defined as h(S, I) = βSI/1 + mI², which incorporates a nonlinear media effect on transmission. Similarly, Barman et al. ¹³ introduced a fractional-order SEIR model featuring a saturated media influence. Their study examined the role of two critical control variables: media-induced awareness and treatment interventions. The media effect was modeled as a saturation term PSM/1 + qSM. These investigations emphasize that rising awareness due to media exposure encourages behavioral changes, such as increased social distancing. More recent studies, such as those by Hu et al.,^14,2 have further quantified the impact of media coverage on transmission dynamics within fractional-order COVID-19 models, demonstrating its critical role in shaping epidemic outcomes. Platforms like news media, online networks and social media play a pivotal role in spreading information, which helps individuals adopt preventive practices. Consequently, informed individuals adjust their interactions to minimize exposure to susceptible populations. To capture these dynamics, it is essential to revise the disease incidence functions to account for behaviorally driven modifications in transmission patterns.

Dano et al. ¹⁵ developed a SIR model based on fractional calculus, incorporating a saturated incidence rate expressed as β ^α S(t)I(t)/1 + γ ^α I(t). They applied two control strategies—vaccination and treatment with the aim of reducing disease prevalence across different levels of viral transmissibility. Furthermore, Kumar et al. ¹⁶ formulated an epidemic model that integrates both immunization and treatment constraints to determine optimal control strategies. Their study evaluated the practicality and cost-effectiveness of such control techniques, as well as the impact of the fundamental reproduction number R₀ on intervention outcomes.

This article is structured in the following manner. While recent approaches have also explored hybrid methods like Dynamics Informed Neural Networks,^17–19 our work focuses on a classical analytical and optimal control framework to derive clear public health insights. In “Review of Basic Concepts”, we lay the groundwork by introducing essential definitions and key lemmas related to fractional calculus, which form the theoretical basis for the analysis that follows. “Mathematical Model Formulation” presents the construction of a fractional order compartmental framework to capture the behavior of the infectious disease under investigation. In “Existence of Positive and Bounded solutions”, we investigate the basic qualitative characteristics of the established model, addressing aspects such as the non-negativity of solutions and the bounded nature of the system. These characteristics are essential to confirm that the model retains both biological relevance and mathematical validity.

Moving forward, “Equilibrium Analysis and Reproduction Number” is dedicated to analyzing the presence of equilibrium states and deriving the basic reproduction number, which is essential to comprehending how the disease might spread. In “Analysis of System’s Stability”, we carry out a detailed analysis of the stability of the biologically feasible equilibrium points. The next section “Sensitivity Analysis” explores the sensitivity of the basic reproduction number to different parameters, highlighting those that most significantly affect the transmission of the disease and offering valuable insights for targeted interventions. Also the section “Optimal Control Analysis” is dedicated to formulating the optimal control problem related to the model and employs Pontryagin’s Maximum Principle to derive both the necessary and adequate conditions for optimal performance. In “Cost Effectiveness Analysis”, cost-efficiency analysis of the intervention strategies is conducted to assess their affordability and effectiveness in improving public health outcomes.“Numerical Simulation” includes numerical simulations performed in MATLAB to support and validate our theoretical results, showcasing the behavior of the model under various scenarios. Lastly, in Conclusion, we discusses their relevance and provides final thoughts to conclude the study.

Review of basic concepts

Here, we present crucial terms and principles from fractional calculus that underpin the investigation.

Definition 2.1

Consider a function $g:C^k[{\mathrm{t}}_0,\infty )\to \mathbb{R}$ , then the Caputo fractional derivative of order p > 0 ²⁰ is expressed as

${}^{C}D^{\mathrm{p}}g(\mathrm{t})=\frac{1}{\mathrm{\Gamma }(k-\mathrm{p})}{\int }_{{\text{t}}_0}^{\text{t}}\frac{g^{(k)}(\xi )}{{(\mathrm{t}-\xi )}^{\text{p}-k+1}}d\xi ,$

where C ^k [t₀, ∞) represents the set of functions that are continuously differentiable up to k-times on [t₀, ∞). Here, Γ represents the Gamma function and $k\in {\mathbb{Z}}^{+}$ is satisfying k − 1 < p < k, with t > t₀.

Specifically, when p ∈ (0, 1), the Caputo derivative is given as

${}^{C}D^{\mathrm{p}}g(\mathrm{t})=\frac{1}{\mathrm{\Gamma }(1-\mathrm{p})}{\int }_{{\text{t}}_0}^{\text{t}}\frac{g^{\prime }(\xi )}{{(\mathrm{t}-\xi )}^{\text{p}}}d\xi .$

Definition 2.2

Consider a real valued continuous function g(t) defined on [t₀, ∞), the fractional integral of order p > 0 is given as follows ²¹

$I^{\text{p}}g(\mathrm{t})=\frac{1}{\mathrm{\Gamma }(\mathrm{p})}{\int }_{{\text{t}}_0}^{\text{t}}{(\mathrm{t}-\xi )}^{\text{p}-1}g(\xi )d\xi ,$

Definition 2.3

For $z\in \mathbb{C}$ , the Mittag-Leffler function with one and two parameters ²¹ is represented by

$E_{\alpha }(z)=\sum _{r=0}^{\infty }\frac{z^r}{\mathrm{\Gamma }(\alpha r+1)}$

and

$E_{\alpha ,\gamma }(z)=\sum _{r=0}^{\infty }\frac{z^r}{\mathrm{\Gamma }(\alpha r+\gamma )},$

where α, γ > 0.

Lemma 2.1

²² Suppose that h ∈ C¹[a, b] and 0 < p < 1

• (1) For a ≤ t ≤ b, $D_{\text{t}}^{\text{p}}h(\mathrm{t})\ge 0\Rightarrow h(\mathrm{t})$ is monotonic non-decreasing.

• (2) For a ≤ t ≤ b, $D_{\text{t}}^{\text{p}}h(\mathrm{t})\le 0\Rightarrow h(\mathrm{t})$ is monotonic non-increasing.

Lemma 2.2

²⁰ Let Y: [t₀, ∞) be a real valued continuous function and satisfying that

$D_{\text{t}}^{\text{p}}Y(\mathrm{t})+\lambda Y(\mathrm{t})\le \eta ,$

$Y({\mathrm{t}}_0)=Y_0,$

where 0 < p ≤ 1, λ ≠ 0 and $\lambda ,\eta \in \mathbb{R}$ . Then, we obtain

$Y(\mathrm{t})\le \left(Y_0-\frac{\eta }{\lambda }\right)E_{\text{p}}\left(-\lambda {(\mathrm{t}-{\mathrm{t}}_0)}^{\text{p}}\right)+\frac{\eta }{\lambda }.$

Mathematical model formulation

We formulate a fractional order mathematical model comprising five distinct classes to simulate the progression and regulation of an infectious disease. These compartments—Susceptible S(t), Vaccinated V(t), Exposed E(t), Infected I(t), and Recovered R(t). The compartment of Susceptible population S(t) who have not yet been infected but are at risk of contracting it. Individuals who have been immunized through vaccination and are considered protected from infection fall into the compartment V(t). The compartment of exposed individuals E(t) comprises those who have been infected but are in the latent phase and not yet infectious. The infected compartment I(t) represents symptomatic individuals who are disseminating the disease to others. Lastly, the recovered class R(t) represents the individuals who have overcome the infection and are assumed to have developed some level of immunity.

The sum of individuals in all these compartments gives the entire population, denoted by N(t) as follows

$\mathrm{N}(\mathrm{t})=\mathrm{S}(\mathrm{t})+\mathrm{V}(\mathrm{t})+\mathrm{E}(\mathrm{t})+\mathrm{I}(\mathrm{t})+\mathrm{R}(\mathrm{t}).$

To enhance the comprehension of the disease’s progression within the community, we outline some key assumptions and parameters used in constructing the model:

New individuals are added to the susceptible group S(t) at a constant rate of recruitment d. Every individual across all compartments is subject to a natural death rate g₁. A fraction of the susceptible individuals undergoes vaccination at a rate f and they shift to the recovered class R(t) at a rate n.

The transmission of the disease is governed by a Holling type II saturated incidence rate. The saturated incidence rate can capture the effect of reduced contact due to awareness, indirectly relating it to the media effects modeled explicitly in Hu et al. ¹⁴ and Barman et al. ¹³ The exposure rate for susceptible individuals is given by f₁(I) = b ^p S(t) I(t)/(1 + r ^p I(t)), where b represents the interaction rate between susceptible and infected persons and r accounts for the inhibition effects arising from awareness or behavioral changes that limit disease transmission. Individuals in the exposed class E(t) become infected at a progression rate q and the rate of mortality due to the disease among infected individuals is g₂.

A fraction of infected persons may experience natural recovery at a rate h. However, due to waning immunity, a segment of the recovered individuals lose immunity and transition back to the susceptible class at the rate m.

Based on the above assumptions and mechanisms, we construct a system of nonlinear fractional order differential equations to describe the temporal evolution of each compartment and to investigate the disease’s spread and potential control strategies

$\begin{array}{cc}\hfill {}^{C}D_{\text{t}}^{\mathrm{p}}\mathrm{S}(\mathrm{t})& ={\mathrm{d}}^{\mathrm{p}}-{\mathrm{f}}_1(\mathrm{I})-({\mathrm{g}}_{\mathrm{1}}^{\mathrm{p}}+{\mathrm{f}}^{\mathrm{p}})\mathrm{S}+{\mathrm{m}}^{\mathrm{p}}\mathrm{R},\hfill \\ \hfill {}^{C}D_{\text{t}}^{\mathrm{p}}\mathrm{V}(\mathrm{t})& ={\mathrm{f}}^{\mathrm{p}}\mathrm{S}-({\mathrm{g}}_{\mathrm{1}}^{\mathrm{p}}+{\mathrm{n}}^{\mathrm{p}})\mathrm{V},\hfill \\ \hfill {}^{C}D_{\text{t}}^{\mathrm{p}}\mathrm{E}(\mathrm{t})& ={\mathrm{f}}_1(\mathrm{I})-({\mathrm{g}}_{\mathrm{1}}^{\mathrm{p}}+{\mathrm{q}}^{\mathrm{p}})\mathrm{E},\hfill \\ \hfill {}^{C}D_{\text{t}}^{\mathrm{p}}\mathrm{I}(\mathrm{t})& ={\mathrm{q}}^{\mathrm{p}}\mathrm{E}-({\mathrm{g}}_{\mathrm{1}}^{\mathrm{p}}+{\mathrm{g}}_{\mathrm{2}}^{\mathrm{p}}+{\mathrm{h}}^{\mathrm{p}})\mathrm{I},\hfill \\ \hfill {}^{C}D_{\text{t}}^{\mathrm{p}}\mathrm{R}(\mathrm{t})& ={\mathrm{h}}^{\mathrm{p}}\mathrm{I}+{\mathrm{n}}^{\mathrm{p}}\mathrm{V}-({\mathrm{g}}_{\mathrm{1}}^{\mathrm{p}}+{\mathrm{m}}^{\mathrm{p}})\mathrm{R},\hfill \end{array}$

(1)

where, f₁(I) = b^p S I/(1 + r^p I) under the initial values:

$\mathrm{S}(0)>0,\mathrm{V}(0)\ge 0,\mathrm{E}(0)\ge 0,\mathrm{I}(0)\ge 0,\mathrm{R}(0)\ge 0.$

(2)

The schematic of the suggested model appears in Figure 1, while Table 1 outlines the biological explanations and values of the associated parameters.

Figure 1.

Flowchart of SVEIR model.

Table 1.

The ecological relevance of parameters and their associated values ¹³ .

Parameter	Ecological relevance	Values
d	Recruitment rate	100
b	Transmission rate	0.0045
g1	Natural mortality rate	0.05
g2	Infectious fatality rate	0.02
h	Intrinsic recovery rate	0.3405
q	Transition rate of exposed population	0.522
m	Transition rate of recovered population	0.2857
n	Transition rate of vaccinated population	0.024

In the following section, we examine the properties of the non-negativity and boundedness of the model (1).

Existence of positive and bounded solutions

In infectious disease modeling, it is essential to ensure that solutions remain both positive and bounded to maintain biological relevance. These conditions guarantee that population variables stay within realistic and meaningful limits.

Theorem 4.1

If S(0) > 0, V(0) ≥ 0, E(0) ≥ 0, I(0) ≥ 0, R(0) ≥ 0, then the solutions of model (1) (S(t), V(t), E(t), I(t), R(t)) are positive for p ∈ (0, 1], ∀t ≥ 0.

Proof

We have

${}^{C}D_t^{\mathrm{p}}\mathrm{S}{|}_{\text{S}=0}={\mathrm{d}}^{\mathrm{p}}+{\mathrm{m}}^{\mathrm{p}}\mathrm{R}\ge 0,$

${}^{C}D_t^{\mathrm{p}}\mathrm{V}{|}_{\text{V}=0}={\mathrm{f}}^{\mathrm{p}}\mathrm{S}\ge 0,$

${}^{C}D_t^{\mathrm{p}}\mathrm{E}{|}_{\text{E}=0}={\mathrm{f}}_1(\mathrm{I})\ge 0,$

${}^{C}D_t^{\mathrm{p}}\mathrm{I}{|}_{\text{I}=0}={\mathrm{q}}^{\mathrm{p}}\mathrm{E}\ge 0,$

${}^{C}D_t^{\mathrm{p}}\mathrm{R}{|}_{\text{R}=0}={\mathrm{h}}^{\mathrm{p}}\mathrm{I}+{\mathrm{n}}^{\mathrm{p}}\mathrm{V}\ge 0.$

Since these equations are positive for all t ≥ 0. Hence, using Lemma 2.1, S(t), V(t), E(t), I(t), R(t) ≥ 0, ∀t ≥ 0. □

Theorem 4.2

The solutions corresponding to system (1) are uniformly bounded for every t ≥ 0.

Proof

We have

${}^{C}D^{\mathrm{p}}\mathrm{N}{=}_0{}^{C}D^{\mathrm{p}}\mathrm{S}+{}^{C}D^{\mathrm{p}}\mathrm{V}+{}^{C}D^{\mathrm{p}}\mathrm{E}+{}^{C}D^{\mathrm{p}}\mathrm{I}+{}^{C}D^{\mathrm{p}}\mathrm{R}.$

From system (1), this simplifies to

${}^{C}D^{\mathrm{p}}\mathrm{N}={\mathrm{d}}^{\mathrm{p}}-{\mathrm{g}}_1^{\text{p}}\mathrm{N}-{\mathrm{g}}_2^{\text{p}}\mathrm{I},$

${}^{C}D^{\mathrm{p}}\mathrm{N}+{\mathrm{g}}_1^{\text{p}}\mathrm{N}\le {\mathrm{d}}^{\mathrm{p}}.$

Applying Lemma 2.2, we obtain

$\mathrm{N}(\mathrm{t})\le \left(\mathrm{N}(0)-\frac{{\mathrm{d}}^{\mathrm{p}}}{{\mathrm{g}}_1^{\text{p}}}\right)E_{\text{p}}\left(-{\mathrm{g}}_1^{\text{p}}{\mathrm{t}}^{\text{p}}\right)+\frac{{\mathrm{d}}^{\mathrm{p}}}{{\mathrm{g}}_1^{\text{p}}}.$

Taking t → ∞ and conclude that

$\mathrm{N}(\mathrm{t})\to \frac{{\mathrm{d}}^{\mathrm{p}}}{{\mathrm{g}}_1^{\text{p}}}.$

Furthermore, this is the system’s feasible region

${\psi }^{+}=\left\{(\mathrm{S},\mathrm{V},\mathrm{E},\mathrm{I},\mathrm{R})\in {\mathbb{R}}_{+}^5|0\le \mathrm{N}(\mathrm{t})\le {\left(\frac{\mathrm{d}}{{\mathrm{g}}_1}\right)}^{\text{p}}\right\}.$

As a result, the system (1) has bounded solutions.

□

In the following section, we analyze the system by locating its equilibrium states and evaluating the associated basic reproduction number. This analysis is essential for understanding the disease dynamics and threshold conditions for its spread.

Equilibrium analysis and reproduction number

This section focuses on identifying the steady states and calculating the basic reproduction number, which together help in understanding the system’s behavior over time. In the context of infectious disease models, they help assess whether the infection will vanish over time or remain within the population. A disease-free equilibrium represents a stable situation with no active cases, while an endemic steady state signifies the continual presence and propagation of the infection within the community.

Disease-free equilibrium (DFE)

For the disease-free equilibrium (DFE), we take E = 0, I = 0 in model (1) and obtain

${\mathcal{E}}_0=({\mathrm{S}}^0,{\mathrm{V}}^0,\mathrm{0,0},{\mathrm{R}}^0),$

where

${\mathrm{S}}^0=\frac{{\mathrm{d}}^{\text{p}}}{{\mathrm{f}}^{\text{p}}+{\mathrm{g}}_1^{\text{p}}},{\mathrm{V}}^0=\frac{{\mathrm{f}}^{\text{p}}{\mathrm{S}}^0}{{\mathrm{g}}_1^{\text{p}}+{\mathrm{n}}^{\text{p}}}\text{and}{\mathrm{R}}^0=\frac{{\mathrm{n}}^{\text{p}}{\mathrm{V}}^0}{{\mathrm{g}}_1^{\text{p}}+{\mathrm{m}}^{\text{p}}}.$

Fundamental reproduction number

The fundamental reproduction number, represented by ${\mathcal{R}}_0$ , measures the typical quantity of secondary infections produced by one infected person in a community that is fully susceptible. To ascertain ${\mathcal{R}}_0$ , we utilize the next-generation matrix technique. ¹³

According to the model (1) equations, the matrices $\mathcal{F}$ and $\mathcal{V}$ , which stand for the terms for new infections and the compartment transition processes, respectively, are as follows

$\mathcal{F}=\left[\begin{array}{c}\hfill \frac{{\mathrm{b}}^{\mathrm{p}}\mathrm{S}\mathrm{I}}{1+{\mathrm{r}}^{\mathrm{p}}\mathrm{I}}\hfill \\ \hfill 0\hfill \end{array}\right],\mathcal{V}=\left[\begin{array}{c}\hfill ({\mathrm{g}}_1^{\text{p}}+{\mathrm{q}}^{\mathrm{p}})\mathrm{E}\hfill \\ \hfill -{\mathrm{q}}^{\mathrm{p}}\mathrm{E}+({\mathrm{g}}_1^{\text{p}}+{\mathrm{g}}_2^{\text{p}}+{\mathrm{h}}^{\mathrm{p}})\mathrm{I}\hfill \end{array}\right]$

At the DFE point ${\mathcal{E}}_0=({\mathrm{S}}^0,{\mathrm{V}}^0,\mathrm{0,0},{\mathrm{R}}^0)$ , we then continue to obtain the Jacobian matrices of $\mathcal{F}$ and $\mathcal{V}$

${\mathcal{F}}_1=\left[\begin{array}{cc}\hfill 0\hfill & \hfill {\mathrm{b}}^{\mathrm{p}}{\mathrm{S}}^0\hfill \\ \hfill 0\hfill & \hfill 0\hfill \end{array}\right],{\mathcal{V}}_1=\left[\begin{array}{cc}\hfill {\mathrm{g}}_1^{\text{p}}+{\mathrm{q}}^{\mathrm{p}}\hfill & \hfill 0\hfill \\ \hfill -{\mathrm{q}}^{\mathrm{p}}\hfill & \hfill {\mathrm{g}}_1^{\text{p}}+{\mathrm{g}}_2^{\text{p}}+{\mathrm{h}}^{\mathrm{p}}\hfill \end{array}\right].$

Next, we find

${\mathcal{F}}_1{\mathcal{V}}_1^{-1}=\left[\begin{array}{cc}\hfill \frac{{\mathrm{b}}^{\mathrm{p}}{\mathrm{q}}^{\mathrm{p}}{\mathrm{S}}^0}{({\mathrm{g}}_1^{\text{p}}+{\mathrm{q}}^{\mathrm{p}})({\mathrm{g}}_1^{\text{p}}+{\mathrm{g}}_2^{\text{p}}+{\mathrm{h}}^{\text{p}})}\hfill & \hfill \frac{{\mathrm{b}}^{\mathrm{p}}{\mathrm{S}}^0}{({\mathrm{g}}_1^{\text{p}}+{\mathrm{g}}_2^{\text{p}}+{\mathrm{h}}^{\text{p}})}\hfill \\ \hfill 0\hfill & \hfill 0\hfill \end{array}\right].$

Consequently, ${\mathcal{R}}_0=\rho ({\mathcal{F}}_1{\mathcal{V}}_1^{-1})$ yields the fundamental reproduction number for model (1), where $\rho ({\mathcal{F}}_1{\mathcal{V}}_1^{-1})$ represents the spectral radius of the matrix product ${\mathcal{F}}_1{\mathcal{V}}_1^{-1}$ , which can be computed as follows

${\mathcal{R}}_0=\frac{{\mathrm{b}}^{\mathrm{p}}{\mathrm{q}}^{\mathrm{p}}{\mathrm{S}}^0}{({\mathrm{g}}_1^{\text{p}}+{\mathrm{q}}^{\mathrm{p}})({\mathrm{g}}_1^{\text{p}}+{\mathrm{g}}_2^{\text{p}}+{\mathrm{h}}^{\text{p}})}.$

(3)

Endemic equilibrium (EE) point

To ascertain the EE point of the system (1) ${\mathcal{E}}_1=({\mathrm{S}}^{*},{\mathrm{V}}^{*},{\mathrm{E}}^{*},{\mathrm{I}}^{*},{\mathrm{R}}^{*})$ , we have

${\mathrm{d}}^{\mathrm{p}}-{\mathrm{f}}_1({\mathrm{I}}^{*})-({\mathrm{g}}_{\mathrm{1}}^{\mathrm{p}}+{\mathrm{f}}^{\mathrm{p}}){\mathrm{S}}^{*}+{\mathrm{m}}^{\mathrm{p}}{\mathrm{R}}^{*}=0,$

(4)

${\mathrm{f}}^{\mathrm{p}}{\mathrm{S}}^{*}-({\mathrm{g}}_{\mathrm{1}}^{\mathrm{p}}+{\mathrm{n}}^{\mathrm{p}}){\mathrm{V}}^{*}=0,$

(5)

${\mathrm{f}}_1({\mathrm{I}}^{*})-({\mathrm{g}}_{\mathrm{1}}^{\mathrm{p}}+{\mathrm{q}}^{\mathrm{p}}){\mathrm{E}}^{*}=0,$

(6)

${\mathrm{q}}^{\mathrm{p}}{\mathrm{E}}^{*}-({\mathrm{g}}_{\mathrm{1}}^{\mathrm{p}}+{\mathrm{g}}_{\mathrm{2}}^{\mathrm{p}}+{\mathrm{h}}^{\mathrm{p}}){\mathrm{I}}^{*}=0,$

(7)

${\mathrm{h}}^{\mathrm{p}}{\mathrm{I}}^{*}+{\mathrm{n}}^{\mathrm{p}}{\mathrm{V}}^{*}-({\mathrm{g}}_{\mathrm{1}}^{\mathrm{p}}+{\mathrm{m}}^{\mathrm{p}}){\mathrm{R}}^{*}=0,$

(8)

where

${\mathrm{f}}_1({\mathrm{I}}^{*})=\frac{{\mathrm{b}}^{\mathrm{p}}{\mathrm{S}}^{*}{\mathrm{I}}^{*}}{1+{\mathrm{r}}^{\mathrm{p}}{\mathrm{I}}^{*}}.$

On solving the above equations, we obtain

${\mathrm{S}}^{*}=\frac{{\mathrm{d}}^{\text{p}}+\frac{{\mathrm{m}}^{\text{p}}{\mathrm{h}}^{\text{p}}{\mathrm{I}}^{*}}{{\mathrm{g}}_1^{\text{p}}+{\mathrm{m}}^{\text{p}}}}{M},{\mathrm{V}}^{*}=\frac{{\mathrm{f}}^{\text{p}}{\mathrm{S}}^{*}}{{\mathrm{g}}_1^{\text{p}}+{\mathrm{n}}^{\text{p}}},$

${\mathrm{E}}^{*}=\frac{{\mathrm{b}}^{\text{p}}{\mathrm{I}}^{*}{\mathrm{S}}^{*}}{(1+{\mathrm{r}}^{\text{p}}{\mathrm{I}}^{*})({\mathrm{g}}_1^{\text{p}}+{\mathrm{q}}^{\text{p}})},{\mathrm{R}}^{*}=\frac{{\mathrm{h}}^{\text{p}}{\mathrm{I}}^{*}+{\mathrm{n}}^{\text{p}}{\mathrm{V}}^{*}}{{\mathrm{g}}_1^{\text{p}}+{\mathrm{m}}^{\text{p}}},$

where

$M=\left({\mathrm{g}}_1^{\text{p}}+{\mathrm{f}}^{\mathrm{p}}+\frac{{{\mathrm{b}}^{\mathrm{p}}\mathrm{I}}^{*}}{1+{{\mathrm{r}}^{\mathrm{p}}\mathrm{I}}^{*}}\right)\left(1-\frac{{\mathrm{m}}^{\mathrm{p}}{\mathrm{n}}^{\mathrm{p}}{\mathrm{f}}^{\mathrm{p}}}{({\mathrm{g}}_1^{\text{p}}+{\mathrm{n}}^{\mathrm{p}})({\mathrm{g}}_1^{\text{p}}+{\mathrm{m}}^{\mathrm{p}})}\right).$

Substituting S*, V*, E*, I* and R* in (7), we get

$B{\mathrm{I}}^{*2}+C{\mathrm{I}}^{*}+D=0,$

where

$B=A_1{\mathrm{m}}^{\mathrm{p}}{\mathrm{h}}^{\mathrm{p}}-A_2{\mathrm{r}}^{\mathrm{p}}\left(({\mathrm{g}}_1^{\text{p}}+{\mathrm{f}}^{\mathrm{p}}){\mathrm{r}}^{\mathrm{p}}+{\mathrm{b}}^{\mathrm{p}}\right),C=A_1\left({\mathrm{d}}^{\mathrm{p}}{\mathrm{r}}^{\mathrm{p}}({\mathrm{g}}_1^{\text{p}}+{\mathrm{m}}^{\mathrm{p}})+{\mathrm{m}}^{\mathrm{p}}{\mathrm{h}}^{\mathrm{p}}\right)-A_2\left({\mathrm{b}}^{\mathrm{p}}+2{\mathrm{r}}^{\mathrm{p}}({\mathrm{g}}_1^{\text{p}}+{\mathrm{f}}^{\mathrm{p}})\right),$

$D=A_1{\mathrm{d}}^{\text{p}}({\mathrm{g}}_1^{\text{p}}+{\mathrm{m}}^{\text{p}})-A_2({\mathrm{g}}_1^{\text{p}}+{\mathrm{f}}^{\text{p}}),A_1={\mathrm{q}}^{\mathrm{p}}{\mathrm{b}}^{\mathrm{p}}({\mathrm{g}}_1^{\text{p}}+{\mathrm{n}}^{\mathrm{p}}),$

$A_2=({\mathrm{g}}_1^{\text{p}}+{\mathrm{g}}_2^{\text{p}}+{\mathrm{h}}^{\mathrm{p}})({\mathrm{g}}_1^{\text{p}}+{\mathrm{q}}^{\mathrm{p}})\left(({\mathrm{g}}_1^{\text{p}}+{\mathrm{n}}^{\mathrm{p}})({\mathrm{g}}_1^{\text{p}}+{\mathrm{m}}^{\mathrm{p}})-{\mathrm{m}}^{\mathrm{p}}{\mathrm{n}}^{\mathrm{p}}{\mathrm{f}}^{\mathrm{p}}\right).$

Analysis of System’s stability

The current section focuses on analyzing the stability behavior of model (1) around the disease-free and endemic equilibria.

Theorem 6.1

When ${\mathcal{R}}_0<1$ , the disease-free equilibrium point ${\mathcal{E}}_0$ =(S⁰, V⁰, 0, 0, R⁰) is locally asymptotically stable. ¹⁵

Proof

The Jacobian matrix J, corresponding to model (1) at the disease-free equilibrium ${\mathcal{E}}_0$ is defined as follows

$J=\left[\begin{array}{ccccc}\hfill -({\mathrm{g}}_1^{\text{p}}+{\mathrm{f}}^{\mathrm{p}})\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -{\mathrm{b}}^{\mathrm{p}}{\mathrm{S}}^0\hfill & \hfill {\mathrm{m}}^{\mathrm{p}}\hfill \\ \hfill {\mathrm{f}}^{\mathrm{p}}\hfill & \hfill -({\mathrm{g}}_1^{\text{p}}+{\mathrm{n}}^{\mathrm{p}})\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill -({\mathrm{g}}_1^{\text{p}}+{\mathrm{q}}^{\mathrm{p}})\hfill & \hfill {\mathrm{b}}^{\mathrm{p}}{\mathrm{S}}^0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill {\mathrm{q}}^{\mathrm{p}}\hfill & \hfill -({\mathrm{g}}_1^{\text{p}}+{\mathrm{g}}_2^{\text{p}}+{\mathrm{h}}^{\mathrm{p}})\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {\mathrm{n}}^{\mathrm{p}}\hfill & \hfill 0\hfill & \hfill {\mathrm{h}}^{\mathrm{p}}\hfill & \hfill -({\mathrm{g}}_1^{\text{p}}+{\mathrm{m}}^{\mathrm{p}})\hfill \end{array}\right]$

The eigenvalues of the matrix J are given by

${\lambda }_1=-{\mathrm{g}}_1^{\text{p}},$

${\lambda }_{\mathrm{2,3}}=-\left({\mathrm{g}}_1^{\text{p}}+\frac{{\mathrm{g}}_2^{\text{p}}}{2}+\frac{{\mathrm{h}}^{\mathrm{p}}}{2}+\frac{{\mathrm{q}}^{\mathrm{p}}}{2}\right)\stackrel{+}{-}\frac{P}{2},$

${\lambda }_{\mathrm{4,5}}=-\left({\mathrm{g}}_1^{\text{p}}+\frac{{\mathrm{n}}^{\mathrm{p}}}{2}+\frac{{\mathrm{f}}^{\mathrm{p}}}{2}+\frac{{\mathrm{m}}^{\mathrm{p}}}{2}\right)\stackrel{+}{-}\frac{1}{2}{({\mathrm{n}}^{2\mathrm{p}}+{\mathrm{f}}^{2\mathrm{p}}+{\mathrm{m}}^{2\mathrm{p}}-2{\mathrm{f}}^{\mathrm{p}}{\mathrm{m}}^{\mathrm{p}}-2{\mathrm{m}}^{\mathrm{p}}{\mathrm{n}}^{\mathrm{p}}-2{\mathrm{f}}^{\mathrm{p}}{\mathrm{n}}^{\mathrm{p}})}^{\frac{1}{2}},$

where

$P={\left(\frac{({\mathrm{g}}_1^{\text{p}}+{\mathrm{f}}^{\mathrm{p}})({\mathrm{g}}_2^{2\mathrm{p}}+{\mathrm{h}}^{2\mathrm{p}}+{\mathrm{q}}^{2\mathrm{p}})+2{\mathrm{q}}^{\mathrm{p}}(2{\mathrm{b}}^{\mathrm{p}}{\mathrm{d}}^{\mathrm{p}}-{\mathrm{g}}_1^{\text{p}}{\mathrm{g}}_2^{\text{p}}-{\mathrm{g}}_1^{\text{p}}{\mathrm{h}}^{\mathrm{p}})+2{\mathrm{f}}^{\mathrm{p}}({\mathrm{g}}_2^{\text{p}}{\mathrm{h}}^{\mathrm{p}}-{\mathrm{g}}_2^{\text{p}}{\mathrm{q}}^{\mathrm{p}}-{\mathrm{h}}^{\mathrm{p}}{\mathrm{q}}^{\mathrm{p}})}{{\mathrm{g}}_1^{\text{p}}+{\mathrm{f}}^{\mathrm{p}}}\right)}^{\frac{1}{2}}.$

Clearly, eigenvalues λ₁, λ₃, and λ₅ are always negative. Further, λ₂ is negative if $\left({\mathrm{g}}_1^{\text{p}}+{\mathrm{g}}_2^{\text{p}}/2+{\mathrm{h}}^{\mathrm{p}}/2+{\mathrm{q}}^{\mathrm{p}}/2\right)>P$ and λ₄ is negative if $\left({\mathrm{g}}_1^{\text{p}}+{\mathrm{n}}^{\mathrm{p}}/2+{\mathrm{f}}^{\mathrm{p}}/2+{\mathrm{m}}^{\mathrm{p}}/2\right)>1/2{({\mathrm{n}}^{2\mathrm{p}}+{\mathrm{f}}^{2\mathrm{p}}+{\mathrm{m}}^{2\mathrm{p}}-2{\mathrm{f}}^{\mathrm{p}}{\mathrm{m}}^{\mathrm{p}}-2{\mathrm{m}}^{\mathrm{p}}{\mathrm{n}}^{\mathrm{p}}-2{\mathrm{f}}^{\mathrm{p}}{\mathrm{n}}^{\mathrm{p}})}^{1/2}$ . Then, the DFE point is locally asymptotically stable when ${\mathcal{R}}_0<1$ .

□

Theorem 6.2

The endemic equilibrium point ${\mathcal{E}}_1$ =(S*, V*, E*, I*, R*) of model (1) is locally asymptotically stable, if ${\mathcal{R}}_0>1$ . ¹⁵

Proof

At the EE point ${\mathcal{E}}_1$ , the Jacobian matrix J* for model (1) is defined as follows

$J^{*}=\left[\begin{array}{ccccc}\hfill -b_1-a_{11}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -b_2\hfill & \hfill a_{15}\hfill \\ \hfill a_{21}\hfill & \hfill a_{22}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill b_1\hfill & \hfill 0\hfill & \hfill a_{33}\hfill & \hfill b_2\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill a_{43}\hfill & \hfill a_{44}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill a_{52}\hfill & \hfill 0\hfill & \hfill a_{54}\hfill & \hfill a_{55}\hfill \end{array}\right]$

where

$b_1=\frac{{\mathrm{b}}^{\text{p}}\mathrm{I}}{1+{\mathrm{r}}^{\text{p}}\mathrm{I}},b_2=\frac{{\mathrm{b}}^{\mathrm{p}}\mathrm{S}}{{(1+r^{p}I)}^2},a_{11}={\mathrm{f}}^{\text{p}}+{\mathrm{g}}_1^{\text{p}},a_{15}={\mathrm{m}}^{\text{p}},a_{21}={\mathrm{f}}^{\text{p}},a_{22}=-{\mathrm{g}}_1^{\text{p}}-{\mathrm{n}}^{\text{p}},$

$a_{33}=-{\mathrm{g}}_1^{\text{p}}-{\mathrm{q}}^{\text{p}},a_{43}={\mathrm{q}}^{\mathrm{p}},a_{44}=-{\mathrm{g}}_1^{\text{p}}-{\mathrm{g}}_2^{\text{p}}-{\mathrm{h}}^{\text{p}},a_{52}={\mathrm{n}}^{\mathrm{p}},a_{54}={\mathrm{h}}^{\mathrm{p}},a_{55}=-{\mathrm{g}}_1^{\text{p}}-{\mathrm{m}}^{\text{p}}.$

The characteristic polynomial of matrix J* is given by

${\lambda }^5+c_1{\lambda }^4+c_2{\lambda }^3+c_3{\lambda }^2+c_4\lambda +c_5=0,$

where

$\begin{array}{ll}\hfill c_1& =a_{11}-a_{22}-a_{33}-a_{44}-a_{55}+b_1,\hfill \\ \hfill c_2& =a_{22}{a}_{33}+a_{22}{a}_{44}+a_{22}{a}_{55}+a_{33}{a}_{44}+a_{33}{a}_{55}+a_{44}{a}_{55}-a_{11}{a}_{22}-a_{11}{a}_{33}-a_{11}{a}_{44}-a_{11}{a}_{55}\hfill \\ \hfill & -a_{22}{b}_1-a_{33}{b}_1-a_{44}{b}_1-a_{55}{b}_1-a_{43}{b}_2,\hfill \\ \hfill c_3& =a_{11}{a}_{22}{a}_{33}+a_{11}{a}_{22}{a}_{44}+a_{11}{a}_{22}{a}_{55}+a_{11}{a}_{33}{a}_{44}+a_{11}{a}_{33}{a}_{55}+a_{11}{a}_{44}{a}_{55}-a_{15}{a}_{21}{a}_{52}\hfill \\ \hfill & -a_{22}{a}_{33}{a}_{44}-a_{22}{a}_{33}{a}_{55}-a_{22}{a}_{44}{a}_{55}-a_{33}{a}_{44}{a}_{55}-a_{11}{a}_{43}{b}_2+a_{22}{a}_{33}{b}_1+a_{22}{a}_{43}{b}_2\hfill \\ \hfill & +a_{22}{a}_{44}{b}_1+a_{22}{a}_{55}{b}_1+a_{33}{a}_{44}{b}_1+a_{33}{a}_{55}{b}_1+a_{43}{a}_{55}{b}_2+a_{44}{a}_{55}{b}_1,\hfill \\ \hfill c_4& =a_{15}{a}_{21}{a}_{33}{a}_{52}+a_{15}{a}_{21}{a}_{44}{a}_{52}-a_{15}{a}_{43}{a}_{54}{b}_1-a_{11}{a}_{22}{a}_{33}{a}_{55}-a_{11}{a}_{22}{a}_{33}{a}_{44}-a_{11}{a}_{22}{a}_{44}{a}_{55}\hfill \\ \hfill & -a_{11}{a}_{33}{a}_{44}{a}_{55}+a_{22}{a}_{33}{a}_{44}{a}_{55}+a_{11}{a}_{22}{a}_{43}{b}_2+a_{11}{a}_{43}{a}_{55}{b}_2-a_{22}{a}_{33}{a}_{44}{b}_1-a_{22}{a}_{33}{a}_{55}{b}_1\hfill \\ \hfill & -a_{22}{a}_{43}{a}_{55}{b}_2-a_{22}{a}_{44}{a}_{55}{b}_1-a_{33}{a}_{44}{a}_{55}{b}_1,\hfill \\ \hfill c_5& =a_{11}{a}_{22}{a}_{33}{a}_{44}{a}_{55}-a_{15}{a}_{21}{a}_{33}{a}_{44}{a}_{52}-a_{11}{a}_{22}{a}_{43}{a}_{55}{b}_2+a_{15}{a}_{21}{a}_{43}{a}_{52}{b}_2+a_{15}{a}_{22}{a}_{43}{a}_{54}{b}_1\hfill \\ \hfill & +a_{22}{a}_{33}{a}_{44}{a}_{55}{b}_1.\hfill \end{array}$

Now, according to the Routh–Hurwitz stability criterion, ²³ the roots of the characteristic equation have a negative real part if

$c_1{c}_2-c_3>0,c_5(-c_4{c}_1^2+c_2{c}_1{c}_3+c_5{c}_1-c_3^2)>0.$

Hence, the endemic equilibrium will be locally asymptotically stable.

□

The next section illustrates how the basic reproduction number responds to variations in the model’s parameters through a sensitivity investigation. This aids in determining which parameters have a major effect on the dynamics of the transmission.

Sensitivity analysis

This section examines the sensitivity of the basic reproduction number ${\mathcal{R}}_0$ to changes in the model parameters, aiming to pinpoint the most influential factors driving disease spread. The sensitivity of reproduction number to a parameter p is determined using the normalized forward sensitivity index, ¹⁵ defined as

${\mathrm{\Gamma }}_{{\mathcal{R}}_0}^p=\frac{\partial {\mathcal{R}}_0}{\partial p}\times \frac{p}{{\mathcal{R}}_0}.$

If the sensitivity index is positive, then an increase in parameter p leads to a rise in ${\mathcal{R}}_0$ ; conversely, a negative index suggests that ${\mathcal{R}}_0$ decreases as p increases. This analysis is crucial for determining which parameters should be targeted for effective disease control and intervention strategies. Now, we compute sensitivity indices of parameters of the system (1) to the basic reproduction number, as presented below

${\mathrm{\Gamma }}_{{\mathcal{R}}_0}^{\text{d}}=1,{\mathrm{\Gamma }}_{{\mathcal{R}}_0}^{\text{b}}=1,{\mathrm{\Gamma }}_{{\mathcal{R}}_0}^{\text{m}}=0,$

${\mathrm{\Gamma }}_{{\mathcal{R}}_0}^{\text{n}}=0,{\mathrm{\Gamma }}_{{\mathcal{R}}_0}^{\text{r}}=0,$

${\mathrm{\Gamma }}_{{\mathcal{R}}_0}^{\text{q}}=1-\frac{\mathrm{q}}{{\mathrm{g}}_1+\mathrm{q}},{\mathrm{\Gamma }}_{{\mathcal{R}}_0}^{\text{f}}=-\frac{\mathrm{f}}{{\mathrm{g}}_1+\mathrm{f}},$

${\mathrm{\Gamma }}_{{\mathcal{R}}_0}^{\text{h}}=-\frac{\mathrm{h}}{{\mathrm{g}}_1+{\mathrm{g}}_2+\mathrm{h}},{\mathrm{\Gamma }}_{{\mathcal{R}}_0}^{{\text{g}}_2}=-\frac{{\mathrm{g}}_2}{{\mathrm{g}}_1+{\mathrm{g}}_2+\mathrm{h}},$

${\mathrm{\Gamma }}_{{\mathcal{R}}_0}^{{\text{g}}_1}=-\frac{{\mathrm{g}}_1\left(2\mathrm{f}{\mathrm{g}}_1+\mathrm{f}{\mathrm{g}}_2+\mathrm{f}\mathrm{h}+2{\mathrm{g}}_1{\mathrm{g}}_2+2{\mathrm{g}}_1\mathrm{h}+\mathrm{f}\mathrm{q}+2{\mathrm{g}}_1\mathrm{q}+{\mathrm{g}}_2\mathrm{q}+\mathrm{h}\mathrm{q}+3{\mathrm{g}}_1^2\right)}{(\mathrm{f}+{\mathrm{g}}_1)({\mathrm{g}}_1+\mathrm{q})({\mathrm{g}}_1+{\mathrm{g}}_2+\mathrm{h})}.$

Based on the parameter values outlined in Table 1, the corresponding sensitivity indices are displayed in Table 2. Parameters d and b possess positive values of sensitivity indices, demonstrating a directly proportional impact on ${\mathcal{R}}_0$ . In particular, each has a sensitivity index of 1, indicating that a 100% change in either parameter leads to an equivalent 100% change in ${\mathcal{R}}_0$ . Moreover, increasing q by 100% leads to an 83% increase in ${\mathcal{R}}_0$ , indicating a significant positive influence.

Table 2.

Sensitivity analysis of model parameters for p = 1.

Parameter	d	b	g1	g2	f	h	q
Sensitivity index	1	1	−1.3671	−0.1922	−0.9468	−0.3272	0.8333

On the other hand, parameters with negative sensitivity indices, namely g₁, g₂, f, and h exert an inverse effect on ${\mathcal{R}}_0$ . For instance, a 100% increase in f and h results in a reduction of ${\mathcal{R}}_0$ by 94% and 32%, respectively. Similarly, decreasing g₁ and g₂ by 10% leads to an increase in ${\mathcal{R}}_0$ by approximately 13% and 3.2%, respectively.

These relationships are visually summarized in Figure 2, which illustrates how each parameter impacts the basic reproduction number ${\mathcal{R}}_0$ .

Figure 2.

Forward sensitivity indices for reproduction number $({\mathcal{R}}_0)$ .

Optimal control analysis

The main purpose of analyzing the patterns of infection spread and disease progression is to develop effective strategies for controlling disease spread and ultimately eradicating the infection. To accomplish this, we use Pontryagin’s Maximum Principle (PMP), a powerful mathematical tool in optimal control theory. PMP helps determine the optimal control functions that minimize both disease prevalence and intervention costs.

In this analysis, we introduce two time-dependent variables, namely, vaccine effort u₁(t) and treatment effort u₂(t). These controls are implemented under the assumption that intervention measures should be maximized when necessary to curb disease transmission while ensuring that the associated costs remain minimal.

Our objective is to strike a balance between two competing factors:

(1) Decreasing the number of individuals infected within the community.

(2) Limiting monetary expenses and societal burdens associated with vaccination and treatment interventions.

Calculation of total cost

Here, we want to evaluate the total expenditure associated with applying control strategies, which needs to be kept low to maintain efficient disease management.

Cost associated with the disease burden

The economic impact of the infection, including loss of productivity and increased healthcare needs, is represented as

${\int }_0^T\left({\mathrm{X}}_1\mathrm{E}(\mathrm{t})+{\mathrm{X}}_2\mathrm{I}(\mathrm{t})\right)d\mathrm{t}.$

Cost of vaccination effort

The expense incurred in implementing vaccination programs, including vaccine procurement, distribution, and administration, is given by

${\int }_0^T{\mathrm{X}}_3{\mathrm{u}}_1^2(\mathrm{t})d\mathrm{t}.$

Cost of treatment effort

Providing treatment to infected individuals involves hospitalization, medical care, and pharmaceutical expenses. The total cost associated with treatment interventions is

${\int }_0^T{\mathrm{X}}_4{\mathrm{u}}_2^2(\mathrm{t})d\mathrm{t}.$

To formulate an effective control strategy, we define the objective functional as detailed below

$\mathcal{J}({\mathrm{u}}_1(\mathrm{t}),{\mathrm{u}}_2(\mathrm{t}))={\int }_0^T\left(({\mathrm{X}}_1\mathrm{E}(\mathrm{t})+{\mathrm{X}}_2\mathrm{I}(\mathrm{t}))+\frac{1}{2}\left({\mathrm{X}}_3{\mathrm{u}}_1^2(\mathrm{t})+{\mathrm{X}}_4{\mathrm{u}}_2^2(\mathrm{t})\right)\right)d\mathrm{t}.$

(9)

subjected to

$\begin{array}{ll}\hfill {}^{C}D_t^{\mathrm{p}}\mathrm{S}(\mathrm{t})& ={\mathrm{d}}^{\mathrm{p}}-{\mathrm{f}}_1(\mathrm{I})-({\mathrm{g}}_{\mathrm{1}}^{\mathrm{p}}+{\mathrm{u}}_1)\mathrm{S}+{\mathrm{m}}^{\mathrm{p}}\mathrm{R},\hfill \\ \hfill {}^{C}D_t^{\mathrm{p}}\mathrm{V}(\mathrm{t})& ={\mathrm{u}}_1\mathrm{S}-({\mathrm{g}}_{\mathrm{1}}^{\mathrm{p}}+{\mathrm{n}}^{\mathrm{p}})\mathrm{V},\hfill \\ \hfill {}^{C}D_t^{\mathrm{p}}\mathrm{E}(\mathrm{t})& ={\mathrm{f}}_1(\mathrm{I})-({\mathrm{g}}_{\mathrm{1}}^{\mathrm{p}}+{\mathrm{q}}^{\mathrm{p}})\mathrm{E},\hfill \\ \hfill {}^{C}D_t^{\mathrm{p}}\mathrm{I}(\mathrm{t})& ={\mathrm{q}}^{\mathrm{p}}\mathrm{E}-({\mathrm{g}}_{\mathrm{1}}^{\mathrm{p}}+{\mathrm{g}}_{\mathrm{2}}^{\mathrm{p}})\mathrm{I}-({\mathrm{h}}^{\mathrm{p}}+b{\mathrm{u}}_2)\mathrm{I},\hfill \\ \hfill {}^{C}D_t^{\mathrm{p}}\mathrm{R}(\mathrm{t})& =({\mathrm{h}}^{\mathrm{p}}+b{\mathrm{u}}_2)\mathrm{I}+{\mathrm{n}}^{\mathrm{p}}\mathrm{V}-({\mathrm{g}}_{\mathrm{1}}^{\mathrm{p}}+{\mathrm{m}}^{\mathrm{p}})\mathrm{R},\hfill \end{array}$

(10)

where f₁(I) = b^p S I/(1 + r^p I) under the initial values

$\mathrm{S}(0)>0,\mathrm{V}(0)\ge 0,\mathrm{E}(0)\ge 0,\mathrm{I}(0)\ge 0,\mathrm{R}(0)\ge 0.$

With the implementation of treatment control, a certain fraction of the infected individuals recover and move to the recovered compartment. Specifically, a fraction bu₂ of the infected individuals receive effective treatment and are cured, where b represents the effectiveness of the treatment strategy.

The optimal control problem aims to reduce the disease’s effects while considering the associated costs of vaccination and treatment. The objective functional for this optimal control problem is provided below

$\mathcal{J}({\mathrm{u}}_1,{\mathrm{u}}_2)=\underset{0\le {\mathrm{u}}_1,{\mathrm{u}}_2\le 1}{\min }{\int }_0^T\left(({\mathrm{X}}_1\mathrm{E}(\mathrm{t})+{\mathrm{X}}_2\mathrm{I}(\mathrm{t}))+\frac{1}{2}\left({\mathrm{X}}_3{\mathrm{u}}_1^2(\mathrm{t})+{\mathrm{X}}_4{\mathrm{u}}_2^2(\mathrm{t})\right)\right)d\mathrm{t}.$

(11)

Here, the parameters X _i (i = 1, 2, 3, 4) are positive weights assigned to appropriately scale and prioritize the individual components of the objective function. In infectious disease control,^15,24 quadratic costs are used for the control measures, where $1/2{\mathrm{X}}_3{\mathrm{u}}_1^2$ represents the cost of implementing the vaccination program and $1/2{\mathrm{X}}_4{\mathrm{u}}_2^2$ accounts for the treatment cost. The control pair ${\mathrm{u}}^{*}=({\mathrm{u}}_1^{*},{\mathrm{u}}_2^{*})$ satisfies

$\mathcal{J}({\mathrm{u}}_1^{*},{\mathrm{u}}_2^{*})=\min \left(\mathcal{J}({\mathrm{u}}_1,{\mathrm{u}}_2):({\mathrm{u}}_1,{\mathrm{u}}_2)\in \mathrm{U}\right),$

(12)

where U denotes the set of admissible controls given as

$\mathrm{U}=\{({\mathrm{u}}_1,{\mathrm{u}}_2):{\mathrm{u}}_1\text{\hspace{0.17em}and\hspace{0.17em}}{\mathrm{u}}_2\text{\hspace{0.17em}are\hspace{0.17em}measurable,\hspace{0.17em}}0\le {\mathrm{u}}_1,{\mathrm{u}}_2\le 1\}.$

(13)

The impact of implementing control functions u₁ and u₂ on the reproduction number

This section focuses on evaluating how the control functions u₁ and u₂ affect the basic reproduction number, represented by ${\mathcal{R}}_{0c}$ , for the system (10). To compute ${\mathcal{R}}_{0c}$ , we adopt the next-generation matrix approach as presented in 25. Consequently, we derive the Jacobian matrices $\mathcal{F}$ and $\mathcal{V}$ , as follows

$\mathcal{F}=\left[\begin{array}{cc}\hfill 0\hfill & \hfill \frac{{\mathrm{b}}^{\mathrm{p}}{\mathrm{d}}^{\mathrm{p}}}{{\mathrm{g}}_1^{\text{p}}+{\mathrm{f}}^{\mathrm{p}}}\hfill \\ \hfill 0\hfill & \hfill 0\hfill \end{array}\right],\mathcal{V}=\left[\begin{array}{cc}\hfill {\mathrm{g}}_1^{\text{p}}+{\mathrm{q}}^{\mathrm{p}}\hfill & \hfill 0\hfill \\ \hfill -{\mathrm{q}}^{\mathrm{p}}\hfill & \hfill {\mathrm{g}}_1^{\text{p}}+{\mathrm{g}}_2^{\text{p}}+{\mathrm{h}}^{\mathrm{p}}+b{\mathrm{u}}_2\hfill \end{array}\right].$

Therefore

${\mathcal{F}\mathcal{V}}^{-1}=\left[\begin{array}{cc}\hfill \frac{{\mathrm{b}}^{\mathrm{p}}{\mathrm{d}}^{\mathrm{p}}{\mathrm{q}}^{\mathrm{p}}}{({\mathrm{g}}_1^{\text{p}}+{\mathrm{f}}^{\mathrm{p}})({\mathrm{g}}_1^{\text{p}}+{\mathrm{q}}^{\mathrm{p}})({\mathrm{g}}_1^{\text{p}}+{\mathrm{g}}_2^{\text{p}}+{\mathrm{h}}^{\mathrm{p}}+b{\mathrm{u}}_2)}\hfill & \hfill \frac{{\mathrm{b}}^{\mathrm{p}}{\mathrm{d}}^{\mathrm{p}}{\mathrm{q}}^{\mathrm{p}}}{({\mathrm{g}}_1^{\text{p}}+{\mathrm{f}}^{\mathrm{p}})({\mathrm{g}}_1^{\text{p}}+{\mathrm{g}}_2^{\text{p}}+{\mathrm{h}}^{\mathrm{p}}+b{\mathrm{u}}_2)}\hfill \\ \hfill 0\hfill & \hfill 0\hfill \end{array}\right].$

Hence, R_0c of the system (10) is given as

${\mathcal{R}}_{0c}=\frac{{\mathrm{b}}^{\mathrm{p}}{\mathrm{d}}^{\mathrm{p}}{\mathrm{q}}^{\mathrm{p}}}{({\mathrm{g}}_1^{\text{p}}+{\mathrm{f}}^{\mathrm{p}})({\mathrm{g}}_1^{\text{p}}+{\mathrm{q}}^{\mathrm{p}})({\mathrm{g}}_1^{\text{p}}+{\mathrm{g}}_2^{\text{p}}+{\mathrm{h}}^{\mathrm{p}}+b{\mathrm{u}}_2)}.$

(14)

Figures 3(a)–(c) illustrate how ${\mathcal{R}}_{0c}$ changes when the control measures u₁ and u₂ are applied at different values of p. On observing, we see that the value of ${\mathcal{R}}_{0c}$ is least when both controls u₁ and u₂ are applied. In Figure 3(a), when no control is applied, the value of ${\mathcal{R}}_{0c}$ is 1.46,722 and when both controls are applied the value of ${\mathcal{R}}_{0c}$ is 0.31,967. Similarly, we see this variation in Figures 3(b) and (c). Also, we can observe that fractional order significantly impacts the values of ${\mathcal{R}}_{0c}$ . When full control is applied, the value of ${\mathcal{R}}_{0c}$ is decreased as the fractional order decreases.

Figure 3.

Impact of optimal control measures u₁ and u₂ on basic reproduction number R_0c at (a) p = 0.8, (b) p = 0.9 and (c) p = 1.

Optimal control problem

Here, we construct the necessary conditions for optimal control by utilizing Pontryagin’s Maximum Principle ²⁶ to the Hamiltonian function. We define u₁(t) = u₁ and u₂(t) = u₂ for simplicity’s sake. The Hamiltonian function is defined as follows

$\begin{array}{ll}\hfill \mathcal{H}& ={\mathrm{X}}_1\mathrm{E}+{\mathrm{X}}_2\mathrm{I}+\frac{1}{2}\left({\mathrm{X}}_3{\mathrm{u}}_1^2+{\mathrm{X}}_4{\mathrm{u}}_2^2\right)+{\lambda }_1{}_0{}^{C}D_t^{\mathrm{p}}\mathrm{S}(\mathrm{t})+{\lambda }_2{}_0{}^{C}D_t^{\mathrm{p}}\mathrm{V}(\mathrm{t})\hfill \\ \hfill & +{\lambda }_3{}_0{}^{C}D_t^{\mathrm{p}}\mathrm{E}(\mathrm{t})+{\lambda }_4{}_0{}^{C}D_t^{\mathrm{p}}\mathrm{I}(\mathrm{t})+{\lambda }_5{}_0{}^{C}D_t^{\mathrm{p}}\mathrm{R}(\mathrm{t}).\hfill \end{array}$

Here, λ _i (i = 1, 2, 3, 4, 5) represent the adjoint variables corresponding to the state variables of system (10).

Now, from (10) we obtain

$\begin{array}{ll}\hfill \mathcal{H}& ={\mathrm{X}}_1\mathrm{E}+{\mathrm{X}}_2\mathrm{I}+\frac{1}{2}\left({\mathrm{X}}_3{\mathrm{u}}_1^2+{\mathrm{X}}_4{\mathrm{u}}_2^2\right)\hfill \\ \hfill & +{\lambda }_1\left({\mathrm{d}}^{\mathrm{p}}-{\mathrm{f}}_1(\mathrm{I})-({\mathrm{g}}_{\mathrm{1}}^{\mathrm{p}}+{\mathrm{u}}_1)\mathrm{S}+{\mathrm{m}}^{\mathrm{p}}\mathrm{R}\right)\hfill \\ \hfill & +{\lambda }_2\left({\mathrm{u}}_1\mathrm{S}-({\mathrm{g}}_{\mathrm{1}}^{\mathrm{p}}+{\mathrm{n}}^{\mathrm{p}})\mathrm{V}\right)\hfill \\ \hfill & +{\lambda }_3\left({\mathrm{f}}_1(\mathrm{I})-({\mathrm{g}}_{\mathrm{1}}^{\mathrm{p}}+{\mathrm{q}}^{\mathrm{p}})\mathrm{E}\right)\hfill \\ \hfill & +{\lambda }_4\left({\mathrm{q}}^{\mathrm{p}}\mathrm{E}-({\mathrm{g}}_{\mathrm{1}}^{\mathrm{p}}+{\mathrm{g}}_{\mathrm{2}}^{\mathrm{p}})\mathrm{I}-({\mathrm{h}}^{\mathrm{p}}+b{\mathrm{u}}_2)\mathrm{I}\right)\hfill \\ \hfill & +{\lambda }_5\left(({\mathrm{h}}^{\mathrm{p}}+b{\mathrm{u}}_2)\mathrm{I}+{\mathrm{n}}^{\mathrm{p}}\mathrm{V}-({\mathrm{g}}_{\mathrm{1}}^{\mathrm{p}}+{\mathrm{m}}^{\mathrm{p}})\mathrm{R}\right).\hfill \end{array}$

(15)

Theorem 8.1

Assume that (S*, V*, E*, I*, R*) represents the steady-state solution of system (10) under the optimal control function $({\mathrm{u}}_1^{*},{\mathrm{u}}_2^{*})$ , then there exist corresponding adjoint variables ${\lambda }_i\in \mathbb{R}$ such that

$\begin{array}{ll}\hfill {}^{C}D_t^{\mathrm{p}}{\lambda }_1(\mathrm{t})& =({\lambda }_1-{\lambda }_3)\frac{{\mathrm{b}}^{\mathrm{p}}\mathrm{I}}{(1+{\mathrm{r}}^{\mathrm{p}}\mathrm{I})}+({\lambda }_1-{\lambda }_2){\mathrm{u}}_1+{\lambda }_1{\mathrm{g}}_1^{\text{p}},\hfill \\ \hfill {}^{C}D_t^{\mathrm{p}}{\lambda }_2(\mathrm{t})& =({\lambda }_2-{\lambda }_5){\mathrm{n}}^{\text{p}}+{\lambda }_2{\mathrm{g}}_1^{\text{p}},\hfill \\ \hfill {}^{C}D_t^{\mathrm{p}}{\lambda }_3(\mathrm{t})& =({\lambda }_3-{\lambda }_4)q^p+{\lambda }_3{\mathrm{g}}_1^{\text{p}}-{\mathrm{X}}_1,\hfill \\ \hfill {}^{C}D_t^{\mathrm{p}}{\lambda }_4(\mathrm{t})& =({\lambda }_1-{\lambda }_3)\frac{{\mathrm{b}}^{\mathrm{p}}\mathrm{S}}{{(1+{\mathrm{r}}^{\mathrm{p}}\mathrm{I})}^2}+{\lambda }_4({\mathrm{g}}_1^{\text{p}}+{\mathrm{g}}_2^{\text{p}})+({\lambda }_4-{\lambda }_5)({\mathrm{h}}^{\mathrm{p}}+b{\mathrm{u}}_2)-{\mathrm{X}}_2,\hfill \\ \hfill {}^{C}D_t^{\mathrm{p}}{\lambda }_5(\mathrm{t})& =({\lambda }_5-{\lambda }_1){\mathrm{m}}^{\text{p}}+{\lambda }_5{\mathrm{g}}_1^{\text{p}},\hfill \end{array}$

(16)

along the transversality conditions^27,13 λ _i (T) = 0, for i = 1, 2, 3, 4, 5.

Proof

Assume that ${\mathrm{u}}_1^{*},{\mathrm{u}}_2^{*}$ are the optimal control functions and S*, V*, E*, I*, R* denote the associated steady-state values of the system variables. Based on Pontryagin’s Maximum Principle, the presence of adjoint variables λ _i for i = 1, 2, 3, 4, 5 is established through the following conditions

${}^{C}D_t^{\mathrm{p}}{\lambda }_1=-\frac{\partial \mathcal{H}}{\partial \mathrm{S}},{}^{C}D_t^{\mathrm{p}}{\lambda }_2=-\frac{\partial \mathcal{H}}{\partial \mathrm{V}},{}^{C}D_t^{\mathrm{p}}{\lambda }_3=-\frac{\partial \mathcal{H}}{\partial \mathrm{E}},{}^{C}D_t^{\mathrm{p}}{\lambda }_4=-\frac{\partial \mathcal{H}}{\partial \mathrm{I}},{}^{C}D_t^{\mathrm{p}}{\lambda }_5=-\frac{\partial \mathcal{H}}{\partial \mathrm{R}},$

subject to the transversality condition and the Hamiltonian function (15).

□

Theorem 8.2

The optimal control pair $({\mathrm{u}}_1^{*},{\mathrm{u}}_2^{*})$ , which achieves the minimum value of the objective functional $\mathcal{J}$ under the system constraints (10) and within the control set U, is identified as

$\left\{\begin{array}{c}{\mathrm{u}}_1^{*}=\min \left(1,\max \left(0,\frac{({\lambda }_1-{\lambda }_2)\mathrm{S}}{{\mathrm{X}}_3}\right)\right),\hfill \\ {\mathrm{u}}_2^{*}=\min \left(1,\max \left(0,\frac{({\lambda }_4-{\lambda }_5)b\mathrm{I}}{{\mathrm{X}}_4}\right)\right).\hfill \end{array}\right.$

(17)

Proof

Now using the transversality conditions λ _i (T) = 0, for i = 1, 2, 3, 4, 5 and the optimality conditions, we have

$\begin{array}{ll}\hfill \frac{\partial \mathcal{H}}{\partial {\mathrm{u}}_1}& ={\mathrm{X}}_3{\mathrm{u}}_1-({\lambda }_1-{\lambda }_2)\mathrm{S}=0,\hfill \\ \hfill \frac{\partial \mathcal{H}}{\partial {\mathrm{u}}_2}& ={\mathrm{X}}_4{\mathrm{u}}_2-({\lambda }_4-{\lambda }_5)b\mathrm{I}=0.\hfill \end{array}$

Solving for u₁ and u₂, we obtain

${\mathrm{u}}_1^{*}=\frac{({\lambda }_1-{\lambda }_2)\mathrm{S}}{{\mathrm{X}}_3}={\tilde{\mathrm{u}}}_1,{\mathrm{u}}_2^{*}=\frac{({\lambda }_4-{\lambda }_5)b\mathrm{I}}{{\mathrm{X}}_4}={\tilde{\mathrm{u}}}_2,$

(18)

with the following conditions

${\mathrm{u}}_1^{*}=\left\{\begin{array}{cc}0,\hfill & \text{if\hspace{0.17em}}{\tilde{\mathrm{u}}}_1\le 0,\hfill \\ {\tilde{\mathrm{u}}}_1,\hfill & \text{if\hspace{0.17em}}0<{\tilde{\mathrm{u}}}_1<1,\hfill \\ 1,\hfill & \text{if\hspace{0.17em}}{\tilde{\mathrm{u}}}_1\ge 1.\hfill \end{array}\right.{\mathrm{u}}_2^{*}=\left\{\begin{array}{cc}0,\hfill & \text{if\hspace{0.17em}}{\tilde{\mathrm{u}}}_2\le 0,\hfill \\ {\tilde{\mathrm{u}}}_2,\hfill & \text{if\hspace{0.17em}}0<{\tilde{\mathrm{u}}}_2<1,\hfill \\ 1,\hfill & \text{if\hspace{0.17em}}{\tilde{\mathrm{u}}}_2\ge 1.\hfill \end{array}\right.$

□

Optimal problem model

The optimal problem model corresponding to the minimum objective function $\mathcal{J}$ is expressed as follows

$\begin{array}{ll}\hfill {}^{C}D_t^{\mathrm{p}}\mathrm{S}(\mathrm{t})& ={\mathrm{d}}^{\mathrm{p}}-{\mathrm{f}}_1(\mathrm{I})-({\mathrm{g}}_{\mathrm{1}}^{\mathrm{p}}+{\mathrm{u}}_1)\mathrm{S}+{\mathrm{m}}^{\mathrm{p}}\mathrm{R},\hfill \\ \hfill {}^{C}D_t^{\mathrm{p}}\mathrm{V}(\mathrm{t})& ={\mathrm{u}}_1\mathrm{S}-({\mathrm{g}}_{\mathrm{1}}^{\mathrm{p}}+{\mathrm{n}}^{\mathrm{p}})\mathrm{V},\hfill \\ \hfill {}^{C}D_t^{\mathrm{p}}\mathrm{E}(\mathrm{t})& ={\mathrm{f}}_1(\mathrm{I})-({\mathrm{g}}_{\mathrm{1}}^{\mathrm{p}}+{\mathrm{q}}^{\mathrm{p}})\mathrm{E},\hfill \\ \hfill {}^{C}D_t^{\mathrm{p}}\mathrm{I}(\mathrm{t})& ={\mathrm{q}}^{\mathrm{p}}\mathrm{E}-({\mathrm{g}}_{\mathrm{1}}^{\mathrm{p}}+{\mathrm{g}}_{\mathrm{2}}^{\mathrm{p}})\mathrm{I}-({\mathrm{h}}^{\mathrm{p}}+b{\mathrm{u}}_2)\mathrm{I},\hfill \\ \hfill {}^{C}D_t^{\mathrm{p}}\mathrm{R}(\mathrm{t})& =({\mathrm{h}}^{\mathrm{p}}+b{\mathrm{u}}_2)\mathrm{I}+{\mathrm{n}}^{\mathrm{p}}\mathrm{V}-({\mathrm{g}}_{\mathrm{1}}^{\mathrm{p}}+{\mathrm{m}}^{\mathrm{p}})\mathrm{R},\hfill \end{array}$

(19)

where f₁(I) = b^p S I/1 + r^p I under the initial values

$\mathrm{S}(0)>0,\mathrm{V}(0)\ge 0,\mathrm{E}(0)\ge 0,\mathrm{I}(0)\ge 0,\mathrm{R}(0)\ge 0.$

The corresponding adjoint model is given by

$\begin{array}{ll}\hfill {}^{C}D_t^{\mathrm{p}}{\lambda }_1(\mathrm{t})& =({\lambda }_1-{\lambda }_3)\frac{{\mathrm{b}}^{\mathrm{p}}\mathrm{I}}{(1+{\mathrm{r}}^{\mathrm{p}}\mathrm{I})}+({\lambda }_1-{\lambda }_2){\mathrm{u}}_1+{\lambda }_1{\mathrm{g}}_1^{\text{p}},\hfill \\ \hfill {}^{C}D_t^{\mathrm{p}}{\lambda }_2(\mathrm{t})& =({\lambda }_2-{\lambda }_5){\mathrm{n}}^{\text{p}}+{\lambda }_2{\mathrm{g}}_1^{\text{p}},\hfill \\ \hfill {}^{C}D_t^{\mathrm{p}}{\lambda }_3(\mathrm{t})& =({\lambda }_3-{\lambda }_4)q^p+{\lambda }_3{\mathrm{g}}_1^{\text{p}}-{\mathrm{X}}_1,\hfill \\ \hfill {}^{C}D_t^{\mathrm{p}}{\lambda }_4(\mathrm{t})& =({\lambda }_1-{\lambda }_3)\frac{{\mathrm{b}}^{\mathrm{p}}\mathrm{S}}{{(1+{\mathrm{r}}^{\mathrm{p}}\mathrm{I})}^2}+{\lambda }_4({\mathrm{g}}_1^{\text{p}}+{\mathrm{g}}_2^{\text{p}})+({\lambda }_4-{\lambda }_5)({\mathrm{h}}^{\mathrm{p}}+b{\mathrm{u}}_2)-{\mathrm{X}}_2,\hfill \\ \hfill {}^{C}D_t^{\mathrm{p}}{\lambda }_5(\mathrm{t})& =({\lambda }_5-{\lambda }_1){\mathrm{m}}^{\text{p}}+{\lambda }_5{\mathrm{g}}_1^{\text{p}},\hfill \end{array}$

(20)

under the transversality conditions λ _i (T) = 0, for i = 1, 2, 3, 4, 5. The control pair $({\mathrm{u}}_1^{*},{\mathrm{u}}_2^{*})$ is obtained from equation (18).

Cost effectiveness analysis

Monitoring and managing the outbreak and diffusion of transmissible diseases within a population demands considerable financial and logistical resources. Hence, identifying and implementing cost-efficient strategies for disease prevention is of paramount importance. Additionally, factors such as lifestyle and public awareness significantly influence the transmission dynamics of an infection. Consequently, countries might implement different approaches to efficiently handle and mitigate the disease. Drawing from the outcomes of the optimal control simulations, it becomes essential to identify and choose the most economically viable intervention strategy from among the available combinations.

To assess cost-effectiveness, two approaches are employed:

Average cost-effectiveness ratio (ACER)

ACER measures how cost-efficient an intervention is by comparing the overall expenses associated with the intervention relative to the number of infections it prevents. It is given by below formula

$\text{ACER}=\frac{\text{Total\hspace{0.17em}expenses\hspace{0.17em}of\hspace{0.17em}implementing\hspace{0.17em}an\hspace{0.17em}intervention}}{\text{Total\hspace{0.17em}infections\hspace{0.17em}prevented\hspace{0.17em}by\hspace{0.17em}intervention}}.$

The overall cost required for carrying out a specific intervention strategy is given by

$\frac{1}{2}{\int }_0^T\left({\mathrm{X}}_3{\mathrm{u}}_1^2(\mathrm{t})+{\mathrm{X}}_4{\mathrm{u}}_2^2(\mathrm{t})\right)d\mathrm{t}.$

(21)

Incremental cost-effectiveness ratio (ICER)

ICER is a metric used to evaluate the cost-effectiveness of two competing intervention strategies by comparing their differences in cost and health outcomes. For two strategies, A and B, the ICER is expressed as

$\text{ICER}=\frac{\text{Change\hspace{0.17em}in\hspace{0.17em}total\hspace{0.17em}cost\hspace{0.17em}between\hspace{0.17em}strategies\hspace{0.17em}}A\text{\hspace{0.17em}and\hspace{0.17em}}B}{\text{Change\hspace{0.17em}in\hspace{0.17em}total\hspace{0.17em}infections\hspace{0.17em}prevented\hspace{0.17em}by\hspace{0.17em}strategies\hspace{0.17em}}A\text{\hspace{0.17em}and\hspace{0.17em}}B}.$

(22)

The change in the number of infections avoided by each technique is captured by the denominator, while the numerator shows the difference in disease control costs between the two interventions. To demonstrate the effects of various optimal control strategies on infectious disease management, we analyzed three distinct combinations of the control variables u₁ and u₂:

• Strategy A: Implementation of only the vaccination measure (u₁ ≠ 0, u₂ = 0).

• Strategy B: Implementation of only the treatment-focused strategy (u₂ ≠ 0, u₁ = 0).

• Strategy C: A combined approach incorporating both vaccination and treatment controls (u₁ ≠ 0, u₂ ≠ 0).

To implement ICER, we perform a simulation of the system with three different intervention approaches. The strategies are ordered in increasing sequence based on the total number of infections prevented, as determined from the simulation results of the optimal control problem. In this case, Strategy B mitigates the least number of infections, followed by Strategy A, while Strategy C achieves the highest reduction in infections, as shown in Table 3. For each strategy, the ICER values calculated using (22), we obtain

$ICER(B)=\frac{4.5334-0}{47.4094-0}\approx 0.0956,ICER(A)=\frac{5.6643-4.5334}{379.2726-47.4094}\approx 0.00340,$

$ICER(C)=\frac{10.1977-5.6643}{421.9749-379.2726}\approx 0.10616.$

Table 3.

The total infections prevented, the overall cost of the control strategies and the respective ACER and ICER values.

Strategy	Infection averted	Total cost ($)	ACER values	ICER values
Strategy B	47.4094	4.5354	0.0956	0.0956
Strategy A	379.2726	5.6643	0.01493	0.00340
Strategy C	421.9749	10.1977	0.02417	0.10616

While Strategy C prevents the highest number of infections, its incremental cost per additional infection averted over Strategy A is significantly higher. Therefore, from a pure cost-effectiveness standpoint, Strategy A is the most efficient use of resources. This finding aligns with the ACER results and underscores the high value of vaccination as a primary intervention.

Next, we evaluate Strategy A against Strategy B by recalculating the ICER

$ICER(B)=\frac{4.5334-0}{47.4094-0}\approx 0.0956,ICER(A)=\frac{5.6643-4.5334}{379.2726-47.4094}\approx 0.00340.$

The comparison between the two strategies shows that Strategy A is more cost-effective compared to Strategy B, as Strategy B is less efficient and more expensive than Strategy A. Consequently, in this study, Strategy A is regarded as the most economical infection control measure. This conclusion aligns with the ACER method results obtained earlier. Furthermore, the total infections prevented, the overall cost of the control strategies and the respective ACER and ICER values are shown in Figures 4(a)–(d).

Figure 4.

Strategic comparison focusing on infections prevented, overall cost, ACER, and ICER metrics.

Numerical method

To obtain approximate solutions of the fractional–order system, we employ a fractional Euler scheme ²⁸ constructed for differential equations involving the Caputo derivative.

Consider the fractional differential equation (FDE) of Caputo type

${}^{C}D_{0,t}^{\alpha }y(t)=f(t,y(t)),t\in (0,T],T>0,$

(23)

subject to the initial conditions

$y^{(j)}(0)=y_0^{(j)},j=\mathrm{0,1},\dots ,n-1,$

(24)

where 0 < α < 1 and n = ⌈α⌉.

Solving equation (23) is often converted to solving the following equivalent Volterra integral equation

$y(t)=\sum _{j=0}^{n-1}\frac{t^j}{j!}{y}^{(j)}(0)+\frac{1}{\mathrm{\Gamma }(\alpha )}{\int }_0^t{(t-s)}^{\alpha -1}f(s,y(s))ds.$

(25)

Let the interval [0, T] be divided uniformly as

$t_j=jh,h=\frac{T}{N},j=\mathrm{0,1},\dots ,N.$

Consider the following integral term

$I_{k+1}={\int }_0^{t_{k+1}}{(t_{k+1}-s)}^{\alpha -1}g(s)ds,k=\mathrm{0,1},\dots ,N-1,$

A common strategy is to approximate this integral by

$I_{k+1}\approx {\int }_0^{t_{k+1}}{(t_{k+1}-s)}^{\alpha -1}{\tilde{g}}_{k+1}(s)ds,$

where ${\tilde{g}}_{k+1}(s)$ is an approximation of g(s) on [0, t_k+1]. Different choices of ${\tilde{g}}_{k+1}(s)$ lead to different numerical methods.

Fractional forward Euler approximation

If ${\tilde{g}}_{k+1}(s)$ is chosen piecewise constant as

${\tilde{g}}_{k+1}(s){|}_{[t_j,t_{j+1})}=g(t_j),0\le j\le k,$

then the resulting discretization yields the fractional forward Euler method

$y_{k+1}=\sum _{j=0}^{n-1}\frac{t_{k+1}^j}{j!}{y}_0^{(j)}+\frac{h^{\alpha }}{\mathrm{\Gamma }(\alpha +1)}\sum _{j=0}^k{b}_{j,k+1}f(t_j,y_j),$

(26)

where the coefficients b_j,k+1 are given by

$b_{j,k+1}={(k-j+1)}^{\alpha }-{(k-j)}^{\alpha }.$

(27)

Fractional backward Euler approximation

If instead we take

${\tilde{g}}_{k+1}(s){|}_{(t_j,t_{j+1}]}=g(t_{j+1}),0\le j\le k,$

we obtain the fractional backward Euler method

$y_{k+1}=\sum _{j=0}^{n-1}\frac{t_{k+1}^j}{j!}{y}_0^{(j)}+\frac{h^{\alpha }}{\mathrm{\Gamma }(\alpha +1)}\sum _{j=0}^k{b}_{j,k+1}f(t_{j+1},y_{j+1}),$

(28)

where the coefficients b_j,k+1 retain the definition in (27).

The stability properties of the fractional Euler scheme have been rigorously analyzed in the literature. In particular, Theorems 3.1 and 3.2 in 28 provide conditions under which the method yields stable numerical approximations.

Forward-backward sweep method for fractional-order optimal control problem

In order to obtain the optimal control for the proposed fractional-order model, we apply a forward–backward sweep method, ²⁶ which iteratively updates the state, adjoint and control variables within a unified computational framework.

We consider a fractional-order optimal control problem

$J(u)=\underset{u(\cdot )}{\min }{\int }_0^{T}L(x(t),u(t),t)dt,$

(29)

subject to

${}^{C}D_{0,t}^{\alpha }x(t)=f(x(t),u(t),t),0<t\le T,$

(30)

with initial value x(0) = x₀. Here, x(t) denotes the state variable and u(t) represents the control measure.

According to Pontryagin’s Maximum Principle, ²⁹ an optimal pair (x*, u*) must satisfy the given Hamiltonian system.

$\begin{array}{ll}\hfill {}^{C}D_{0,t}^{\alpha }x(t)& =f(x(t),u(t),t),\hfill \end{array}$

(31)

$\begin{array}{ll}\hfill {}^{C}D_{t,T}^{\alpha }\lambda (t)& =-\frac{\partial H}{\partial x}(x(t),u(t),\lambda (t),t),\hfill \end{array}$

(32)

$\begin{array}{ll}\hfill 0& =\frac{\partial H}{\partial u}(x(t),u(t),\lambda (t),t),\hfill \end{array}$

(33)

where the Hamiltonian is expressed as

$H(x,u,\lambda ,t)=L(x,u,t)+\lambda f(x,u,t),$

with x(0) = x₀ and the transversality condition λ(T) = 0.

To compute the optimal solution numerically, the fractional state and adjoint equations are discretized using the fractional Euler approximations described earlier. A general outline of the numerical algorithm is summarized below.

Numerical simulation

This section is dedicated to performing a computational simulation of the system defined by (1) and solving the corresponding optimal control problem given in (19).

Numerical simulation of model (1)

Upon computing basic reproduction number using q = 0.01, h = 0.03,405 ¹³ and considering the parameters r = 0.3 and f = 0.89 based on the literature^27,20,13 as biologically reasonable and suitable for the dynamics of the proposed model and utilizing the other parameters from Table 1, we find that ${\mathcal{R}}_0=0.76682<1$ , suggesting that the infection will ultimately be eradicated. The Jacobian matrix J at this equilibrium has eigenvalues λ₁ = −0.9279, λ₂ = −0.0500, λ₃ = −0.3718, λ₄ = −0.1582, and λ₅ = −0.0058, confirming the local asymptotic stability of ${\mathcal{E}}_0$ .

Considering the initial values S(0) = 50, V(0) = 20, E(0) = 15, I(0) = 10, and R(0) = 5, Figure 5 demonstrates that when ${\mathcal{R}}_0<1$ , the disease gradually disappears from the population. Moreover, as the fractional order decreases, the infected population approaches zero more rapidly.

Figure 5.

Dynamic behavior of different compartments for DFE at ${\mathcal{R}}_0=0.76682<1$ .

Next, by taking b = 0.075, f = 0.70 and the additional parameters are obtained from Table 1. We compute ${\mathcal{R}}_0=2.2231$ and the eigenvalues of J* are λ₁ = −0.8401, λ₂ = −0.6739, λ₃ = −0.4073, λ₄ = −0.500, and λ₅ = −0.1869, confirming that ${\mathcal{E}}_1$ is locally asymptotically stable. Moreover, Figure 6 shows that disease persistence in the population occurs when R₀ > 1. The number of infectious and exposed individuals increases when the fractional order decreases.

Figure 6.

Time series behavior of population densities for EE at ${\mathcal{R}}_0=2.2231>1$ .

Numerical simulation of the optimal control problem

This portion of the study involves numerically simulating the optimal control system (10) and the objective functional (9). The optimality system (19)–(20), along with the optimality conditions (18) and the control set (13), is solved using Sweep method with the Fractional Euler forward-backward method over [0,100] days. The state equations are computed with the forward Fractional Euler method starting from an initial guess, while the adjoint system is tackled using the backward Fractional Euler method. Additionally, we assess the influence of the control measures u₁ and u₂ in decreasing the infected population. We consider the values of weight parameters X₁ = 3, X₂ = 2, X₃ = 4, and X₄ = 9. ²⁷

Figure 7 displays the optimal control profiles over time and illustrates that both control measures follow a similar temporal pattern, applying maximum effort initially and gradually reducing their intensity as the system stabilizes. To assess the effectiveness of different control interventions, we consider the control variables u₁ = 0.2 and u₂ = 0.7 and simulated the model under multiple strategy scenarios, as illustrated in Figures 8–11. Figure 8 illustrates the behavior of state variables after implementing vaccination control. The results show a notable decline in the exposed and infected populations, while the recovered population increases rapidly. Furthermore, with a decrease in the fractional order, the count of exposed and infected individuals also lessens, underscoring the essential contribution of curbing disease propagation and mitigating its overall impact.

Figure 7.

Comparison of control profiles u₁ and u₂.

Figure 8.

Susceptible, Vaccinated, Exposed, Infected, and Recovered populations with and without vaccination control.

Figure 9.

Susceptible, Vaccinated, Exposed, Infected, and Recovered populations with and without treatment control.

Figure 10.

Susceptible, Vaccinated, Exposed, Infected, and Recovered populations without and with control using a combination of vaccination and treatment.

Figure 11.

Susceptible, Vaccinated, Exposed, Infected, and Recovered populations without and with different control strategies at p = 1.

Figure 9 demonstrates the impact of treatment control on the state variables. The findings show a gradual decrease in the exposed and infected populations over time, along with a rise in both the susceptible and recovered groups. The fractional order is key in determining the dynamics of all compartments, affecting the speed of disease spread and recovery. A lower fractional order further accelerates the decline in infections, highlighting the memory effect in fractional order systems. These findings emphasize the effectiveness of treatment strategies in limiting disease transmission, enhancing recovery rates and ultimately contributing to better epidemic control.

Figure 10 depicts the total influence of applying both vaccination and treatment interventions on the various population groups. The results demonstrate that implementing both interventions simultaneously is more effective than applying either intervention individually. The integrated approach notably decreases the exposed and infected populations while considerably boosting the recovered population. Moreover, as the fractional order diminishes, the infection burden further declines, indicating the influence of memory effects in fractional order models. The synergy between vaccination and treatment enhances disease mitigation by not only preventing new infections but also accelerating recovery, thereby reducing the overall disease burden.

Figure 11 depicts the dynamics of susceptible, vaccinated, exposed, infected, and recovered populations under different control strategies when p = 1. As shown in Figure 11(a), the susceptible population is highest in strategy B, followed by strategies A and C. Likewise, Figure 11(b) shows that strategies C and A result in a swift rise in the vaccinated population. Additionally, Figures 11(c) and (d) demonstrate that strategy C is highly effective in lowering the number of exposed and infected individuals. Lastly, Figure 11(e) shows that strategy A results in the highest number of recovered individuals compared to strategies B and C.

Conclusion

In this research, we proposed a fractional-order SVEIR model to investigate the dynamics and control of infectious disease transmission. The inclusion of fractional derivatives allows the model to account for memory effects and hereditary properties, providing a more realistic representation of the underlying biological processes compared to classical models. We derived the basic reproduction number ${\mathcal{R}}_0$ and analyzed the existence and local stability of both disease-free and endemic equilibrium points. Sensitivity analysis was performed to identify the most influential parameters affecting disease spread, offering valuable insights into effective control measures.

To minimize the disease’s impact, an optimal control problem was formulated incorporating time-dependent vaccination and treatment strategies. Using Pontryagin’s Maximum Principle, the necessary conditions for optimality were obtained. Numerical simulations were conducted to illustrate the effectiveness of these control strategies in reducing the number of infections. Additionally, a cost-effectiveness analysis was carried out to compare the control strategies. Our cost-effectiveness conclusion, which strongly favors vaccination, is supported by the dynamics observed in other complex models, such as those incorporating nucleic acid testing and protective awareness, where prevention-oriented controls often yield superior efficiency. ³⁰

The fractional order p ∈ (0, 1) determines the strength of memory included in the epidemic dynamics. For p = 1, the system behaves as a classical model without memory, whereas for p < 1, the present evolution depends on the accumulated history of infection and recovery. Biologically, this accounts for delayed immune responses, waning immunity and behavior shaped by past disease levels. As p decreases, memory effects intensify, leading to slower transitions and smoother epidemic curves. Hence, fractional models can more realistically capture history-dependent disease processes. ³¹ Overall, this study emphasizes the potential of fractional-order models in epidemiological research and underscores the importance of integrating optimal control and economic evaluation to inform disease intervention policies.