Endemic coexistence and competition of virus variants under partial cross-immunity

1.   Introduction

The ongoing threat of infectious diseases, particularly those capable of producing multiple viral variants, has underscored the need for advanced mathematical models to understand and predict the behavior of such diseases. New virus variants, as observed in the COVID-19 pandemic, are prone to exhibit significant differences in transmissibility, immune evasion capacity, and cross-immunity between variants. These differences are difficult for public health and epidemic modeling, considering that understanding the interaction of multiple virus variants in a population is essential in forecasting the trajectory of an epidemic and guiding intervention measures. Mathematical modeling has also been applied to investigate co-infections and the dynamics of their interactions, e.g., the dynamics of interaction between HTLV-2 and HIV-1 ^[1]. Here, an attempt is made to contribute towards such knowledge by developing and studying an SIR (susceptible-infected-recovered) model framework incorporating two competing virus variants with distinct epidemiological characteristics.

Compartmental models like the SIR model, developed by Kermack and McKendrick ^[2], have been fundamental to epidemiology in dividing populations into susceptible, infected, and recovered compartments for simpler modeling of disease dynamics. Though effective in handling single-strain epidemics, standard SIR models are not flexible enough to capture complexity from interacting variants, hence making extensions for incorporating multi-variant dynamics a source of worry for scientists. The initial development of multi-strain or multi-variant models was necessitated by the reality that most infectious diseases, such as influenza and dengue fever, occur in more than one serotype or strain with varying epidemiological features. Multi-variant models attempt to capture these processes by introducing unique compartments or sub-compartments for each variant, each with varying transmission rates, recovery rates, and potential interactions, such as partial cross-immunity between variants. Initial pioneering work by Andreasen et al. ^[3] on multi-strain influenza models entailed competing strains coexisting in the presence of cross-immunity, and they found that transmissibility and degree of cross-immunity have a significant role in determining if variant strains can coexist or if one strain comes to dominate the other.

Additional studies, such as those by Kamo and Sasaki ^[4], expanded these models by introducing seasonal forcing and analyzing how cross-immunity can lead to synchronized and chaotic oscillations in multi-strain epidemics. Garba and Gumel ^[5] further developed the understanding of influenza transmission dynamics by examining the effects of cross-immunity and backward bifurcation, underscoring the complex nature of multi-strain interactions. Similarly, Sneppen et al. ^[6] proposed a simplified model in which hosts can transmit only the most recent infection, creating a "spreading of immunity" effect that captures the diversity and duration of strain coexistence within a population.

New research has developed these early tenets to include immune escape features and viral mutations, particularly for quickly evolving conditions like COVID-19. Mathematical models incorporating partial cross-immunity and asymmetric temporary durations of immunity have shown that minor differences in such factors can yield significant competitive disparity between strains ^[7,8]. For example, Ogura and Preciado's adaptive susceptible-infected-susceptible model demonstrates how individuals can sever connections to infected nodes to avoid disease spread, offering a new network-based epidemic control view ^[9]. Otunga ^[10] further developed this approach by modeling COVID-19 infections with focus on the Delta and Omicron variants, investigating how vaccination and recovery dynamics influence variant spread and control. Burbano Lombana et al. ^[11] employed a new approach by incorporating human behavior and memory impacts in models of simultaneous epidemic strains, which is important for the impact of interventions on multiple strains during COVID-19 epidemics. Olumoyin and Khaliq ^[12] employed a data-driven deep learning model for COVID-19 variant dynamics, which adopted methods to improve short-term prediction based on time-varying transmission rates. Bessonov et al. ^[13] also contributed to this understanding by considering immune responses, i.e., cytotoxic T lymphocytes and neutralizing antibodies, and how they can lead to the emergence and competition of respiratory virus variants and how immune escape mechanisms are established. The models considered in ^[7] show that early COVID-19 variants were more prone to the dramatic takeover while newer strains can potentially co-exist.

