EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS 2025, Vol. 18, Issue 4, Article Number 6254 ISSN 1307-5543 – ejpam.com Published by New York Business Global Mathematical Modeling of SEIRV Epidemic Model With Delay and Optimal Control of Holling Type II Shumaila Irum1, Anwar Zeb1,∗, Ahmed A. Mohsen2,3, Thoraya N. Alharthi4, Ilyas Khan5,6,7, Osama Oqilat8, Wei Sin Koh9 1 Department of Mathematics, COMSATS University Islamabad, Abbottabad Campus, Pakistan 2 Department of Mathematics, Open Education College, Iraq 3 Department of Mathematics, College of Education for Pure Science (Ibn Al-Haitham), University of Baghdad, Baghdad, Iraq 4 Department of Mathematics, College of Science, University of Bisha, P.O. Box 551, Bisha 61922, Saudi Arabia 5 Department of Mathematical Sciences, Saveetha School of Engineering, SIMATS, Chennai, Tamil Nadu, India 6 Hourani Center for Applied Scientific Research, Al-Ahliyya Amman University, Amman, Jordan 7 Department of Mathematics, College of Science Al-Zulfi, Majmaah University, Al-Majmaah 11952, Saudi Arabia 8 Hourani Center for Applied Scientific Research, Department of Basic Sciences, Faculty of Arts and Science, Al-Ahliyya Amman University, Amman, Jordan 9 INTI International University, Persiaran Perdana BBN Putra Nilai, 71800 Nilai, Negeri Sembilan, Malaysia Abstract. This work proposes a mathematical model of the SEIRV epidemic framework incorporating time delays and optimal control strategies based on Holling Type II functional responses. The SEIRV model, which includes compartments for Susceptible, Exposed, Infectious, Recovered, and Vaccinated in- dividuals, is extended to account for delays in the transmission and vaccination processes. The stability and behavior of the model are analyzed using differential equations and delay differential equations. Op- timal control technique is applied to minimize the spread of infection and optimize vaccination strategies, considering resource limitations and practical constraints. Numerical simulations demonstrate the effec- tiveness of the proposed control strategies in reducing infection rates and achieving disease eradication. The findings contribute to the understanding of epidemic dynamics and provide valuable information for public health policy and intervention planning. 2020 Mathematics Subject Classifications: 92D30, 49N90, 34K20 Key Words and Phrases: Diseases, seirv epidemic model, time delay, optimal control ∗Corresponding author. DOI: https://doi.org/10.29020/nybg.ejpam.v18i4.6254 Email addresses: shumaila@cuiatd.edu.pk (S. Irum), anwar@cuiatd.edu.pk (A. Zeb), aamuhseen@gmail.com (A. A. Mohsen), talhrthe@ub.edu.sa (T. N. Alharthi), i.said@mu.edu.sa (I. Khan), weisin.koh@newinti.edu.my (W. S. Koh) https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0) S. Irum et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6254 2 of 23 1. Introduction Throughout the last several decades, infectious diseases have posed significant chal- lenges to public health worldwide. Despite advances in medical science, diseases such as malaria and tuberculosis continue to affect millions. Although numerous studies have been conducted on prevention and treatment, a significant knowledge gap persists on the mechanisms of disease spread under varying environmental conditions. This research aims to develop a comprehensive model for predicting infectious disease outbreaks. The study tells us about the influence of climate variations on disease transmission. And also tells the role of human behaviors in spreading infections. This research aims to contribute to the development of more effective disease control strategies. However, this study is limited to data from specific regions and may not be universally applicable. This paper argues that incorporating environmental factors into disease prediction models significantly improves their accuracy. Mathematical modeling has proven indispensable in understanding and containing in- fectious diseases. The SEIRV model, which includes compartments for susceptible, ex- posed, infectious, recovered, and vaccinated individuals, provides a comprehensive frame- work for studying disease dynamics. These models were also used during the COVID-19 pandemic [1–6]. The classical SEIR (Susceptible–Exposed–Infectious–Recovered) model has been widely used to study the dynamics of infectious diseases. However, the emergence of vaccination as a pivotal public health strategy necessitates the inclusion of a Vaccinated (V) compartment, resulting in the extended SEIRV framework. The motivation for this work stems from the need to understand how vaccination policies influence disease spread, especially in realistic settings that involve time delays, control strategies, and nonlinear incidence rates. The novelty of this study lies in incorporating a time-delayed transmission dynamic and an optimal control approach in the SEIRV model, enabling the analysis of how delays in exposure or treatment impact epidemic control. Furthermore, the model accommodates a Holling Type II incidence function, which better captures the saturation effects in disease transmission, reflecting the limitations of the behavioral and medical response. This extended SEIRV model is applicable to the evaluation and design of ef- fective vaccination campaigns and public health interventions during outbreaks such as COVID-19, measles, or influenza. It can assist policymakers in determining optimal re- source allocation for vaccination and treatment, timing of intervention strategies, and understanding the long-term impact of delayed responses. The model is particularly rele- vant for guiding decision making in regions with constrained healthcare infrastructure and variable public response [7, 8]. However, existing models often overlook the impact of de- lays and optimal control mechanisms. This study aims to develop and analyze a delayed SEIRV epidemic model incorporating optimal control with a Holling-type II functional response. By incorporating delays, we seek to understand their effect on disease spread and identify optimal control strategies to mitigate outbreaks. This study addressed how do delays influence the transmission dynamics of the disease and what are the most ef- fective control strategies under different scenarios. The scope of this research is restricted to theoretical analysis and simulations based on specific parameter values. This paper S. Irum et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6254 3 of 23 argues that incorporating delays and optimal control into the SEIRV model significantly enhances its predictive accuracy and effectiveness in disease management. Mathematical modeling using the SEIRV epidemic framework (especially with its com- plexities like delays and the quest for optimal control) is far more than an academic exer- cise. It stands as a guiding light in the global fight against infectious diseases. In exploring this model’s depth, we are not merely working with equations or computations but we are developing powerful tools to forecast, manage, and ultimately overcome epidemics. Diving into the SEIRV epidemic model with delay and Holling type delay optimal con- trol, we are not just working with numbers and equations we are forging advanced tools for public health. This model is not merely theoretical it is a blueprint for predicting and controlling outbreaks, potentially saving countless lives. Stay motivated, as our contri- butions could well be the key to unlocking the mysteries of infectious disease dynamics. This model helps in understanding the spread and control of infectious diseases through the incorporation of various factors such as susceptible, exposed, infected, recovered, and vaccinated individuals. The delay aspect accounts for the incubation period of the disease, making the model more realistic [9, 10]. After reading the paper [11], we were inspired to extend the research. We observed that incorporating time delay and optimal control had not yet been explored in this context, so we introduced these elements as novel contributions. By applying optimal control theory, the model identifies the most effective strategies to curb the spread of infectious diseases while accounting for resource limitations. This approach is essential for effective public health planning and intervention [12]. To better capture the dynamics between the disease and the population, the model utilizes the Holling Type II functional response, which offers a more realistic depiction of how infections progress and impact the population [13]. The insights gained from this model can aid policymakers in determining optimal strategies for controlling epidemics, such as vaccination, quarantine, and other public health measures [12, 14]. This framework also lays the groundwork for future epidemiological research, enabling scientists to develop innovative approaches for disease prevention and health improvement [9]. Ultimately, studying this model empowers researchers to contribute significantly to efforts aimed at controlling and preventing the spread of infectious diseases, thereby saving lives and enhancing public health outcomes [12]. This paper is structured as follows: In Section 2, we introduce the mathematical model with time delay and its solution. In Section 3, we formulate the mathematical epidemic model SEIRV with optimal control and its solution are find. In Section 4, we provide a numerical method alongside the relevant simulation results. In summery, Section 5 compiles the final conclusion and discussion. 2. Formulation of SEIRV Model with Time Delay In this paper, we consider the model presented in [11], which was discussed with- out delay and optimal control. In their model, they have divided the total population N(t) at time t > 0 into five compartments. These are S(t) susceptible populations that are at risk of infection but have not been exposed to the infected yet, exposed E(t) de- S. Irum et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6254 4 of 23 notes the exposed population that has come into contact with the infected people but is not infectious. These people are in a latent state, I(t) refers to infectious people that spread the infection to other people, R(t) indicates recovered people who are no longer infected by the infection. ,and V (t) stands for the vaccinated population. Therefore N(t)=S(t)+E(t)+I(t)+R(t)+V (t). Each compartment presents a unique scenarios.For example the population of S(t) compartment increased by recruitment rate ∆ and passed away naturally at a rate of µ, vaccination at a rate of γ2 and transmission at a rate of β with term βSI 1+αI , where α being the saturated rate.The population of E(t) compartment passed away naturally at a rate of µ, with latency period σ . The population of I(t) compartment decreased by natural death at a rate of µ, with speed of recovery rate γ1 . The population of the R(t) compartment decreased of natural causes at a rate of µ, and vaccinated at a rate of γ3. Similarly, the population of the V (t) compartment passed away naturally at a rate of µ. However, delay and optimal control play a crucial role in the dynamics of epidemic models. Therefore, in this study, we investigate a system of differential equations incorporating time delays to capture these essential effects more accurately. Here τ > 0 represents the response time of disease, which follows a similar argument in [15, 16] in the model presented in [11]. The model is as follows for the SEIRV framework. dS(t) dt = ∆− βS(t)I(t) 1 + αI(t) − µS(t)− γ2S(t), dE(t) dt = βS(t)I(t) 1 + αI(t) − σE(t− τ)− µE(t), dI(t) dt = σE(t− τ)− γ1I(t)− µI(t), dR(t) dt = γ1I(t)− γ3R(t)− µR(t), dV (t) dt = γ2S(t) + γ3R(t)− µV (t). (1) The initial conditions of the above system as follows:{ S(0) = S(θ) > 0, E(0) = E(θ) ≥ 0, I(0) = I(θ) ≥ 0, R(0) = R(θ) ≥ 0, V (0) = V (θ) ≥ 0, (2) where, θ ∈ [−τ, 0] is assumed. The proposed model has two equilibrium points that are disease-free equilibrium point: E0(S0.E0, I0, R0, V0) = E0( Λ µ+γ2 , 0, 0, 0, γ2Λ µ(µ+γ2) ), and endemic equilibrium point: E∗(S∗.E∗, I∗, R∗, V∗) = E∗( (σ+µ)(γ1+µ)(1+αI∗) βσ , (γ1+µ)(γ2+µ)(R0−1) σ(β+α(µ+γ2)) , (µ+γ2)(R0−1) β+α(µ+γ2) , γ1(R0−1) β+α(µ+γ2) , γ2(S∗+R∗) µ ). Reproductive number is give as: R0 = βΛσ (σ+µ)(γ1+µ)(γ2+µ) . On the same way, the first three equations of system (1) are independent of fourth and fifth variables R and V , enabling us to analyze the system by examining the subsystem, S. Irum et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6254 5 of 23 as dS(t) dt = ∆− βS(t)I(t) 1 + αI(t) − µS(t)− γ2S(t), dE(t) dt = βS(t)I(t) 1 + αI(t) − σE(t− τ)− µE(t), dI(t) dt = σE(t− τ)− γ1I(t)− µI(t). (3) Theorem 1. The solutions of system (3) remain bounded and non-negative for all time t ≥ 0, provided that the initial conditions are non-negative. That is, every solution remains in a positively invariant set Ω ⊂ R3 +, and (S(t), E(t), I(t)) ∈ R3 + for all t ≥ 0. Proof. Suppose N(t)=S(t)+E(t)+I(t), then dN dt = d dt (S(t)) + d dt (E(t)) + d dt (I(t)), dN dt = ∆− µ (S(t) + E(t) + I(t))− γ2S(t)− γ1I(t), dN(t) dt = ∆− µN(t)− γ2S(t)− γ1I(t) ≤ ∆− µN(t), d dt ( Neµt ) ≤ ∆eµt, implies that N(t) ≤ ∆ µ + Be−µ(t). Thus, 0 ≤ N(t) ≤ ∆ µ as t→ ∞. Consequently, the dynamic of the model can be analyzed within the region: Ω = {(S(t), E(t), I(t)) ∈ R3 + : N ≤ ∆ µ }, then the set Ω is positive invariant with respect system (3) and hence it is closed. 2.1. Linearization To linearized the system (3), we apply Taylor series. For simplicity, we represent this process as follows: dS dt = f(S(t), E(t), I(t)), where f(S(t), E(t), I(t)) = ∆− βS(t)I(t) 1 + αI(t) − µS(t)− γ2S(t) dE dt = g(S(t), E(t), I(t)), S. Irum et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6254 6 of 23 where g(S(t), E(t), I(t)) = βS(t)I(t) 1 + αI(t) − σE(t− τ)− µE(t) dI dt = h(S(t), E(t), I(t)) where h(S(t), E(t), I(t)) = σE(t− τ)− γ1I(t)− µI(t). We define S − S∗ = u, E − E∗ = v, I − I∗ = w. Now expansion of model (3) around the point (S∗, E∗, I∗) is given by: u̇ = u ( − βI∗ 1 + αI∗ − µ− γ2 ) + v(0) + w ( − βS∗ (1 + αI∗)2 ) , similarly v̇ = u ( βI∗ 1 + αI∗ − ) + v(−µ− σe−λτ ) + w ( βS∗ (1 + αI∗)2 ) , ẇ = u(0) + v(−σ)e−λτ + w(−µ− γ1). 2.2. Matrix Representation of the Linearized System The linearized system for (3) is expressed in matrix form as follows: u̇ v̇ ẇ  =  ∂f ∂S ∂f ∂E ∂f ∂I ∂g ∂S ∂g ∂E ∂g ∂I ∂h ∂S ∂h ∂E ∂g ∂I  u v w  . Alternatively, in simplified notation: Ẋ = JX, where J represents the Jacobian matrix, composed of the derivatives of each term on the right- hand side of the model’s equations, namely (S,E, I). So, with linearization the system obtains the following form:  du dt = a11u+ a13w, dv dt = a21u+ ua22v + b22u ( t− τ ) + a23w, dw dt = b32v ( t− τ ) + a33w, (4) where a11 = − βI∗ 1 + αI∗ − µ− γ2, a13 = − βS∗ (1 + αI∗)2 , a21 = βI∗ 1 + αI∗ − µ− γ2, S. Irum et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6254 7 of 23 a22 = −µ, a23 = βS∗ (1 + αI∗)2 , a33 = −γ1 − µ, b22 = −σ b32 = σ. Thus, the Jacobian of system (4) is J = a11 0 a13 a21 a22 + b22e −λτ a23 0 b32e −λτ a33  . For eigen values, solve |λI − J | = 0,∣∣∣∣∣∣ λ− a11 0 −a13 −a21 λ− (a22 + b22e −λτ ) −a23 0 −b32e−λτ λ− a33 ∣∣∣∣∣∣ = 0. After solving the determinant we acquire the characteristic equation: λ3 + p1λ 2 + p2λ+ p3 + ( q1λ 2 + q2λ+ q3 ) e−λτ = 0. (5) Here p1 = −a11 + a22 + a33, p2 = a11a22 + a11a33 + a22a33, p3 = −a11a22a33, q1 = −b22, q2 = a33b22 − a23b32 + a11b22, q3 = a13a23b32 + a13a22b32d6 − a11b22a33 When τ = 0, then λ3 + (p1 + q1)λ 2 + (p2 + q2)λ+ (p3 + q3) = 0. (6) Let suppose: D2 = p1 + q1, D1 = p2 + q2, D0 = p3 + q3. Equation (6) become: λ3 +D2λ 2 +D1λ+D0 = 0. (7) By using the Hurwitz criteria, E∗ exhibits local asymptotic stability at τ = 0, under the following conditions: S. Irum et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6254 8 of 23 det1 = D2 > 0, (8) det2 = ( D2 1 D0 D1 ) > 0, (9) det3 = D2 1 0 D0 D1 D2 0 0 D0  > 0. (10) Since the model parameters are all positive, then it is always asymptotically stable, which shows all roots lie on the left half plane. In the same way, when τ > 0, we consider the characteristics equation provided in eq (5) and, setting λ = ιω, we have λ3 + p1λ 2 + p2λ+ p3 + ( q1λ 2 + q2λ+ q3 ) e−ιωτ = 0. (11) By separated the real and imaginary parts of above equation we have (q3 − q1ω 2) cos(ωτ) + q2ω sin(ωτ) = p1ω 2 − p3, (12) q2ω cos(ωτ)− (q3 − q1ω 2) sin(ωτ) = ω3 − p2ω (13) which leads to ω6 + ( p21 + 2p2 − q21 ) ω4 + (p22 − 2p1p3 − q22 + 2q1q3)ω 2 + (p23 + q23) = 0. (14) If t = ω2, then above equation becomes t3 + s1t 2 + s2t+ s3 = 0. (15) where  s1 = ( p21 + 2p2 − q21 ) , s2 = (p22 − 2p1p3 − q22 + 2q1q3), s3 = (p23 + q23), Thus by substituting s1, s2, s3 in above equation, we have negative roots which implies that it is asymptotically stable. 3. Epidemic Model SEIRV-Formulation with optimal control In this section, we examine control strategy of model model presented in [11]. The population of E(t) reduced with control u1 as quarantine. The population of I(t) compartment decreased with control u2 as treatment . The population of the R(t) compartment is increased by control u1 of the exposed class and cotrol u2 of infected class. Thus the optimal control model can be expressed as: S. Irum et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6254 9 of 23 dS(t) dt = ∆− βS(t)I(t) 1 + αI(t) − µS(t)− γ2S(t), dE(t) dt = βS(t)I(t) 1 + αI(t) − σE(t)− µE(t)− u1E(t), dI(t) dt = σE(t)− γ1I(t)− µI(t)− u2I(t), dR(t) dt = γ1I(t)− γ3R(t)− µR(t) + u1E(t) + u2I(t), dV (t) dt = γ2S(t) + γ3R(t)− µV (t). (16) Our aim to decrease the number of Exposed and Infected individuals while increase the count of recovering throughout the epidemic. Mathematically, for a fixed terminal time ti, the problem is to minimize the objective function, which is given as: J(u1, u2) = ∫ T 0 [ A1S(t) +A2E(t) +A3I(t)−A4R(t) +A5V (t) + A6 2 u21(t) + A7 2 u22 ] dt. (17) In the above expression, Ai ≥ 0 (for i = 1..7) represents the scaling factors that adjust the contribution of each term in the functional. The objective is to find the optimal control pairs u∗1,u ∗ 2 such that: J(u∗1, u∗2) = min{J(u1, u2) : u1, u2 ∈ U}. (18) In this setting, U denotes the pair of feasible controls, given by: U = {u = ui(t) : 0 ≤ ui(t) ≤ umax i ≤ 1, t ∈ [0, T ], i = 1, 2}. (19) Here u is Lebesgue measurable. The flow chart of the model (16) is presented in Fig. (1). S. Irum et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6254 10 of 23 Figure 1: Flow chart of model (16 ) Now, we first find the kernels of control model. For this purpose, we introduce the functions Ki where i ∈ {1, . . . , 5}, while assuming L[0, 1]. We describe these functions in the following manner: K1(t, S) = ∆− βS(t)I(t) 1 + αI(t) − (γ2 + µ)S(t), K2(t, E) = βS(t)I(t) 1 + αI(t) − (σ + µ)E(t)− u1E(t), K3(t, I) = σE(t)− (γ1 + µ)I(t)− u2I(t) K4(t, R) = γ1I(t)− γ3R(t)− µR(t) + u1E(t) + u2I(t), K5(t, V ) = γ2S(t) + γ3R(t)− µV (t) Theorem 2. Given the assumptions stated above, the kernels Ki satisfy the Lipschitz condition and are contractions if and only if the Lipschitz constants ψj < 1 where j ∈ {1, . . . , 5}, for each i ∈ {1, . . . , 5}. Proof. For proof of this theorem we are employing certain conditions on the kernels: Lipschitz Condition for K1(S): ∥K1(S)−K1(S ∗)∥ = ∥∥∥∥∆− βSI 1 + αI − (γ2 + µ)S − ( ∆− βS∗I 1 + αI − (γ2 + µ)S∗ ) ∥, ≤ ∥∥∥∥ βI 1 + αI + (γ2 + µ) ∥∥∥∥ ∥S − S∗∥, = ψ1∥S − S∗∥, S. Irum et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6254 11 of 23 where ψ1 = ∥∥∥∥ βI 1 + αI + (γ2 + µ) ∥∥∥∥ . Lipschitz Condition for K2(E): ∥K2(E)−K2(E ∗)∥ = ∥∥∥∥ βSI 1 + αI − (σ + µ)E − u1E − [ βSI 1 + αI − (σ + µ)E∗ − u1E ∗ ] ∥, ≤ ∥∥∥∥ βSI 1 + αI + (σ + µ) + u1 ∥∥∥∥ ∥E − E∗∥, = ψ2∥E − E∗∥, where ψ2 = ∥∥∥∥ βSI 1 + αI + (σ + µ) + u1 ∥∥∥∥ . Lipschitz Condition for K3(I): ∥K3(I)−K3(I ∗)∥ = ∥σE − (γ1 + µ)I − u2I − [σE − (γ1 + µ)I∗ − u2I ∗]∥ , ≤ ∥(γ1 + µ) + u2∥ ∥I − I∗∥, = ψ3∥I − I∗∥, where ψ3 = ∥(γ1 + µ) + u2∥ . Lipschitz Condition for K4(R): ∥K4(R)−K4(R ∗)∥ = ∥γ1I − µR+ u1E + u2I − [γ1I − µR∗ + u1E + u2I]∥ , ≤ ∥µ∥ ∥R−R∗∥, = ψ3∥R−R∗∥, where ψ4 = ∥µ∥ . Lipschitz Condition for K5(V ): ∥K5(V )−K5(V ∗)∥ = ∥γ2S − γ3R− µV − [γ2S − γ3R− µV ∗]∥ , ≤ ∥µ∥ ∥V − V ∗∥, = ψ5∥V − V ∗∥, where ψ5 = ∥µ∥ . Therefore, the analysis of the Lipschitz conditions demonstrates that the kernels satisfy the criteria for being contractions. S. Irum et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6254 12 of 23 3.1. Existence of an optimal control pair This section presents a finding on whether an optimal control exists within the model we proposed. Theorem 3. Consider the control problem with the system (16). There exists an optimal control pair u∗i ∈ U such that: J(u∗i ) = min{J(ui) : ui ∈ U, i = 1, 2}. (20) Proof. To establish the existence of an optimal control, we utilize a result by Fleming and Rishel [17]. The following conditions are be checked: (i) The set of controls and the corresponding state variables must be nonempty. The existence of solution for system (16) is proved using result by Lukes [18]. (ii) The set U of control is convex and closed. (iii) The right-hand side of the system (16) is bounded by a linear function of the state and control variables. (iv) The integrand of the objective functional is convex with respect to U. (v) There exist constants m1,m2 > 0 and ρ > 1 such that the integrand L(S,E, I, u1, u2) of the objective functional satisfies: L(S,E, I, u1, u2) ≥ m2 +m1(|u|2)ρ/2. Therefore, using these conditions, we confirm that an optimal control exists. 3.2. Characterization of the Optimal Control For deriving the necessary conditions for a control to be optimal, Pontryagin’s Maximum Principle is applied. Based on this principle we transform our system into a minimization problem for the Hamiltonian H, which is given by: H = A1S(t) +A2E(t) +A3I(t)−A4R(t) +A5V (t) + A6 2 u21(t) + A7 2 u22 + 5∑ i=1 λiki. The right-hand terms of differential equations are represented by ki corresponding to the i-th state variable. Theorem 4. Given an optimal control u∗ = (u∗1, u ∗ 2) ∈ U and corresponding variables S, E, I, R, and V , there exist adjoint variables λ1, ..., λ4, and λ5 that must fulfill the following adjoint equations: λ̇1 = −A1 + λ1 ( βI 1 + αI + (γ2 + µ) ) − λ2 ( βI 1 + αI ) − λ5γ2, λ̇2 = −A2 + λ2 (σ + µ+ u1) + σλ3 − u1λ4, λ̇3 = −A3 + λ1 ( βS (1 + αI)2 ) − λ2 ( βS (1 + αI)2 ) + λ3(γ1 + µ+ u2)− λ4(γ1 + u2) S. Irum et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6254 13 of 23 λ̇4 = A4 + µλ4 + γ2λ5, λ̇5 = A5 + µλ5, following the transversality conditions: λ1(T ) = λ2(T ) = λ3(T ) = λ4(T ) = λ5(T ) = 0. Additionally, the optimal control u∗1(t) and u ∗ 2(t) is characterized by: u∗1(t) = min ( 1,max ( 0, (λ4 − λ2) A1 E )) , u∗2(t) = min ( 1,max ( 0, (λ4 − λ3) A2 I )) . Using Pontryagin’s maximum principle adjoint equations and transversality conditions are derived. Specifically: (i) Adjoint Equations: The adjoint equations take the following form: λ̇1 = −∂H ∂S , λ̇2 = −∂H ∂E , λ̇3 = −∂H ∂I , λ̇4 = −∂H ∂R , λ̇5 = −∂H ∂V . (ii) Optimal Control: The optimal control u∗1 and u∗2 can be determined from the optimality condition, which involves partially differentiating the Hamiltonian in terms of ui to zero: ∂H ∂u1 = 0. ∂H ∂u2 = 0. This yields: −A1u ∗ 1 + (λ4 − λ2)(E) = 0. −A2u ∗ 2 + (λ4 − λ3)(I) = 0. S. Irum et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6254 14 of 23 Given the bounds on the control ui specified by the set U , we express the optimal control u∗1 and u∗2 in the following form: u∗1(t) = min ( 1,max ( 0, (λ4 − λ2) A1 E )) , u∗2(t) = min ( 1,max ( 0, (λ4 − λ3) A2 I )) . 4. Numerical Simulation To discuss the dynamics of disease spread and study the impact of both time delay and control strategies, we performed a numerical simulation of the solution trajectories of the model described in model (3), where classes (S,I,E) are independent of ( R, V) and the adjoint model in Theorem (4.2). The parameter values model simulations are shown in Table 1. Table 1. Parameter values of model (3). Parameters Value ∆ 0.87 β 0.01 α 0.1 µ 0.1 γ2 = γ3 0.1 γ1 0.01 σ 0.2 The simulation of the model (3) illustrated in Figure 1, shows the dynamic behavior of the model (3) that approaches the extinction of the disease. S. Irum et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6254 15 of 23 0 0.2 0.2 5 I( t) E(t) 0 S(t) 0.4 -0.2 0 0 500 1000 1500 2000 t 0 2 4 6 S (t ) 0 500 1000 1500 2000 t 0 0.05 0.1 0.15 0.2 E (t ) 0 500 1000 1500 2000 t 0 0.1 0.2 0.3 0.4 I( t) Stable point Figure 2: The trajectory of the system (3) using the parameters set in Table 1. Now, we show the effect of increasing the infection rate on the dynamics of the model (3), by keeping the values of all parameters in Table 1 except the value of β we change it from 0.01 to 0.1 and τ < 18.5, we get the solution of model 1, which tends to the endemic point; see Figure 2. 0 2 2 4 I( t) E(t) 1 S(t) 4 2 0 0 0 500 1000 1500 2000 t 0 1 2 3 4 S (t ) 0 500 1000 1500 2000 t 0 0.5 1 1.5 2 E (t ) 0 500 1000 1500 2000 t 0 1 2 3 I( t) Stable point Figure 3: The trajectory of the system (3) using the parameters set in Table 1 with β = 0.1. S. Irum et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6254 16 of 23 In addition, we also need to know the impact of time delay on the solution of the model (3). In Figure 3, if we choose the value of the parameter τ = τ0 = 18.5. The figure tells us loss of stability of the endemic point and change of the system’s behavior to the periodic state (occurring Hopf bifurcation). However, both figures (4) and (5) tell us the high periodic state when we increase the time delay such as τ = 20, 25, respectively. 0 4 4 I( t) E(t) 2 S(t) 5 2 0 0 0 500 1000 1500 2000 t 0 1 2 3 4 S (t ) 0 500 1000 1500 2000 t 0 1 2 3 E (t ) 0 500 1000 1500 2000 t 0 2 4 6 I( t) Figure 4: The trajectory of the model (3) using the parameters set in Table 1 with τ = 18.5. S. Irum et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6254 17 of 23 0 4 5 4 I( t) E(t) 2 S(t) 10 2 0 0 0 500 1000 1500 2000 t 0 1 2 3 4 S (t ) 0 500 1000 1500 2000 t 0 1 2 3 4 E (t ) 0 500 1000 1500 2000 t 0 2 4 6 I( t) Figure 5: The trajectory of the model (3) using the parameters set in Table 1 with τ = 20. 0 4 5 4 I( t) E(t) 2 S(t) 10 2 0 0 0 500 1000 1500 2000 t 0 1 2 3 4 S (t ) 0 500 1000 1500 2000 t 0 1 2 3 4 E (t ) 0 500 1000 1500 2000 t 0 2 4 6 8 I( t) Figure 6: The trajectory of the system (3) using the parameters set in Table 1 with τ = 25. For the simulation analysis to the each control strategies are quarantine and treatment and represented by the parameters u1 and u2 respectively. The four scenarios for numerical simulation of the adjoint model are given below: (i) Scenario A: There is no quarantine and no treatment available for the infected (u1 = u2 = 0). S. Irum et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6254 18 of 23 (ii) Scenario B: Only quarantine is implemented (u1 ̸= 0, u2 = 0). (iii) Scenario C: Only treatment for the infected is provided (u1 = 0, u2 ̸= 0). (iv) Scenario D: A combination of both quarantine and treatment is implemented (u1 ̸= 0, u2 ̸= 0). It is clear that, figure 6, shows that the impact of both scenarios (A) and (B) and this figure tell us when we increasing the quarantine rate (i.e. u1 ≥ 0.5) with keeping the other parameters in table 1 with β = 0.1 and τ < 18.5, we get the dynamical behavior of system (3) tend towards disease extinction point. 0 100 200 300 400 Days 0 2 4 6 S µ 1 =µ 2 =0 µ 1 =0.1,µ 2 =0 µ 1 =0.2,µ 2 =0 µ 1 =0.5,µ 2 =0 0 100 200 300 400 Days 0 0.5 1 E 0 100 200 300 400 Days 0 0.5 1 1.5 2 I 0 1 1 5 I E 0.5 S 2 0 0 Figure 7: The trajectory of the model (3) with µ2 = 0 and different values of µ1. Thus, figure 7, presents that the impact of both scenarios (A) and (C) and this figure tell us when we increasing the treatment rate (i.e. u2 ≥ 0.2) with keeping the other parameters in table 1 with β = 0.1 and τ < 18.5, we get the dynamical behavior of model (3) tend towards disease extinction point. S. Irum et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6254 19 of 23 0 100 200 300 400 Days 0 2 4 6 S µ 1 =0,µ 2 =0 µ 1 =0,µ 2 =0.1 µ 1 =0,µ 2 =0.2 µ 1 =0,µ 2 =0.5 0 100 200 300 400 Days 0 0.5 1 E 0 100 200 300 400 Days 0 0.5 1 1.5 2 I 0 1 1 5 I E 0.5 S 2 0 0 Figure 8: The trajectory of the system (3) with µ1 = 0 and different values of µ2. Finally, as depicted in Figure 7, scenario (D) leads to a more rapid progression toward dis- ease extinction compared to scenarios (B) and (C), highlighting the enhanced effectiveness of the interventions assumed in this scenario. 0 100 200 300 400 Days 0 2 4 6 S µ 1 =0,µ 2 =0 µ 1 =0.1,µ 2 =0.1 µ 1 =0.2,µ 2 =0.2 0 100 200 300 400 Days 0 0.5 1 E 0 100 200 300 400 Days 0 0.5 1 1.5 2 I 0 1 1 6 I E 0.5 4 S 2 2 0 0 Figure 9: The trajectory of the system (3) with µ1 = 0.2 and µ2 = 0.2. S. Irum et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6254 20 of 23 0 5 10 15 20 25 30 0 100 200 300 400 500 600 700 800 900 Time (days) P o p u la tio n SEIRV Model with Control Susceptible Exposed Infected Recovered Vaccinated Figure 10: Control model (16). The Figure 10, demonstrates the impact of optimal control strategies on disease progression. The interventions reduce the infection peak and overall number of cases, showcasing the effective- ness of control measures like vaccination, treatment, or public awareness. 0 5 10 15 20 25 30 0 100 200 300 400 500 600 700 800 Time (days) P o p u la tio n Simple SEIRV Model Susceptible Exposed Infected Recovered Vaccinated Figure 11: Without delay and control for model (16). S. Irum et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6254 21 of 23 The Figure 11, illustrates the baseline epidemic dynamics without any delay or control in- terventions. The disease spreads naturally according to the model’s parameters, showing typical infection peaks and declines. 0 5 10 15 20 25 30 0 100 200 300 400 500 600 700 800 Time (days) P o p u la tio n SEIRV Model with Delay Susceptible Exposed Infected Recovered Veccinated Figure 12: Without delay for model (16). The Figure 12, reflects the influence of time delays in the epidemic model, such as delayed response to infection or incubation periods. The delay affects the timing and height of the epidemic peak, potentially causing oscillatory or more prolonged outbreaks. S. Irum et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6254 22 of 23 5. Conclusion and Discussion In recent paper, we formulate an epidemic model with time delay, consisting of five state vari- ables and nine parameters. Based on the analytical results, two equilibrium points are identified: disease-free equilibrium and endemic equilibrium. The results show that the system exhibits stable behavior in disease-free equilibrium for all time delay values, given the parameter values presented in Table 1 (see Figure 2). However, when the infection rate increases, the system transitions to a stable endemic equilibrium, provided that the time delay is below 18.5 (see figure 3). This stability is lost when the time delay exceeds this threshold (see Figures 4-6). Furthermore, control variables are incorporated into the model, namely quarantine and treatment. Pontryagin’s Maximum Prin- ciple is applied to solve the optimal control problem. Numerical simulations with specified weights indicate that the optimal strategy to control the disease spread is the combination of both control measures (see figure 7-9). Further direction is to apply the optimal control over a delayed system. References [1] Elham Taghizadeh and Ali Mohammad-Djafari. Seir modeling, simulation, parameter esti- mation, and their application for covid-19 epidemic prediction. Phys. Sci. Forum, 5(1):18, 2022. [2] R. Paulo and J.P.Z. Zingano. A matlab code to compute reproduction numbers with appli- cations to the covid-19 outbreak, 2020. arXiv:2006.13752. [3] S. He, Y. Peng, and K. Sun. Seir modeling of the covid-19. Nonlinear Dyn, 101:1667–1680, 2020. [4] J.P. Gleeson, T.B. Murphy, J.D. Brien, and D.J. Sullivan. A population-level seir model for covid-19 scenarios, 2020. Technical Note of the Irish Epidemiological Modelling Advisory Group to NPHET, Department of Health, Government of Ireland. [5] Ahmed A. Mohsen, Hassan F. AL-Husseiny, and Raid Kamel Naji. The dynamics of corona virus pandemic disease model in the existence of a curfew strategy. Journal of Interdisciplinary Mathematics, 2022. [6] Hassan F. AL-Husseiny, Nidhal F. Ali, and Ahmed A. Mohsen. The effect of epidemic disease outbreaks on the dynamic behavior of a prey-predator model with holling type ii functional response. Commun. Math. Biol. Neurosci., 2021. Article ID 72. [7] A. Acharya, S. Paul, M.A. Biswas, A. Mahata, S. Mukherjee, and B. Roy. Study of seirv epidemic model in imprecise environment. In Advances in Intelligent System and Computing, ICMMCS, pages 371–380. 2023. [8] S. Paul, A. Acharya, M.A. Biswas, A. Mahata, S. Mukherjee, P.C. Mali, and B. Roy. The scenario of covid 19 pandemic in brazil using seir epidemic model. In Advances in Intelligent System and Computing, ICMMCS, pages 419–426. 2023. [9] Ping Yan and Shengqiang Li. Seir epidimic model with delay. ANZIAM J, 48:119–134, 2006. [10] Ahmed A. Mohsen and Raid Kamel Naji. Stability and bifurcation of a delay cancer model in the polluted environment. Advances in Systems Science and Applications, 22(3):1–17, 2022. [11] Sajal and Fahad Mostafa. Seirv-holling type-ii: Implementing holling type functional response in covid-19 disease modeling. SSRN. https://ssrn.com/abstract=4953603. [12] Animesh, Subrata Paul, Supriya Mukherjee, Meghadri Das, and Banamali Roy. Dynamics of caputo fractional order seirv epidimic model with optimal control and stability analysis. International Journal of Applied and Computational Mathematics, 8:28, 2022. [13] Sarah A. Al-Sheikh. Modeling and analysis of seir epidimic model with limited resource of S. Irum et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6254 23 of 23 treatment. Global Journal of Science Frontier Research Mathematics and Decision Sciences. Print ISSN: 0975-5896. [14] Ahmed A. Mohsen and Raid Kamel Naji. Dynamical analysis within-host and between- host for hiv/aids with the application of optimal control strategy. Iraqi Journal of Science, 61(5):1173–1189, 2020. [15] Oluwatosin Babasola, Oshinubi Kayode, Olumuyiwa James Peter, Faithful Chiagoziem On- wuegbuche, and Festus Abiodun Oguntolu. Time-delayed modelling of the covid-19 dynamics with a convex incidence rate. Informatics in Medicine Unlocked, 35:101124, 2022. [16] S. Tipsri and W. Chinviriyasit. The effect of time delay on the dynamics of an seir model with nonlinear incidence. Chaos Solitons Fractals, 75:153–172, 2015. [17] W.H. Fleming and R.W. Rishel. Deterministic and Stochastic Optimal Control. Springer, New York, NY, USA, 1975. [18] D.L. Lukes. Differential Equations: Classical to Controlled. Academic Press, New York, 1982.