Dynamics and Control Study of the New H1N1 Influenza with Two Roots of Infection: Impact of Optimal Vaccination and Treatment

1. Introduction

The H1N1 virus produces symptoms similar to those caused by other influenza viruses [1]. Influenza viruses are broadly classified into four types: A, B, C, and D. Among these, types A and B are primarily responsible for human flu infections, while types C and D cause only mild illnesses and are typically restricted to animals. Seasonal influenza affects approximately one in six individuals each year, leading to 3 to 5 million cases of severe illness and an estimated 290,000 to 650,000 respiratory-related deaths globally [2,3]. The H1N1 pandemic virus spreads through person-to-person contact, airborne particles, coughing, sneezing, and contaminated surfaces (fomites) [4]. Most symptoms appear rapidly and may include: Fever, however sometimes, sore muscles, sweats and chills, Cough, throat pain, runny or congested nose. Vomiting and feeling ill in one’s stomach are prevalent in youngsters but less so in adults [1]. Seasonal influenza spreads rapidly in crowded environments such as mass gatherings (World Health Organization, 2018) [5]. To reduce transmission, individuals should practice regular hand hygiene and cover their mouth and nose when coughing or sneezing. The World Health Organisation (WHO) proclaimed the H1N1 influenza epidemic on 11 June 2009. In that year, the virus is estimated to have killed 284,400 people globally. In August 2010, the WHO pronounced the epidemic to be finished. However, one of the viruses that cause seasonal flu is the H1N1 virus from the epidemic [6].

The initial line of defense against influenza viruses is often vaccination [7]. Despite the fact that the H1N1 2009 pandemic vaccination was made accessible in the U.S. in October 2009, there are major shortages in the number of vaccinations accessible. The full vaccine production process takes at least six months to complete [8]. Non-pharmaceutical measures, such as school closings and thermal scans at airports, were put into place during the recent pandemic (H1N1) 2009 epidemic to prevent the spread of the illness [9]. Using facemasks is another non-pharmaceutical method. Many people used facemasks during the SARS outbreak in 2003 to lessen their risk of getting sick. $76\%$ of Hong Kong people claimed to have used masks during the SARS outbreak in 2003 [10]. Infectious disease spread In order to determine how well facemasks can stop the transmission of sickness, particularly during the pandemic (H1N1) 2009, mathematical models of the spread of infectious disease might be helpful. The primary influenza activity period is divided into two categories based on when it started: the primary influenza activity period starts after October, and the primary influenza activity period starts after April. The WHO recommends administering the seasonal influenza vaccine before the start of the primary period of increased influenza activity [11,12].

Vaccination is the most effective method of intervention for halting the spread of diseases [13]. The WHO has been working with scientists and officials on a worldwide scale for more than 50 years to offer a coordinated approach to the development of influenza vaccines, including their effective use and distribution, manufacture, testing, and regulatory oversight [14]. In addition to reducing illness incidence, it also lessens the social and financial costs to society. Every year, in February and September, the World Health Organisation (WHO) holds technical consultations to suggest viruses for inclusion in seasonal influenza vaccinations for the northern and southern hemispheres [15].

Mathematical modeling is a helpful tool for investigating transmission patterns and assessing the outcomes of intervention initiatives. The mathematical models of influenza outbreaks have been studied and evaluated by several academics(since influenza is an air-borne disease). Tracht et al.[16] studied the effectiveness of facemasks in controlling the spread of the H1N1 influenza virus. Tchuenche et al.[17] studied the public views and reactions to infectious illnesses and the impact of media coverage. A studied H1N1 virus disease dynamics and sensitivity by Krishnapriya et.al.[18]. Kanyiri et al.[19] studied a Mathematical analysis of influenza transmission and viral control drug resistance. On the development of an age-dependent vaccination strategy for the 2009 A/H1N1 influenza epidemic in the Republic of Korea studied by Kim et al. [20]. Dalal et al.[21] evaluated the process of creating the genosensor and its effectiveness in identifying the H1N1 virus.

The role of children in the spread of pandemic influenza and the importance of immunization in reducing the occurrence of influenza. studied by Ratre et al.[22]. Baba et al.[23] studied A mathematical model that analyses the coexistence of two influenza strains and examines their method of transmission and treatment options is presented. Alharbi et al.[24] studied the transmission dynamics of influenza strains among pilgrims from different hemispheres and the effectiveness of vaccination in preventing infections. Ackerman et al.[25] studied mathematical models to identify the differences in immune processes between infections with different strains of the virus.

The numerical findings of our study are obtained using MATLAB’s built-in functions since we are more interested in examining the qualitative dynamical behaviors of the model under examination than in the precision, rate of convergence, etc. of the generated numerical solutions. To reduce vaccination and maximize profit, we have also created an optimum control problem for the system.

Sensitivity analysis in mathematical modeling of H1N1 helps identify which parameters most influence the model’s outcomes, allowing researchers to refine predictions and optimize control strategies. Studies have used techniques like Partial Rank Correlation Coefficient (PRCC) and Latin Hypercube Sampling to assess the impact of different variables on disease transmission [26]. One approach involves analyzing the basic reproduction number to determine how changes in parameters affect the spread of the virus. Another study explores bifurcation analysis, examining stability conditions and control measures like vaccination and treatment [27]. These methods help policymakers design effective interventions to mitigate outbreaks [28].

In this research, As a result, we suggested epidemiological models for H1N1 seasonal influenza viruses that take into account the virus’ antigenic alterations. The models contained details on modifications to the HA proteins’ amino acid sequences at epitope locations in the models for transmission. To do this, we first estimated the rate of time-varying antigenic change for each influenza subtype using the sequencing data. Finally, we showed that viral antigenic alterations may significantly alter the dynamics of seasonal influenza transmission at the population level. Thus, we proposed epidemiological models for the H1N1 seasonal influenza viruses. So, for the H1N1 seasonal influenza viruses, we suggested epidemiological models.

The study [29] examines the impact of contact with symptomatic and asymptomatic individuals on infectious disease spread using the SEIR model. It reveals that relative infectiousness and asymptomatic cases significantly influence disease transmission, emphasizing the need for isolation measures. [30] involves studying the spread of swine flu, differentiating between symptomatic and asymptomatic individuals, and applying optimal control strategies to manage the infection effectively. It typically involves mathematical modeling, using differential equations to analyze disease dynamics and control interventions. A recent research [27] explore the bifurcation analysis of influenza A (H1N1) models with treatment and vaccination models, using numerical simulations and mathematical theories. They investigate transcritical, Hopf, and backward bifurcations to understand disease transmission dynamics and the basic reproduction number’s role in disease persistence or death.

Sensitivity analysis is a crucial tool in disease modeling, assessing how variations in model parameters affect predictions [31]. It helps identify key factors affecting disease spread and control strategies, improves model accuracy, enhances decision-making, optimizes resource allocation, and supports uncertainty quantification [32]. It is widely used in infectious disease modeling, epidemiology, and biomedical research to refine strategies and minimize risks. Sensitivity analysis is particularly useful for evaluating the impact of different assumptions on disease forecasts, improving preparedness [31].

There are several studies that explore sensitivity analysis and optimal control strategies for H1N1 influenza. One study focuses on mathematical modeling of H1N1 transmission and control, using real-world data from Mexico, Italy, and South Africa to determine critical illness factors and forecast trends [33]. Another research paper introduces a Susceptible-Exposed-Infectious-Quarantined-Recovered (SEIQR) model, incorporating quarantine as a key intervention and applying optimal control theory to balance epidemiological impact with cost-effectiveness [26]. Additionally, a study examines an influenza model with vaccination and treatment, applying Pontryagin’s Maximum Principle to identify the best strategies for disease control [34].

The structure of the article is as follows. We develop the mathematical model in part 2, which is the following part. Section 3 has an analysis of fundamental qualities such as well-posedness, nonnegativity, boundedness, etc. In Section 4, equilibrium and stability analysis are covered. The optimum control issue is presented in Section 5. We presented the numerical simulations in Section 6. Finally, Section 7’s discussion brings the article to a close.

Table 1. List of recent influenza pandemics throughout the past century [19,35,36].