The addition of cross-immunity and reinfection to such models has developed a more nuanced understanding of multi-variant epidemics. Cross-immunity occurs when immunity to a single strain cross-protects against another, a particularly relevant idea for viruses that are extremely mutable and lead to variants that can escape immune responses to previous infections. Reich et al. ^[14] employed a novel framework to examine cross-immunity among dengue serotypes and showed that cross-protection duration and intensity were of key importance to disease incidence and emphasized long-term immunization. Ferguson et al. ^[15] and Andreasen ^[16] studied the impact of moderate levels of cross-immunity on interaction dynamics and concluded that variants can become dominant temporarily before being replaced by new strains with greater immune evasion capacity, a phenomenon seen in influenza where strains become dominant in successive seasons due to loss of immunity.

Further work has extended our understanding on how cross-immunity contributes to multi-variant behavior. Sachak-Patwa et al. ^[17] demonstrated that the integration of cross-immunity into models of influenza significantly enhances epidemic prediction, with the potential for further progress in public health policy. Atienza-Diez and Seoane ^[18] explored different types of cross-immunity, including sterilizing and attenuating, and concluded that higher levels of cross-immunity led to less frequent, reduced outbreaks. This shift in epidemic patterns emphasizes the effect of cross-immunity on the herd immunity levels and danger of future outbreaks. Chung and Lui ^[19] employed cross-immunity in a two-strain influenza model to demonstrate that alterations in cross-immunity levels may have significant effects on system stability and periodicity of outbreaks. For studies of COVID-19 variants, Niu et al. ^[20] developed a compartmental model with cross-immunity and heterogeneity of transmissibility for investigating stability and competition dynamics of variants.

A number of studies provide broader perspectives on pathogen dynamics within host and ecological systems. For instance, Ojosnegros et al. ^[21] investigated the competition-colonization trade-offs among viral strains, finding that shifts in virulence distribution support a stable coexistence of low-virulence strains, highlighting viral evolution's complex interplay. Seabloom et al. ^[22] reviewed ecological interactions in pathogen spread, emphasizing the impact of species composition and community structure on infection dynamics, particularly in co-infection and host-pathogen relationships. Gjini et al. ^[23] focused on pneumococcus serotypes and the role of competition in pathogen dynamics, highlighting how direct competition shapes strain coexistence and challenges for vaccination strategies. Ackleh et al. ^[24] also contributed a reaction-diffusion model addressing competitive exclusion and conditions for pathogen strain coexistence. Amador et al. ^[25] analyzed stochastic interactions between antibiotic-sensitive and antibiotic-resistant bacterial strains, applying an extreme values approach to understand epidemic severity in hospital settings. Finally, Jover et al. ^[26] investigated host-phage dynamics in bacterial communities, finding that trade-offs in infection networks can influence the evolutionary and community structure of microbial ecosystems.

While some researchers have made substantial progress in studying the dynamics of multi-variant systems, key knowledge gaps remain about what ensures the persistence, proliferation, or disappearance of competing molecular variants. Current models often crudely parameterize cross-immunity, treating it as either complete or non-existent, while empirical data suggest that cross-immunity is only partial, and depends on factors such as genetic relatedness between variants and immune responses of the host. Third, several of these assume that transmissibility and recovery rates of competing variants are constant, ignoring evolutionary pressures that will favor the emergence of more transmissible or immune-evasive variants.

We aim to bridge these gaps by introducing an SIR model that incorporates two virus variants with distinct transmission and recovery rates, alongside a flexible framework for variable cross-immunity. Our aim is to explain how these factors interact to determine the competitive dynamics between variants. Using stability analysis and numerical simulations, we determine the equilibria behavior of variants — coexistence, dominance, and eradication — in the population. The biological significance of the model lies in its ability to simulate the dynamics of dual-variant epidemics, providing insights into how factors like transmissibility, immune evasion, and partial immunity influence viral spread and informing public health strategies.

2.   Mathematical model formulation

The proposed model aims to simulate the dynamics of an epidemic caused by two simultaneously circulating variants of a virus, which differ in terms of transmissibility, immune evasion, and recovery rates. The model divides the population into five compartments: Susceptible ($S$), infected by variant 1 ($I_1$), infected by variant 2 ($I_2$), and recovered individuals from variants 1 and 2 ($R_1$ and $R_2$, respectively).

