Dynamic analysis and optimal control of co-infection system under different outbreak times of mutant strains

1.   Introduction

Infectious diseases are widespread globally, posing significant challenges to human health and socioeconomic development [1]. In recent years, severe outbreaks of certain infectious diseases have garnered global attention. For instance, the severe acute coronavirus disease 2019 (COVID-19) [2] has transmitted worldwide since its initial detection in December 2019, resulting in a large number of infections and deaths. As of March 2024, there have been over 774 million people infected with COVID-19, with above 7 million deaths reported by the World Health Organization (WHO). Apart from COVID-19, other infectious diseases including malaria, tuberculosis, and HIV/AIDS also pose significant threats to human health, leading to millions of deaths annually [3]. Additionally, outbreaks of vaccine-preventable diseases like measles, pertussis, and influenza can result in large-scale epidemics, casualties, and economic losses [4].

In recent years, infectious disease models have flourished [5], [6], and the emergence of various infectious disease vaccines has alleviated the burden of epidemic control. The history of mathematical epidemiology models can be traced back to Bernoulli [7]. Sean T. McQuade et al. considered social distancing, testing, contact tracing, and isolation to establish an SEIR model [8], forecasting and controlling the epidemic in New Jersey and offering constructive suggestions for epidemic prevention. Zhen Wang and colleagues revealed the extensive impact of SIR and related models on vaccination and epidemic control [9], Jiaxing Chen et al. proposed a susceptible-infectious-quarantined-recovered-susceptible (SIQRS) model with quarantine and studied its evolution on simple complexes [10], Raul Nistal et al. introduced the concept of custom next-generation matrices for an innovative extended SEIR model (SI(n)R model) and explored the model's stability [11], Qun Liu and his team established a stochastic multi-group SIQR epidemic model with standard incidence, providing sufficient conditions for the disease's eradication [12]. Besides, there are an abundance of mathematical models for dual-strain diseases in the literature [13], [14]. Most research on the dynamics of two-virus infections considers cross-immunity and coinfection, while others describe the roles of competition between two competing strains with the characteristics such as cross-immunity and incorporate control parameters.

Research on complex networks has significantly propelled the advancement of epidemiology, with the effectiveness of complex networks in controlling the spread of infectious diseases lying in their provision of a more detailed and comprehensive perspective for understanding and addressing public health challenges. In recent years, the study of complex networks in epidemiology has yielded a wealth of research findings. J. Chen et al. proposed a composite effective degree Markov chain approach (CEDMA) to describe discrete-time epidemic dynamics on higher-order networks [15]. This study effectively predicts the critical points of continuous phase transitions and confirms the existence of discontinuous phase transitions in the susceptible-infected-susceptible (SIS) process. Y. Wu et al. proposed a complex network two-strain epidemic model with imperfect vaccination and immune decline, verified the global stability of each strain dominance equilibrium point, and further derived the critical value[16]. Their study is helpful to understand the dynamic behavior of multi-strain epidemics from the perspective of complex networks. Contrary to most papers, the evolution of a highly pathogenic and genetically unstable virus through mutations is considered, representing a less explored but highly relevant direction. For instance, since the outbreak of COVID-19, various mutant strains including Alpha, Beta, Gamma, and Omicron [17] have emerged, each with different epidemiological characteristics. How will these mutant strains compete and evolve compared to the original strain, affecting the trajectory of the epidemic? What impact does the uncertainty in the timing of mutant strain emergence have on epidemic control? Under the background of the COVID-19 pandemic, studying the mutation models of the virus is of great practical significance for epidemic control and prediction.

Current research typically assumes that the variant strain or competing strain and the original strain emerge simultaneously, and recommendations for epidemic prevention and control are based on this assumption. This is clearly an idealized scenario, with limited applicability in practice. To make the research more relevant to real-world situations and address this gap, we considered a more general scenario in which the variant strain and the original strain do not emerge at the same time. Through theoretical derivation and numerical simulation, we found that the timing of the emergence of the two strains does not affect the final outcome of the epidemic. After discussing the epidemic trends under different emergence times of the variant strain, we further explored the control measures for the epidemic based on different emergence times of the variant strain. By establishing an optimal control system and conducting numerical simulations, we demonstrated the nonlinear relationship between the emergence time of the variant strain and the cost of control.

The main contributions of this work can be summarized as follows:

The structure of the paper is shown below: Section 1 offers an introduction, highlighting the hazards of infectious diseases and the current research status of infectious disease dynamics. Section 2 establishes a mathematical model for unpredictable virus mutations, describing the dynamic process of disease transmission. Section 3 proposes optimal control strategies for COVID-19 and its mutant strains, constructing a control system. Section 4 validates and evaluates the control strategies under different outbreak times for mutant strains using numerical simulation methods, considering control costs. Finally, the last section summarizes the research findings, offers suggestions for infectious disease prevention, and outlines future research directions.

2.   The inter-strain model and basic reproduction number $\mathcal R_{0}$

2.1.   The co-infection model of inter-strain

The COVID-19 virus undoubtedly represents the most prominent class of infectious diseases currently. The ongoing COVID-19 pandemic is still a main public health critical situation significantly disrupting healthcare systems and society as a whole. Since its initial detection at the end of 2019, the incidence of this lethal disease has continuously increased, with a persistently high global mortality rate [18]. Although the dynamics and administration of this disease are researched extensively, healthcare systems still meet multiple barriers to decreasing its threat to human health.

The variation of various viruses has increased the difficulty of controlling disease cases, which is one of the major challenges. Recently discovered variants serve as typical examples, as they exhibit high transmissibility in both vaccinated and unvaccinated populations, despite presenting milder symptoms [19]. Furthermore, other infectious diseases with analogous symptoms pose greater challenges to the diagnosis of infected individuals, further complicating the task of decreasing the risk of COVID-19 for healthcare mechanisms. Influenza, malaria, and tuberculosis are a few instances of such infectious diseases [20].

It is noteworthy that COVID-19 is also impacting other diseases that threaten human life despite the current attention to the ongoing pandemic. Based on a study in [21], the prevalence of COVID-19 has increased the risk of developing drug-resistant strains and has also led to delays in diagnosing and treating patients. Additionally, it has increased the risk of intra-household transmission during lockdown periods, thereby impeding the control of infectious diseases. Furthermore, some studies have indicated that patients concurrently infected with COVID-19 and other infectious diseases have a higher mortality rate [20].

To investigate the internal competition and evolutionary relationship between two strains of the same virus, this paper establishes a model for dual-strain cross-infection. The population in the environment is divided into seven compartments based on epidemiological characteristics: susceptible individuals (S), individuals infected with Strain 1 (I₁), individuals infected with Strain 2 (I₂), individuals co-infected with Strain 1 and Strain 2 (I₁₂), individuals who have recovered from Strain 1 (R₁), individuals who have recovered from Strain 2 (R₂), and individuals who have recovered from cross-infection (R₁₂), with temporary immunity for the recovered individuals. The total population at time t can be represented as:

where N(t) represents the total population at time t, and S(t), I₁(t), I₂(t), I₁₂(t), R₁(t), R₂(t), and R₁₂(t) represent the population in every compartment at time t.

To make the study more realistic, the following assumptions are made:

Let β₁ and β₂ represent the infection rates of the two strains toward susceptible individuals per unit time. The probability of a patient infected with Strain 1 (or Strain 2) being reinfected with Strain 2 (or Strain 1) is denoted as β₂ (or β₁). The probability of susceptible individuals coming into contact with patients infected with Strain 1 and Strain 2 is given by I₁/N and I₂/N, respectively. Therefore, the incremental number of new infections for Strain 1 and Strain 2 at time t, denoted as ΔI₁ and ΔI₂, respectively, can be calculated as β₁I₁S/N and β₂I₂S/N. According to the previous assumptions, susceptible individuals do not directly become cross-infected, but instead become infected with one strain before becoming infected with the other strain implying the derivation of the cross-infected population from the population that was initially infected with only one strain. The effective contact rate between individuals infected with Strain 1 and Strain 2 is I₂/N. Similarly, the effective contact rate between individuals infected with Strain 2 and individuals infected with Strain 1 is I₁/N. Consequently, the incremental number of cross-infected individuals at time t, denoted as ΔI₁₂, is provided by (β₁ + β₂)I₁I₂/N. Parameter definitions for System (2) are presented in Table 1, and the flowchart is explained in Figure 1.

Then the dynamic model for the dual-strain system can be established as follows:

The initial conditions satisfy the following inequalities:

For the System (2) with non-negative initial values, its solutions are non-negative and bounded. The proof is straightforward and can be found in [24]. Positivity is great for the existence of biologically feasible solutions, while boundedness indicates finite solutions. Next, it is proven that the solutions of System (2) are contained within a bounded area Ω.

Lemma 1. The set $\Omega=\left\{\left(S, \,I_{1},\, I_{2}, \,I_{12},\, R_{1},\,R_{2},\,R_{12}\right) \in \mathbb{R}_{+}^{7}: N \leq {\Lambda}/{d}\right\}$ is a positively invariant and attractive set for System (2).

Proof. By summing all the equations in System (2), we acquire

Applying the standard comparison principle [25], we have

If N(0)≤Λ/d, then N(t)≤Λ/d. Therefore, Ω means a positively invariant set for System (2). However, when N(0)≥Λ/d, the solutions of the model will enter Ω in a finite time, or as t→∞, we have N(t)→Λ/d.

Hence, this region is attractive for all solutions to the system. Therefore, the set Ω means a positively invariant and attractive set for the system.

The set Ω means a positively invariant and attractive set for System (2). The model is well-posed in both epidemiology and mathematics.

2.2.   The basic reproduction number $\mathcal {R}_0$

The disease-free and endemic equilibrium constitute two critical special points within the field of epidemiology. The analysis presented in this study commences from an exploration of these two pivotal equilibria. Obviously, the disease-free equilibrium of System (2) is M0=(S0,0,0,0,0,0,0), where S0=Λ/d. The next-generation matrix method [26] is a commonly used approach to calculate the basic reproduction number, $\mathcal {R}_0$. Next, this method will be applied in conjunction with the disease-free equilibrium M₀ to calculate the basic reproduction number of System (2).

Theorem 1. The basic reproduction number is $\mathcal{R}_{0} =\max{\left \{\mathcal{R}_{1},\mathcal{R}_{2} \right \} }$, where $\mathcal R_{1} = {\beta_{1}}/({\gamma_{1}+d+d_{1}})$, $\mathcal R_{2} = {\beta_{2}}/({\gamma_{2}+d+d_{2}})$.

Proof. Utilizing the next generation matrix method [26], the related next generation matrix is obtained as

The individual-to-next-generation transition matrix is given by

By taking partial derivatives of the row vectors of matrices $\mathcal{F}$ and $\mathcal{V}$ in terms of I₁ and I₂, and replacing the values at the disease-free equilibrium, the matrices F and V are obtained as

Therefore, we have

The spectral radius or maximum eigenvalue of FV⁻¹ provides the basic reproduction number as

Upon substituting the values at the disease-free equilibrium, i.e., S₀ = N₀, the basic reproduction number is acquired as

Let $\mathcal R_{1} = {\beta_{1}}/({\gamma_{1}+d+d_{1}})$ and $\mathcal R_{2} = {\beta_{2}}/({\gamma_{2}+d+d_{2}})$. Hence, the basic reproduction number is given by $R_{0}=\max{\left \{ \mathcal{R}_{1},\mathcal{R}_{2} \right \}}$.

The basic reproduction number has served as a threshold for determining the existence of disease-free equilibrium or endemic equilibrium in a system. Then the association between the basic reproduction number and the steadiness of the disease-free equilibrium will be explored.

Theorem 2. If $\mathcal{R}_0 > 1$, the disease-free equilibrium M₀ is unstable. If $\mathcal{R}_0 < 1$, the disease-free equilibrium M₀ is locally asymptotically steady.

Proof. The Jacobian matrix J(M₀) of the system at equilibrium point M₀ is given by

The eigenvalues of J(M₀) are $\lambda_{1}=-d,\lambda_{2}=(\gamma_{1}+d+d_{1})(\mathcal{R}_{1}-1),\lambda_{3}=(\gamma_{2}+d+d_{2})(\mathcal{R}_{2}-1),\lambda_{4}=-\gamma_{12}-d-d_{12},\lambda_{5}=-m_{1}-d,\lambda_{6}=-m_{2}-d,$ and λ12=−m12−d. If $\mathcal{R}_{0}<1$, that is, $\mathcal{R}_{1},\mathcal{R}_{2}<1$, then λ2,λ3<0. According to Lyapunov's first method [27], it is concluded that the disease-free equilibrium M₀ is locally asymptotically stable. If $\mathcal{R}_{0}>1$, the disease-free equilibrium M₀ is unstable.

