In this work, we study a fractional two-strain epidemic SIIR model with a general incidence rate for both strains, incorporating treatment for each strain. The population is divided into four compartments representing susceptible individuals, individuals infected with the first strain, individuals infected with the second strain, and recovered individuals. The proposed model is formulated using a system of four fractional-order differential equations in the sense of Caputo, which allows us to capture memory effects and hereditary properties inherent in the transmission dynamics of infectious diseases.

The mathematical analysis begins with the proof of the existence and uniqueness of positive solutions, ensuring the well-posedness and biological feasibility of the model. Positivity and boundedness of solutions are also established. Several equilibrium points are derived, including the disease-free equilibrium and endemic equilibria associated with each strain. The global stability of these equilibria is investigated by constructing suitable Lyapunov functions and applying fractional-order stability theory.

Numerical simulations are performed to validate the theoretical findings and to illustrate the influence of the fractional-order parameter on the qualitative behavior of the solutions and the convergence speed toward equilibrium states. In addition, a sensitivity analysis is conducted to examine the role of treatment efficiency in reducing the prevalence of infected individuals for both strains. The results highlight the importance of effective treatment strategies and fractional dynamics in controlling multi-strain epidemic outbreaks.