We assume that both strains of the virus spread separately, but their interaction through partial cross-immunity complicates the dynamics of the infection. The recovered population is assumed to have partial immunity to the other strain, i.e., individuals recovered from one strain can still be infected by the other but with reduced susceptibility. This provides a form of immunity that diminishes the risk of severe outcomes of reinfection, but not the risk of further transmission. We also assume that death from disease is negligible for both infectious strains, so it is not modeled. This enables the transmission dynamics, immunity, and recovery dynamics to be highlighted.

The model is constructed under several assumptions to simplify the complexities of real-world epidemic scenarios. First, we assume that the total population remains constant, with births ($\Lambda$) and natural death ($\mu$) rates included in the system. The transmission of each variant depends on the contact rate between susceptible individuals and infected individuals, with different transmission rates for each variant ($\beta_1$ and $\beta_2$). Furthermore, the recovery rates for each variant differ ($\gamma_1$ and $\gamma_2$), as do the rates of immunity loss for those recovered from either variant ($\delta$).

The model, whose flow diagram is given in Figure 1, is represented by a system of ordinary differential equations (ODEs) that describe the time evolution of the compartments. These equations incorporate the interactions between the compartments, including the effects of partial immunity due to prior recovery, as well as the dynamics of infection and recovery for each variant. The differential equations are as follows:

Here, the variables represent the population sizes at time $t$ for each compartment: $S(t)$ for susceptible individuals, $I_1(t)$, and $I_2(t)$ for individuals infected with variants 1 and 2, respectively, and $R_1(t)$ and $R_2(t)$ for those recovered from variants 1 and 2, respectively. The parameters govern the interactions between these compartments, including the transmission rates for each variant ($\beta_1$ and $\beta_2$), the recovery rates ($\gamma_1$ and $\gamma_2$), the degree of cross-immunity between the variants ($\alpha_1$ and $\alpha_2$), and the rate of immunity loss ($\delta$).

This model provides a framework for understanding the dynamics of dual-variant epidemics, incorporating key features such as partial immunity and the potential for reinfection, which are crucial for modeling real-world viral outbreaks where multiple variants of a pathogen are in circulation. Through numerical simulations of the system of ODEs, the model can offer insights into the spread of both variants, the impact of different intervention strategies, and the long-term dynamics of dual-variant epidemics.

3.   Equilibria and stability analysis

The analysis of equilibria and stability forms the basis of all that is known about the dynamics of mathematical models in epidemiology. By determining the steady states and examining their stability, we can comprehend the potential for disease occurrence, eradication, and persistence. In this section, we rigorously explore the equilibria of the model, starting with the disease-free equilibrium, where there is no infection, and progressing towards endemic conditions to encompass the coexistence or dominance of infectious agents. With the help of the next-generation matrix, Jacobian analysis, and bifurcation analysis, we aim to explore thresholds and conditions to switch between these states. To find the steady states, we need to solve the system of equations

3.1.   Disease-free equilibrium and $R_0$

We see that disease-free equilibrium (DFE) is $(S, I_1, I_2, R_1, R_2) = (\Lambda/\mu, 0, 0, 0, 0).$ To compute the basic reproduction number $R_0$, we use the next-generation matrix method. The dynamics of the infectious compartments are given by the equation

where $\mathbf{X} = \begin{bmatrix} I_1 \\ I_2 \end{bmatrix}$, $\mathcal{F}$ represents the new infections, and $\mathcal{V}$ represents the transitions out of the infectious compartments. These are given as

The Jacobian matrices at the disease-free equilibrium (DFE) are

The next-generation matrix is given by $FV^{-1}$, where

The basic reproduction number $R_0$ is the spectral radius of $FV^{-1}$, which in this case is the maximum of the diagonal entries:

Theorem 1. The system has basic reproduction number

and DFE is locally stable if and only if $R_0 < 1.$

