Mathematical modelling of infectious diseases has evolved over more than two centuries, beginning with early formal attempts by Daniel Bernoulli, who analysed smallpox dynamics in 1766 [1, 2], and followed by William Farr’s statistical studies of cholera outbreaks in 1852 [3, 2]. Later, Ronald Ross developed a threshold-based compartmental framework for malaria transmission between 1908 and 1911 [4, 5], providing a conceptual foundation for modern infectious-disease modelling. This foundation was significantly strengthened in 1927 by Kermack and McKendrick, who introduced the classical SIR model and the epidemic threshold theorem [6, 7, 8]. Subsequent extensions, including SEIR-type formulations incorporating latent periods [9, 10], and the formalization of the modern use of the basic reproduction number by Dietz, Anderson, and May [11], shaped contemporary epidemiological theory.

Despite major advances in the 21st century, infectious diseases continue to impose substantial global burdens. The recent COVID-19 pandemic illustrated the dificulty of managing multi-strain diseases driven by emergent variants [12]. Similar mutational mechanisms are observed in diseases such as tuberculosis, dengue, and HIV [13, 14, 15]. Modelling these systems requires accounting for multi-strain structure and incidence mechanisms that incorporate infectiousness, behavioural responses, public awareness, and healthcare availability.

Diferent strains often vary in both their transmissibility and their perceived severity, leading to behavioural heterogeneity among susceptible individuals. As a result, distinct incidence functions may more realistically capture strain-specific transmission. The literature documents several incidence functions, including bilinear [16], saturated [17, 18], Beddington–DeAngelis [19, 20], Crowley–Martin [21, 22], and non-monotone structures [23, 24, 25, 26].

Researchers have therefore developed multi-strain SEIR models. Some studies [27, 28, 29] assume bilinear incidence rate for all strains, while others [30, 31] incorporate non-monotone incidence functions. A recent study [32] considers general incidence forms for all strains and derives necessary conditions for equilibria and global stability under certain structural assumptions. However, the epidemiological interpretation of these assumptions is not clearly articulated, and suficient conditions for the existence of equilibria are not explicitly established. Moreover, all of these studies efectively impose homogeneous behavioural responses across strains, which limits their ability to capture realistic behaviour-driven diferences in transmission. Despite these advances, there remains no general framework that provides behaviourally interpretable, mathematically rigorous conditions for the existence and stability of two-strain endemic equilibria under heterogeneous incidence structures.

Behavioural heterogeneity plays an important role in multi-strain outbreaks. Populations often react more cautiously to newly emerged or mutated strains because of uncertainty and perceived risk, whereas older, better-understood strains tend to elicit weaker precautionary responses. This motivates modelling the original strain using bilinear incidence and the mutated strain using non-monotone incidence. Such a formulation captures realistic, behaviour-driven diferences in contact patterns. A model of this type is analysed in [33]. For such multi-strain diseases, understanding strain coexistence requires analysing the co-dynamics before either strain is eradicated. Thus, the study of the two-strain endemic equilibrium is a foundational step.

The authors in [33] show that a two-strain (co-existence) endemic equilibrium exists and is globally stable if 1 \end{document} ]]>. However, because the strains difer in transmission mechanisms and behavioural responses, and because the condition does not involve any behavioural parameter, this symmetry requirement is rarely satisfied in realistic settings. Motivated by this, we re-examine the co-dynamics of a two-strain SEIR model incorporating behavioural heterogeneity and derive general, epidemiologically interpretable conditions for the existence and stability of a two-strain endemic equilibrium.

In this article, we study a two-strain SEIR model in which one strain follows bilinear incidence and the other a non-monotone incidence rate reflecting behavioural efect. The key contributions of this manuscript are: (i) the development of mathematically rigorous, necessary and suficient conditions for the existence of the two-strain endemic equilibrium; (ii) the formulation of global stability criteria expressed in terms of the basic reproduction numbers and the behavioural parameter; and (iii) the epidemiological interpretation of these analytical conditions. Finally, all analytical results are supported by numerical simulations to verify consistency between theory and model behaviour.

We consider a two-strain SEIR model, as illustrated in Figure 1, with a bilinear incidence rate for strain 1 and a non-monotone incidence rate for strain 2 to reflect behavioural changes among susceptible individuals. The model consists of six compartments: susceptible , exposed to strain-1 and strain-2 ， infectious with strain-1 and strain-2 , and recovered , described in Table 1; following the structure in [33] with a key modification.

