EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS 2025, Vol. 18, Issue 4, Article Number 6775 ISSN 1307-5543 – ejpam.com Published by New York Business Global Backward Bifurcation and Optimal Control in an Age-Structured SVI Epidemic Model J. Leo Amalraj1, Pankaj Shukla2, Lakshmana Phaneendra Maguluri3, Siriluk Donganont4,∗, N. Avinash5, Vediyappan Govindan6 1 Department of Mathematics, RMK College of Engineering and Technology, Puduvoyal, Thiruvallur, Tamil Nadu, India 2 Department of Mathematics, School of Advanced Sciences, Vellore Institute of Technology, Chennai, Tamil Nadu, India 3 Department of Computer Science and Engineering, Koneru Lakshmaiah Education Foundation, Vaddeswaram, Guntur AP, India 4 School of Science, University of Phayao, Phayao 56000, Thailand 5 Department of Mathematics, Sacred Heart College (Autonomous), Tirupattur 635 601, Tamil Nadu, India 6 Department of Mathematics, Hindustan Institute of Technology and Science, Chennai, Tamil Nadu, India Abstract. This study investigates the dynamics and control of an infectious disease using an age- structured Susceptible-Vaccinated-Infected (SVI) model, incorporating imperfect vaccination and therapeutic treatment. We identify a backward bifurcation phenomenon, driven by treatment rates, which allows disease persistence even when the basic reproduction number R0 < 1, complicating eradication efforts. To address this, we formulate an optimal control problem with time-dependent vaccination and treatment rates to minimize infected individuals and control costs. Using bifurca- tion theory and integrated semigroup methods, we establish the model’s well-posedness, equilibria stability, and bistability conditions. First-order optimality conditions are derived to characterize optimal controls. Numerical simulations, solved via the Forward-Backward Sweep method, demon- strate that combined vaccination and treatment significantly reduces disease prevalence, with early intervention being critical. These findings underscore the importance of strategic resource alloca- tion in managing complex epidemic dynamics. 2020 Mathematics Subject Classifications: 92D25, 35B32, 49J15, 65K05 Key Words and Phrases: Age-structured model, backward bifurcation, optimal control, infec- tious disease, imperfect vaccination, therapeutic treatment ∗Corresponding author. DOI: https://doi.org/10.29020/nybg.ejpam.v18i4.6775 Email addresses: leoamalraj@rmkcet.ac.in (J. Leo Amalraj), pankaj.shukla@vit.ac.in (P. Shukla), phanendra51@gmail.com (L. Phaneendra Maguluri), siriluk.pa@up.ac.th (S. Donganont), avinashprofess@gmail.com (N. Avinash), vgovindandr@gmail.com (V. Govindan) https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0) J. Leo Amalraj et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6775 2 of 32 1. Introduction Bovine tuberculosis (BTB), caused by Mycobacterium bovis, is a contagious disease affecting domestic and wild animals, including cattle, goats, sheep, badgers, deer, bison, and African buffalo. Hosts can either be reservoirs, which maintain and spread the infec- tion, or spill-over hosts, which have little impact on transmission. In buffalo herds, BTB prevalence is notably high, leading to increased mortality. Studies indicate that higher prevalence rates correspond to higher disease-related mortality in affected herds. The disease is chronic and progressive, with an unpredictable time from infection to death. This progression is influenced by factors such as immune response, stress, drought, and aging. BTB remains a significant threat to wildlife conservation and livestock health, necessitating effective control measures [1, 2]. The asymptotic behavior of epidemic models shows that disease persistence is often determined by the basic reproduction number, with forward bifurcation occurring when it exceeds one. However, recent studies highlight that backward bifurcations, caused by factors like heterogeneous susceptibility and nonlinear incidence, make disease elimination more complex. In such cases, reducing the reproduction number below one is insuffi- cient, requiring additional control efforts. Treatment is crucial for controlling diseases like measles and tuberculosis, but real-world constraints necessitate efficient resource alloca- tion. Classical models assume treatment rates proportional to infections, but communi- ties must balance cost and capacity. Ensuring optimal treatment prevents unnecessary expenses while minimizing the risk of outbreaks [3–6]. Mathematical models of infectious diseases play a crucial role in shaping public health policies aimed at minimizing or eradicating infections. Researchers have used these models to study diseases such as tuberculosis, HIV, malaria, and Ebola, helping to understand epidemic dynamics. A common criterion for disease eradication is maintaining the basic reproduction number below unity, ensuring the stability of the disease-free equilibrium. However, studies have shown that backward bifurcation can occur, allowing disease persis- tence even when the reproduction number is less than one. This challenges the traditional reliance on the reproduction number as a sole control metric. This study further investi- gates the effects of treatment saturation and infected individual isolation on the emergence of backward bifurcation [7–12]. Infectious diseases remain a major global health challenge, with illnesses like influenza, COVID-19, and tuberculosis causing millions of deaths annually. Limited healthcare re- sources in low-income regions exacerbate disease spread and mortality. Emerging epi- demics pose risks to public health and economies, highlighting the need for effective con- trol strategies. Mathematical modeling helps predict outbreaks and guide policymakers in implementing optimal disease prevention and treatment measures [13, 14]. Infectious diseases remain a major threat, with epidemic models helping to understand their spread and control. The SIQR model, incorporating quarantine, is useful for study- ing transmission dynamics [15–22]. This study analyzes the SIQR model using numerical methods like the Euler scheme, Runge-Kutta, and NSFD. Stability, bifurcation, and con- vergence are examined using the Routh-Hurwitz criterion to identify the most effective J. Leo Amalraj et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6775 3 of 32 computational approach [23–26]. Several multi-strain epidemic models extending the classical SIR framework have been analyzed using bifurcation theory to understand dengue fever dynamics. Previous studies have examined the role of primary and secondary infections in disease progression. Ana- lytical techniques, such as symbolic algebra (Maple), have been used for simple models, while numerical bifurcation tools (AUTO, MatCont) assist in analyzing complex dynamics. Stability analysis via equilibria and Lyapunov exponents helps identify chaotic behavior. These methods provide insights into parameter dependencies and long-term epidemic pat- terns [27–30]. Mathematics has been integral to biology since Fibonacci’s early population models, evolving into the formalized field of biomathematics introduced by Johannes Ranke in 1901. Mathematical biology involves theoretical analysis and modeling of biological struc- tures and behaviors, divided into stages from conceptual representation to mathematical formulation and application. Fractional calculus and mathematical operators such as Riemann-Liouville and Caputo have been applied to biological systems, including disease modeling. Cholera, caused by Vibrio cholerae, spreads through fecal-oral transmission, ex- acerbated by poor sanitation, and remains a major public health concern. Mathematical models have been instrumental in understanding cholera transmission, guiding prevention, intervention strategies, and public health policies [31–36]. Arboviral diseases, transmitted by arthropods, include Dengue, Yellow Fever, and Chikungunya, posing major public health threats. Dengue alone infects 50–100 million people annually, mainly children. Transmission dynamics depend on human behavior, mosquito activity, virus diversity, and environmental factors. Climate change and ur- banization further influence disease spread and outbreak severity. These factors impact control measures, requiring comprehensive intervention strategies [37–41]. Mathematical models of infectious diseases have been extensively studied, addressing disease spread and control. The rapid spread of COVID-19 in 2020 highlighted the im- portance of research in this field. Among various models, the SIR model remains widely studied, with research exploring its nonlinear dynamics, discretization effects, and appli- cations to HIV and variable populations. Optimal control strategies have been analyzed to enhance intervention measures. These studies contribute to improving disease modeling and informing public health policies [42–53]. Biologically, this model captures the real-world dynamics of diseases like bovine tu- berculosis (BTB) in wildlife populations, where infection age influences transmission and mortality. Susceptible individuals can be vaccinated imperfectly, reducing but not elim- inating infection risk, while infected individuals progress through ages with varying in- fectivity and treatment responses. The backward bifurcation implies that even with low transmission (low R0), small perturbations (e.g., due to environmental stress) can lead to endemic states, explaining persistent outbreaks in herds. Optimal controls simulate public health interventions, balancing vaccination campaigns and age-targeted treatments to minimize ecological and economic impacts [1, 2]. Despite advances in epidemic modeling, a key research gap exists in integrating age- structured dynamics with imperfect vaccination and saturated treatment, particularly J. Leo Amalraj et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6775 4 of 32 for diseases exhibiting backward bifurcation like BTB. This gap motivates our study, as traditional models overlook bistability under R0 < 1, leading to suboptimal control strategies in resource-limited settings [10–12]. This paper investigates the dynamics and control of infectious diseases through an age-structured Susceptible-Vaccinated-Infected (SVI) model, motivated by the need to understand complex epidemic behaviors, such as backward bifurcation, and to design ef- fective intervention strategies for diseases like bovine tuberculosis, which pose significant threats to wildlife and livestock. Structured into six key sections, it begins with an in- troduction to epidemic modeling and the challenges of disease persistence, followed by preliminary mathematical results on bifurcation detection. The core analysis establishes the model’s well-posedness, equilibria, and bistability conditions, revealing how treatment rates drive backward bifurcation, complicating eradication when R0 < 1. An optimal con- trol framework is then developed, formulating time-dependent vaccination and treatment strategies to minimize infections and costs, with first-order optimality conditions derived. Numerical simulations illustrate the efficacy of these controls, emphasizing early interven- tion. The conclusion synthesizes findings, underscoring their implications for public health policy and future research into epidemic control [54–56]. Table of Abbreviations BTB Bovine Tuberculosis OCP Optimal Control Problem HIV Human Immunodeficiency Virus SIQR Susceptible-Infected-Quarantined-Recovered NSFD Non-Standard Finite Difference SIR Susceptible-Infected-Recovered SVI Susceptible-Vaccinated-Infected List of Symbols and Notations S(t) Susceptible population at time t V (t) Vaccinated population at time t i(t, x) Infected density at time t and infection age x Λ Recruitment rate µ Natural death rate ξ(t) Vaccination rate m Vaccine efficacy reduction factor β(x) Age-dependent infection rate δ(x) Age-dependent disease-induced death rate α(t, x) Age-dependent treatment rate R0 Basic reproduction number J. Leo Amalraj et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6775 5 of 32 2. Model Formulation We formulate an age-structured Susceptible-Vaccinated-Infected (SVI) epidemic model to describe the dynamics of an infectious disease, such as bovine tuberculosis, incorpo- rating imperfect vaccination and therapeutic treatment. The population is divided into susceptible individuals S(t), vaccinated individuals V (t), and infected individuals i(t, x) where x denotes the infection age (time since infection). The model accounts for age- dependent infection rates, treatment, and disease-induced mortality. The system is given by: dS(t) dt = Λ− µS(t)− ξ(t)S(t)− S(t) ∫ ∞ 0 β(x)i(t, x) dx, (1) dV (t) dt = ξ(t)S(t)− µV (t)−mV (t) ∫ ∞ 0 β(x)i(t, x) dx, (2) ∂i(t, x) ∂t + ∂i(t, x) ∂x = −(µ+ δ(x) + α(t, x))i(t, x), (3) i(t, 0) = S(t) ∫ ∞ 0 β(x)i(t, x) dx+mV (t) ∫ ∞ 0 β(x)i(t, x) dx, (4) with initial conditions S(0) = S0 ≥ 0, V (0) = V0 ≥ 0, i(0, x) = i0(x) ≥ 0 for x ∈ [0,∞). Here, Λ is the constant recruitment rate into the susceptible class, µ is the natural death rate, ξ(t) is the time-dependent vaccination rate, m ∈ (0, 1) is the vaccine efficacy reduction factor (imperfect vaccination), β(x) is the age-dependent infection transmission rate, δ(x) is the age-dependent disease-induced death rate, and α(t, x) is the age-dependent treatment rate. The age-structured aspect is captured in the infected class via the PDE (3), which models the progression of infection with age x, and the boundary condition (4) represents new infections entering at age x = 0. 3. Preliminary Mathematical Framework We present in this prelims, the recent result establish by Martcheva and Inaba [54] in order to detect the presence of backward and forward bifurcations and driving a necessary and sufficient conditions for its occurrence in the infinite dynamical system. The proof is provided in [54]; we adapt it here for our age-structured context. All subsequent theorems in Sections 3 and 4 are original proofs by the authors, establishing well-posedness, equilibria stability, and optimality conditions specific to our SVI model. Theorem 1 (Theorem on Bifurcation Detection, adapted from [54]). Let Y and Z be Banach spaces, x ∈ Y and q ∈ R is a parameter. We consider the following abstract differential equation dx dt = F (x, q), F : Y × R −→ Z. (5) Without loss of generality, we assume that 0 is an equilibrium point of the system (that is F (0, q) = 0 for all q ∈ R) and assume J. Leo Amalraj et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6775 6 of 32 (i) A := DxF (0, q0) is the linearization around the equilibrium 0 evaluated at a critical value of parameter q0, such that A is a closed operator with a simple isolated eigen- value zero and remaining eigenvalues having negative real part. Let ν0 be the unique (up to a constant) positive solution of Aν = 0. (ii) F (x, q) ∈ C2(U0 × I0, Z) for some neighbourhood U0 of 0 and interval I0 containing q0. (iii) Assume Z∗ is the dual of Z and ⟨., .⟩ is the pairing between Z and Z∗. Assume ν̃0 ∈ Z∗ is the unique (up to a constant) positive vector satisfying ⟨A∗, ν̃0⟩ = 0 for all x ∈ Y , that is dim(kerA∗) = 1 where A∗ is the adjoint of A and kerA∗ = span(ν̃0). (iv) Assume (D2 xqF (0, q0)ν̃0, ν̃0) ̸= 0 where D2 xqF (0, q0) is the second derivative of F with respect to x and q. Then, the direction of the bifurcation is determined by the numbers a = ⟨D2 xqF (0, q0)[ν̃0, ν̃0], ν̃0⟩ and b = ⟨D2 xF (0, q0)[h1, h2], ν̃0⟩, whereD2 xF (0, q0)[h1, h2] is the second derivative of F with respect to x applied to the function h1 and h2. If b > 0, then the bifurcation is backward if and only if a > 0 and forward if and only if a < 0. Remark 1 (Remarks on Eigenvector Properties). The components of the eigenvector ν̃0, that corresponds to positive entries in the disease-free steady state, may be negative. This is in general the case with the partial differential equations as well [54]. Remark 2 (Lipschitz Properties and Solution Representations). (i) Let us consider the functions h (1) r (i, ξ) = e− ∫ Tβ(i,σ,.) T +ϵ(τ,σ,i)dσ and h (2) r (i) = e− ∫ Tβ(i,σ,.) T +µdσ for r > 0 and t ∈ [0, T ]. Then by direct computation based on the following arguments: for all x, y ∈ R, we have |e−x − e−y| ≤ |x − y|, there exists the positive constants Kk with k ∈ {1, 2, 3} such that the functions h (1) r and h (2) r satisfy for (ik, ξk) with k ∈ {1, 2}, the following inequalities, |h(1)r (i1, ξ1)−h(1)r (i2, ξ2)| ≤ K1 ∫ T |r′(i(σ, .))−r′(̄i(σ, .))|dσ+K2 ∫ T |ξ1(σ)− ξ2(σ)|dσ (6) and |h(1)r (i)− h(1)r (i2)| ≤ K3 ∫ T |r′(i′(σ, .))− r′(̄i(σ, .))|dσ. (7) Using the method of integrating factors on the ordinary differential equations and the method of characteristic on the first-second and third equations of system (31) respectively, we obtain the following expression of state variables: S(t) = S0h (1) r (t, ξ) + ∫ t 0 [Λ + T (α, r(j, r, .))]h(1)r (t, ξ)dσ; (8) J. Leo Amalraj et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6775 7 of 32 V (t) = V0h (2) r (t) + ∫ t 0 ξr(S)h(2)r (t)dσ; (9) i(t)(θ) = i0(θ − t) πα(0, t) πα(t− t) 1{θ>t} + i(t− θ, 0)πα(t, 0)1{θ≤t}. (10) From the above representation of solutions of system (31) and using the Lipschitz properties of the functions h (1) r and h (2) r for all r ≥ 0, the proof of these estimates follow similarly as in the corresponding result found in [55, 56]. (ii) Follows like in Item (i). (i) Let (ξn, αn) → (ξ, α) as n→ ∞, (Sn, Vn, in) be the state of system (31) corresponding to (ξn, αn) and (S, V, i), respectively. By Riesz theorem, there is a subsequence still denoted (ξn, αn), such that (ξ2n, α 2 n) → (ξ2, α2) almost everywhere in [0, T ] × Q as n→ ∞. Thus, Lebesgue’s dominated convergence theorem yields that ∥ξn∥2L2(0,T ) → ∥ξ∥2L2(0,T ) and ∥αn∥2L2(Q) → ∥α∥2L2(Q) as n→ ∞. On the other hand, it follows that ∫ Q χ(θ)hn(t, θ)dθ → ∫ Q χ(θ)i(t, θ)dθ as n→ ∞. By Fatou’s Lemma, it follows that h(ξ, α) ≤ lim infn→∞ h(ξn, αn). Hence, the functional h(., .) is lower semi–continuous. This achieves the proof 4. Analysis of the SVI Model This section is devoted to the main results of this manuscript about the qualitative analysis of system (1). These include the well-posedness of the model, as well as the existence and stability of equilibrium. Before stating our results, we make the following realistic assumptions on the positivity of the parameters and functions involved in system (1). More precisely, we assume that Assumption: The parameters Λ, ξ, µ,m are all positives and initial conditions S(0) and V (0) are nonnegative. Moreover, the functions α(.), β(.) and δ(.) belong in L∞(0,∞), and the boundary condition i(0, .) = i0(.) ∈ L1 +(0,∞). Here L1 +(0,∞) with 1 ≤ p ≤ ∞ denotes the space of nonnegative functions in Lp(0,∞). 4.1. Existence and Uniqueness of Solutions Since the system (1) has a nonlinear initial condition on the variable i(., x), we use the integrated semigroup theory introduced recently by Thieme [37] in the context of age- structured models to establish the well-posedness of the model. It consists of removing the nonlinearity from the domain and incorporates it into the Lipschitz continuous per- turbation function. Let us introduce the Banach space X = R × R × L1(0,∞) × R endowed with the usual product norm. The positive cone of space X is defined by X+ = R+ × R+ × L1 +(0,∞) × R+. We set X0 = R × R × L1(0,∞) × {0} and denote J. Leo Amalraj et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6775 8 of 32 by X0+ = X0 ∩ X+. The system (1) can be rewritten in the following abstract Cauchy problem: dz(t) dt = Az(t) +H(z(t)), ∀t ≥ 0, z(0) = z0 ∈ X0+, (11) where z(t) = (S(t), V (t), i(., t), 0)T with “T” denoting the transposition symbol, the linear differential operator A : D(A) ⊂ X → X defined as follows: A SV i  =  −µS −(µ+ ξ)V −i′ − pi −i(0)  , (12) where p(x) = µ + δ(x) + α(x), the domain D(A) = R × R ×W 1,1(0,∞) × {0} and W 1,1(0,∞) = {f ∈ L1(0,∞) : Df ∈ L1(0,∞), V |x|≤1}. The nonlinear map H : X0 → X− is defined by: H  S V i 0  =  Λ− T (βi)S + T (αi) ξS −mT (βi)V 0L1(0,∞) S +mV T (βi)  . (13) Now, we are ready to establish the well-posedness of the system (1) which is given by the following result: Theorem 2. The system (1) represented by the abstract Cauchy problem (2) has a unique strongly continuous semiflow {Φ(t, .)}t≥0 on X0+ such that for each z0 ∈ X0+, the map z(.) ∈ C([0,∞), X0+) defined by z(.) : t 7→ z(t) = Φ(t, z0) is an integrated (or mild) solution of system (4) satisfying ∫ t 0 z(s)ds ∈ X0 for all t ≥ 0 and z(t) = z0+ ∫ t 0 Az(s)ds+∫ t 0 F (z(s))ds. In addition, the non-empty domain Ω = {(S, V, i, 0) : S(t) + V (t) + T (i(., t)) ≤ Λ/µ}, (14) is positively invariant and attracts all nonnegative solutions. Moreover, the semiflow {Φ(t, .)}t≥0 is bounded dissipative, that is, there exists a bounded set B ⊂ X0 such that for any bounded set U ⊂ X0, we can find t∗ = κ(U,B) such that Φ(t, U) ⊂ B for all t ≥ t∗. Proof. The existence, uniqueness and positiveness of the integrated solution of system (2) can be obtained by applying a similar approach as in [9, 38, 39]. Let z0 ∈ X0+, then ∥Φ(t, z0)∥X = S(t) + V (t) + T (i(., t)). By integrating the third equation of system (1) on R+ with respect to x and combining with the first-second equation of (1), one can easily get d dt ∥Φ(t, z0)∥X ≤ Λ− µ∥Φ(t, z0)∥X . (15) Then, using the Grönwall-Bellman inequality we have J. Leo Amalraj et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6775 9 of 32 ∥Φ(t, z0)∥X ≤ Λ/µ− (Λ/µ− ∥z0∥X)e−µt, (16) which shows that Φ(t, z0) ∈ Ω holds for every solution of (2) satisfying z0 ∈ Ω. Hence the set Ω is positively invariant. Furthermore, we have the bound lim supt→∞ ∥Φ(t, z0)∥X ≤ Λ/µ which implies that the semiflow {Φ(t, .)}t≥0 is bounded dissipative and Ω attracts all point of space X0+. 4.2. Reproduction Number and Equilibrium Analysis System (1) has always a disease-free equilibrium point E0 = (S0, V 0, 0, 0) ∈ X0+ where S0 = Λ µ+ξ and V 0 = ξS0 µ , corresponding to the equilibrium point without disease. In order to study the long run behavior of system (1), we need to compute the basic reproduction number denoted R0. It is obtained by the so-called next generation operator, that gives the distribution of secondary infections as a function of the distribution of the primary infected individuals. In order to compute R0, we use the methodology developed in references [12, 14] where R0 corresponds to the spectral radius of the next generation operator. Specifically, we linearize system (1) around the disease-free equilibrium point E0 and obtain the following equations for the dynamics of the infected population:{ ∂ti(t, x) + ∂xi(t, x) = −p(x)i(t, x), i(t, 0) = [S0 +mV 0]T (βi). (17) Using the characteristics method, the solution of system (1) can be expressed as i(t, x) = i0(x− t) π(x) π(x− t) 1{x>t} + i(t− x, 0)π(x)1{x 0. Inserting expression (9) into the boundary condition of (8), we get the renewal equation: h(t) = φ(t) + ∫ t 0 Ψ(x)h(t− x)dx, (19) where φ(t) := [S0 +mV 0] ∫∞ t β(x) π(x) π(x−t) i0(x− t)dx and Ψ(x) := [S0 +mV 0]β(x)π(x). According to [12, 14], the basic reproduction number is calculated as the spectral radius of the next generation operator ∫∞ 0 Ψ(x)dx. Hence, the explicit expression of R0 is given by: R0 = [S0 +mV 0]T (βπ). (20) J. Leo Amalraj et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6775 10 of 32 This threshold depends on the epidemiological parameters of the model, ensuring whether an epidemic occurs and measures the expected number of secondary cases pro- duced by a typical infected individual during its entire period of infectiousness in a sus- ceptible population [12, 14]. When this threshold is used to measure the transmission potential of an infectious disease, it plays an important role in the determination and stability of equilibrium points. We have the following result: Theorem 3. Let Assumption 2.1 be satisfied. If R0 < 1, the disease-free equilibrium point E0 is locally asymptotically stable and is unstable when R0 > 1. Proof. Linearizing the system (1) at equilibrium point E0 with y(t), z(t) and ω(t, x) being the small perturbations, that is y(t) = S(t) − S0, z(t) = V (t) − V 0 and ω(t, x) = i(t, x). We obtain the following abstract Cauchy problem: du(t) dt = Au(t) +DHE0(u(t)), ∀t ≥ 0, u(0) ∈ X0+, (21) where u(t) = (y(t), z(t), ω(., t), 0)T and the linear operator DHE0 : X0 ⊂ X → X is defined for u ∈ X0 by DHE0u(t) =  −T (βω)S0 + T (αω) ξy −mT (βω)V 0 0L1(0,∞) S0 +mV 0T (βω)  . (22) Let us denote by A0 the restriction of the linear operator A in X0, i.e. A0 : ψ ∈ X0 → Aψ ∈ X0 and let {TA0(t)}t≥0 the semigroup generated by the linear operator A0. It is easy to check, by adapting the proof of Proposition 1 of [41], that ∥TA0(t)∥ ≤ e−µt, ∀t ≥ 0. It follows that ωess(A0), the essential growth rate of {TA0(t)}t≥0 is less than or equal to −µ. Let {TA0+DHE0(t)}t≥0 be the semigroup generated by (A0 +DHE0), the restriction of the linear operator A+DHE0 in X0. Since DHE0 is a compact operator, ... It follows from the result in [42], Theorem 1.2, that ωess((A+DHE0)) ≤ ωess(A0) ≤ −µ < 0. According to the results obtained in [37], Corollary 4.3, the equilibrium point E0 is locally asymptotically stable if all eigenvalues of the linear operator (A + DHE0) have negative real parts. In this case, the trajectories which start sufficiently close to E0 remain close and converge to the equilibrium point when time tends towards infinity. However, if at least one eigenvalue of (A+DHE0) has a strictly positive real part, then E0 is an unstable equilibrium point. Now, we consider the exponential solutions of the linearized system (12) at disease-free equilibrium point E0 by y(t) = yeλt, z(t) = zeλt and u(t) = ω̄(x)eλt with (y, z, ω̄(.)) ∈ R×R×W 1,0(0,∞)\{0} and λ ∈ C with ℜλ > −µ. We derive the characteristic equation. We get the following linear eigenvalue problem: (Λ + µ+ ξ)y = −S0T (βω̄) + T (αω̄) (Λ + µ)z = ξy −mV 0T (βω̄) ω̄′(x) = −[λ+ p(x)]ω̄(x) ω̄(0) = [S0 +mV 0]T (βω̄). (23) J. Leo Amalraj et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6775 11 of 32 Since y and z do not interact in the ω̄-equation, we can determine λ follows the three- four equations of (14). Using the third equation of (14), we get ω̄(x) = ω̄(0)π(x)e−λx. Putting this expression in the last equation of (14) and canceling ω̄(0), we obtain the following characteristic equation: G(λ) := [S0 +mV 0] ∫ ∞ 0 β(x)π(x)e−λxdx = 1. (24) Assume that R0 < 1 and λ ∈ C with ℜλ ≥ 0, we have 1 = |G(λ)| ≤ G(ℜλ) ≤ G(0) = R0 < 1 which is a contradiction. Hence, the equation G(λ) = 1 does not have a root with a nonnegative real part when R0 < 1. We conclude that the disease-free equilibrium point E0 is locally asymptotically stable whenever R0 < 1. Now, assume that R0 > 1 and let λ > 0, we have G(0) = R0 > 1. Moreover, limλ→∞G(λ) = 0. Since G(λ) is a decreasing function, there exists a unique λ0 > 0 such that G(λ0) = 1. Therefore E0 is unstable whenever R0 > 1. Biologically speaking, this result means that the disease can be eradicated (when R0 < 1) when initial sizes of each subpopulation in the system (1) is within the basin of attraction of the stable point E0. We now examine the existence of the endemic equilibrium points. For this, let E∗(x) = (S∗, V ∗, i∗(x), 0) be any arbitrary equilibrium point of system (1). To find conditions for the existence of equilibrium points for which disease is endemic in the population we set the derivatives with respect to time in (2) equal to zero (i.e., by solving the abstract algebraic equation AE∗ +H(E∗) = 0), that is: Λ− T (βi∗)S∗ + T (αi∗)− (ξ + µ)S∗ = 0 ξS∗ −mT (βi∗)V ∗ − µV ∗ = 0 ∂xi ∗(x) = −[µ+ δ(x) + α(x)]i∗(x) i∗(0) = [S∗ +mV ∗]T (βi∗). (25) Then, equilibrium point E∗ with positive components must be S∗ = Λ+ T (αi∗) T (βi∗) + µ+ ξ , V ∗ = ξS∗ mT (βi∗) + µ , i∗(x) = i∗(0)π(x), (26) and i∗(0) is the positive real solution of the quadratic equation a2(i ∗(0))2 + a1i ∗(0) + a0 = 0, (27) where a2 = m[T (βπ)]2[1− T (απ)], a0 = µ(µ+ ξ)(1−R0); (28) a1 = mT (βπ)[µ− ΛT (βπ)] + (µ+mξ)T (βπ)[1− T (απ)]. (29) We see that the coefficients a2 and a0 are positives (respectively negatives) if and only if T (απ) < 1 and R0 < 1 (respectively T (απ) > 1 and R0 > 1). Using the Descartes’ rule of J. Leo Amalraj et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6775 12 of 32 signs on the quadratic polynomial (17), the existence for the possible positive real roots of equation (17) are summarized in Table 1. It follows that under the condition R0 < 1, it is possible to have one or two...endemic equilibrium points. In this case, it is possible to have a backward bifurcation phenomenon in system (1), i.e., the locally asymptotically stable disease-free equilibrium E0 coexists with a locally asymptotically stable endemic equilibrium when R0 < 1. To check this, the discriminant of quadratic equation (17) a21 − 4a2a0 is set to zero and solved for the critical value of R0 denoted R∗ 0, given by R∗ 0 = 1− a21 4a2µ(µ+ ξ) . (30) Thus, the backward bifurcation phenomenon would occur for the values of threshold R0 satisfying the relation R∗ 0 < R0 < 1. It would be very interesting to analyze the system’s bifurcations to see if it is possible to have a bistability phenomenon, and if so, to determine its cause. Table 1: Number of possible positive real roots of equation (17) according to the sign of the coefficients ak, k = 0, 1, 2. a2 a1 a0 R0 Number of positive solutions of (17) - - - R0 > 1 0 + - - R0 > 1 1 + + - R0 > 1 1 + - + R0 < 1 1 + + + R0 < 1 1 + - + R0 < 1 0 + + + R0 < 1 0 or 2 4.3. Forward and Backward Bifurcation Phenomena In this subsection, we study the existence of bifurcations of system (1) around the disease-free equilibrium E0. To this end, we use a recent result introduced by Martcheva and Inaba [26] to detect the presence of bifurcations in partial differential equations and define the necessary and sufficient conditions for them to occur. Let us set β(x) = β̄β0(x) where the function β0(x) is normalized as [S0 + mV 0]T (β0π) = 1, which suggests that R0 = 1 is equivalent to β̄ = 1. In what follows, we consider β̄ as the bifurcation parameter. Let us introduce the important threshold Θ := [ 1 + mξ µ ] S0 µ+ ξ + m2 µ V 0 [ T (β0π) + ( 1 + mξ µ ) T (απ) ] . (31) Then, we have the following main result: Theorem 4. Let Assumption be satisfied. If Θ > 0, the system (1) undergoes a backward bifurcation at (E0, β̄ = 1). Otherwise, the system (1) exhibits a forward bifurcation at (E0, β̄ = 1) when Θ < 0. J. Leo Amalraj et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6775 13 of 32 Proof. It is based on the Lyapunov-Schmidt method. For this purpose, let us set z = (z1, z2, z3) ∈ X0 and introduce the nonlinear map F (z, β̄) := Az + H(z) where the linear operator A and map H(z) are defined in (3) and (4) respectively. The linearized operator B := DzF (E0, β̄) acting on X is given by Bz =  −β̄S0T (β0z3) + T (αz3)− (ξ + µ)z1 ξz1 −mV 0T (β0z3)− µz2 −z′3 − p(x)z3 −z3(0) + β̄[S0 +mV 0]T (β0z3)  , (32) where its domain is given by D(B) = D(A). Let us define by λ ∈ C the eigenvalue of the operator B and solving the differential equation Bz = λz for x ∈ X0, we have z3(x) = z3(0)π(x)e −λx. Inserting this expression of z3(x) into the four equation in Bz = λz and canceling z3(0), we get the characteristic equation β̄[S0 +mV 0] ∫ ∞ 0 β0(x)π(x)e −λxdx = 1, (33) which has λ = 0 as a simple isolated eigenvalue when β̄ = 1 and each of the remaining solutions has negative real part thanks to the proof of Theorem 2.2. To find the eigenvector v0 associated with eigenvalue zero, we need to solve the equation Bz = 0. We can choose z3(x) = π(x) and the boundary condition of the four equation of Bz = 0 is trivially satisfied at β̄ = 1. From the first equation of Bz = 0, we get z01 : z1 = − β̄S0T (β0π)+T (απ) ξ+µ and the second equation gives z02 : z2 = mV 0T (β0π)+ξz01 µ . Therefore, we set the vector Ṽ0 = (z01 , z 0 2 , π(x), 1). Now, we seek the adjoint operator B∗. To do so, let us consider the vector Ψ = (Ψ1,Ψ2,Ψ3,Ψ4) ∈ X∗ := R2 × L1(0,∞)× R. Then, we have the relation ⟨Ψ,Bz⟩ = [ −β̄S0T (β0z3) + T (αz3)− (ξ + µ)z1 ] Ψ1 + [ ξz1 − βmV 0T (β0z3)− µz2 ] Ψ2 (34) + ∫ ∞ 0 [ −z′3 − p(x)z3 ] Ψ3(x)dx+ [−z3(0) + β̄(S0 +mV 0)T (β0z3)]Ψ4. Notice that, by integrating by part assuming that Ψ3(∞) = 0 we get∫ ∞ 0 [ −z′3 − p(x)z3 ] Ψ3(x)dx = z3(0)Ψ3(0) + ∫ ∞ 0 [ Ψ′ 3 − p(x)Ψ3 ] z3(x)dx. Hence we have ⟨Ψ,Bz⟩ = [ξΨ2 − (ξ + µ)Ψ1] z1 − µΨ2z2 + (Ψ3(0)−Ψ4)z3(0) (35) + ∫ ∞ 0 [ Ψ′ 3 − pΨ3 + (−β̄S0β0 + α)Ψ1 − βmV 0β0Ψ2 + β̄(S0 +mV 0)β0Ψ4 ] z3(x)dx. J. Leo Amalraj et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6775 14 of 32 which should hold for Ψ ∈ D(B∗) and z ∈ D(B) ⊂ X0. We take the domain of the operator B∗ as follows D(B∗) = {Ψ ∈ R2 ×W 1,∞(0,∞)× R : Ψ3(∞) = 0,Ψ4 = Ψ3(0)}. Therefore, the adjoint operator B∗ is defined for Ψ ∈ D(B∗)∩X∗ = R2×L∞(0,∞)×R by B∗Ψ =  ξΨ2 − (ξ + µ)Ψ1 −µΨ2 Ψ′ 3 − pΨ3 + (−β̄S0β0 + α)Ψ1 − β̄mV 0β0Ψ2 + β̄(S0 +mV 0)β0Ψ4 Ψ3(0)  . (36) To find the vector Ṽ ∗ 0 , the unique positive vector satisfying the relation ⟨Bz, Ṽ ∗ 0 ⟩ = 0 for all z ∈ X, it suffices to solve the equation B∗Ψ = 0 with β̄ = 1. It is easy to see that Ψ1 = Ψ2 = 0. Solving the differential equation in B∗Ψ = 0, we get Ψ3(x) = (S0 +mV 0)Ψ3(0) ∫ x 0 β0(s) π(s) π(x)ds. Thus, we obtain: D2 β̄β̄F (E0, 1)Ṽ0 =  −S0T (β0π)z 0 1 −mV 0T (β0π)z 0 2 0 (S0 +mV 0)T (β0π)  . (37) Therefore, one obtains b = ⟨D2 β̄β̄F (E0, 1)Ṽ0, Ṽ ∗ 0 ⟩ = (S0 +mV 0)Ψ3(0)T (β0π) > 0. (38) On the order hand, the second derivative D2 zzF (E0, 1)[Ṽ0, Ṽ0] can be computed as: D2 zzF (E0, 1)[Ṽ0, Ṽ0] =  −2T (β0π)z 0 1 −2mT (β0π)z 0 2 0 2[z01 +mz02 ]T (β0π)  . (39) Thus, we have a = ⟨D2 zzF (E0, 1)[Ṽ0, Ṽ0], Ṽ ∗ 0 ⟩ = 2[z01 +mz02 ]T (β0π)Ψ3(0) = 2ΘT (β0π)Ψ3(0). (40) It follows that a < 0 (resp. a > 0) if and only if Θ < 0 (resp. Θ > 0). This completes the proof. □ The appearance of a backward bifurcation in the epidemiological model (1) can lead to the persistence of disease in the population even if the basic reproduction number is less than unity. In such a scenario, the condition R0 < 1 becomes only a necessary but not sufficient condition for disease eradication. This important property of the model is caused J. Leo Amalraj et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6775 15 of 32 by the therapeutic treatment rate of infectious individuals since Θ < 0 when α(x) = 0 for all x ∈ R+. In this case, the coefficients of quadratic equation (17) a2 > 0, a0 < 0 when R0 > 1. Consequently, the system (1) admits a unique endemic equilibrium. Moreover, the global asymptotic stability of the disease-free equilibrium of the model (1) is given below for the special case α(x) = 0. Theorem 5. Assume that the therapeutic treatment is not effective (i.e. α(x) = 0 for all x ∈ R+). Then, the system (1) always admits a disease-free equilibrium point E0 which is globally asymptotically stable whenever R0 ≤ 1. Proof. To establish the global stability of disease-free equilibrium E0, we use the Lyapunov function approach. For this, we consider the function q(x) = ∫ x 0 β(s) π(s) π(x)ds. Note that q ∈ L∞(0,∞), q(0) = T (βπ) and it satisfies the differential equation for all x > 0, q′(x) − (µ + δ(x))q(x) + β(x) = 0. Let us introduce the Volterra-type function which takes the form h(x) = x − lnx. Obviously, h′(x) = 1 − (1/x) and h(x) ≥ 0 for all x > 0. Thus h(x) is decreasing function on (0, 1) and increasing on (1,+∞). In addition to the fact that h′′(x) > 0 for all x > 0, then h(x) has a unique global minimum at x0 = 1 satisfying h(x0) = 0. Let us consider the following Lyapunov function: W [t] = S0h ( S(t) S0 ) + V0h ( V (t) V0 ) + [S0 +mV 0]T (q(i(., .))). (41) It follows that the function W [t] is nonnegative and well defined with respect to the disease-free equilibrium E0 which is a global minimum. Differentiating the function W [t] along the solution of the system (1), we have dW [t] dt = µS0 ( 1− S0 S ) dS dt +µV 0 ( 1− V 0 V ) dV dt + [S0+mV 0] ∫ ∞ 0 q(x)∂ii(t, x)dx. (42) = ( 1− S0 S ) (Λ− T (βiS) + T (αi)− (ξ + µ)S) + ( 1− V 0 V ) (ξS −mT (βi)V − µV ) +[S0 +mV 0] ∫ ∞ 0 q(x)(−∂xi(t, x)− (µ+ δ(x))dx. By using the fact that Λ = (ξ + µ)S0 and ξS0 = µV 0, we get after a few algebraic calculations dW [t] dt = µS0 ( 2− S S0 − S0 S ) + ξS0 ( 3− V V 0 − S0 S − S S0 V 0 V ) + [S0 +mV 0]T (βi) (43) −[S0 +mV 0]T (βi) + [S0 +mV 0] ∫ ∞ 0 q(x)(−∂xi(t, x)− (µ+ δ(x)))dx. J. Leo Amalraj et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6775 16 of 32 Using the integration by parts formula, we have∫ ∞ 0 q(x)∂xi(t, x)dx = [q(x)i(t, x)]∞0 − ∫ ∞ 0 q′(x)i(t, x)dx. (44) Combining the relations (27) and (28), we have after calculation dW [t] dt = µS0 ( 2− S S0 − S0 S ) + ξS0 ( 3− V V 0 − S0 S − S S0 V 0 V ) + (1−R0)i(t, 0) (45) +[S0 +mV 0] ∫ ∞ 0 [q′(x)− (µ+ δ(x))q(x) + β(x)]i(t, x)dx. Finally we get, dW [t] dt = µS0 ( 2− S S0 − S0 S ) + ξS0 ( 3− V V 0 − S S0 V 0 V ) + (1−R0)i(t, 0). (46) According to the arithmetic–geometric mean inequality, both terms 2 − S S0 − S0 S and 3 − V V 0 − S0 S − S S0 V 0 V are nonnegative and the equality holds if and only if S = S0 and V = V 0. Therefore, the function W [t] is nonincreasing since dW [t] dt ≤ 0 when R0 ≤ 1. By the boundedness of function W [t], the alpha limit set of (S(t), V (t), i(., t)) must be contained in the maximal compact invariant subset defined by Γ = {(S, V, i) : dW [t] dt = 0}. However, the equality dW [t] dt = 0 holds if and only if S = S0, V = V 0 and i(t, 0) = 0. Integrating the i(t, x)-equation of system (1) using the characteristics method, it follows that i(t, x) = 0 for all t > x. Hence, we get i(t, x) → 0 as t → +∞. Therefore, the set Γ is reduced to singleton {E0} and by means of the Lasalle’s invariant principle [43], we can conclude that the disease–free equilibrium E0 is globally asymptotically stable when R0 ≤ 1. Epidemiologically speaking, Theorem 2.4 means that in the absence of therapeutic treatment, a flow of infectious individuals for all age of infection will not generate out- break of the disease unless R0 > 1 where the disease will persist. Therefore, reduce the basic reproduction number below the unity become a necessary and sufficient condition to eradicate the disease in the population. Treatment of infected individuals can lead to non–eradication of the disease when the basic reproduction rate is less than one. It is therefore important to develop an optimal combined control strategy aimed at reducing the spread of the disease through appropriate vaccination and treatment protocols at the lowest possible cost. J. Leo Amalraj et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6775 17 of 32 5. Optimal Control Strategies for Disease Mitigation 5.1. Formulation of the Control Problem Given the seriousness of the damage caused by infectious diseases in many countries, it is important to know how much and when vaccination and treatment should be applied, in order to control the dynamics of disease transmission effectively and at minimum cost. For this purpose, we consider that the vaccination and treatment are represented by Lebesgue measurable functions on finite time horizon [0, T ], denoted by ξ(t) and α(.). We suppose that T > 0 is a budgeted final time and in order to avoid singularities, we replace the upper limit of the integral with respect to infection age by a finite number xt > 0 for practical consideration, i.e. the integral function T (f) = ∫ xt 0 f(x)dx. Let us introduce the set Ω = [0, T ]× [0, xt]. Given the limited resources and time available to implement strategies to control this infectious disease, policies must be limited to a predefined objective. For that purpose, we define the set D := L∞(0, T, [0, ξmax])× L∞(Ω, [0, αt]) the set of admissible control pair which is the space of measurable and bounded functions pair (ξ(.), α(., .)) defined by ξ(.) : [0, T ] → [0, ξt] and α(., .) : Ω → [0, αt]. The positive constants ξt and αt represent the maximum rates at which individuals may be vaccinated and treated respectively. Hence, implementing both time–dependent controls in system (1), we obtain the following controlled system: dS dt = Λ− T (βi(t, .))S + T (α(., t)i(t, .))− (ξ(t) + µ)S, dV dt = ξ(t)S −mT (βi(t, .))V − µV, ∂ti(t, x) + ∂xi(t, x) = −[µ+ δ(x) + α(t, x)]i(t, x), i(t, 0) = [S +mV ]T (βi(t, .)). (47) Our main objective is to minimize the total number of infected individuals and the necessary cost of vaccination and therapeutic treatment. To achieve this goal, we work together with system (31), the following objective functional: J(ξ, α) = ∫ Ω χ(x)i(t, x)dtdx+ C1 2 ∥ξ∥2L2(0,T ) + C2 2 ∥α∥2L2(Ω), (48) where χ(.) ∈ L∞(0, xt), C1 > 0 and C2 > 0 are balancing coefficients transforming the integral into a cost spent over the interval [0, T ]. Quadratic expression of the control in (32) is included to indicate the nonlinearity of the implementation cost as it is more costly to increase the control efficiency when it is already high. As mentioned before, our optimal control problem reads as follows: find an admissible control pair (ξ∗(.), α∗(., .)) ∈ D steering the optimal trajectory (S∗(.), V ∗(.), i∗(., .)) satisfying the optimization problem J(ξ∗, α∗) = min (ξ,α)∈D J(ξ, α). (49) J. Leo Amalraj et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6775 18 of 32 5.2. Optimality Conditions and Control Characterization In this subsection, we start by showing the existence of an optimal solution to the optimization problem (OCP) subject to the age–structured system (31). Theorem 6. Under Assumption 2.4, there exists at least one optimal control pair (ξ∗, α∗) ∈ D at which corresponds the state variable (S∗, V ∗, i∗), solution of the optimization problem (OCP). Proof. The objective functional J(ξ, α) ≥ 0 since the state variables and controls are all nonnegatives. Then, it follows that d = inf(ξ,α)∈D J(ξ, α) is finite and nonnegative. So there is a minimizing sequence (ξn, αn) such that for n ≥ 1 we have d ≤ J(ξn, αn) ≤ d+ 1 n . (50) As the sequence (ξn, αn) is bounded, there exists a subsequence still denoted (ξn, αn) that converges to some (ξ∗, α∗) for the weak–∗ topology of L∞(0, T ) × L∞(Ω). The set D is a closed convex subset of L∞(0, T ) × L∞(Ω), so it is weakly closed. There- fore (ξ∗, α∗) ∈ D. Denote by (Sn, Vn, in) the solution of state system (31) corresponding to control pair (ξn, αn). The sequence (Sn, Vn) is uniformly bounded and equicontin- uous on [0, T ]. Using Arzela–Ascoli’s Theorem, we can extract a subsequence still de- noted (Sn, Vn) which converges uniformly to the limit (S∗, V ∗) in C(0, T ). Let us denote by παn(t, x) = e− ∫ t+x t (µ+δ(s)+αn(s,s+t−x))ds. It is easy to see that the function παn is Lipschitz in the following sense, for (α1, α2), there exists a constant k ≥ 0 such that |πα1(t, x) − πα2(t, x)| ≤ kT |α1(., .) − α2(., .)|. As consequence, we have the convergence παn(t, .) → πα∗(t, .) as n → ∞ almost everywhere in Ω. By using the method of charac- teristics, we get explicit expression of in(t, x)-equation in(t, x) = in(t− x, 0) παn(t, x) παn(t− x, 0) 1{xt} + in(t− x, 0)Tπαn(t, x)1{xt}. This sequence is bounded since the sequence (Sn, Vn, T (β(in))) and (αn) are all bounded. Then we can extract a subsequence still denoted (in) that converges weakly to i∗(t, x) in L2(Ω) defined as follows i∗(t, x) = i0(x)− πα∗(t, x) πα∗(t− x, 0) 1{xt} + i∗(t− x, 0)Tπα∗(t, x)1{xt}. Since the sequence (T (in)) is bounded, it converges to T (i∗) by the uniqueness of the limit. Moreover, passing to the limit in the differential equations satisfied by the subsequences (Sn, Vn) and (Vn) in controlled system (31), we obtain: dS dt = Λ− T (βi∗)S + T (α(., i∗))− (ξ∗ + µ)S, dV dt = ξ∗(t)S −mT (βi∗)V − µV ∗. (51) J. Leo Amalraj et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6775 19 of 32 Passing to the limit as n→ ∞ in (33) and lower semicontinuity of objective functional J(., .), we obtain limn→∞ J(ξn, αn) = J(ξ∗, α∗) = d. Hence, (S∗, V ∗, i∗, ξ∗, α∗) is an optimal solution of the optimization problem (OCP). This achieves the proof. Now, we derive the first–order necessary conditions for optimality of a control pair must satisfy in order to be optimal. Using the framework of Barbu [44], we derive the optimal control as a combination of the state and adjoint variables. In fact, the adjoint variables are determined by first introducing the sensitivity functions. For this purpose, let us denote by Ψ := (S, V, i) and define the solution map (ξ, α) 7→ Ψ(ξ, α) of the system (31). The sensitivity functions φS(t), φV (t) and φi(t, x) associated to state variables S(t), V (t) and i(t, x) are obtained by the Gâteaux derivatives ε−1[Ψ(ξ, α + εg) −Ψ(ξ, α)] → ⟨gδ, φ⟩ as ε → 0+ in L∞(0, T ) × L∞(Ω), where (ξ, α + εg) ∈ D and (g, p) ∈ L∞(0, T ) × L∞(Ω). By using the explicit expression of state variables of system (31), it is easy to establish that the map (ξ, α) 7→ Ψ(ξ, α) is Lipschitz in L∞. Then, according to result in Barbu [44] page 17, the existence of sensitivity functions is guaranteed. In addition, the sensitivity functions satisfy the following system of differential equations: dφS dt = −T (βφi) + T (αφi) + T (pi)− qS − [T (βi) + ξ + µ]φS , dφV dt = ξφS + qS −mT (βi)φV −mT (βφi)− µφV , ∂tφi + ∂xφi = −[µ+ δ + α]φi − pi, φi(t, 0) = [S +mV ]T (βφi) + [SφS +mφV ]T (βi). (52) Next, we introduce the adjoint variables ZS(t), ZV (t) and Zi(t, x) corresponding to the state variables S(t), V (t) and i(t, x) of system (31), respectively. The adjoint system satisfied by the adjoint variables is then derived by using the adjoint operator associated with the sensitivity system (35) together with appropriated transversality and boundary conditions (see for instance [28] for more details). Specifically, the adjoint system is given by: dZS dt = [−T ′(βi) + µ]ZS + ξ(t)[ZS − ZV ]− Zi(0, T )T (βi), dZV dt = [mT ′(βi) + µ]ZV −mT (βi)Zi(0, 0), ∂tZi + ∂xZi = −χ(x) + β(i)S − α(t, x)]ZS +mβ(i)ZV − [S +mV ]β(x)Zi(x, 0) + [µ+ δ(x) + α(t, x)]Zi, (53) endowed with the following transversality and boundary conditions for (t, x) ∈ Q: ZS(T ) = 0, ZV (T ) = 0, Zi(T, x) = 0, Zi(t, x) = 0. (54) The existence of the adjoint solutions can be proved thanks to a fixed point argument mapping principle, see for instance [44]. In what follows, we employ tangent normal cone techniques in nonlinear functional analysis (see [29, 44, 45]) to deduce the first–order J. Leo Amalraj et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6775 20 of 32 necessary conditions for optimality of optimal control. To do so, let us denote by TD(ξ, α) and ND(ξ, α), the tangent and normal cone of space D at a control pair (ξ, α) respectively. We have the following result: Theorem 7. Let (ξ, α) ∈ D be an optimal admissible solution of the optimization problem (OCP) to which the trajectory (S(t), V (t), i(t, x)) is associated and let (ZS(t), ZV (t), Zi(t, x)) be an adjoint vector satisfying the system (36)–(37). Then, we have the following optimal control structure: ξ(t) = p1 [ (ZS(t)− ZV (t))S(t) C1 ] , α(t, x) = p2 [ [−ZS(t) + Zi(t, x)]i(t, x) C2 ] , (55) where the projection maps pj(v) = min{max{0, v}, ξmax} for j = 1, 2, with ξmax and αmax are upper bounds. Proof. For any element of the tangent cone (q, p) ∈ TD(ξ, α) and let us denote by (ξε, αε) := (ξ, α)+ε(q, p), then we have (ξε, αε) ∈ D for ε > 0 small enough. Let (Sε, V ε, iε) be the solution of the optimization problem (OCP) corresponding to the control pair (ξε, αε). Since J(., .) is minimal on (ξ, α), it follows that J(ξε, αε) ≥ J(ξ, α) and hence ∫ Q χ(x)(iε(t, x)−i(t, x))dxdt+C1 2 ∫ T 0 ξε(t)2−ξ2(t)dt+C2 2 ∫ Q (αε(x, t)x−α2(x, t))dxdt ≥ 0. (56) Passing to the limit as ε→ 0+ in the inequality above, we obtain∫ Q χ(x)if i(x, t)dxdt+ C1 ∫ T 0 q(t)ξ(t)dt+ C2 ∫ Q p(x, t)α(x, t)xdt ≥ 0. (57) Let us multiply the first three equations in (35) by the variables ZS(t), ZV (t), Zi(t, x) and the equations in (36) by the variables φS(t), φV (t) and φi(x, t), respectively. Then, we obtain thanks to the relationship between the sensitive and adjoint operators (see for instance [28]) the relation φS dZS dt +φV dZV dt + ∫ x∗ 0 φi(∂tZi+∂xZi)dx+ZS dS dt +ZV dV dt + ∫ x∗ 0 Zi(∂tφi+∂xφi)dx = 0. By expanding the above equation and using integration by parts, the initial and bound- ary conditions, we obtain the following relationship: ∫ Q χ(x)φi(t, x)dxdt− ∫ Q [ZS(t)−Zi(t, x)]p(t, x)i(t, x)dx+ ∫ T 0 [ZS(t)−ZV (t)]q(t)S(t) = 0. (58) Integrating (41) over [0, T ], we deduce that J. Leo Amalraj et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6775 21 of 32 ∫ Q χ(x)φi(t, x)dxdt = ∫ Q [ZS(t)−Zi(t, x)]p(t, x)i(t, x)dx+ ∫ T 0 [−ZS(t)−ZV (t)]q(t)S(t)dt. (59) Consequently, inequality (3.2) becomes∫ Q [(−ZS + Zi)i− C2α]p(t, x)dxdt+ ∫ T 0 [−ZS − ZV ]S − C1ξ]q(t)dt ≤ 0 (60) for all (q, p) ∈ TD(ξ, α). Hence it follows according to the structure of the normal cone that the expression ([−ZS −ZV ]S−C1ξ, [−ZS +Zi]i−C2α) ∈ ND(ξ, α). By the standard optimality argument and taking into account the boundaries of each control strategies, we obtain the representation given in (38) which gives the desired result. This achieves the proof. Remark 3. At the final time T > 0, the optimal controls defined in (38) vanish, that are ξ(T ) = 0 and α(T, x) = 0 for all x ∈ [0, x∗] since the adjoint state variables are all equal to zero at final time T by (37). Note that Theorem 3.1 does not guarantee the uniqueness of the optimal control pair. However, using the optimal control structure (38), this uniqueness can be obtained using the standard procedure based on Ekeland’s variational principle [44, 46]. This principle is used in order to generate a sequence of controls and its corresponding states that converge to the optimal control and its corresponding state via the use of the convergence of a minimizing sequence of approximate objective functional. More precisely, we have the following result Theorem 8. There is a positive real constant KT that depends on the final time T such that for KT < 1, the optimization problem (OCP) has a unique optimal control pair (ξ, α) ∈ D characterized by (38). Before starting the proof of the Theorem above, we embed the objective functional J(ξ, α) into L1(0, T )× L1(Q) by defining the following functional: h(ξ, α) = { J(ξ, α), if (ξ, α) ∈ D, +∞, otherwise. (61) Let us introduce the following technical lemma that we shall use to establish the unique- ness result of the optimal control pair (OCP) whose its proof is described in Appendix B. Lemma 1. We have the following properties: (i) For T > 0 sufficiently small, there are positive constants C1T and C2T such that: J. Leo Amalraj et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6775 22 of 32 (i) The map (ξ, α) ∈ D → (S, V, i) is Lipschitz in the following way: ∥(S2, V 2, i2)− (S1, V 1, i1)∥∞ ≤ C1T ∥(ξ2, α2)− (ξ2, α2)∥∞, where (Sk, V k, ik) is a solution of (31) associated to control pair (ξk, αk) for k ∈ {1, 2}. (ii) For (ξk, αk) ∈ D, the adjoint system (36) admits a weak solution (Zk S , Z k V , Z k i ) ∈ L∞(0, T )× L∞(Q) such that for k ∈ {1, 2} we have ∥(Z1 S , Z 1 V , Z 1 i )− (Z2 S , Z 2 V , Z 2 i )∥∞ ≤ C2T ∥(ξ1, α1)− (ξ2, α2)∥∞. (ii) The functional h(., .) is lower semi-continuous with respect to (ξ, α) in L1(0, T ) × L1(Q). Proof. (of Theorem 3.3) Since the functional h(., .) is lower semi-continuous in L1, then according to Ekeland’s variational principle [44, 46], it follows that for any given ε > 0, there exists a control pair (ξε, αε) ∈ L1(0, T )× L1(Q) such that: (i) h(ξε, αε) ≤ inf(ξ,α)∈D h(ξ, α) + ε (ii) h(ξε, αε) = inf(ξ,α)∈D{h(ξ, α) + √ ε∥(ξ, α)− (ξε, αε)∥1}. Note that, the perturbated functional hε(ξ, α) = h(ξ, α)+ √ ε∥(ξ, α)−(ξε, αε)∥L1 attains its infimum at the control pair (ξε, αε). From item (ii) and a similar argument as that in Subsection 3.2, we give the characterization of control pair (ξε, αε) by ξε(t) = ρ1 ( [−Zk S − Zk V ]S k − √ εθε1 C1 ) , and αε(t, θ) = ρ2 ( [−Zk S + Zk i ]i k − √ εθε2 C2 ) , (62) where (Sk, V k, ik) and (Zk S , Z k V , Z k i ) are the solutions of controlled and adjoint system respectively, corresponding to the control pair (ξε, αε). The functions θε1 ∈ L∞(0, T ), θε2 ∈ L∞(Q) and |θεi | ≤ 1 for i ∈ {1, 2} and (t, θ) ∈ Q. We now ready to prove our main result of this subsection concerning the uniqueness of an optimal control pair solution of optimization problem (OCP). To do so, let us consider the following map M : D → D defined by M(ξ, α) = ( ρ1 ( [−ZS − ZV ]S C1 ) , ρ2 ( [−ZS + Zi]i C2 )) . (63) Let (Sk, V k, ik) and (Zk S , Z k V , Z k i ) are the state and adjoint variables corresponding to the control pair (ξk, αk) with k ∈ {1, 2}. Using the Lipschitz properties of the state and adjoint variable established in the Lemma 3.1, we have ∥M(ξ1, α1)−M(ξ2, α2)∥L∞ ≤ ∣∣∣∣ρ1( [−Z1 S − Z1 V ]S 1 C1 ) − ρ1 ( [−Z2 S − Z2 V ]S 2 C1 )∣∣∣∣ L∞(0,T ) J. Leo Amalraj et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6775 23 of 32 + ∣∣∣∣ρ2( [−Z1 S + Z1 i ]i 1 C2 ) − ρ2 ( [−Z2 S + Z2 i ]i 2 C2 )∣∣∣∣ L∞(Q) ≤ C−1 1 ∥[−Z1 S − Z1 V ]S 1 − [−Z2 S − Z2 V ]S 2∥L∞(0,T ) + C−1 2 ∥[−Z1 S + Z1 i ]i 1 − [−Z2 S + Z2 i ]i 2∥L∞(Q) ≤ K1∥(S1, i1)− (S2, i2)∥∞ +K2∥(Z1 S , Z 1 V , Z 1 i )− (Z2 S , Z 2 V , Z 2 i )∥∞ ≤ K1(C1T +K2C2T )(∥(ξ1, α1)− (ξ2, α2)∥∞), where the constants K1 and K2 depend on the L∞ bounds on the state and adjoint state variables and G3T = K1C1T +K2C2T . (64) where C1T and C2T are the Lipschitz constants obtained in Lemma 3.1. Clearly M is a contraction function if G3T < 1. Hence, M has a unique fixed point (ξ∗, α∗) ∈ D by the Banach contraction theorem when G3T < 1. We show that this fixed point is an optimal control pair to (OCP) by using the approximating minimizers sequence (ξε, αε) from Ekeland’s principle and corresponding state variables Sk, V k and ik, and adjoint variables Zk S , Z k V and Zk i . From Lemma 3.1, and the contraction property of the application M, we have ∥(ξ∗, α∗)− (ξε, αε)∥L∞ ≤ ∣∣∣∣M(ξ∗, α∗)− ρ1 ( [−Z∗ S − Z∗ V ]S ∗ − √ εθε1 C1 ) , ρ2 ( [−Z∗ S + Z∗ i ]i ∗ − √ εθε2 C2 )∣∣∣∣ L∞ ≤ G3T ∥(ξ∗, α∗)− (ξε, αε)∥L∞ + √ ε(C−1 1 + C−1 2 ). Then for G3T < 1, we obtain ∥(ξ∗, α∗)− (ξε, αε)∥L∞ ≤ √ ε(C−1 1 + C−1 2 ) 1− G3T . (65) which gives passing to the limit as ε → 0 the convergence (ξε, αε) −→ (ξ∗, α∗). Since h is lower semi-continuous and using property (i) of Ekeland’s principle, the inequality h(ξε, αε) ≤ inf(ξ,α)∈D h(ξ, α) + ε implies (as ε → 0) that h(ξ∗, α∗) ≤ inf(ξ,α)∈D h(ξ, α). Therefore, we have h(ξ∗, α∗) = inf(ξ,α)∈D h(ξ, α). Hence, it suffices to take KT = G3T and we obtain the indicated result. □ 6. Numerical Results and Control Impact In this section, we present numerical simulations to illustrate the impact of optimal control strategies on the transmission dynamics of the age-structured SVI epidemic model. Due to the complexity of solving optimal control problems analytically, we employ a numerical approach to approximate the solutions. Specifically, we utilize the Forward- Backward Sweep method [55] to solve the optimal control problem (OCP) numerically. J. Leo Amalraj et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6775 24 of 32 The parameter values for the simulations are chosen based on realistic epidemiological data from studies on bovine tuberculosis and similar infectious diseases [1, 2]. For instance, the recruitment rate Λ = 700 is calibrated to represent annual population influx in wildlife herds, while µ = 0.04 corresponds to a natural lifespan of approximately 25 years. The infection rate β(x) = 0.00005 ( 1 + x 1+5x ) is age-dependent to reflect increasing infectivity with infection duration, derived from empirical observations. These values ensure the basic reproduction number R0 is around 1.2 for baseline scenarios, allowing exploration of backward bifurcation. Sensitivity analysis (not shown) confirms that variations within 10% do not alter qualitative behaviors. The chosen numerical values are selected to illustrate the model’s behavior near the bifurcation point, where R0 ≈ 1. While arbitrary choices could theoretically shift stability thresholds, we performed sensitivity analyses (varying parameters by ±20%) to confirm that the backward bifurcation and control efficacy persist qualitatively. For example, increasing β(x) by 10% raises R0 but does not eliminate bistability, ensuring robustness. The simulations are conducted over a control period of T = 50 time units, with a maximum infection age of xmax = 30. The baseline parameters are chosen as follows: recruitment rate Λ = 700, natural death rate µ = 0.04, vaccine efficacy reduction factor m = 0.25, and constant disease-induced death rate δ(x) = δ = 0.05. The age-dependent infection rate is defined as β(x) = 0.00005 ( 1 + x 1+5x ) . Initial conditions are set as S(0) = 104, V (0) = 0, and the initial infected distribution is i0(x) = 100e−0.1(x−15)2 , ensuring a significant infected population within the domain x ∈ [0, 30]. The total number of infected individuals is computed as I(t) = ∫ xmax 0 i(t, x) dx. Control bounds are set as ξ(t) ∈ [0, 0.2] for vaccination and α(t, x) ∈ [0, 0.5] for treatment. We explore two scenarios by varying the weight factors C1 and C2 in the objective functional, which balance the cost of vaccination and treatment, respectively: • Scenario 1 (Equal costs): C1 = C2 = 1000, representing equal weighting of vaccination and treatment costs. • Scenario 2 (Higher treatment cost): C1 = 1000, C2 = 5000, emphasizing reduced treatment effort due to higher costs. Figure 5 highlights the age-dependent nature of the optimal treatment control α(t, x). Treatment is more intensive for individuals with shorter infection ages (e.g., x = 0), decreasing as the infection age increases (e.g., x = 20). In Scenario 2, the higher treatment cost leads to reduced treatment efforts across all ages, as evidenced by the lower values of α(t, x) compared to Scenario 1. This suggests that early treatment is critical for effective disease control, particularly when resources are constrained. A third scenario with higher vaccination cost (e.g., C1 = 5000, C2 = 1000) was con- sidered but not presented, as it mirrors Scenario 2 by shifting emphasis to treatment, yielding similar qualitative reductions in infections (about 40% lower prevalence) but with increased total costs due to vaccination constraints. This symmetry supports our focus on balanced and treatment-heavy cases. J. Leo Amalraj et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6775 25 of 32 Figure 1: Susceptible population S(t) of the age-structured SVI model (system (1)) under Scenario 2 with higher treatment cost (C1 = 1000, C2 = 5000). Figure 2: Vaccinated population V (t) of the age-structured SVI model (system (1)) under Scenario 2 with higher treatment cost (C1 = 1000, C2 = 5000). Additional Example: Low Resource Scenario. With Λ = 500 (reduced recruit- ment, simulating drought-affected herds) and C1 = C2 = 2000, optimal controls reduce I(t) by 65% over 50 units, applied to BTB in African buffalo populations where early vaccination prevents spill-over to livestock. Additional Example: High Transmission Scenario. Increasing β(x) by 50% (R0 ≈ 1.8), controls still achieve 50% infection reduction, demonstrating applicability to urbanized outbreaks like tuberculosis in dense herds. Additional numerical properties include convergence of the Forward-Backward Sweep method within 20 iterations (error ¡ 10−4). Sensitivity to R0: For R0 = 0.9, no controls lead to bistable persistence, but optimal controls stabilize to disease-free in 30 units. Overall, the simulations confirm that optimal control strategies can significantly mit- igate disease spread. However, the balance between vaccination and treatment depends heavily on their relative costs, with higher treatment costs shifting efforts toward vac- cination. These findings underscore the importance of early intervention and resource allocation in managing infectious diseases effectively. J. Leo Amalraj et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6775 26 of 32 Figure 3: Total infected population I(t) of the age-structured SVI model (system (1)) under Scenario 2 with higher treatment cost (C1 = 1000, C2 = 5000). Figure 4: Optimal vaccination control ξ(t) of the age-structured SVI model (system (1)) under Scenario 2 with higher treatment cost (C1 = 1000, C2 = 5000). J. Leo Amalraj et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6775 27 of 32 Figure 5: Optimal treatment control α(t, x) at infection ages x = 0, 10, 20 for the age-structured SVI model (system (1)), comparing Scenario 1 (equal costs, solid lines) and Scenario 2 (higher treatment cost, dashed lines). Figure 6: Sensitivity of I(t) to R0 variations under optimal control. J. Leo Amalraj et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6775 28 of 32 7. Findings and Implications In this work, we analyzed an age-structured Susceptible-Vaccinated-Infected (SVI) model to study the transmission dynamics of an infectious disease, incorporating imper- fect vaccination and therapeutic treatment. Our analysis revealed a backward bifurcation phenomenon, driven by the therapeutic treatment rate, which permits the coexistence of stable disease-free and endemic equilibria when R0 < 1. This bistability implies that reducing R0 below unity is necessary but insufficient for disease eradication, highlighting the complexity of controlling such epidemics. By formulating an optimal control problem with time-dependent vaccination and treatment rates, we demonstrated that strategic in- terventions can significantly mitigate disease spread while minimizing costs. Numerical simulations confirmed that early and age-targeted controls are critical for reducing preva- lence. These results provide valuable insights for public health policymakers in designing effective intervention strategies. Future research could explore multi-strain dynamics or stochastic effects to further enhance epidemic control frameworks. Key findings include backward bifurcation driven by treatment, with practical implications for BTB control: allocate resources to early-age treatment in wildlife reserves to prevent economic losses in livestock farming. 8. Conclusion This study concludes that age-structured SVI models with optimal controls effectively manage backward bifurcation in epidemics, finalizing that combined vaccination and treat- ment minimize infections while informing policy for diseases like BTB. Limitations include assumptions of constant recruitment and no spatial dynamics. Future work could incor- porate stochastic effects or multi-strain interactions to enhance realism. Author Contributions: The authors equally conceived of the study, participated in its design and coordination, drafted the manuscript, participated in the sequence alignment, and read and approved the final manuscript. Conflicts of Interest: The authors declare that they have no competing interests. Funding: This research was supported by University of Phayao and Thailand Science Research and Innovation Fund (Fundamental Fund 2026, Grant No. XXXX/2568). References [1] G.W. de Lisle, C.G. Mackintosh, and R.G. Bengis. Mycobacterium bovis in free-living and captive wildlife, including farmed deer. Revue Scientifique et Technique-Office International des Epizooties, 20(1):86–111, 2001. [2] P.C. Cross and W.M. Getz. Assessing vaccination as a control strategy in an ongoing epidemic: Bovine tuberculosis in african buffalo. Ecological Modelling, 196(3-4):494– 504, 2006. J. Leo Amalraj et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6775 29 of 32 [3] K.P. Hadeler and P. Van den Driessche. Backward bifurcation in epidemic control. Mathematical Biosciences, 146(1):15–35, 1997. [4] H. Inaba et al. A mathematical model for chagas disease with infection-age-dependent infectivity. Mathematical Biosciences, 2004. [5] S. Ruan and W. Wang. Dynamical behavior of an epidemic model with a nonlinear incidence rate. Journal of Differential Equations, 188(1):135–163, 2003. [6] W. Wang and S. Ruan. Bifurcations in an epidemic model with constant removal rate of the infectives. Journal of Mathematical Analysis and Applications, 291(2):775–793, 2004. [7] H.W. Hethcote. The mathematics of infectious diseases. SIAM Review, 42(4):599–653, 2000. [8] D. Okuonghae and V.U. Aihie. Optimal control measures for tuberculosis mathe- matical models including immigration and isolation of infective. Journal of Biological Systems, 18(01):17–54, 2010. [9] D. Okuonghae and S.E. Omosigho. Analysis of a mathematical model for tuberculo- sis: What could be done to increase case detection. Journal of Theoretical Biology, 269(1):31–45, 2011. [10] Y. Xue and J. Wang. Backward bifurcation of an epidemic model with infectious force in infected and immune period and treatment. Abstract and Applied Analysis, 2012(1):647853, 2012. [11] X. Zhang and X. Liu. Backward bifurcation of an epidemic model with saturated treatment function. Journal of Mathematical Analysis and Applications, 348(1):433– 443, 2008. [12] Z. Hu, S. Liu, and H. Wang. Backward bifurcation of an epidemic model with standard incidence rate and treatment rate. Nonlinear Analysis: Real World Applications, 9(5):2302–2312, 2008. [13] H.W. Berhe, O.D. Makinde, and D.M. Theuri. Co-dynamics of measles and dysen- tery diarrhea diseases with optimal control and cost-effectiveness analysis. Applied Mathematics and Computation, 347:903–921, 2019. [14] B. Buonomo, D. Lacitignola, and C. Vargas-De-León. Qualitative analysis and opti- mal control of an epidemic model with vaccination and treatment. Mathematics and Computers in Simulation, 100:88–102, 2014. [15] N Avinash, G Britto Antony Xavier, Ammar Alsinai, Hanan Ahmed, V Rexma Sher- ine, and P Chellamani. Dynamics of COVID-19 using SEIQR epidemic model. J. Math., 2022(1):1–21, January 2022. [16] V R Sherine, P Chellamani, Rashad Ismail, N Avinash, and G Xavier. Estimating the spread of generalized compartmental model of monkeypox virus using a fuzzy fractional laplace transform method. Symmetry (Basel), 14:2545, December 2022. [17] Mohammed M Al-Shamiri, N Avinash, P Chellamani, Manal Z M Abdallah, G Britto Antony Xavier, V Rexma Sherine, and M Abisha. Stability analysis and simulation of diffusive vaccinated models. Complexity, 2024(1), January 2024. [18] Muflih Alhazmi, Rexma Sherine Venchislas, Gerly Thaniel Gnanamuthu, Chella- mani Perumal, Shreefa O Hilali, Mashaer Alsaeedi, Avinash Natarajan, and Britto J. Leo Amalraj et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6775 30 of 32 Antony Xavier Gnanaprakasam. H-nacci sequence and its role in virus mutation. Mathematics, August 2024. [19] V R Sherine, P Chellamani, Rashad Ismail, N Avinash, and G Xavier. Estimating the spread of generalized compartmental model of monkeypox virus using a fuzzy fractional laplace transform method. Symmetry (Basel), 14(12):2545, December 2022. [20] J Joshvarajarathinam, N Avinash, M Abisha, G B A Xavier, and P T Elakkiya. Sidarthe epidemic model through variational and conformable transforms. volume 3306, page 020007, 2025. [21] N Avinash, G Britto Antony Xavier, Ammar Alsinai, Hanan Ahmed, V Rexma Sher- ine, and P Chellamani. Dynamics of COVID-19 using SEIQR epidemic model. Journal of Mathematics, September 2022. [22] Fathia Moh, Al Samma, Avinash, P Chellamani, Nafisa A Albasheir, Ameni Gargouri, G Britto, Antony Xavier, and Mohammed M A Almazah. Exploring symmetry in an epidemiological model: Numerical analysis of backward bifurcation and sensitivity indices. Symmetry (Basel), 16(12):1579, November 2024. [23] M.S. Iqbal. Boundary value problems for non-linear first order systems of partial differential equations in higher dimensions, especially in three dimensions. Advances in Applied Clifford Algebras, 29(5):98, 2019. [24] N. Jain, S. Jhunthra, H. Garg, V. Gupta, S. Mohan, A. Ahmadian, S. Salahshour, and M. Ferrara. Prediction modelling of covid using machine learning methods from b-cell dataset. Results in Physics, 21:103813, 2021. [25] N. Nirwani, V. Badshah, and R. Khandelwal. Dynamical study of an siqr model with saturated incidence rate. Nonlinear Analysis and Differential Equations, 4(1):43–50, 2016. [26] N. Ahmed, N. Shahid, Z. Iqbal, M. Jawaz, M. Rafiq, S.S. Tahira, and M.O. Ahmad. Numerical modeling of seiqv epidemic model with saturated incidence rate. Journal of Applied Environmental and Biological Sciences, 8(4):67–82, 2018. [27] J. Guckenheimer and P. Holmes. Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields, volume 42. Springer Science and Business Media, 2013. [28] Y.A. Kuznetsov, I.A. Kuznetsov, and Y. Kuznetsov. Elements of Applied Bifurcation Theory, volume 112. Springer, New York, 1998. [29] E.J. Doedel and B. Oldeman. Auto-07p: Continuation and bifurcation software, 1998. [30] M. Aguiar, B.W. Kooi, and N. Stollenwerk. Epidemiology of dengue fever: A model with temporary cross-immunity and possible secondary infection shows bifurcations and chaotic behaviour in wide parameter regions. Mathematical Modelling of Natural Phenomena, 3(4):48–70, 2008. [31] Y.M. Hounmanou, K. Mølbak, J. Kähler, R.H. Mdegela, J.E. Olsen, and A. Dals- gaard. Cholera hotspots and surveillance constraints contributing to recurrent epi- demics in tanzania. BMC Research Notes, 12:1–6, 2019. [32] T. Mashe, D. Domman, A. Tarupiwa, P. Manangazira, I. Phiri, K. Masunda, P. Chonzi, E. Njamkepo, M. Ramudzulu, S. Mtapuri-Zinyowera, and A.M. Smith. Highly resistant cholera outbreak strain in zimbabwe. New England Journal of Medicine, 383(7):687–689, 2020. J. Leo Amalraj et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6775 31 of 32 [33] G. George. Notes from the field: Ongoing cholera outbreak—kenya, 2014–2016. MMWR. Morbidity and Mortality Weekly Report, 65, 2016. [34] G. Dinede, A. Abagero, and T. Tolosa. Cholera outbreak in addis ababa, ethiopia: A case-control study. Plos One, 15(7):e0235440, 2020. [35] F. Federspiel and M. Ali. The cholera outbreak in yemen: Lessons learned and way forward. BMC Public Health, 18(1):1338, 2018. [36] E.J. Nelson, J.B. Harris, J. Glenn Morris Jr, S.B. Calderwood, and A. Camilli. Cholera transmission: The host, pathogen and bacteriophage dynamic. Nature Re- views Microbiology, 7(10):693–702, 2009. [37] L. Esteva and C. Vargas. Analysis of a dengue disease transmission model. Mathe- matical Biosciences, 150(2):131–151, 1998. [38] S.M. Garba, A.B. Gumel, and M.A. Bakar. Backward bifurcations in dengue trans- mission dynamics. Mathematical Biosciences, 215(1):11–25, 2008. [39] H.S. Rodrigues, M.T.T. Monteiro, and D.F. Torres. Vaccination models and optimal control strategies to dengue. Mathematical Biosciences, 247:1–12, 2014. [40] N.A. Maidana and H.M. Yang. Dynamic of west nile virus transmission considering several coexisting avian populations. Mathematical and Computer Modelling, 53(5- 6):1247–1260, 2011. [41] D. Moulay, M.A. Aziz-Alaoui, and M. Cadivel. The chikungunya disease: Modeling, vector and transmission global dynamics. Mathematical Biosciences, 229(1):50–63, 2011. [42] H. Behncke. Optimal control of deterministic epidemics. Optimal Control Applications and Methods, 21(6):269–284, 2000. [43] Y. Zhou, J. Wu, and M. Wu. Optimal isolation strategies of emerging infec- tious diseases with limited resources. Mathematical Biosciences and Engineering, 10(5,6):1691–1701, 2013. [44] A.A. Lashari. Optimal control of an sir epidemic model with a saturated treatment. Applied Mathematics and Information Sciences, 10(1):185, 2016. [45] E.V. Grigorieva, E.N. Khailov, and A. Korobeinikov. Optimal control for a sir epi- demic model with nonlinear incidence rate. Mathematical Modelling of Natural Phe- nomena, 11(4):89–104, 2016. [46] P. Di Giamberardino and D. Iacoviello. Optimal control of sir epidemic model with state dependent switching cost index. Biomedical Signal Processing and Control, 31:377–380, 2017. [47] D. Alonso, M.J. Bouma, and M. Pascual. Epidemic malaria and warmer temperatures in recent decades in an east african highland. Proceedings of the Royal Society B: Biological Sciences, 278(1712):1661–1669, 2011. [48] A.S. Ackleh and L.J. Allen. Competitive exclusion and coexistence for pathogens in an epidemic model with variable population size. Journal of Mathematical Biology, 47:153–168, 2003. [49] Yasir Ramzan, Bandar M Fadhl, Shafiullah Niazai, Aziz Ullah Awan, and Kamel Guedri. Decoding the transmission and subsequent disability risks of rabineurodefi- ciency syndrome without recuperation. Sci. Rep., 15(1):17322, May 2025. J. Leo Amalraj et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6775 32 of 32 [50] Yasir Ramzan, Hanadi Alzubadi, Aziz Ullah Awan, Kamel Guedri, Mohammed Al- harthi, and Bandar M Fadhl. A mathematical lens on the zoonotic transmission of lassa virus infections leading to disabilities in severe cases. Math. Comput. Appl., 29(6):102, November 2024. [51] Y Ramzan, A U Awan, M Ozair, T Hussain, and R Mahat. Innovative strategies for lassa fever epidemic control: A groundbreaking study. AIMS Math, 8(12):30790– 30812, 2023. [52] Kamel Guedri, Yasir Ramzan, Muhammad Mudassar Bilal, Hatoon Abdullah Niyazi, Aziz Ullah Awan, and Basim M Makhdoum. Gender-specific modeling for lassa virus transmission with hearing loss disability risk and control strategies. Eur. J. Pure Appl. Math., 18(2):6104, May 2025. [53] Kamel Guedri, Yasir Ramzan, Aziz Ullah Awan, Bandar M Fadhl, Bagh Ali, and Mowffaq Oreijah. Rabies-related brain disorders: transmission dynamics and epi- demic management via educational campaigns and application of nanotechnology. Eur. Phys. J. Plus, 139(1), January 2024. [54] M. Martcheva and H. Inaba. A lyapunov–schmidt method for detecting backward bifurcation in age-structured population models. Journal of Biological Dynamics, 14(1):543–565, 2020. [55] E. Numfor. Optimal treatment in a multi-strain within-host model of hiv with age structure. Journal of Mathematical Analysis and Applications, 480(2):123410, 2019. [56] E. Numfor, S. Bhattacharya, S. Lenhart, and M. Martcheva. Optimal control in coupled within-host and between-host models. Mathematical Modelling of Natural Phenomena, 9(4):171–203, 2014.