As $R_0$ crosses the critical value of 1, the system undergoes a bifurcation, transitioning from a stable disease-free state to a scenario where one or both viruses can persist in the population. This bifurcation marks a change in the qualitative dynamics of the system, highlighting the pivotal role of $R_0$ in determining the system's behavior. For convenience, we fix the following notation

so that $R_0 = \max(R_0^{(1)}, R_0^{(2)}).$

3.2.   Endemic equilibrium

An endemic equilibrium occurs when at least one of the infectious compartments is non-zero.

3.2.1.   Single virus dominance equilibria

Due to the symmetry of the model with respect to the two virus components, we focus on one of the two cases for a single-virus equilibrium in the population. Specifically, we consider the endemic equilibrium where $I_1 = 0$ and $I_2 \neq 0$. The situation where $I_2 = 0$ and $I_1 \neq 0$ is analogous and can be treated in a similar manner.

Lemma 1. If $R_0^{(2)} > 1$, the system admits a single virus endemic equilibrium (SVEE) given by

Proof. When $I_1 = 0$ and $I_2 \ne 0$ we have $R_1 = 0, S = \frac{\gamma_2+\mu}{\beta_2}$ and the system of Eq (3.1) simplifies to

Solving the second for $R_2$ we get $R_2 = \frac{\gamma_2}{\mu+\delta_2}I_2$ and substituting this expression for $R_2$ into the first equation

Solving for $I_2$ we get the desired result.

Lemma 2. The endemic equilibrium from the above lemma is locally asymptotically stable if and only if

In particular, if $R_0^{(1)} \ge R_0^{(2)}$ the endemic equilibrium is unstable.

Proof. Linearizing the system (2.1) around the single virus endemic equilibrium (3.2) gives

for

where

Let us consider the linearization around the endemic equilibrium as above. In particular, $I_1^* = R_1^* = 0$ and $A_{33} = 0$ and $S^*, I_2^*, R_2^*$ are as in Lemma 1, giving

We proceed by calculating the eigenvalues.

We see that

is an eigenvalue and the remaining eigenvalues are determined from

Another eigenvalue is $A_{44} = -\delta-\mu - \alpha_2 \beta_2 I_2^* < 0$, leaving

To compute the eigenvalues of the matrix we get

Noting $A_{55} = -\delta_2-\mu$ we arrive at

It follows that the last three eigenvalues are negative. Except for $A_{22} = \beta_1 S^* + \alpha_1 \beta_1 R_2^* -\gamma_1 - \mu$ all other eigenvalues have negative real part. Hence, the stability depends on $A_{22}$, which proves the first assertion of the lemma. As for the second assertion, we note that when $\beta_1 S^* + \alpha_1 \beta_1 R_2^* -\gamma_1 - \mu < 0,$ it follows that

We conclude that the endemic equilibrium becomes unstable if $R_0^{(1)} \ge R_0^{(2)}.$

3.2.2.   Coexistence equilibrium

In this section, we study the endemic equilibrium when both viruses coexist, i.e., when $I_1 \neq 0$ and $I_2 \neq 0$. Under these conditions, the system of equations simplifies, enabling us to express key variables such as the susceptible population S in terms of a cubic equation. To simplify the formulas, we fix the following notations:

Proposition 1. If system (2.1) admits a coexistence equilibrium ($I_1, I_2 \ne 0$) then the number of susceptible $S$ at the equilibrium satisfies the cubic polynomial

with coefficients

Moreover, $S, I_1, I_2$ at the equilibrium must satisfy

It remains unclear whether the necessary conditions stated in the proposition are also sufficient. While it is straightforward to verify that the coefficients $c_0$ and $c_3$ are both positive, ensuring that at least one root of the cubic equation is negative, the possibility of the existence of two distinct coexistence equilibria is not established. This ambiguity warrants further investigation. Our numerical analysis in § 5 demonstrates the existence of a coexistence equilibrium under certain conditions. In addition to establishing existence, conducting a stability analysis is equally crucial for a comprehensive understanding. Again, our numerical analysis in § 5 shows that the existence of coexistence equilibrium does not necessarily imply asymptotic stability.