Pandemic	Year	Strain	People perished
Spanish flu	$1918–1920$	H1N1	$40–100$ million
Asian flu	$1957–1958$	H2N2	$1–2$ million
Hong Kong flu	$1968–1970$	H3N2	$0.5–2$ million
Swine flu	$2009–2010$	H1N1	Up to $\mathrm{575,000}$

Table 1. List of recent influenza pandemics throughout the past century [19,35,36].

Pandemic	Year	Strain	People perished
Spanish flu	$1918–1920$	H1N1	$40–100$ million
Asian flu	$1957–1958$	H2N2	$1–2$ million
Hong Kong flu	$1968–1970$	H3N2	$0.5–2$ million
Swine flu	$2009–2010$	H1N1	Up to $\mathrm{575,000}$

2. Derivation of the Mathematical Model

To capture the transmission dynamics of the Influenza virus, the total human population at any time t, denoted by $N\left(t\right)$, is stratified into six epidemiological compartments, each representing a distinct health status or disease progression phase:

(i) Susceptible population ($S\left(t\right)$): Individuals who have not yet contracted the virus and are at risk of infection through contact with infectious persons.

(ii) Exposed population ($E\left(t\right)$): Individuals who have been exposed to the virus and are in the incubation period. Although they do not show symptoms, they remain capable of transmitting the virus, thus contributing to its spread in the population.

(iii) Symptomatic infected population ($I_S\left(t\right)$): Individuals who have developed clinical symptoms of influenza and are actively contagious. This group represents the primary source of detectable and severe transmission cases.

(iv) Asymptomatic infected population ($I_A\left(t\right)$): Individuals who are infected but do not exhibit clinical symptoms. Despite their asymptomatic status, they play a critical role in the silent transmission of the virus.

(v) Recovered population ($R\left(t\right)$): Individuals who have overcome the infection and acquired immunity. They no longer participate in the transmission chain.

(vi) Vaccinated population ($V\left(t\right)$): Individuals who have received a vaccine and possess varying levels of protection depending on vaccine efficacy.

It is important to note that the exposed ($E\left(t\right)$), symptomatic infected ($I_S\left(t\right)$), and asymptomatic infected ($I_A\left(t\right)$) compartments collectively represent the infected individuals, differentiated based on the presence or absence of clinical symptoms and the stage of infection. This stratification allows for a more nuanced depiction of viral transmission, particularly in accounting for asymptomatic carriers and pre-symptomatic transmission—a key challenge in managing influenza outbreaks.

To describe the progression and spread of influenza within the stratified population, and to incorporate the influence of media-driven public health interventions, we introduce the following deterministic system of nonlinear ordinary differential equations. This system models the time-dependent interactions among the compartments under the influence of media awareness, which can alter behavioral patterns such as adherence to preventive measures, rate of vaccination, or reduction in contact frequency with infected individuals.

$$\begin{array}{ccc}\hfill {\displaystyle \frac{dS}{dt}}& =& \lambda -{\displaystyle \frac{\beta S(I_S+I_A)}{N}}-(\mu +a)S+\theta V,\hfill \\ \hfill {\displaystyle \frac{dE}{dt}}& =& {\displaystyle \frac{\beta S(I_S+I_A)}{N}}-(\kappa +\mu )E,\hfill \\ \hfill {\displaystyle \frac{d{I}_S}{dt}}& =& \rho \kappa E+b{I}_A-({\sigma }_2+\mu )I_S,\hfill \\ \hfill {\displaystyle \frac{d{I}_A}{dt}}& =& (1-\rho )\kappa E-(b+{\sigma }_1+\mu )I_A,\hfill \\ \hfill {\displaystyle \frac{dV}{dt}}& =& aS-(\theta +\delta +\mu )V,\hfill \\ \hfill {\displaystyle \frac{dR}{dt}}& =& \delta V+{\rho }_1{I}_A+{\sigma }_2{I}_S-\mu R,\hfill \end{array}$$

(1)

with initial values

$$\begin{array}{c}\hfill S\left(0\right)>0,I_S\left(0\right)>0,I_A\left(0\right)>0,V\left(0\right)>0,R\left(0\right)>0.\end{array}$$

(2)

In the formulation of the epidemiological model described by equation (1), each parameter carries specific biological significance. These parameters, which collectively define the dynamics of population transitions under viral infection and vaccination strategies, are elucidated as follows:

- $\lambda$ denotes the constant production rate of viral agents within the infected host population. This parameter captures the continual generation of infectious particles that contribute to disease propagation.

- $\beta$ represents the infection transmission rate, quantifying the probability of successful disease spread upon contact between susceptible and infectious individuals.

- The natural death rate of individuals within the model is denoted by $\mu$, encompassing mortality unrelated to the disease or vaccination process.

- The parameter a characterizes the rate of vaccination administration across the population, effectively indicating the proportion of susceptible individuals receiving immunization over time.

- Vaccination efficacy is described by $\theta$, which is assumed to be low; this reflects the limited protective effect conferred by the vaccine, potentially due to weak immune response or pathogen resistance. - The rate of progression from the exposed class to the infectious class is captured by $\kappa$, indicating the transition speed from latent infection to active disease.

- $\rho$ corresponds to the proportion of the population subject to disease transmission and progression, influencing the overall dynamics of the epidemic spread.

- The parameters ${\sigma }_1$ and ${\sigma }_2$ represent distinct recovery rates from infection, potentially differentiating recovery outcomes based on treatment status or disease severity.

- The recovery rate of vaccinated individuals who may still undergo infection is denoted by $\delta$, incorporating post-vaccination outcomes within the model framework.

- The symbol b appears to be unassigned in the current context; clarification or specification is required to complete its interpretation.

The numerical values and descriptions of these parameters are comprehensively listed in Table 2, which supports the analytical and simulation aspects of the model.

Figure 1. Diagrammatic representation of the model (1).

Figure 1. Diagrammatic representation of the model (1).

3. Positive Invariance and Boundedness of Solutions

In order to analyze the behavior of the influenza virus, we present the following characteristics.

3.1. Positive Invariance

We have the following theorem on the nonnegative of solutions.

Theorem 1.

Given the initial conditions in (2), all solutions to system (1) for all $t>0,$ are non-negative.

Proof. 

Let

$$\begin{array}{c}\hfill T_1=sup\{t>0:S\left(t\right)>0,E\left(t\right)>0,I_S\left(t\right)>0,I_A\left(t\right)>0,V\left(t\right)>0,R\left(t\right)>0\}.\end{array}$$

Using initial condition (2), we have $T_1>0$. Also, $S,E,I_S,I_A,V,andR$ are all equal to zero at $T_1$ if $T_1<\infty .$

From the first equation of (1), we can write

$$\begin{array}{ccc}\hfill {\displaystyle \frac{dS}{dt}}& =& \lambda -{\displaystyle \frac{\beta S(I_S+I_A)}{N}}-(\mu +a)S+\theta V.\hfill \end{array}$$

That is

$$\begin{array}{ccc}\hfill {\displaystyle \frac{dS}{dt}}+\left[{\displaystyle \frac{\beta (I_S+I_A)}{N}}+(\mu +a)\right]S& =& (\lambda +\theta V)\hfill \end{array}$$

Thus,

$$\begin{array}{ccc}\hfill {\displaystyle \frac{dS}{dt}}\left\{S\left(t\right)exp\left[{\int }_0^t\left({\displaystyle \frac{\beta (I_S+I_A)}{N}}+(\mu +a)\right)du\right]\right\}& =& (\lambda +\theta V)exp\left[{\int }_0^t\left({\displaystyle \frac{\beta (I_S+I_A)}{N}}+(\mu +a)\right)du\right]\hfill \end{array}$$

Hence,

$$\begin{array}{ccc}& & S\left(t\right)exp\left\{\left[{\int }_0^t\left({\displaystyle \frac{\beta (I_S+I_A)}{N}}+(\mu +a)\right)du\right]\right\}\hfill \\ & & =S\left(0\right)+{\int }_0^t\left\{(\lambda +\theta V)exp\left[{\int }_0^t\left({\displaystyle \frac{\beta (I_S+I_A)}{N}}+(\mu +a)\right)du\right]\right\}dz.\hfill \end{array}$$

So that,