Theorem 3. If $\mathcal{R}_{0}<1$, the disease-free equilibrium point M₀ is a globally asymptotically stable equilibrium point.

Proof. We construct the Lyapunov function as

Since I1,I2,I12>0, it follows that H(S,I1,I2,I12)>0, and H(S,I1,I2)=0 holds if and only if I1=I2=0.

Next, we prove that dH/dt≤0:

where H ˙=0 holds if and only if I1=I2=0. Therefore, if $\mathcal{R}_{0}<1$, by LaSalle's Invariant Principle [28], it follows that the equilibrium point M₀ is globally asymptotically stable in Ω.

Theorem 4. The following conclusions hold at the local equilibrium points of System (2):

(1) The equilibrium point M1=(S*,I1*,0,0,R1*,0,0) exists if and only if $\mathcal R_{1}>1>\mathcal R_{2}$.

(2) The equilibrium point M2=(S*,0,I2*,0,0,R2*,0) exists if and only if $\mathcal R_{2}>1>\mathcal R_{1}$.

Proof. (1) Set α1=β1I1/N and α2=β2I2/N. Then, System (2) can be rewritten as

Substitute the equilibrium point M₁ into the above equation to get

Then we obtain

where q0=(γ1+d+d1)(m1+d)/γ1, q1=α1*+d.

At this point, the total number of individuals in the population is

Expressions for S*,I1*,R1*, and N* in terms of α1* have now been obtained.

Since α1*=β1I1*/N*, substituting the values of I1* and N* yields

Simplifying System (19) yields the following quadratic equation

where c₂, c₁, and c₀ are the coefficients of the quadratic, linear, and constant terms of the variable α1*, respectively.

It can be seen that c₂ > 0, and when $\mathcal R_{1}>1$, c₀ < 0, i.e., c₀c₂ < 0. The discriminant of System (21) is thus Δ=c12−4c0c2>0, suggesting that the equation has two solutions. Moreover, c₀/c₂ < 0, indicating that the two solutions of System (20) have opposite signs. Thus, when $\mathcal R_{1}>1$, the system has a distinct local endemic equilibrium M1=(S*,I1*,0,0,R1*,0,0).

(2) The proof for $\mathcal R_{2}>1$ is similar and equivalent to M2=(S*,0,I2*,0,0,R2*,0).

When the basic reproduction number $\mathcal R_{0}$ is below 1, the system will eventually reach a disease-free equilibrium where no infected individuals exist in the environment. However, if a pathogen persists in the environment, its basic reproduction number must be above 1.

2.3.   Bifurcation analysis

Next, we study the bifurcation phenomena of System (2) using the centre manifold theory [29]. For the sake of clarity, we define (S(t),I1(t),I2(t),I12(t),R1(t),R2(t),R12(t))=(x1,x2,x3,x4,x5,x6,x7)=x, and then System (2) can be expressed as

where ζ1=x1x2/∑i=17xi, ζ2=x1x3/∑i=17xi, and ζ3=x2x3/∑i=17xi.

The Jacobian matrix of System (2) at the disease-free equilibrium point M₀ can be obtained as System (12). Considering when $\mathcal{R}_0=\mathcal{R}_1={\beta_1}/(\gamma_1+d+d_1)$, choose β1=β1* as a bifurcation coefficient. Through solving $\mathcal{R}_1=1$ for β1*, β1*=β1/(γ1+d+d1) is obtained. Thus, the Jacobian matrix can be expressed as

It is then obtained that the left eigenvector w and the right eigenvector v of the Jacobian matrix Jβ1* are equal, with w=v=[0,β2−(γ2+d+d2),(γ12+d+d12),(m1+d),(m2+d),(m12+d),−d]. To determine the direction of the bifurcation, the signs of the two bifurcation coefficients $\mathcal{A}$ and $\mathcal{B}$ will next be calculated.

By substituting the specific values of the parameters from Table 2, it is not difficult to obtain that $\mathcal{A} < 0$ and $\mathcal{B} > 0$. According to Theorem 4.1 in [30], this excludes the probability of a backward bifurcation at $\mathcal{R}_0 = 1$ in System (2), meaning that a stable disease-free equilibrium and a stable endemic equilibrium cannot coexist simultaneously. Therefore, at the threshold coefficient $\mathcal{R}_1 = 1$, a forward or transcritical bifurcation exists as shown in Figure 2. Consequently, the endemic equilibrium point is distinct, and an appropriate Lyapunov function can be constructed to certify its global asymptotic stability [31].