Proof. When $I_1 \ne 0$ and $I_2 \ne 0$, the system of equations simplifies to

Solving for $R_1$ and $R_2$ we obtain

Since $R_1, R_2 \ge 0$ we obtain the last assertion in the proposition. Letting $x = \beta_1 S I_1, y = \beta_2 S I_2$ and substituting for $\alpha_1 R_2$ and $\alpha_2 R_1$, we get

This gives

Solving for $x$ and $y$:

Substituting this into

and solving for $S$ yields cubic equation:

Simplifying this yields the desired cubic polynomial.

4.   Global stability of the disease-free equilibrium

In Theorem 1, we establish the local stability of the disease-free equilibrium (DFE). In this section, we establish the global stability of DFE. To achieve this, we construct a suitable Lyapunov function and analyze the system dynamics within the positively invariant set $\Omega$. The set $\Omega$ ensures that all solutions remain biologically feasible, and using standard techniques, we demonstrate that when $R_0 < 1$, the DFE is asymptotically stable within $\Omega$. This result highlights the conditions under which the disease can be eradicated from the population, building on approaches seen in similar analyses ^[27].

We define

Thus, we have the following.

Lemma 3. The set $\Omega$ given in (4.1) is positively invariant, that is any solution to (2.1) with initial values in $\Omega$ stays in $\Omega$ for $t \ge 0.$

This lemma follows from standard arguments, as those in ^[28].

Theorem 2. When $R_0 < 1$, the DFE $(S, I_1, I_2, R_1, R_2) = (\Lambda/\mu, 0, 0, 0, 0)$ is asymptotically stable on $\Omega$.

Proof. Let us consider the Lyapunov function given by

where the coefficients $a, b > 0$ are sufficiently small to be determined later. We now show that $V$ is positive definitive and $\frac{dV}{dt}$ is negative definite on $\Omega.$

We note that $V > 0$ on $\Omega$ except at DFE. Hence, it is positive definite. To show negative definiteness, let us write $V = G+F$ where

We have

Since $\alpha_1, \alpha_2 \le 1$ we see that $S+\alpha_1 R_2 \le S+ R_2 \le \Lambda/\mu$ on $\Omega$. Similarly, $S+\alpha_2 R_1 \le \Lambda/\mu$. Thus,

Next, we consider the time derivative of $G:$

Using $\Lambda - \mu S \ge 0$ on $\Omega$ and denoting by $H$ the negative terms

we get

We are given that $R_0 < 1$. We may pick $a > 0$ sufficiently small so that

With this choice of $a$ we have

Combining (4.2) and (4.4) we arrive at

Since $R_0 < 1$, we see from (4.3) that $dV/dt < 0$ on $\Omega$ except at DFE. Hence, it follows from Lyapunov stability theorems, see e.g., [29, Theorem 5.2] and [30, Theorem 4.1], that DFE is asymptotically stable on $\Omega$.

5.   Model calibration and numerical analysis

For calibration, we use dynamics of Delta and Omicron variants of COVID-19. The competitive dynamics of virus variants, such as Delta and Omicron, highlight the pivotal role of partial cross-immunity in shaping epidemiological outcomes ^[31]. The Delta variant, first identified in India in late 2020, became the dominant strain worldwide by mid-2021 due to its high transmissibility and virulence. Omicron, detected in South Africa and Botswana in November 2021, rapidly outpaced the Delta strain to become the dominant variant by December 2021 in many regions, including the United States. This rapid displacement was attributed to the reproduction number of Omicron, which was approximately 3.5 times higher than that of Delta, coupled with enhanced immune evasion capabilities.