$$\begin{array}{ccc}\hfill S\left(t\right)& =& S\left(0\right)exp\left[-{\int }_0^t\left({\displaystyle \frac{\beta (I_S+I_A)}{N}}+(\mu +a)\right)du\right]\hfill \\ & & +exp\left[-{\int }_0^t\left({\displaystyle \frac{\beta (I_S+I_A)}{N}}+(\mu +a)\right)du\right]\hfill \\ & & \times {\int }_0^t\left\{(\lambda +\theta V)exp\left[\left({\displaystyle \frac{\beta (I_S+I_A)}{N}}+(\mu +a)\right)du\right]\right\}dz>0\hfill \end{array}$$

Again, from system (1), we get

$$\begin{array}{ccc}\hfill {\displaystyle \frac{dE}{dt}}& =& {\displaystyle \frac{\beta S(I_S+I_A)}{N}}-(\kappa +\mu )E.\hfill \end{array}$$

That is

$$\begin{array}{ccc}\hfill {\displaystyle \frac{dE}{dt}}+(\kappa +\mu )E& =& {\displaystyle \frac{\beta S(I_S+I_A)}{N}}\hfill \end{array}$$

Thus,

$$\begin{array}{ccc}\hfill {\displaystyle \frac{dE}{dt}}\left\{E\left(t\right)exp\left[{\int }_0^t(\kappa +\mu )du\right]\right\}& =& {\displaystyle \frac{\beta S(I_S+I_A)}{N}}exp\left[{\int }_0^t(\kappa +\mu )du\right]\hfill \end{array}$$

Hence,

$$\begin{array}{ccc}\hfill E\left(t\right)exp\left[{\int }_0^t(\kappa +\mu )du\right]-E\left(0\right)& =& {\int }_0^t\left\{{\displaystyle \frac{\beta S(I_S+I_A)}{N}}exp\left[{\int }_0^t(\kappa +\mu )du\right]\right\}dz\hfill \end{array}$$

So that,

$$\begin{array}{ccc}\hfill E\left(t\right)& =& E\left(0\right)exp\left[-{\int }_0^t(\kappa +\mu )du\right]+exp\left[-{\int }_0^t(\kappa +\mu )du\right]\hfill \\ & & \times {\int }_0^t\left\{{\displaystyle \frac{\beta S(I_S+I_A)}{N}}exp\left[{\int }_0^t(\kappa +\mu )du\right]\right\}dz>0\hfill \end{array}$$

Similarly, from (1), we get

$$\begin{array}{ccc}\hfill I_S\left(t\right)& =& I_S\left(0\right)exp[-{\int }_0^t({\sigma }_2+\mu )du]+exp[-{\int }_0^t({\sigma }_2+\mu )du]\hfill \\ & & \times {\int }_0^t\{(\rho \kappa E+b{I}_A)exp\left[{\int }_0^t({\sigma }_2+\mu )du\right]\}dz>0,\hfill \end{array}$$

$$\begin{array}{ccc}\hfill I_A\left(t\right)& =& I_A\left(0\right)exp\left[-{\int }_0^t(b+{\sigma }_1+\mu )du\right]+exp\left[-{\int }_0^t(b+{\sigma }_1+\mu )du\right]\hfill \\ & & \times {\int }_0^t\left\{(1-\rho )\kappa Eexp\left[(b+{\sigma }_1+\mu )du\right]\right\}dz>0,\hfill \end{array}$$

$$\begin{array}{ccc}\hfill V\left(t\right)& =& V\left(0\right)exp\left[-{\int }_0^t(\theta +\delta +\mu )du\right]+exp[-{\int }_0^t(\theta +\delta +\mu )du]\hfill \\ & & \times {\int }_0^t\{aSexp\left[{\int }_0^t(\theta +\delta +\mu )du\right]\}dz>0,\hfill \end{array}$$

and

$$\begin{array}{ccc}\hfill R\left(t\right)& =& R\left(0\right)exp\left[-{\int }_0^t\left(\mu R\right)t]+exp[-{\int }_0^t\left(\mu R\right)du\right]\hfill \\ & & \times {\int }_0^t\{(\delta V+{\rho }_1{I}_A+{\sigma }_2{I}_S)exp\left[{\int }_0^t\left(\mu R\right)du\right]\}dz>0.\hfill \end{array}$$

Following the same procedure, we can show that for all$t>0$, we get $S\left(t\right)>0,E\left(t\right)>0,I_S\left(t\right)>0,I_A\left(t\right)>0,V\left(t\right)>0$ and $R\left(t\right)>0$. Thus, we conclude that all trajectories of system (1) remain positive in ${\mathbb{R}}_{+}^6$ for any all $t>0$. Hence, the proof is completed. □

3.2. Boundedness

We now begin with the theorem that guarantees that the solutions to system (1) are constrained to have nonnegative beginning values.

Theorem 2.

All solutions of system (1) are bounded in the region Ψ, defined by

$$\begin{array}{ccc}\hfill \Psi & =& \left\{(S,E,I_S,I_A,V,R)\in R^5:0\le S+E+I_S+I_A+V+R\le {\displaystyle \frac{\lambda }{\mu }}\right\}.\hfill \end{array}$$

Proof. 

Here, $N=S+E+I_S+I_A+V+R$. Thus, using system (1), we can write

$$\begin{array}{ccc}\hfill {\displaystyle \frac{dN}{dt}}& =& {\displaystyle \frac{dS}{dt}}+{\displaystyle \frac{dE}{dt}}+{\displaystyle \frac{d{I}_S}{dt}}+{\displaystyle \frac{d{I}_A}{dt}}+{\displaystyle \frac{dV}{dt}}+{\displaystyle \frac{dR}{dt}}\hfill \\ & \le & \lambda -\mu N.\hfill \end{array}$$

This gives, ${\displaystyle \frac{dN}{dt}}+\mu N\le \lambda$.

Using the theory of functional differential equation, $N\left(t\right)\le {\displaystyle \frac{\lambda }{\mu }}(1-e^{-\mu t})+N_0{e}^{-\mu t}.$

Thus, we have

$$\begin{array}{c}\hfill \underset{t\to 0}{lim\; sup}N\left(t\right)={\displaystyle \frac{\lambda }{\mu }}\end{array}$$

(3)

That means $\Psi$ is an invariant set for the system (1) and each solutions that starts in $\Psi$, always remains in $\Psi$.

□

4. Basic Reproduction Number, Equilibria and Stability Analysis

This section contains the derivation of basic reproduction number, existence of steady states and their stability analysis.

4.1. The Basic Reproduction Number

To calculate $R_0$, we use next generation method [38]. We consider two vector $\mathcal{F}$ and $\mathcal{G}$ defined for the system (1).

$$\begin{array}{ccc}\hfill \mathcal{F}& =& \left(\begin{array}{c}{\displaystyle \frac{\beta S(I_S+I_A)}{N}}\\ 0\\ 0\end{array}\right),\mathcal{G}=\left(\begin{array}{c}(\kappa +\mu )E\\ -\rho \kappa E-b{I}_A+(\mu +{\sigma }_2)I_S\\ -(1-\rho )\kappa E+(b+{\sigma }_1+\mu )I_A\end{array}\right).\hfill \end{array}$$

(4)

Next, we are setting that the entry-wise non-negative new infection matrix is F, and G be the non-singular Metzler matrix [39] define the transitions of influenza infection between the infection compartments and the matrices, which are given as follows:

$$\begin{array}{ccc}\hfill F& =& \left(\begin{array}{ccc}0& {\displaystyle \frac{\beta S}{N}}& {\displaystyle \frac{\beta S}{N}}\\ 0& 0& 0\\ 0& 0& 0\end{array}\right),G=\left(\begin{array}{ccc}(\kappa +\mu )& 0& 0\\ -\rho \kappa & (\mu +{\sigma }_2)& -b\\ (\rho -1)\kappa & 0& (b+{\sigma }_1+\mu )\end{array}\right),\hfill \end{array}$$