3.   Control system

Infectious diseases have had a severe impact on the economies, societies, and healthcare systems of countries worldwide, posing significant challenges to human society and public health. To discuss the feasibility of achieving complete control of the COVID-19 virus, optimal control methods [24], [27], [32] are employed, and the dynamics are simulated using System (2). Based on the discussion of the basic reproduction numbers $\mathcal R_{0}, \mathcal R_{1}, \mathcal R_{2}$ and equilibrium points, the transmission rate (β) and the treatment recovery rate (γ) are chosen as control parameters. However, due to the considerable difficulties associated with implementing control measures for disease fatality rates (d₁, d₂, d₁₂) and immune degradation rates (m₁, m₂, m₁₂) in the field of public health, incorporation of control over these two parameters is not considered.

Specific optimal control strategies are presented next, and the cost-effectiveness of prevention and control meathods is analyzed on the basis of the timing of new strain outbreaks.

The specific control flow chart is given by Figure 3. Therefore, considering the aforementioned two public health methods, the original System (2) is modified into the control mechanism below

Up to this point, some measures have been introduced to intervene in the epidemic trend in System (2), namely government measures for epidemic control and prevention of infectious diseases, an increase in healthcare workers, and the improvement of medical facilities, resulting in the optimal control System (25). The control system builds upon the original framework by incorporating the impact of epidemic prevention measures on the disease's trajectory. It integrates control elements as coefficients within the transmission and recovery rates to mirror the real-time variations in the virus's spread and healing rates, with the degree of change directly tied to the strength of the preventive measures. It becomes clear that achieving a transmission rate of 0 and a recovery rate of 1 signifies the epitome of epidemic prevention, yet this is often challenging to realize. In practice, the cost of such measures must be considered. We are faced with the task of balancing the costs against epidemic management, striving to reduce the cumulative number of infections to the greatest extent possible while incurring the minimal cost. Next, the control objective function for System (25) will be constructed.

Let the control time be T₀, and the control variables u₁ and u₂ be bounded and Lebesgue measurable on the interval [0, 1]. Define the set U={(u1,u2)|0≤ui≤1,t∈[0,T0],i=1,2}, then the objective function is represented as

where A_i (i = 1, 2, 3) stands for the cost incurred by treating infected people by Strain 1, Strain 2, and cross-infections. B_i (i = 1, 2) represents the cost incurred by population control measures and medical inputs per unit time. Optimal control aims at minimizing the number of infected people for Strain 1, Strain 2, and cross-infections while minimizing the control cost. The optimal control variables u1* and u2* are sought that satisfy

subject to the controlled System (25). It is easy to verify that the system is compact in the control and state space and that the objective function is convex. Based on Theorem 4.1 in Chapter 3 of Flemming and Rishel [40], we can ensure the existence of an optimal control solution for this control system.

To solve the optimal control Eqs (25)–(27), the Hamiltonian operator $\mathcal{H}$ is defined as

where λ_i (i = 1, 2, ..., 7) are undetermined coefficients. According to Pontryagin's maximum principle [41], [42], the optimal solutions u1* and u2* of System (27) must satisfy ${\partial \mathcal{H}(t, x, u, \lambda)}/{\partial u} = 0$. Setting ${\partial \mathcal{H}}/{\partial u_{i}} = 0$, (i = 1, 2), we have

Solving these equations yields

Then u1* and u2* are obtained as

According to Pontryagin's maximum principle [41], [42], the control functions u₁(t), u₂(t), and the state variable x(t) of System (27) must satisfy

In System (2), the state variable x(t) represents (S(t),I1(t),I2(t),R1(t),R2(t),R12(t)). Substituting the Hamiltonian operator Syetem (28) into Syetem (32), the coefficient λ is obtained as

In the control System (25), since the state variables do not have terminal values, the boundary condition at time T₀ is

Combining Syetems (2), (31), and (33), the optimal control set U that satisfies Syetem (27) is determined.