Table 1 presents the parameter values used for model calibration. Unless otherwise specified, we use these values in our numerical analysis. Using these values, the basic reproduction numbers for the Delta and Omicron variants are calculated as $R_0^{(1)} = \frac{\beta_1 \Lambda}{\mu(\gamma_1+\mu)} = 3.80$ and $R_0^{(2)} = \frac{\beta_2 \Lambda}{\mu(\gamma_2+\mu)} = 8.29$, respectively. These values fall within the ranges previously reported in the literature (see, e.g., ^[35,36]). We note that letting $\mu = \Lambda$ we assume the constant population model with a total population density equal to 1.

Table 2 outlines the initial population values used in the numerical simulations. These baseline conditions provide a consistent framework for analyzing the infection dynamics and steady-state behavior across the scenarios presented in this study.

Figure 2 illustrates the dynamics of the model through 100 simulations, where the parameters are randomly selected within a $\pm20\%$ range as specified in Table 1.

In numerical analysis, we are interested in studying the endemic prevalence $(I_1 + I_2)(t)$ and the infection ratio $(I_1/(I_1+I_2))(t)$, particularly their final values.

5.1.   Prevalence and infection ratio

While Lemma 2 provides criteria for local stability, it does not offer insights into global dynamics. Therefore, we conduct numerical simulations to better illustrate the global dynamics, the prevalence, and their relationship to the stability of the single-virus dominance equilibrium.

In Figure 3, we present three representative examples of the infection dynamics over 100,000 days for different transmission rates ($\beta_1$ and $\beta_2$). These figures illustrate the behavior of the two-virus system, highlighting how variations in transmission parameters influence the temporal evolution of the infected populations ($I_1(t)$ and $I_2(t)$). The results show distinct patterns of dominance, coexistence, or decline of the viruses, providing a visual foundation for the numerical analysis summarized in the following table. Looking at the last 1000 days of the simulation, we observe periodic behavior, indicating that the infection may not necessarily converge to any of the equilibria.

Table 3 presents numerical results for varying values of $\beta_1$ and $\beta_2$. The table includes the basic reproduction numbers $R_0$, $R_0^{(1)}$, and $R_0^{(2)}$, as well as the corresponding endemic prevalences and infection ratios after $100k$ days. The values of $I_1^*$ and $I_2^*$ represent the endemic equilibria for the single-virus cases where $I_1 = 0$ and $I_2 = 0$, respectively. In addition, the coexistence equilibrium $(I_{1c}^*, I_{2c}^*)$ is shown, reflecting the steady-state values of both infected populations when both viruses coexist. The notation, (u) and (s) indicate whether the equilibrium is locally asymptotically unstable or stable, respectively. We note that when the infection ratio (IR) lies strictly between 0 and 1, both virus variants coexist in the population. From the table, it is evident that as $R_0$ decreases, the infection ratio (IR) approaches 1.0, indicating a higher prevalence of one virus over the other or a dominance of the coexistence equilibrium. The bifurcation of the system occurs as $R_0$ passes critical thresholds, transitioning between stable and unstable equilibria, with the corresponding changes in the infection dynamics. The entries marked as "DNE" (Does Not Exist) suggest that certain conditions lead to a scenario where no valid equilibrium can be reached. These cases typically occur when the basic reproduction number $R_0$ is too low to sustain either virus at endemic levels. Endemic Prevalence (EP) represents the steady-state infection levels of each variant in the population. It reflects the long-term dynamics of the virus after the system has reached equilibrium, providing insights into how the infection stabilizes over time under different conditions.

For instance, when $R_0 = 9.65$, $I_1^*$ does not exist (denoted by DNE), while $I_2^*$ has a small value, and the coexistence equilibrium is nonzero. This suggests that the system tends to focus on the second virus as the more dominant strain when the first virus does not establish a steady state. As $R_0$ decreases to values near 1, the table indicates that the system shifts toward a state where both viruses coexist with low prevalences, resulting in a low infection ratio (IR), and eventually both viruses fail to persist at all when $R_0$ falls below 1.

In the next, we consider prevalence and infection ratio for varying cross immunities $\alpha_1$ and $\alpha_2$ (Figure 4). This is done by considering the 10,000 day iteration of the system and the variablees are calculated by taking the average of the last 100 days from simulations. The left panel shows that higher cross-immunity tends to reduce total prevalence, while the right panel indicates that infection ratio dynamics are more sensitive to variations in $\alpha_1$, suggesting potential asymmetric effects of cross-immunity on variant competition.

