In this paper, we suggest a fractional two-strain SVLIR epidemic model incorporating non-monotonic incidence rates, vaccination and quarantine strategies. The model consists of seven ordinary differential equations that describe the interactions between the susceptible, vaccinated, exposed, infected, and removed individuals. We begin by establishing several theorems related to the existence, positivity, and boundedness of the model’s solutions. The analysis reveals four equilibrium points: a disease-free equilibrium, an endemic equilibrium for strain 1, an endemic equilibrium for strain 2, and an endemic equilibrium with respect to both strains. The global stability of these equilibria has been demonstrated using some suitable Lyapunov functions using Lyapunov’s method and LaSalle’s invariance principle. The basic reproduction number \(R_0\) is calculated using the next-generation method and depends on the basic reproduction numbers of strain 1 ( \(R_0^1\) ) and strain 2 ( \(R_0^2\) ). We have demonstrated that the disease dies out when the value of the basic reproduction number is less than one. Additionally, we have observed that the global stability of the endemic steady states depends on the strain basic reproduction number and on the strain inhibitory effect reproduction number. It was remarked also that the strain with higher value of basic reproduction number dominates the other strain. Numerical simulation is provided to validate our theoretical results about the global stability of equilibria and to illustrate the impact of the fractional derivative order parameter and quarantine parameters on infection eradication. A numerical comparison with COVID-19 vaccination data is also included. Finally, we conducted a sensitivity analysis to identify the model parameters with the greatest potential to influence infection dynamics. We notice that our suggested model has some limitations and does not allow predicting the long-term dynamics in other reproduction numbers scenarios.