In the established optimal control system, two key public health measures have been introduced to intervene in the development trend of the COVID-19 epidemic. First, u₁(t) represents government measures for epidemic control and prevention of infectious diseases. This includes limiting population mobility and canceling large-scale gatherings, effectively reducing the effective contact rate between susceptible and infected individuals. The influence of u₁(t) on the transmission rate β is significant as it directly changes the disease's capacity of spreading within the population. Second, u₂(t) encompasses measures to increase healthcare workers and improve medical facilities. This reflects an enhancement in the healthcare system's capability to manage the pandemic, directly affecting the recovery rate γ. By growing the value of u₂(t), the rate that infected individuals recover from the disease is elevated, thereby decreasing the number of active infections in the population. Building upon System (2), System (25) aims to minimize the total number of infections in the population at the lowest possible social control cost. The objective function J(u₁, u₂) quantifies the trade-off between reducing infections and the costs associated with implementing control measures. Here, A_i represents the cost of treating infections, and B_i represents the social cost per unit time of population control measures and medical investments. In the process of minimizing J(u₁, u₂), the optimal control variables u1* and u2* are sought to achieve an ideal balance between controlling the epidemic and maintaining manageable costs. Next, the optimal control solution for Syetem (25) will be found through numerical simulation, and the different evolutionary trends of the epidemic before and after optimal control will be analyzed.

4.   Numerical simulations

The evolution of the virus varies depending on the appearance time of the new strain. In this section, the evolution of two strains under different outbreak times of the new strain is discussed, along with the evolution after the implementation of control measures. Let Strain 1 represent the novel coronavirus, and Strain 2 represent the mutated strain of the novel coronavirus. The initial values of Syetems (25) and Syetem (2) are assumed to be set as S(0) = 10⁶, I₁(0) = 10, I₂(0) = 10, I₁₂(0) = 0, R₁(0) = 0, R₂(0) = 0, and R₁₂(0) = 0. That is, at the beginning of the emergence of the two strains, there were only ten infected people, and there was no evolution of cross-infection at this time. Considering the controlled System (25), the adjoint Eq (33), the optimal control set Syetem (31), and the initial conditions identical to the System (2), the terminal condition is obtained from Syetem (34). The weight values are set as follows: A₁ = A₂ = A₃ = 1, B₁ = B₂ = 0.5 × 10⁷. For the balance between the control variables u₁ and u₂, which are much smaller than the number of infected individuals, larger balance coefficients (cost coefficients) B_i are needed, where i = 1, 2 are chosen. The remaining cost coefficients A_i, where i = 1, 2, 3, are set to 1. The control model is simulated using the backward difference method, and Table 2 shows all the coefficients taking part in the system. Consequently, the basic reproduction numbers for the novel coronavirus and the mutated strain can be obtained as $\mathcal R_{1}=2.1667$ and $\mathcal R_{2}=3.5435$, respectively.

4.1.   Impact of different outbreak times of the new strain on the bistrain model

In this section, the influence of different outbreak times of the new strain on the evolution of two strains is investigated using the bistrain cross-infection Model (2). According to epidemiological theory, in the case of a basic reproduction number below 1, the disease is expected to be controlled in the population, while a value greater than 1 suggests the potential persistence of the disease. In an environment where two different strains coexist, with one strain having a basic reproduction number greater than 1 and the other strain having a number less than 1, the strain with a higher reproduction number is more likely to persist in the population. However, in certain cases with potential superinfection dynamics, the replacement of the disease may depend on the outcome of competition among hosts. Therefore, the dynamics of the model under various outbreak times of the new strain will be studied.

According to Theorem 3, from an epidemiological perspective, when the basic reproduction number of a disease is below 1, it will be eliminated from a population. Therefore, when two diseases compete in a population, regardless of the time interval between the outbreaks of the two strains, the strain with a threshold quantity less than 1 will be eradicated, while if the threshold quantity is greater than 1, it allows the invasion of the competing disease. Then the evolution between the two strains under local equilibrium states ($\mathcal R_{1}>1, \mathcal R_{2}>1$) of Strain 2 with different outbreak times is investigated.

As shown in Figure 4, the competition between Strain 2 and Strain 1 was simulated within 100 days after the outbreak of Strain 1 at different time points. When both strains have a basic reproduction number, $\mathcal R_{1}$ and $\mathcal R_{2}$, greater than 1, regardless of the time point at which Strain 2 outbreaks, Strain 1 will eventually disappear from the environment while Strain 2 becomes the dominant strain and persists in the environment in the long term, according to Figure 4(b),(c). With the disappearance of Strain 1, cross-infection cases also gradually diminish from the environment, as illustrated in Figure 4(d). From Figure 4(a), it is observed that regardless of the outbreak time of Strain 2, the susceptible population in the environment stabilizes at the same level after the disappearance of Strain 1.