To enrich the practical applicability of our results, we note the implications for public health strategies of the multi-variant nature of the epidemics. Indeed, as our analysis points out, multiple endemic equilibria are possible depending on certain factors like cross-immunities and transmission rates. If reaching the DFE is not immediate, strategies might be developed that would drive the dynamics toward the equilibrium with the lowest magnitude of prevalence. This could be achieved by interventions targeting a selective reduction in the dominance of variants of higher transmissibility or immune evasion that, in the end, reduce overall disease burden.

6.   Sensitivity analysis

We analyze the sensitivity of two key outcomes in a multi-strain infectious disease model—total prevalence and infection ratio—to various epidemiological parameters. Sensitivity analysis is a crucial tool in understanding how the variability in model outputs can be attributed to different input variables ^[37]. Using Latin Hypercube Sampling (LHS), 2000 parameter sets are generated for 8 parameters $(\beta_1, \beta_2, \gamma_1, \gamma_2, \alpha_1, \alpha_2, \delta_1, \delta_2)$ across ranges spanning from 0.5 to 1.5 times the base values, as provided in Table 1. The model, solved numerically using ODEs, simulates the dynamics of susceptible, infected, and recovered populations, calculating average steady-state values of the two dependent variables.

Sensitivity analysis is conducted by calculating the Partial Rank Correlation Coefficient (PRCC) using Spearman's rank correlation method, and the results are visualized in Figure 5. The bar charts show the strength and direction of correlations for each parameter, with statistically significant results $(p < 0.05)$ marked by red asterisks. Parameters such as transmission rates $(\beta_1, \beta_2)$ strongly influence Total Infected, while interaction terms $(\alpha_1, \alpha_2)$ significantly impact infection ratio. These findings highlight the most influential parameters, providing insights into disease dynamics and intervention strategies.

Additionally, convergence of the PRCC results was assessed to ensure that the sensitivity analysis stabilized with increasing sample size. Figure 6 displays the convergence plots for total infected (a) and infection ratio (b), demonstrating the stabilization of partial rank correlation coefficients as the number of simulations increases. These plots confirm the reliability of the sensitivity analysis, with results becoming more consistent as the simulation count approaches 1000.

For completeness, scatter plots illustrating the relationships between parameters and dependent variables are provided in Figures 7 and 8. These plots offer a detailed visual representation of how each epidemiological parameter influences the model outcomes—endemic prevalence and infection ratio—across simulation scenarios. By plotting the parameter values against the corresponding model outputs, these scatter plots enable a clearer understanding of the magnitude and direction of the parameter effects, complementing the sensitivity analysis presented earlier.

Figure 5(a) demonstrates the sensitivity of total prevalence to various epidemiological parameters using Partial Rank Correlation Coefficients (PRCC). The transmission rates ($\beta_1$ and $\beta_2$) exhibit the strongest positive correlations, indicating that higher transmission rates significantly increase the total prevalence of infection. Recovery rates ($\gamma_1$ and $\gamma_2$) are negatively correlated with prevalence, underscoring that faster recovery leads to reduced infection levels. Waning immunity rates ($\delta_1$ and $\delta_2$) show weaker positive correlations, suggesting that quicker loss of immunity slightly amplifies the spread by increasing susceptibility to reinfection. These findings highlight the pivotal role of transmission and recovery dynamics in shaping the overall disease burden.

Figure 5(b) reveals the sensitivity of the infection ratio to key model parameters, focusing on the distribution of infections between the two virus variants. Cross-immunity coefficients ($\alpha_1$ and $\alpha_2$) emerge as the most influential factors, with significant correlations that reflect their role in modulating variant competition and dominance. The transmission rates ($\beta_1$ and $\beta_2$) also influence the infection ratio, demonstrating how disparities in transmissibility between variants can shift the balance of prevalence.