$$\begin{array}{ccc}\hfill G^{-1}& =& \left(\begin{array}{ccc}{\displaystyle \frac{1}{(\kappa +\mu )}}& 0& 0\\ {\displaystyle \frac{\kappa (\rho {\sigma }_1+\mu \rho +b)}{(\kappa +\mu )(\mu +{\sigma }_2)(b+{\sigma }_1+\mu )}}& {\displaystyle \frac{1}{(\mu +{\sigma }_2)}}& {\displaystyle \frac{b}{(\mu +{\sigma }_2)(b+{\sigma }_1+\mu )}}\\ {\displaystyle \frac{(1-\rho )\kappa }{(\kappa +\mu )(b+{\sigma }_1+\mu )}}& 0& {\displaystyle \frac{1}{(b+{\sigma }_1+\mu )}}\end{array}\right).\hfill \end{array}$$

Therefore, we have,

$$\begin{array}{ccc}\hfill F{G}^{-1}& =& \left[\begin{array}{ccc}0& {\displaystyle \frac{\beta S}{N}}& {\displaystyle \frac{\beta S}{N}}\\ 0& 0& 0\\ 0& 0& 0\end{array}\right.\times \left[\begin{array}{ccc}{\displaystyle \frac{1}{(\kappa +\mu )}}& 0& 0\\ {\displaystyle \frac{\kappa (\rho {\sigma }_1+\mu \rho +b)}{(\kappa +\mu )(\mu +{\sigma }_2)(b+{\sigma }_1+\mu )}}& {\displaystyle \frac{1}{(\mu +{\sigma }_2)}}& {\displaystyle \frac{b}{(\mu +{\sigma }_2)(b+{\sigma }_1+\mu )}}\\ {\displaystyle \frac{(1-\rho )\kappa }{(\kappa +\mu )(b+{\sigma }_1+\mu )}}& 0& {\displaystyle \frac{1}{(b+{\sigma }_1+\mu )}}\end{array}\right.,\hfill \end{array}$$

(5)

$$\begin{array}{ccc}& =& \left[\begin{array}{ccc}{\displaystyle \frac{\beta S^0\kappa [(b+\mu +{\sigma }_2)+\rho ({\sigma }_1-{\sigma }_2)]}{N(\kappa +\mu )(\mu +{\sigma }_2)(b+{\sigma }_1+\mu )}}& {\displaystyle \frac{\beta S^0}{N(\mu +{\sigma }_2)}}& {\displaystyle \frac{\beta S^0(b+{\sigma }_2+\mu )}{N(\mu +{\sigma }_2)(b+{\sigma }_1+\mu )}}\\ 0& 0& 0\\ 0& 0& 0\end{array}\right].\hfill \end{array}$$

(6)

With the help of Next generation matrix, we have $R_0$ the largest eigenvalue of $F{G}^{-1}$ at $E_0$

$$\begin{array}{ccc}\hfill R_0& =& {\displaystyle \frac{\beta S^0\kappa [(b+\mu +{\sigma }_2)+\rho ({\sigma }_1-{\sigma }_2)]}{N_0(\kappa +\mu )(\mu +{\sigma }_2)(b+{\sigma }_1+\mu )}}\hfill \end{array}$$

Theorem 3.

Model (1) has a basic reproduction number given in (6) at $E_0$. For $R_0>1$, the system (1) switches to endemic state.

4.2. Equilibrium Points

There are two equilibria of the system:

(i). The infection-free equilibrium is

$$\begin{array}{c}\hfill E_0\left({\displaystyle \frac{\lambda (\theta +\delta +\mu )}{(\mu +a)\left(\right(\theta +\delta +\mu )-a\theta )}},0,0,0,{\displaystyle \frac{\lambda a}{(\mu +a)\left(\right(\theta +\delta +\mu )-a\theta )}},{\displaystyle \frac{a\delta \lambda }{\mu \left(\right(\theta +\delta +\mu )-a\theta )}}\right)\end{array}$$

(ii). The disease state equilibrium is $E_1(S^{*},E^{*},I_S^{*},I_A^{*},V^{*},R^{*}),$ where

$$\begin{array}{ccc}\hfill S^{*}& =& {\displaystyle \frac{(\kappa +\mu )({\sigma }_2+\mu )(b+\mu +{\sigma }_1)N}{[(b+{\sigma }_2)+\mu (1+\sigma )+\rho ({\sigma }_1-{\sigma }_2-\mu )]\beta \kappa }},\hfill \\ \hfill I_S^{*}& =& {\displaystyle \frac{(b+\rho {\sigma }_1+\rho \mu )\kappa E^{*}}{({\sigma }_2+\mu )(b+\mu +{\sigma }_1)}},I_A^{*}={\displaystyle \frac{(1-\rho )\kappa E^{*}}{(b+\mu +{\sigma }_1)}},\hfill \\ \hfill V^{*}& =& {\displaystyle \frac{(\kappa +\mu )({\sigma }_2+\mu )(b+\mu +{\sigma }_1)Na}{(\theta +\delta +\mu )[(b+{\sigma }_2)+\mu (1+\sigma )+\rho ({\sigma }_1-{\sigma }_2-\mu )]\beta \kappa }},\hfill \\ \hfill R^{*}& =& {\displaystyle \frac{\left[(\kappa +\mu ){(\mu +{\sigma }_2)}^2{(b+\mu +{\sigma }_1)}^2\delta Na\right]+\left[(1-\rho )({\sigma }_2+\mu )(b+\mu +{\sigma }_2+\rho {\sigma }_1-\rho {\sigma }_2)\beta {\rho }_1\kappa E^{*}\right]}{(b+\mu +{\sigma }_2+\rho {\sigma }_1-\rho {\sigma }_2)(b+{\sigma }_1+\mu )({\sigma }_2+\mu )\beta }}\hfill \\ & & +{\displaystyle \frac{[b+\rho {\sigma }_1+\rho \mu )(b+\mu +{\sigma }_2+\rho {\sigma }_1-\rho {\sigma }_2)\beta ]\kappa {\sigma }_2{E}^{*}}{(b+\mu +{\sigma }_2+\rho {\sigma }_1-\rho {\sigma }_2)(b+{\sigma }_1+\mu )({\sigma }_2+\mu )\beta }},\hfill \end{array}$$

and $E^{*}$ is the positive root of the following equation:

$$\lambda N-\beta S^{*}(I_S^{*}+I_A^{*})-(\mu +a)S^{*}N+\theta V^{*}N=0.$$

The calculation is complected. Thus we find a positive value of $E^{*}$ using numerical calculations.

4.3. Local Stability

Theorem 4.

The system is stable at $E_0$ when $R_0<1$, and unstable when $R_0>1$.

Proof. 

At $E_0$, we have,

$$\begin{array}{ccc}\hfill J_{E_0}& =& \left(\begin{array}{cccccc}-(\mu +\theta )& 0& -{\displaystyle \frac{\beta S^0}{N}}& -{\displaystyle \frac{\beta S^0}{N}}& \theta & 0\\ 0& -(\kappa +\mu )& {\displaystyle \frac{\beta S^0}{N}}& {\displaystyle \frac{\beta S^0}{N}}& 0& 0\\ 0& \rho \kappa & -(\mu +{\sigma }_2)& b& 0& 0\\ 0& (1-\rho )\kappa & 0& -(b+{\sigma }_1+\mu )& 0& 0\\ a& 0& 0& 0& -(\theta +\delta +\mu )& 0\\ 0& 0& {\sigma }_1& {\sigma }_2& \delta & -\mu \end{array}\right).\hfill \end{array}$$

(7)

Here the eigenvalues are $\left[\xi =-\mu \right]$ and the other five eigenvalues are calculated from

$$\begin{array}{ccc}\hfill \left|\begin{array}{ccccc}-(\mu +\theta )-\xi & 0& -{\displaystyle \frac{\beta S^0}{N}}& -{\displaystyle \frac{\beta S^0}{N}}& \theta \\ 0& -(\kappa +\mu )-\xi & {\displaystyle \frac{\beta S^0}{N}}& {\displaystyle \frac{\beta S^0}{N}}& 0\\ 0& \rho \kappa & -(\mu +{\sigma }_2)-\xi & b& 0\\ 0& (1-\rho )\kappa & 0& -(b+{\sigma }_1+\mu )-\xi & 0\\ a& 0& 0& 0& -(\theta +\delta +\mu )-\xi \end{array}\right|& =& 0,\hfill \end{array}$$

and the characteristic equation for (7) is

$$\begin{array}{ccc}\hfill {\xi }^5+B_1{\xi }^4+B_2{\xi }^3+B_3{\xi }^2+B_4\xi +B_5& =& 0,\hfill \end{array}$$

(8)

where,

$$\begin{array}{ccc}\hfill B_1& =& -(b_1+b_4),B_2=(b_2+b_5+b_6+b_1{b}_4),\hfill \\ \hfill B_3& =& -(b_3+b_1{b}_5+b_1{b}_6+b_2{b}_4),B_4=(b_2{b}_5+b_2{b}_6+b_3{b}_4),\hfill \\ \hfill B_5& =& -b_3(b_5+b_6),\hfill \end{array}$$

and where,

$$\begin{array}{ccc}\hfill b_1& =& (a_{22}+a_{33}+a_{44}),\hfill \\ \hfill b_2& =& a_{33}{a}_{44}+(a_{22}{a}_{44}-a_{24}{a}_{42})+(a_{22}{a}_{33}-a_{32}{a}_{23}),\hfill \\ \hfill b_3& =& a_{22}{a}_{33}{a}_{44}-a_{23}(a_{32}{a}_{44}-a_{42}{a}_{34})-a_{24}{a}_{42}{a}_{33}.\hfill \end{array}$$

where,

$$\begin{array}{ccc}\hfill a_{22}& =& -(\kappa +\mu ),a_{33}=-(\mu +{\sigma }_2),a_{44}=-(b+{\sigma }_1+\mu ),\hfill \\ \hfill a_{24}& =& a_{23}={\displaystyle \frac{\beta S}{N}},a_{42}=(1-\rho )\kappa ,a_{32}=\rho \kappa ,a_{34}=b.\hfill \end{array}$$

For $A_0\ge 0$ and $A_1>0$. If $A_0>0$, then roots must be negative. For $A_2>0$, we get the threshold condition for $E_0$. We get the condition ${\mu }_1+2{\displaystyle \frac{p{\widehat{H}}_S}{H_T}}<\nu <{\displaystyle \frac{{\mu }_2{\mu }_3}{{\beta }_1{\widehat{H}}_S}}$ which implies that $R_0<1$ which results the removal of infection. □

4.4. Local Equilibrium Stability of the Endemic System

Proof. 

At $E^{*}$, we have the charismatic equation

$$\begin{array}{ccc}\hfill |J_{E^{*}}-\xi |& =& \left(\begin{array}{ccccc}{{\rm Y}}_{11}-\xi & 0& {{\rm Y}}_{13}& {{\rm Y}}_{14}& {{\rm Y}}_{15}\\ {{\rm Y}}_{21}& {{\rm Y}}_{22}-\xi & {{\rm Y}}_{23}& {{\rm Y}}_{24}& 0\\ 0& {{\rm Y}}_{32}& {{\rm Y}}_{33}-\xi & {{\rm Y}}_{34}& 0\\ 0& {{\rm Y}}_{42}& 0& {{\rm Y}}_{44}-\xi & 0\\ {{\rm Y}}_{51}& 0& 0& 0& {{\rm Y}}_{55}-\xi \end{array}\right).\hfill \end{array}$$

(9)

where,

$$\begin{array}{ccc}\hfill {{\rm Y}}_{11}& =& -{\displaystyle \frac{\beta (I_S+I_A)}{N}}-(\mu +\theta ),{{\rm Y}}_{13}=-{\displaystyle \frac{\beta S}{\left(N\right)}},{{\rm Y}}_{14}=-{\displaystyle \frac{\beta S}{\left(N\right)}},\hfill \\ \hfill {{\rm Y}}_{15}& =& \theta ,{{\rm Y}}_{21}={\displaystyle \frac{\beta (I_S+I_A)}{N}},{{\rm Y}}_{22}=-(\kappa +\mu ),{{\rm Y}}_{23}={\displaystyle \frac{\beta S}{N}},\hfill \\ \hfill {{\rm Y}}_{24}& =& {\displaystyle \frac{\beta S}{N}},{{\rm Y}}_{32}=\rho \kappa ,{{\rm Y}}_{33}=-(\mu +{\sigma }_2),{{\rm Y}}_{34}=b,\hfill \\ \hfill {{\rm Y}}_{42}& =& (1-\rho )\kappa ,{{\rm Y}}_{44}=-(b+{\sigma }_1+\mu ),{{\rm Y}}_{51}=a,{{\rm Y}}_{55}=-(\theta +\delta +\mu ).\hfill \end{array}$$

At $E^{*}$, the characteristic equation is

$$\begin{array}{ccc}\hfill {\xi }^5+{\Phi }_1{\xi }^4+{\Phi }_2{\xi }^3+{\Phi }_3{\xi }^2+{\Phi }_4\xi +{\Phi }_5\xi & =& 0.\hfill \end{array}$$

(10)

where,

$$\begin{array}{ccc}& & {\Phi }_1=-(b_1+b_4),{\Phi }_2=(b_2+b_5+b_6+b_1{b}_4)\hfill \\ & & {\Phi }_3=-(b_3+b_1{b}_5+b_1{b}_6+b_2{b}_4),{\Phi }_4=(b_2{b}_5+b_2{b}_6+b_3{b}_4)\hfill \\ & & {\Phi }_5=-(b_3{b}_5+b_3{b}_6)\hfill \end{array}$$

$b_1=$ Sum of the minor of order one along the main diagonal of $\left[J_{E^{*}}\right]$

$b_2=$ Sum of the minors of order two of the diagonal element of $\left[J_{E^{*}}\right]$

$b_3=$Sum of the minors of order two of the diagonal element of$\left[J_{E^{*}}\right]$, etc.

□

5. Parameter Sensitivity Analysis

Sensitivity and identifiability studies were conducted by computing the derivative of the clinical score $\Omega ={[S,E,I_S,I_A,V,R]}^T$ with respect to the vector parameter at each data point. The sensitivities are dimensionless because of their scaling with respect to the values of the variables and parameters. The following is a representation of the system (1) sensitivity functions with respect to any parameter p:

$$\begin{array}{ccc}\hfill {\Omega }_i,\omega & =& {\displaystyle \frac{\partial {\Omega }_i\left(t\right)}{\partial \omega }},i=1,2,...6\hfill \end{array}$$

(11)

Due to these sensitivities, the least important parameter for model output may be found and changed without being used for calibration. The sensitivity functions of DDEs can be found using a variety of techniques [18,40]. However, to determine the sensitivity functions of the system (1), we shall employ the so-called direct technique because of its simplicity.

According to the parameter $\lambda$, the relevant sensitivity system (1), given below

Similarly, the sensitivity of the other parameter can be studied.

6. Optimal Control Problem Formulation

In this section, we have considered two non-medical controls, like awareness through governance and self-learning, to control the disease. Here, our main focus is to minimize the number of infected and exposed people. For this purpose, we have introduced control strategies $(u_v,u_i)$.

• Control $u_v$ represents the role of exposed class control.
• Control $u_i$ plays a crucial role in controlling the infected class.

Thus, our optimal control problem becomes

$$\begin{array}{ccc}\hfill {\displaystyle \frac{dS}{dt}}& =& \lambda -{\displaystyle \frac{\beta S(I_S+I_A)}{N}}-(\mu +u_{v}a)S+\theta V,\hfill \\ \hfill {\displaystyle \frac{dE}{dt}}& =& {\displaystyle \frac{\beta S(I_S+I_A)}{N}}-(u_i\kappa +\mu )E,\hfill \\ \hfill {\displaystyle \frac{d{I}_S}{dt}}& =& u_i\rho \kappa E+b{I}_A-({\sigma }_2+\mu )I_S,\hfill \\ \hfill {\displaystyle \frac{d{I}_A}{dt}}& =& u_i(1-\rho )\kappa E-(b+{\sigma }_1+\mu )I_A,\hfill \\ \hfill {\displaystyle \frac{dV}{dt}}& =& u_{v}aS-(\theta +\delta +\mu )V,\hfill \\ \hfill {\displaystyle \frac{dR}{dt}}& =& \delta V+{\rho }_1{I}_A+{\sigma }_2{I}_S-\mu R.\hfill \end{array}$$

(12)

$$\begin{array}{ccc}\hfill J(u_v,u_i)& =& {\int }_0^T[P_{1}E+P_2{I}_S+P_3{I}_A+Q_1{u}_v^2+Q_2{u}_i^2]dt\hfill \end{array}$$

(13)

In the objection for ${\int }_0^T{P}_{1}Edt$ is the cost induced by exposed clans, ${\int }_0^T{P}_2{I}_{S}dt$ is the cost induced by the infected clam, ${\int }_0^T{P}_3{I}_{A}dt$ is the cost induced by any phase infected class and ${\int }_0^T{Q}_1{u}_v^{2}dt$ on the cost induced by optimal stage, ${\int }_0^T{Q}_2{u}_i^{2}dt$ as the costs at stage.

Here the optimal controls $u_v^{*},u_i^{*}\in \mathcal{U}$ where

$$\begin{array}{ccc}\hfill J(u_v^{*},u_i^{*})& =& minJ(u_v,u_i),\hfill \end{array}$$

where $\mathcal{U}=\{(u_v,u_i)$ such that $u_v,u_i$ is Lebesgue measure on $[0,1]$ with $0\le u_v\le 1,0\le u_i\le 1$}

6.1. Existence of an Optimal Control

Have we the existence of an optimal control model (12)

Theorem 5.

$\mathcal{U}=(u_v,u_i)$ and $\mathcal{M}=(S,E,I_S,I_A,V,R)$ exist for the stake initial value problem (12) and the (13) minimizes $J(u_v,u_i)$ over $\mathcal{U}$

Proof. 

To verify the existence condition, we use the study of family and Rishel’s theorem under certain conditions.

• Systems (12) and (13) with the control fraction V are non-empty.
• The state system should be a linear function dependent on time and time state variables.
• The integral J in (13) is convex on $\mathcal{U}$ and the the exi $\alpha >1$ and there exist $\chi$ a positive numbers and $c_1$ and $c_2$ such that

$$\begin{array}{ccc}\hfill J(t,\mathcal{M},u)& \ge & c_1{|u_v,u_i|}^{\chi }-c_2\hfill \end{array}$$

For every admissible control $\mathcal{U}$, there is a unique solution satisfying the Lipsehiiz condition. Also, we know that the total population is bounded, and it belongs to ${\displaystyle \frac{\delta }{\mu }}$. Also, all the state variables are bounded. Thus, we have a partial derivation of the state system (12) with respect to the state variable system bondedness, which satisfies condition $\left(i\right)$. We can also verify condition $\left(ii\right)$ by checking the liner dependence of the state equation and control $(u_v,u_i)$, for verification function condition (iii). We note that linear and quarantine functions are convex [41]. To prove the bondedness on J, we have $Q_1{u}_v^2\le Q_1$ and $Q_2{u}_i^2\le Q_2$ an $0\le u_v\le 1$ and $0\le u_i\le 1$.

It follows,

$$\begin{array}{ccc}\hfill J(t,\mathcal{M},u)& =& P_{1}E+P_2{I}_S+P_3{I}_A+Q_1{u}_v^2+Q_2{u}_i^2\hfill \\ & \ge & (Q_1{u}_v^2+Q_2{u}_1^2-u_1^2)\hfill \end{array}$$

(14)

from (14) we get

$$\begin{array}{ccc}\hfill J(t,\mathcal{M},u)& \ge & min(Q_1,Q_2)(u_v^2+u_i^2)-u_i\hfill \\ & =& min\left(Q_1{Q}_2\right)\left|\right|u_v,u_i{\left|\right|}^2-u_i.\hfill \end{array}$$

(15)

Thus,

$$\begin{array}{ccc}\hfill J(t,\mathcal{M},u)& \ge & c_1\left|\right|(u_v,u_i){\left|\right|}^{\chi }-c_2,\hfill \end{array}$$

where $c_1=min(u_v,u_i),c_2=u_i$ and $\chi =2$. □

Theorem 6.

When $R_0<1$, $E_0$ is globally asymptotically stable.

Proof. 

Let $L={\psi }_{1}E+{\psi }_2{I}_S+{\psi }_3{I}_A$ be a Lyapunov function I in which we only consider the exposed and infected class. with positive contract ${\psi }_1,{\psi }_2,{\psi }_3$. Then,

$$\begin{array}{ccc}\hfill {\displaystyle \frac{dL}{dt}}& =& {\psi }_1{\displaystyle \frac{dE}{dt}}+{\psi }_2{\displaystyle \frac{d{I}_S}{dt}}+{\psi }_3{\displaystyle \frac{d{I}_A}{dt}}\hfill \\ \hfill & =& {\psi }_1\left[{\displaystyle \frac{\beta S(I_S+I_A)}{N}}-(\kappa +\mu )E\right]\hfill \\ & & +{\psi }_2\left[\rho \kappa E+b{I}_A-({\sigma }_2+\mu )I_S\right]\hfill \\ & & +{\psi }_3\left[(1-\rho )\kappa E-(b+{\sigma }_1+\mu )I_A\right]\hfill \\ & & \le {\psi }_1\left[{\displaystyle \frac{\beta \lambda (I_S+I_A)}{\mu N}}-(\kappa +\mu )E\right]\hfill \\ & & +{\psi }_2\left[\rho \kappa E+b{I}_A-({\sigma }_2+\mu )I_S\right]\hfill \\ & & +{\psi }_3\left[(1-\rho )\kappa E-(b+{\sigma }_1+\mu )I_A\right]\hfill \\ & & =\left[{\psi }_3(1-\rho )\kappa -{\psi }_1(\kappa +\mu )+{\psi }_2\rho \kappa \right]E\hfill \\ & & +\left[{\psi }_1{\displaystyle \frac{\beta \lambda }{\mu N}}-\psi ({\sigma }_2+\mu )\right]I_S\hfill \\ & & +\left[{\psi }_1{\displaystyle \frac{\beta \lambda }{\mu N}}+{\psi }_{2}b-{\psi }_3(b+{\sigma }_2+\mu )\right]I_A\hfill \end{array}$$

(16)

□

6.2. Characterization of the Optimal Control Pair

Let $u^{*}\left(t\right)$ represent the amount of medicine as control input. The essential criteria for an optimum control and related states to meet Pontryagin’s Maximum Principle are represented by the cost function (13) for the system (12). We employ Pontryagin’s maximal principle [42] to choose the best controls for $u_v^{*}\left(t\right)$ and $u_i^{*}\left(t\right)$. By using this technique, we create the Hamiltonian function L with regard to $(u_v\left(t\right),u_i\left(t\right))$, which transforms the system (12) along with (13).

According to the control system provided by (12), we determine the problem’s optimal settings. To define the Lagrangian ($\mathcal{L}$), the Hamiltonian associated with the control problem generated in this way is also included.

$$\begin{array}{ccc}\hfill \mathcal{L}(u_v,u_i)& =& [P_{1}E+P_2{I}_S+P_3{I}_A+Q_1{u}_v^2+Q_2{u}_i^2]dt\hfill \end{array}$$

(17)

The Hamiltonian is built as follows by combining the state variables and adjoint variables:

$$\begin{array}{ccc}\hfill \mathcal{H}& =& P_{1}E+P_2{I}_S+P_3{I}_A+Q_1{u}_v^2+Q_2{u}_i^2\hfill \\ & & +{\xi }_1\left[\lambda -{\displaystyle \frac{\beta S(I_S+I_A)}{N}}-(\mu +u_{v}a)S+\theta V\right],\hfill \\ & & +{\xi }_2\left[{\displaystyle \frac{\beta S(I_S+I_A)}{N}}-(u_i\kappa +\mu )E\right],\hfill \\ & & +{\xi }_3\left[u_i\rho \kappa E+b{I}_A-({\sigma }_2+\mu )I_S\right],\hfill \\ & & +{\xi }_4\left[u_i(1-\rho )\kappa E-(b+{\sigma }_1+\mu )I_A\right],\hfill \\ & & +{\xi }_5\left[u_{v}aS-(\theta +\delta +\mu )V\right],\hfill \\ & & +{\xi }_6\left[\delta V+{\rho }_1{I}_A+{\sigma }_2{I}_S-\mu R\right].\hfill \end{array}$$

(18)

Here, ${\xi }_i,i=1,2,3,4,5,6$

$$\begin{array}{ccc}\hfill {\displaystyle \frac{\partial \mathcal{H}}{\partial S}}& =& -{\xi }_1\left[{\displaystyle \frac{\beta (I_S+I_A)}{N}}+(\mu +u_{v}a)\right]+{\xi }_2\left[{\displaystyle \frac{\beta (I_S+I_A)}{N}}\right]+{\xi }_5{u}_{v}a,\hfill \\ \hfill {\displaystyle \frac{\partial \mathcal{H}}{\partial E}}& =& P_1-{\xi }_2(u_i\kappa +\mu )+{\xi }_3{u}_i\rho \kappa +{\xi }_4{u}_i\kappa (1-\rho ),\hfill \\ \hfill {\displaystyle \frac{\partial \mathcal{H}}{\partial I_S}}& =& P_2-{\xi }_1\left[{\displaystyle \frac{\beta S}{N}}\right]+{\xi }_2\left[{\displaystyle \frac{\beta S}{N}}\right]-{\xi }_3({\sigma }_2+\mu )+{\xi }_6{\sigma }_2,\hfill \\ \hfill {\displaystyle \frac{\partial \mathcal{H}}{\partial I_A}}& =& P_3-{\xi }_1\left[{\displaystyle \frac{\beta S}{N}}\right]+{\xi }_2\left[{\displaystyle \frac{\beta S}{N}}\right]+{\xi }_{3}b-{\xi }_4\left(b+{\sigma }_1+\mu \right)+{\xi }_6{\rho }_1,\hfill \\ \hfill {\displaystyle \frac{\partial \mathcal{H}}{\partial V}}& =& {\xi }_1\theta -{\xi }_5(\theta +\delta +\mu )+{\xi }_6\delta ,\hfill \\ \hfill {\displaystyle \frac{\partial \mathcal{H}}{\partial R}}& =& -{\xi }_6\mu ,\hfill \\ \hfill {\displaystyle \frac{\partial \mathcal{H}}{\partial u_v}}& =& 2Q_1{u}_v-aS({\xi }_1-{\xi }_5),\hfill \\ \hfill {\displaystyle \frac{\partial \mathcal{H}}{\partial u_i}}& =& 2Q_2{u}_i-{\xi }_2\kappa E+{\xi }_3\rho \kappa E+{\xi }_4\kappa E(1-\rho ).\hfill \end{array}$$

(19)

The adjoint system to be estimated for the control input $(u_v\left(t\right),u_i\left(t\right))$ associated with the model state variables $S,E,I_S,I_A,V,R$ is represented as For control input, $(u_v\left(t\right),u_i\left(t\right))$ along with $S,E,I_S,I_A,V,R$ is represented as

$$\begin{array}{ccc}\hfill {\displaystyle \frac{d{\xi }_1}{dt}}& =& -\left[-{\xi }_1\left[{\displaystyle \frac{\beta (I_S+I_A)}{N}}+(\mu +u_{v}a)\right]+{\xi }_2\left[{\displaystyle \frac{\beta (I_S+I_A)}{N}}\right]+{\xi }_5{u}_{v}a\right],\hfill \\ \hfill {\displaystyle \frac{d{\xi }_2}{dt}}& =& -\left[P_1-{\xi }_2(u_i\kappa +\mu )+{\xi }_3{u}_i\rho \kappa +{\xi }_4\kappa u_i(1-\rho )\right],\hfill \\ \hfill {\displaystyle \frac{d{\xi }_3}{dt}}& =& -\left[P_2-{\xi }_1\left({\displaystyle \frac{\beta S}{N}}\right)+{\xi }_2\left({\displaystyle \frac{\beta S}{N}}\right)-{\xi }_3({\sigma }_2+\mu )+{\xi }_6{\sigma }_2\right],\hfill \\ \hfill {\displaystyle \frac{d{\xi }_4}{dt}}& =& -\left[P_3-{\xi }_1\left({\displaystyle \frac{\beta S}{N}}\right)+{\xi }_2\left({\displaystyle \frac{\beta S}{N}}\right)+{\xi }_{3}b-{\xi }_4\left(b+{\sigma }_1+\mu \right)+{\xi }_6{\rho }_1\right],\hfill \\ \hfill {\displaystyle \frac{d{\xi }_5}{dt}}& =& -\left[{\xi }_1\theta -{\xi }_5(\theta +\delta +\mu )+{\xi }_6\delta \right],\hfill \\ \hfill {\displaystyle \frac{d{\xi }_6}{dt}}& =& {\xi }_6\mu .\hfill \end{array}$$

(20)

Here, the transversality conditions are ${\xi }_i\left(T\right)=0,i=1,2,...6.$ Pontryagin’s Maximum Principle [42] has allowed us to ascertain

$$\begin{array}{c}\hfill {\displaystyle \frac{\partial \mathcal{H}}{\partial u^{*}\left(t\right)}}=0.\end{array}$$

(21)

Then,

$$\begin{array}{ccc}\hfill {\displaystyle \frac{\partial \mathcal{H}}{\partial u_v}}& =& 2Q_1{u}_v-aS({\xi }_1-{\xi }_5),\hfill \\ \hfill {\displaystyle \frac{\partial \mathcal{H}}{\partial u_i}}& =& 2Q_2{u}_i-{\xi }_2\kappa E+{\xi }_3\rho \kappa E+{\xi }_4\kappa E(1-\rho ).\hfill \end{array}$$

(22)

Solving(22) for $u_v^{*}\left(t\right)$ and $u_i^{*}\left(t\right)$,

$$\begin{array}{ccc}\hfill u_v^{*}\left(t\right)& =& {\displaystyle \frac{aS({\xi }_1-{\xi }_5)}{2Q_1}},\hfill \\ \hfill u_i^{*}\left(t\right)& =& {\displaystyle \frac{\kappa E[{\xi }_2-{\xi }_3\rho -{\xi }_4(1-\rho )]}{2Q_2}}.\hfill \end{array}$$

(23)

The standard form is

$$\begin{array}{ccc}\hfill u_v^{*}\left(t\right)& =& \left\{\begin{array}{cc}0,\hfill & {\displaystyle \frac{aS({\xi }_1-{\xi }_5)}{2Q_1}}<0,\hfill \\ {\displaystyle \frac{aS({\xi }_1-{\xi }_5)}{2Q_1}},\hfill & 0<{\displaystyle \frac{aS({\xi }_1-{\xi }_5)}{2Q_1}}<1,\hfill \\ 1,\hfill & {\displaystyle \frac{aS({\xi }_1-{\xi }_5)}{2Q_1}}>1,\hfill \end{array}\right.\hfill \end{array}$$

(24)

$$\begin{array}{ccc}\hfill u_i^{*}\left(t\right)& =& \left\{\begin{array}{cc}0,\hfill & {\displaystyle \frac{\kappa E[{\xi }_2-{\xi }_3\rho -{\xi }_4(1-\rho )]}{2Q_2}}<0,\hfill \\ {\displaystyle \frac{\kappa E[{\xi }_2-{\xi }_3\rho -{\xi }_4(1-\rho )]}{2Q_2}},\hfill & 0<{\displaystyle \frac{\kappa E[{\xi }_2-{\xi }_3\rho -{\xi }_4(1-\rho )]}{2Q_2}}<1,\hfill \\ 1,\hfill & {\displaystyle \frac{\kappa E[{\xi }_2-{\xi }_3\rho -{\xi }_4(1-\rho )]}{2Q_2}}>1.\hfill \end{array}\right.\hfill \end{array}$$

(25)

Thus we can write $u_v^{*}\left(t\right)$ as

$$\begin{array}{ccc}\hfill u_v^{*}\left(t\right)& =& max\left(min\left(1,{\displaystyle \frac{aS({\xi }_1-{\xi }_5)}{2Q_1}}\right),0\right).\hfill \end{array}$$

(26)

$$\begin{array}{ccc}\hfill u_i^{*}\left(t\right)& =& max\left(min\left(1,{\displaystyle \frac{\kappa E[{\xi }_2-{\xi }_3\rho -{\xi }_4(1-\rho )]}{2Q_2}}\right),0\right).\hfill \end{array}$$

(27)

Thus, the optimal system becomes

$$\begin{array}{ccc}\hfill {\displaystyle \frac{dS}{dt}}& =& \lambda -{\displaystyle \frac{\beta S(I_S+I_A)}{N}}-(\mu +u_{v}a)S+\theta V,\hfill \\ \hfill {\displaystyle \frac{dE}{dt}}& =& {\displaystyle \frac{\beta S(I_S+I_A)}{N}}-(u_i\kappa +\mu )E,\hfill \\ \hfill {\displaystyle \frac{d{I}_S}{dt}}& =& u_i\rho \kappa E+b{I}_A-({\sigma }_2+\mu )I_S,\hfill \\ \hfill {\displaystyle \frac{d{I}_A}{dt}}& =& u_i(1-\rho )\kappa E-(b+{\sigma }_1+\mu )I_A,\hfill \\ \hfill {\displaystyle \frac{dV}{dt}}& =& u_{v}aS-(\theta +\delta +\mu )V,\hfill \\ \hfill {\displaystyle \frac{dR}{dt}}& =& \delta V+{\rho }_1{I}_A+{\sigma }_2{I}_S-\mu R,\hfill \\ \hfill {\displaystyle \frac{d{\xi }_1}{dt}}& =& -\left[-{\xi }_1\left[{\displaystyle \frac{\beta (I_S+I_A)}{N}}+(\mu +u_{v}a)\right]+{\xi }_2\left[{\displaystyle \frac{\beta (I_S+I_A)}{N}}\right]+{\xi }_5{u}_{v}a\right],\hfill \\ \hfill {\displaystyle \frac{d{\xi }_2}{dt}}& =& -\left[P_1-{\xi }_2(u_i\kappa +\mu )+{\xi }_3{u}_i\rho \kappa +{\xi }_4{u}_i\kappa (1-\rho )\right],\hfill \\ \hfill {\displaystyle \frac{d{\xi }_3}{dt}}& =& -\left[P_2-{\xi }_1\left({\displaystyle \frac{\beta S}{N}}\right)+{\xi }_2\left({\displaystyle \frac{\beta S}{N}}\right)-{\xi }_3({\sigma }_2+\mu )+{\xi }_6{\sigma }_2\right],\hfill \\ \hfill {\displaystyle \frac{d{\xi }_4}{dt}}& =& -\left[P_3-{\xi }_1\left({\displaystyle \frac{\beta S}{N}}\right)+{\xi }_2\left({\displaystyle \frac{\beta S}{N}}\right)+{\xi }_{3}b-{\xi }_4\left(b+{\sigma }_1+\mu \right)+{\xi }_6{\rho }_1\right],\hfill \\ \hfill {\displaystyle \frac{d{\xi }_5}{dt}}& =& -\left[{\xi }_1\theta -{\xi }_5(\theta +\delta +\mu )+{\xi }_6\delta \right],\hfill \\ \hfill {\displaystyle \frac{d{\xi }_6}{dt}}& =& {\xi }_6\mu ,\hfill \end{array}$$

(28)

along with the transversality conditions: ${\xi }_i\left(T\right)=0,i=1,2,...,6$.

7. Numerical Simulations

In this section, we analyze our model and the optimal system graphically using Matlab simulations. All the parameters are taken from Table 2 and some parameters are varied for studying their impacts on the system dynamics.

As analysed in section 5, we plot the sensitivity of model parameters in Figure 2 and also in Table 3. All parameters are sensitive for the study of the H1N1 dynamics. Some parameters, namely $\lambda ,\beta ,\kappa$ etc., are positively sensitive and some parameters are negatively sensitive to the disease.

Figure 3 shows that when $R_0<1$, the system attains its disease-free state. Also, when $R_0>1$, the system moves to a stable endemic steady state. Figure 4 shows the forward bifurcation at the basic reproduction $R_0$. We have varied the infection rate $\beta$, and it is observed that when $R_0<1$, the system becomes disease-free, and it shows its stability around $E_0$. It switches to an endemic state when $R_0>1$ and the system attains stability around $E^{*}$.

Figure 5 shows the stability region for the disease-free state and endemic state with respect to parameters $\beta$ and $\kappa$. From this figure, we reveal that with the increasing values of $\beta$ and $\kappa$, the system switches its stability from $E_0$ to $E^{*}$ that is from disease-free to endemic state. Nonlinear stability of the endemic system is shown in Figure 6. All phase trajectories converge to the same endemic steady state $E^{*}$ showing the global stability of $E^{*}$.

Impact of $\kappa$ is analyzing using the Figure 7 by plotting the equilibrium values of the model population with respect to $\kappa$. Similarly, in Figure 8, the impact of vaccination is shown. Infection is reduces as the vaccinated populations are increased.

Also, we consider the control system with control inputs $u_v\left(t\right)$ and $u_i\left(t\right)$ for the systems (12) and (13). The time-dependent control functions, along with the state system and adjoint system, are analyzed numerically. Here, the optimal control problem is bounded by $[0,T]$. Here, we consider $T=60$ at which the treatment is stopped. Optimal system is plotted in Figure 9. This figure shows the impacts of controls on the disease. In Figure 10, the optimal trajectories show the effective role of control strategies. Also, we have clearly observed that a low weight factor gives a better result in increasing the susceptible population along with class improvement.

8. Discussion and Conclusion

The H1N1 virus is a major contributor to influenza infections, primarily affecting the upper respiratory tract. Severe manifestations are typically observed in individuals with underlying chronic conditions. To better understand and mitigate the spread of this virus, a mathematical model has been developed and extended to incorporate time-dependent control inputs. The model framework includes the establishment of existence conditions and a thorough stability analysis based on the basic reproduction number, $R_0$.

A comprehensive sensitivity analysis reveals that all model parameters significantly influence the dynamics of H1N1 transmission. Notably, parameters such as the transmission rate ($\beta$), recruitment rate ($\lambda$), and contact rate ($\kappa$) exhibit positive sensitivity, indicating that increases in these parameters tend to amplify the spread of infection (Figure 2).

We illustrated the critical role of the basic reproduction number $R_0$ in determining system behavior (Figure 3). When $R_0<1$, the model stabilizes at a disease-free equilibrium ($E_0$), whereas $R_0>1$ leads to a stable endemic equilibrium ($E^{*}$). This transition is further substantiated by the forward bifurcation (Figure 4), which highlights a distinct shift in system dynamics as $R_0$ crosses the threshold value unity. We have seen that when endemic system is feasible, it is global stabile (Figure 6). The impact of the contact rate $\kappa$ is examined which shows that higher values of $\kappa$ correspond to increased equilibrium infection levels (Figure 7). This finding emphasizes the importance of managing contact rates to control disease spread. Additionally, the positive effect of vaccination is studied where an increase in vaccinated individuals leads to a marked reduction in infection prevalence (Figure 8).

For the analysis of the control model, Pontryagin’s Maximum Principle has been employed. The corresponding Hamiltonian function has been constructed, and the adjoint equations, along with the optimal control variables $u_1^{*}$ and $u_2^{*}$, have been derived. To enhance disease mitigation, the model incorporates a control system governed by time-dependent inputs $u_v\left(t\right)$ (vaccination) and $u_i\left(t\right)$ (treatment), analyzed over a finite time horizon $T=60$. We have seen that the optimal control strategies are effective. The simulations reveal that control interventions significantly reduce infection levels. Moreover, a lower weight factor in the cost function yields improved outcomes, including an increase in the susceptible population and enhanced class dynamics, suggesting that more aggressive control measures are beneficial (Figure 9 and Figure 10).

In summary, this study presents a rigorous analysis of the H1N1 transmission model, integrating sensitivity analysis, bifurcation theory, and optimal control strategies. The basic reproduction number $R_0$ serves as a pivotal threshold for system stability. Parameters such as $\beta$, $\lambda$, and $\kappa$ play critical roles in shaping disease dynamics. Vaccination and contact reduction are effective in curbing infection spread. Optimal control strategies, particularly those with lower cost weights, substantially improve population health outcomes.

In conclusion, the results underscore the importance of targeted interventions and strategic parameter management in controlling infectious diseases like H1N1. The integration of mathematical modeling and optimal control provides a valuable framework for public health decision-making and epidemic response planning.