Stability and bifurcation dynamics of a fractional-order two-strain COVID-19 model

We develop and analyse a Caputo fractional-order SEIVR epidemic model for COVID-19 incorporating detected and undetected infectious compartments, saturated treatment reflecting healthcare-capacity constraints, waning vaccine immunity, and memory effects governed by fractional order α ∈ (0,1] α ∈ ( 0 , 1 ] . Positivity and uniform boundedness of all solutions are established, and the unique positively invariant region Ω Ω is identified. The basic reproduction number R₀ = Rd + Rᵤ R 0 = R d + R u is derived via the next-generation matrix method. Local asymptotic stability of the disease-free equilibrium (DFE) is proved using the Matignon criterion when R₀<1 R 0 < 1 , with instability when R₀>1 R 0 > 1 . Global asymptotic stability of both the DFE and the unique endemic equilibrium is established via Volterra-type Lyapunov functions adapted to Caputo systems. A centre-manifold backward bifurcation analysis reveals that a stable endemic equilibrium can coexist with the stable DFE, requiring the vaccination rate to exceed a critical threshold φ c≈ 0.0043\,day⁻¹ ϕ c ≈ 0.0043 day - 1 for disease elimination. Three additional contributions are made: (i) a closed-form vaccination reproduction number Rᵥ(φ ) R v ( ϕ ) with a herd-immunity threshold of $$45.9\,%$$ 45.9 % at baseline R₀=1.85 R 0 = 1.85 ; (ii) a quantitative comparison of fractional versus integer-order stability regions, showing the Matignon angle threshold expands the stable manifold by up to $$30\,%$$ 30 % at α =0.7 α = 0.7 ; and (iii) a monotone sensitivity analysis confirming that φ c ϕ c strictly increases with κ κ , formalising treatment saturation as the key driver of bifurcation severity.