These results highlight the importance of sensitivity analysis in designing public health interventions that control multi-strain infectious diseases. Key parameters such as transmission rates $\beta_1, \beta_2$, cross-immunity coefficients $\alpha_1, \alpha_2$, and immunity waning rates $\delta_1, \delta_2$ strongly influence the robustness of model outcomes. It is seen from the above discussions that parameters, such as the transmission rates $\beta_1$, $\beta_2$ and the cross-immunity coefficients $\alpha_1$, $\alpha_2$, have a profound effect on both the total prevalence and infection ratios to necessitate site-specific interventions adaptive to the change in epidemic conditions. Besides, the minor but significant estimate of immunity waning rates $\delta_1, \delta_2$ points out the need for real-time data application to revise strategies for long-term effective management of the epidemic. It is hoped that researchers may seek the development of adaptive measures that could better achieve an optimized response with resources based on continued shifts in disease dynamics.

The sensitivity analysis highlights that the dynamics of multivariant epidemics are deeply influenced by changes in key parameters, such as transmission rates and cross-immunity. These findings have contributed to informing public health strategies since they may be used to devise interventions for particular epidemic scenarios. For example, targeting the most influential parameters-such as enhancing cross-immunity or adjusting vaccination strategies-may help mitigate the spread of more transmissible or immune-evasive variants. The discussion on adaptive intervention strategies, active mechanisms of feedback, and control-oriented implications might provide greater depth. Integration of adaptive measures from real-time epidemic data is key to optimally responding to fluctuating transmission rates and evolving variants. This enables dynamic adjustments that could work towards effective control of the disease burden. Understanding these factors-precisely how those factors influence the long-term course of the epidemics-will help stakeholders build up adaptive responses that are more effective by optimizing resource allocation and timely intervention.

7.   Discussion and conclusions

In this paper, we provide a comprehensive analysis of the dynamics of two competing virus variants within an SIR framework, emphasizing the role of partial cross-immunity and other epidemiological factors. By exploring the sensitivity of key outcomes, such as total prevalence and infection ratio, to model parameters, we indetify critical drivers of disease dynamics. Transmission rates emerge as the most significant contributors to total prevalence, while cross-immunity levels play a pivotal role in shaping the infection ratio and the competitive interplay between variants.

Our findings underscore the complexity of multi-variant epidemics, highlighting the intricate balance between factors like recovery rates, waning immunity, and cross-immunity, which are invaluable for public health planning. These insights hold great practical implications for designing interventions to effectively reduce disease burden and manage variant competition. For example, actions to enhance cross-immunity or act on enhanced transmissibility of new variants can stem the tide of more virulent or immune-escape strains. Enhanced cross-immunity through vaccination, enhancing monitoring and surveillance infrastructure for tracking variant spread, and creating effective communications framing vaccination and cross-immunity are essential steps to managing variant threats. Furthermore, adaptive response strategies tailored to various epidemiological contexts will enable more effective interventions against newly emerging variants. Together, these findings provide a more nuanced understanding of viral competition dynamics to help stakeholders better prepare for and respond to multi-variant infectious disease challenges.

While the mathematical model yields robust theoretical predictions, their implications are also pertinent for rapidly evolving diseases in real-world settings. Furthermore, researchers may examine how various population structures, behavioral responses, and intervention strategies influence these dynamics. Overall, this study highlights the necessity of a refined understanding of variant interactions for informing adaptive and evidence-based public health interventions.

Author contributions

Shirali Kadyrov: Conceptualization, Formal analysis, Investigation, Methodology, Writing – original draft, Writing – review & editing. Farkhod Khaydarov: Formal analysis, Investigation, Writing – review & editing Khudoyor Mamayusupov: Formal analysis, Investigation, Writing – review & editing Komil Mustayev: Formal analysis, Investigation, Software, Visualization, Writing – original draft, Writing – review & editing.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The first author acknowledges the support of a grant from the Ministry of Science and Higher Education of the Republic of Kazakhstan within the framework of the project AP19676669.

Conflict of interest

Authors declare no conflict of interest.