To capture behavioural responses to strain–2, we adopt the non-monotone incidence function

with 0 \end{document} ]]>. Since behavioural responses to infection typically intensify in a nonlinear manner as disease signals become more prominent, the quadratic term represents the minimal nonlinear form capable of reflecting this efect. It provides a stronger-than-linear reduction in transmission while avoiding the unnecessary complexity, instability, or steepness introduced by cubic or higher-order exponents. In [33], the parameter is considered in the range [0,1], which accounts for mild to moderate behavioural efects. However, empirical studies and simulations (e.g., [34] and [35]) show that stronger behavioural responses may occur in severe outbreaks. To capture this, we extend the model by allowing 0 \end{document} ]]>, including values greater than 1. We exclude the case , which corresponds to a purely monotonic incidence rate, as our focus is on dynamics impacted by behavioural modifications. Transitions between model compartments are represented by various parameters described in Table 1.

We have chosen this compartmental model over others as it has separate compartments for each strain, which help to understand the dynamics of each strain separately. Also, it is not complicating the dynamics as in models with more compartments. We are considering cross-immunity between strains for simplicity and better understanding of co-existence dynamics [36]. There is no loss of temporary immunity, and so the recovered individuals stay in the recovered compartment .

This section discusses global stability analysis of the two-strain equilibrium of the model using the Lyapunov direct method on . Throughout this section, we take

Proposition 3.1. Consider system (2) and its two-strain endemic equilibrium (5). , then .

Proof. Using (6), we get

As , we get,

Corollary 3.2. Consider system (2) and its two-strain endemic equilibrium (5). If , then . In particular, .

Proof. As , we get

Thus, from Equation (9) of Proposition 3.1, . Again, from Proposition 3.1, .

Theorem 3.3. The two-strain endemic equilibrium of (2) is globally asymptomatically stable in

Proof. Consider Lyapunov function on as follows:

Note that, and 0 \end{document} ]]>, over except . Also, is radically unbounded.

By taking time derivative of we have

Using the endemic equilibrium state

and, by collecting appropriate terms, we get

From (6), we have , then

Using the connection between geometric and arithmetic mean, we have

and by Proposition 3.1, . Thus, . Moreover, . Hence, from Lassale’s principle of invariance [38], is globally asymptotically stable.

Corollary 3.4. The two-strain endemic equilibrium of (2) is globally asymptomatically stable in D if and .

The proof of this corollary follows from Theorem 3.3 and Corollary 3.2.

In this section, we numerically illustrate and validate the analytical results obtained in the previous sections. Simulations are performed using the assumed parameter values listed in Table 2, and the ODE system is solved using a standard fourth-order Runge–Kutta scheme with biologically plausible initial conditions. The objective is to confirm the existence and global stability of the two-strain endemic equilibrium under the conditions derived in Lemma 2.3 and Theorem 3.3 .

Throughout, we focus on verifying the inequalities

which characterize the parameter region where coexistence is possible and globally stable.

This study investigated a two-strain SEIR model incorporating behavioural heterogeneity, motivated by the observation that mutated or newly emergent strains often elicit stronger behavioural responses due to unfamiliarity and perceived risk. To capture this mechanism, strain 1 was modelled with a bilinear incidence function, whereas strain 2 was modelled using a non-monotone incidence function modulated by a behavioural parameter .

We established mathematically rigorous, necessary and suficient conditions for the existence and global stability of a two-strain endemic equilibrium. These conditions explain that, because behavioural avoidance suppresses strain 2, it must possess a suficiently higher reproduction number to overcome this behavioural suppression and maintain coexistence. The stability condition further shows that behavioural avoidance modulated by prevents strain 2 from becoming excessively dominant, thereby supporting long-term coexistence. Moreover, suficiently strong behavioural modification may even eliminate strain 2 altogether.

Numerical simulations corroborated all analytical results, confirming the predicted coexistence range and demonstrating how behavioural responses reshape the long-term competition between strains. Behavioural heterogeneity efectively partitions the susceptible population into groups with diferent transmission pathways, thereby relaxing the classical competitive exclusion principle and allowing coexistence even when one strain is intrinsically more transmissible.

Overall, this work provides mathematically grounded and epidemiologically meaningful insights into how behavioural responses influence multi-strain epidemic dynamics. Future investigations may incorporate adaptive behavioural changes, vaccination, treatment interventions, or stochastic efects to further enhance understanding of behaviour-driven epidemiological processes.