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 \(\alpha \in (0,1]\) . Positivity and uniform boundedness of all solutions are established, and the unique positively invariant region \(\Omega \) is identified. The basic reproduction number \(\mathcal {R}_0 = \mathcal {R}_d + \mathcal {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 \(\mathcal {R}_0<1\) , with instability when \(\mathcal {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 \(\phi _c\approx 0.0043\,\text {day}^{-1}\) for disease elimination. Three additional contributions are made: (i) a closed-form vaccination reproduction number \(\mathcal {R}_v(\phi )\) with a herd-immunity threshold of \(45.9\,\%\) at baseline \(\mathcal {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\,\%\) at \(\alpha =0.7\) ; and (iii) a monotone sensitivity analysis confirming that \(\phi _c\) strictly increases with \(\kappa \) , formalising treatment saturation as the key driver of bifurcation severity.