In summary, in the bistrain cross-infection model (2), the emergence time of Strain 2 does not affect the evolution outcome of the two strains. Strain 1 will eventually disappear while Strain 2 persists in the long term. However, the outbreak time of Strain 2 does affect the timing of large-scale epidemics in the environment. The longer the time interval between the outbreak of Strain 2 and Strain 1, the more frequent the fluctuations in the susceptible population; the greater the cumulative number of infections for Strain 1, the higher the peak of its infection count, the fewer fluctuations for Strain 2, and the smaller the number of cross-infections. This is because the evolution process of the two strains in the environment is essentially a competition for “resources” (susceptible population). If Strain 1 and Strain 2 exist separately in the environment, both strains will persist in the environment due to their basic reproduction numbers being greater than 1. If the two strains outbreak simultaneously, as in the case of t₀ = 0 in Figure 4, Strain 1 will eventually disappear, while Strain 2 will persist in the long term. Now, let Strain 1 “develop” in the environment for a certain period of time. In this case, the later Strain 2 enters the environment, the smaller the competition in the environment, and the higher the peak of Strain 1's infection count. On the other hand, the cross-infection cases are influenced by both the infected population of Strain 1 and Strain 2, as they are generated by the reinfection of individuals previously infected by either strain. The more infected individuals there are for both strains in the environment, the greater the number of cross-infections, which aligns with real-world observations.

4.2.   The impact of different outbreak times of the novel strain on the controlled system

After discussing the relationship between the outbreak time of the novel strain and the competition between the two strains, we now explore the relationship between the controlled System (25), which aims to minimize costs, and the outbreak time of Strain 2. From Figure 5, under the optimal control measures, the controlled System (25) reaches a steady state of local disease equilibrium, where Strain 1 becomes extinct and Strain 2 persists in the environment in the long term. Under the objective of minimizing costs, the extinction of Strain 1 is independent of the outbreak time of Strain 2, as shown in Figure 5(b). In other words, regardless of which day within the 100-day period Strain 2 outbreaks, Strain 1 will still become extinct despite having a basic reproductive number greater than 1, as depicted in Figure 5(b). The extinction of Strain 1 is accompanied by the disappearance of the co-infected population, as according to Figure 5(d). With the implementation of control measures, a great decrease in the number of cross-infected individuals is shown, dropping from a peak of 4.6 × 10⁴ to less than 10,000 individuals. In the controlled system, the peak number of infections caused by Strain 2 is higher compared to the uncontrolled System (2). However, the infection count of Strain 2 is effectively controlled in the long run, exhibiting a clear stabilization trend compared to the uncontrolled system, as demonstrated in Figure 5(c). Moreover, as the time interval between the outbreaks of Strain 2 increases, there is greater fluctuation in the appearance of susceptible individuals, and it takes longer for the system to approach stability, as shown in Figure 5(a). Overall, the controlled System (25) effectively mitigates the infectious disease outbreak compared to the uncontrolled System (2).

4.3.   Comparison between the controlled and uncontrolled systems

After implementing control measures, there is a noticeable change in the susceptible population in the environment compared to before the introduction of control measures, as illustrated in Figures 5(a) and Figure 4(a). This is due to the strong intensity of controls applied in the early stages, causing the susceptible population in the environment to rise at a constant natural growth rate until the virus develops to a point where a downward trend is observed.

The number of infections in the environment under the controlled system is not always lower than that of the original system. In the case of optimal epidemic control costs, regardless of the timing of the maximum susceptible population caused by Strain 2, it is always greater than the original System (2), as shown in Figure 6(a). When Strain 2 outbreaks occur at t₀ < 60, the maximum number of infections caused by Strain 1 in the controlled system is significantly lower than that in the original system. However, when t₀ > 60, the maximum number of Strain 1 infections in the original system is greater than that in the controlled system, according to Figure 6(b). Figure 6(c) shows that the maximum number of infections caused by Strain 2 in the controlled system is always greater than that in the original system. Under the control measures, when t₀ < 70, the maximum number of cross-infected individuals in the controlled system is significantly lower than that in the original system. However, when t₀ > 70, the maximum number of cross-infected individuals is comparable between the controlled and original systems, as shown in Figure 6(d).

Although the epidemic eventually stabilizes under the control measures, it does not necessarily mean that the cumulative number of cases in the environment will decrease. When one strain is reduced due to control measures, the resistance to the development of the other strain in the population is diminished. If the control measures are not sufficiently strong at this point, it may lead to a large-scale outbreak of the other strain in the population, resulting in a larger maximum number of infections compared to the coexistence of both strains. Although the controlled system is not superior to the original system at every moment, when the controlled System (25) stabilizes, the number of susceptible individuals in the environment is higher than that in the original System (2), indicating the effectiveness of the two epidemic control methods in curbing the outbreak.

4.4.   Relationship between control cost and timing of the new strain outbreak

In this section, we will discuss the relationship between the timing of the Strain 2 outbreak and control costs. Based on our previous discussions, we know that the control objective of the controlled System (25) is to minimize the cost. So, what is the relationship between the timing of the two strains' outbreaks and the control cost? As shown in Figure 7, we can observe that when the control duration is set to 350 days, the lowest cost is achieved when the two strains simultaneously appear or when they appear separately with a relatively long interval.

When both strains appear simultaneously in the population, the competition between the viruses is the strongest, and neither strain has fully developed in the environment. At this point, the control cost is relatively low. When the two strains appear within a time interval, the infection peaks caused by the two strains in the environment overlap, which adds pressure to epidemic control. In this case, the control cost is relatively high. However, when the outbreak of Strain 2 is sufficiently delayed, and Strain 1 has already been eliminated or is about to be eliminated during the control process, the pressure on epidemic control targeting a single strain becomes lower, resulting in a natural decrease in cost.

In the comprehensive management and control of the COVID-19 pandemic, timely adjustments in vaccination strategies based on model results are crucial. Prioritizing vaccines against the dominant strain, particularly in resource-limited scenarios, can significantly reduce infection risks. Moreover, adjusting the priority order of vaccination according to the predicted transmission trends of strains is a key strategy in reducing infection rates. Simulation results indicate that implementing control measures at different time points significantly impacts the trajectory of the epidemic. Therefore, public health policies should be dynamically adjusted based on actual infection rates and healthcare resources, such as modifying social distancing and lockdown measures according to the progression of the epidemic. Additionally, enhancing medical resource allocation and raising public awareness about epidemic prevention are key to controlling the spread of the disease. Governments and healthcare institutions should increase investment in medical resources, including increasing the number of healthcare workers, enhancing treatment facilities, and ensuring adequate supplies of medical goods and medicines. Promoting the importance of vaccination, proper personal hygiene practices, and adherence to social distancing guidelines are crucial. Even when the pandemic is under control, continuous monitoring of virus mutations and infection trends is necessary to adjust strategies promptly and prevent a resurgence of the epidemic. Long-term monitoring and assessment ensure the effectiveness of pandemic control measures and reduce the long-term threat to public health.

5.   Conclusion

The extensive spread of infectious diseases can lead to significant socio-economic challenges, with the emergence of variant strains potentially complicating epidemic prevention and control efforts, thereby causing a marked increase in both morbidity and mortality rates. This article examines the effects of unpredictable viral mutations on the model at various points in time, exploring the optimal control costs and epidemic trends under different outbreak scenarios to address the gaps in existing research on strains with varying outbreak timings. Our findings indicate that while the timing gap between two strain outbreaks does not alter the epidemic's ultimate trajectory, it does affect the cost of prevention measures during the control phase. The relationship between the interval of strain outbreaks and control costs is nonlinear, exhibiting a peak. A concentrated outbreak of two strains somewhat reduces prevention costs due to inter-strain competition; a larger interval between outbreaks may result in the mutant strain disappearing or causing a minor outbreak due to already implemented preventive measures against the original strain, thereby lowering costs. Compared to current dynamic models, this study incorporates various disease transmission-related factors, including virus mutation, reinfection due to reduced immunity in recovered individuals, population migration, and cross-infection.

Epidemic development in infectious disease research is influenced by numerous factors. However, this study only addresses a subset of epidemiological factors, leaving several others for independent investigation and further exploration, such as the competition and evolution among different types of infections, relapse and incomplete treatment, and the emergence of resistance. Additionally, determining many of the model's parameters discussed here may be challenging, inaccurate, or costly. Incorporating the uncertainty of these parameters into model development to improve prediction accuracy and applicability is a subject that merits further detailed exploration.

Use of AI tools declaration

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