The emergence of a virus variant: dynamics of a competition model with cross-immunity time-delay validated by wastewater surveillance data for COVID-19

We consider the dynamics of a virus spreading through a population that produces a mutant strain with the ability to infect individuals that were infected with the established strain. Temporary cross-immunity is included using a time delay, but is found to be a harmless delay. We provide some sufficient conditions that guarantee local and global asymptotic stability of the disease-free equilibrium and the two boundary equilibria when the two strains outcompete one another. It is shown that, due to the immune evasion of the emerging strain, the reproduction number of the emerging strain must be significantly lower than that of the established strain for the local stability of the established-strain-only boundary equilibrium. To analyze the unique coexistence equilibrium we apply a quasi steady-state argument to reduce the full model to a two-dimensional one that exhibits a global asymptotically stable established-strain-only equilibrium or global asymptotically stable coexistence equilibrium. Our results indicate that the basic reproduction numbers of both strains govern the overall dynamics, but in nontrivial ways due to the inclusion of cross-immunity. The model is applied to study the emergence of the SARS-CoV-2 Delta variant in the presence of the Alpha variant using wastewater surveillance data from the Deer Island Treatment Plant in Massachusetts, USA.


Introduction
Viruses mutate rapidly, which may impact the clinical presentation of the disease, its epidemiology, the efficacy of therapeutics and vaccinations, or the accuracy of diagnostic tools (World Health Organization, 2022).These mutations, along with selection pressures, may result in new variants (or strains) of a pathogen.After the emergence of SARS-CoV-2 in late 2019 (World Health Organization, 2020), for about 11 months, SARS-CoV-2 genomes experienced a period of relative evolutionary stasis.From late 2020, however, multiple countries began reporting the detection of SARS-CoV-2 variants that seemed to be more efficient at spreading.One of the first variants, reported on December 14, 2020 in the United Kingdom, was identified as the B.1.1.7 variant (later renamed the "Alpha" variant).Others include the B.1.351lineage first detected in South Africa and P.1 from four Brazilian travelers at the Haneda (Tokyo) airport (World Health Organization, 2022; National Institute of Infectious Diseases, Japan, 2021).Since then, the World Health Organization has defined five lineages as variants of concern (Alpha, Beta, Gamma, Delta, and Omicron) (World Health Organization, 2022).These SARS-CoV-2 variants possess sets of mutations that confer increased transmissibility and/or altered antigenicity, which the latter likely evolved in response to the immune profile of the human population having changed from naive to having been immune-imprinted from prior infections.Multiple studies have reported the rapid displacement of the Delta variant by Omicron in both clinically reported data and wastewater surveillance data (Lee et al, 2022;Wu et al, 2020).The most recent Omicron BA.4 and BA.5 lineages have also been demonstrated to resist neutralization by full-dose vaccine serum and have reduced neutralization to BA.1 infections (Tuekprakhon et al, 2022).
COVID-19 is now one of the most widely-monitored diseases in human history, allowing for unprecedented insight into variant emergence and competition.While disease surveillance often relies on clinical case data for monitoring (and genetic sequencing to identify new variants), issues related to reporting delays or the underreporting of cases can lead to inaccurate real-time data.Wastewater surveillance was previously used to detect poliovirus (Pöyry et al, 1988), enteroviruses (Gantzer et al, 1998), and illicit drug use (Daughton and Jones-Lepp, 2001); however, it was recently that it came to the forefront by helping fight against the COVID-19 pandemic.The rationale for SARS-CoV-2 detection in wastewater relies on the viral shedding mostly in feces and urine from infected individuals, which gives an alternative approach to recognizing viral presence and penetration in the community (Peccia et al, 2020;Medema et al, 2020;Ahmed et al, 2020;Fall et al, 2022).Quantification of viral concentrations in wastewater thus offers a complementary approach to understanding disease prevalence and predicting viral transmission by integrating with epidemiological modeling, while avoiding the same pitfalls associated with only considering clinical data.
Mathematical models have been used extensively in the study of disease dynamics with applications to the COVID-19 pandemic.Wastewater-based surveillance has increasingly been used in conjunction with mathematical and statistical models.McMahan et al (2021) used an SEIR model to mechanistically relate COVID cases and wastewater data.Phan et al (2023) also used a standard SEIR framework, with the addition of a viral compartment, to estimate the prevalence of COVID-19 using wastewater data; results indicated that true prevalence was approximately 8.6 times higher than reported cases, consistent with previous studies (see Phan et al (2023) and the references therein).Naturally, the SEIR model may be extended to include heterogeneity in the viral shedding other compartments (such as those hospitalized or asymptomatic) as done by Nourbakhsh et al (2022).
Other studies have focused on variant emergence and competition between multiple strains or diseases.A recent study by Miller et al (2022) used a stochastic agent-based model in an attempt to forecast the emergence of SARS-CoV-2 variants without having to previously identify a variant.The authors found that mutations are proportional to the number of transmission events and the the fitness gradient of a strain may provide insight on its persistence (Miller et al, 2022).Fudolig and Howard (2020) presented a modified SIR model with vaccination (in the form of a system of ordinary differential equations) to investigate two-strain dynamics and its local stability properties.Here, individuals infected with the established strain are immediately susceptible to infection by a new strain.The authors determined that the two strains can coexist if the reproduction number of the emerging strain is lower than that of the established strain (Fudolig and Howard, 2020).A general multi-strain model by Arruda and colleagues (Arruda et al, 2021) uses an SEIR-type model for each viral strain and uses an optimal control approach.The authors account for mitigation strategies through the inclusion of a modification terms that can reduce the contact rate of each strain, and individuals infected with a strain will have waning immunity to that same strain.However, the model does not consider cross-immunity between the strains (Arruda et al, 2021).Gonzalez-Parra et al ( 2021) developed a two-strain model of COVID-19 by extending the standard SEIR formulation to include asymptomatic transmission and hospitalization.The study found that the introduction of a slightly more transmissible strain can become dominant in the population (Gonzalez-Parra et al, 2021).These models may also include a time delay to account for various biological phenomena.For example, Rihan et al (2020) developed a delayed stochastic SIR model with cross-immunity, where a time delay was incorporated to adjust for the incubation period of a disease and stochasticity was used to determine the effect of randomness on parameters.
In this paper, we present a four-dimensional modified SIR model to study disease dynamics when two strains are circulating in a population.A time delay is incorporated to account for temporary cross-immunity induced by infection with an established (or dominant) strain.This paper is organized as follows: in section 2, the model is formulated and the equilibria of the full system are analyzed.Interestingly, we find that the time delay does not influence the stability of equilibria and is hence a harmless delay (Gopalsamy, 1983;Driver, 1972).In section 3, we introduce the transient model to study global stability of the coexistence equilibrium, and bifurcation curves are shown.Finally, the model is calibrated using wastewater data and the results are studied using a sensitivity analysis in section 4.

The general model
In this section, we introduced our mathematical model that incorporates two competing virus strains and conduct basic model analysis.
We consider a population-level virus competition model using a compartmental framework.We let S (t), I 1 (t), I 2 (t) and R l (t) be the individuals that are susceptible to both virus strains, infectious with strain 1, infectious with strain 2 and recovered from strain 1 but susceptible to strain 2 at time t, respectively.Let τ be the time it takes for an individual infected with strain 1 to become susceptible to infection by strain 2. We introduce the following two-strain virus competition model with temporary cross-immunity: An ODE version of this model without demography, independently developed, was used recently to describe the evolutionary dynamics of SARS-CoV-2 on the population level Boyle et al (2022).As a practical convention, all parameters in our model are positive.The birth rate of susceptible individuals is constant at rate a. Susceptible individuals die naturally at rate d.Infected individuals with strain 1 or strain 2 die at rate d 1 or d 2 , respectively.To investigate how disease-induced death influences virus strain competition, we make the distinction that d 1 and d 2 are disease-induced death rates, while d is the natural death rate.In practice, d 1 ≥ d and d 2 ≥ d.In system (1), susceptible individuals become infectious when they come into contact with infectious individuals from either strain at rates β 1 and β 2 , respectively.Individuals infected with strain 1 recover at rate γ 1 and enter the R l compartment where they are immune to strain 1, but become susceptible to strain 2 at rate β 2 after τ days has passed.Infectious individuals with strain 2, recover at rate γ 2 .We note that and differentiating with respect to t we have We assume that the transition rates from S to I 1 , S to I 2 and R l to I 2 follow the classical mass action law and all other transition rates are proportional to the compartment being left or entered.Figure 1 shows a summarizing schematic of the model transitions.We note that using standard incidence for the disease transmission rates would make more biological sense, since it shouldn't matter how many people have the disease around you, only how many you come into contact with.Lastly, initial histories for system (1) are prescribed by: where φ i (u), i = 1, 2, 3, 4 are bounded, continuous and nonnegative functions for u ∈ [−τ, 0).

Non-negativity and boundedness
We notice that the vector-valued function (1) and its derivative exist and are continuous.Therefore there exists a unique noncontinuable solution defined on some interval [−τ, s) where s > 0 (Kuang, 1993;Smith, 2011).Our first step is to show that the model produces solutions that are biologically plausible.We prove this with the following two propositions.We first show that if solutions start nonnegative, then they will stay nonnegative on [0, s).After that, we show solutions remain bounded for all time, which then implies s = +∞ by Theorem 3.2 and Remark 3.3 in (Smith, 2011).
Proposition 1 Solutions to system (1) that start nonnegative stay nonnegative.
We see that where This implies that Hence, solutions with nonnegative initial conditions will remain nonnegative.
Throughout the rest of this paper, we assume that S (0) > 0, Proposition 2 Solutions to system (1) are bounded.
The fact that solutions are bounded for all t ≥ 0 implies that s = +∞.

Analysis of equilibria
In order to gain a global understanding of the dynamics of system (1), we study the existence, number and stability of its equilibria.For infectious disease models, the dynamics can usually be characterized using the basic reproduction number (Delamater et al, 2019).By the next generation matrix method in (Driessche and Watmough, 2002) we find the basic reproduction numbers for strain 1 and strain 2 to be R i = aβ i d(d i + γ i ) for i = 1, 2. The full system exhibits four biologically relevant steady states: a disease-free steady state (E 0 ), two steady states where either strain 1 outcompetes strain 2 (E 1 ) or strain 2 outcompetes strain 1 (E 2 ), and a coexistence steady state (E c ).They take the following forms: where The emergence of a virus variant We see that E 0 always exists, E 1 exists when R 1 > 1, and E 2 exists when R 2 > 1.From the form of E c we see that it exists when dd 1 is always true, we can equivalently say dd 1 .We summarize the above discussion in the following proposition.
Proposition 3 The following are true for system (1).
3. The coexistence equilibrium, E c , exists exactly when then we see that f (1) = We have the following result for E 0 .
Proposition 4 E 0 is locally asymptotically stable when Proof The Jacobian matrix evaluated at E 0 is .
The corresponding eigenvalues are Therefore, E 0 is locally asymptotically stable whenever max{R In addition to local stability, we have the following global stability result for E 0 .
Proof Observe that Hence, the region {S (t) : 0 ≤ S (t) ≤ a/d} is positively invariant and attracting.
Since lim t→∞ I 1 (t) = 0, we see that for any ε > 0, there is a and a similar argument can be used to show that R l is eventually bounded by ε/d and hence R l → 0 as ε → 0 and t → ∞.
Finally, the assumption that Thus I 2 (t) is exponentially decreasing.By a similar argument as above, lim t→∞ I 2 (t) = 0. We have shown that lim t→∞ I 1 (t) = 0, lim t→∞ I 2 (t) = 0, and lim t→∞ R l (t) = 0. Thus we obtain the limiting equation dS dt = a − dS which implies that lim t→∞ S (t) = a/d.This concludes the proof.
The conditions that govern the global stability of E 0 make good biological sense.
We have the following result for E 1 .
, then E 1 is unstable.The emergence of a virus variant Proof The Jacobian matrix evaluated at E 1 is .
The corresponding characteristic polynomial factors to Therefore, the corresponding roots are and the roots to the quadratic equation Consequently, λ 1 < 0 and λ 2 < 0 by assumption (1).In addition, by the Routh-Hurwitz stability criterion for quadratic equations (Brauer and Castillo-Chavez, 2012), g(λ) has roots with negative real parts since R 1 > 1.Thus all eigenvalues have negative real part.Lastly, we see that . This concludes the proof.
This previous proposition shows that due to the immune evasion of the emerging strain, the reproduction number of the emerging strain must be significantly lower than that of the established strain for it to competitively exclude the emerging strain.In addition to local asymptotic stability, we have the following result for global stability of E 1 .
, then E 1 is globally asymptotically stable.
Proof Consider the sum S (t) + I 1 (t).Observe that where α = min{d, γ 1 + d 1 }.In practice, we assume that d 1 ≥ d and so α = d.Recall that S (0) + I 1 (0) ≤ a/d.We see that S (t) + I 1 (t) ≤ a/d; hence both S (t) and I 1 (t) are bounded above.Furthermore, because the arithmetic mean is greater than or equal to the geometric Hence we obtain, and therefore, lim sup t→∞ I 1 ≤ we see that there is small constant ε 0 > 0 such that Let ε > 0 and ε 0 = γ 1 d ε.Thus there exists t ε > 0 such that lim sup t→∞ I 1 < and we obtain Therefore, for any ε 0 there exists Hence there exists t 2 > such that for t > t 2 , R l < D + 2ε 0 .In addition, since R 2 < 1, there exists Our previous result ( * ) implies that, for this ε 1 , there exists t 1 > 0 such that S (t) < a/d + ε 1 for t > t 1 .We are now ready to control I 2 .We have the following, Thus, for t > t 1 , I 2 (t) is exponentially decreasing, implying that lim inf t→∞ I 2 (t) ≤ 0. However, since I 2 (t) is non-negative, lim sup t→∞ I 2 (t) ≥ 0. Hence lim t→∞ I 2 (t) = 0.
Observe that once I 2 goes to zero, R l does not impact the dynamics of the model, allowing us to consider the behavior of the resulting two-dimensional system: It's easy to see that this system has a positive equilibrium point, which is globally asymptotically stable.Finally, considering the limiting profile of R l we obtain, lim t→∞ R l (t) = γ 1 β 1 (R 1 − 1) .Therefore, all trajectories of system (1) tend to E 1 .
We have the following result for E 2 . .
The corresponding characteristic polynomial factors to Therefore, the corresponding roots are and the roots to the quadratic equation Consequently, λ 1 < 0 by assumption (1) and λ 2 < 0. In addition, by the Routh-Hurwitz stability criterion for quadratic equations (Brauer and Castillo-Chavez, 2012), g(λ) has roots with negative real parts since R 2 > 1.Thus all eigenvalues have negative real part.Lastly, we see that E 2 is unstable if and only if This concludes the proof.
Remark 2 From the above local stability results we see that for τ ≥ 0, the equilibria E 0 , E 1 , and E 2 do not undergo a delay-induced stability switch.For this reason, we call the delay, τ, a harmless delay (Gopalsamy, 1983;Driver, 1972).We summarize this formally in the next theorem.
Proposition 9 For τ ≥ 0, the equilibria E 0 , E 1 , and E 2 do not undergo a delay-induced stability switch.
Proof We see that by Propositions 4, 6, and 8, τ does not appear in any of the characteristic polynomials and therefore does not influence stability.
and therefore I 1 (t) < I 1 (0)e (β 1 a d +ε −(γ 1 +d 1 ) )t, and so I 1 (t) → 0 as t → ∞.Since I 1 → 0, for ε 1 > 0 there exists t 1 such that γ which implies that lim sup t→∞ R l (t) ≤ ε 1 d .Letting ε 1 → 0 we obtain lim sup t→∞ R l (t) ≤ 0. In addition, since R l ≥ 0 we have lim inf t→∞ R l ≥ 0. Therefore, lim t→∞ R l (t) = 0 and we obtain the 2 dimensional limiting system: Then the derivative with respect to time is given by We have used the steady state relationships Thus all solutions of system (10) tend to , by the Lyapunov-LaSalle Theorem.This shows that all solutions to system (1) tend to E 2 if R 1 < 1 and R 2 > 1.
Figure 3 shows the general stability and existence regions of the equilibrium points of system (1) in the R 1 R 2 plane.Bifurcations occur when crossing from one region to another.To generate these diagrams we parameterize the curve (7) by either β 1 , d 1 and γ 1 .We illustrate the bifurcation from E 0 to E 2 in Figure 2 (panel (a)), the bifurcations from E 1 to E c to E 2 in (panel (b) and (c)), and the bifurcations from E 2 to E c to E 1 in (panel (d)).
We see that the competitive exclusion principle holds for either strain as long as the conditions of either Proposition 6 or Proposition 8 hold (Gause, 1934;Bremermann and Thieme, 1989).However, conditions for one strain to competitively exclude the other are different between the two strains because of temporary crossimmunity.This also suggests that temporary cross-immunity is a mechanism for coexistence of two competing virus strains.

The transient model
To study the stability of the coexistence steady state we make the assumption that the susceptible population is at equilibrium, S (t) = S M and remove the differential equation for S .Furthermore, if I 1 0, then S M = γ 1 +d 1 β 1 .Therefore, dI 1 dt = 0 and we may remove the equation for dI 1 dt , but assume I 1 (0) > 0. We have the following 2 dimensional system of differential equations:  2 for the bifurcation diagram.We note that the curve given by ( 7) can only be plotted as a function of β 1 , d 1 , γ 1 , a and d.We only show plots for the first three since the latter two have similar geometries.3), ( 4), ( 5) and ( 6)).The rest of the parameter values were fixed at γ 1 = 0.25, γ 2 = 0.2, d = .05,d 1 = 0.15, d 2 = 0.15, a=10 and τ = 0.
The reproduction numbers for the transient model are: The emergence of a virus variant , respectively.(13)

Boundedness and positivity
We prove basic positivity and boundedness of solutions for (12).However, we note that if R1 R2 < 1, I 2 (t) becomes unbounded.
Proposition 11 Solutions to system (12) that start positive, remain positive for all time.
Proposition 12 If R1 R2 > 1, then solutions of system (12) are bounded from above. Therefore, Thus we have I 2 + R l ≤ B. Since I 2 > 0 and R l > 0, we have that both I 1 and R l are bounded above.
Lastly, we see that when R1 R2 < 1, then solutions for I 2 are unbounded.
Proof We have R1

R2
< 1, then This implies that I 2 (t) is unbounded for all t > 0.
An interesting implication of Proposition 13 is that if I 1 is not as infectious relative to strain 2, then it cannot control the spread of strain 2 and ultimately strain 2 becomes unbounded in the transient model.

Equilibria of the transient system
For our analysis we would like to have bounded solutions.For this to hold, by Proposition 12 we must have that R1 R2 > 1.Therefore, for the remainder of this section we assume that R1 R2 > 1.Assuming that I 1 (0) > 0, we find two equilibria: the coexistence equilibria, U c , and another equilibrium where the first strain exists, U 1 .They take the following form: We note that U c is dependent on I 1 (0) and is biologically relevant exactly when We note that U c can exist even when both reproduction numbers are less than 1.We have the following theorem on the stability of U c .The emergence of a virus variant Proposition 14 If U c exists, then it is asymptotically stable.
Proof U c is a positive steady state if and only if The Jacobian matrix at U c is .
We find that the trace is Therefore, both eigenvalues have negative real part and U c is locally asymptotically stable whenever it exists.
We have the following theorem on the stability of U 1 .
2. U 1 is unstable when Proof The Jacobian matrix at U 1 is We find that the eigenvalues are We see that λ 2 < 0 exactly when − 1 < 0 and unstable when  A phase portrait of the solution trajectory to the coexistence steady state is shown in Figure 5. Furthermore, it can be shown that U 1 is globally asymptotically stable under certain conditions.
Theorem 17 If β a d < b, then all solutions tend to U 1 .
Proof We prove this result by contradiction.Recall that if β a d < b, then the transient system ( 21) does not attain a positive steady state.Assume that β a d < b and lim t→∞ (x, y) (0, a/d).Furthermore, observe that lim sup t→∞ y(t) ≤ a/d.Thus for ε > 0, there exists a t * > 0 such that y(t) < a/d + ε for t > t * .With this claim, we see that, for t > t * ,  21) has no positive steady state.In other words, the claim yields lim t→∞ x(t) = 0.
Furthermore, since y is bounded, the above result implies that for any ε 1 < a, there exists a t 1 > t * such that βxy < ε 1 for t > t 1 .Therefore Letting ε 1 → 0, we see that lim inf t→∞ y(t) ≥ a/d.As well, our claim indicates that lim sup t→∞ y(t) ≤ a/d.Hence lim t→∞ y(t) = a/d.
In the following, we prove our claim.The proof is divided into three cases: 1. y(0) ≤ a/d; 2. y(0) > a/d and there exists a t 2 > 0 such that y(t) > a/d for t ∈ [0, t 2 ) and y(t 1 ) = a/d; 3. y(t) > a/d for all t > 0. The emergence of a virus variant We consider case 1.We have Hence y(t) < a/d for t > 0, and our claim is true.Consider the second case.From case 1, we see that y(t) < a/d for t > t 1 and again our claim is true.
Finally, consider case 3. Here, dy/dt < 0 and there is a y c ≥ a/d such that lim  From Propositions 14 and 16 we see that both virus strains can coexist as long as the original strain has a higher reproduction number than strain 2 and However, we may solve for R2 in terms of R1 to generate a bifurcation curve between coexistence and competitive exclusion, Figure 6 shows the R2 R1 −bifurcation plane where equation ( 22) is parameterized by For the two strains to coexist together, strain 1 needs to have a higher basic reproduction number than strain 2. However, it can't be too high relative to strain 2 or it will force strain 2 to extinction.The unbounded region corresponds to Proposition 13.In general, the model suggests that viruses which mutate into strains that are slightly less infectious are more likely to coexist together.On the other hand, viruses that mutate into strains that are sufficiently less infectious relative to the original strain, will out-compete the mutated strain.4 Numerical results

Data fitting
The system (1) is validated by fitting to wastewater data from October 1, 2020 to May 13, 2021 obtained from the Deer Island Treatment Plant in Massachusetts (Xiao et al, competitive exclusion of I 2 by I 1 .We do not plot S or I 1 since they are held constant at S M and I 1 (0), respectively.
Table 1 Established results and open questions

Conditions
Results or question See theorem 7 E 1 is globally stable 5.
See Fudolig and Howard (2020) Local stability of E c 7.
Open Global stability of E c 8.
Open Global stability of E 2 with R 1 > 1. 10. Open Influence of τ on the stability of E c . 1.
Inequality (19) U c is globally stable 2022).This plant serves approximately 2.3 million people in the greater Boston area (Xiao et al, 2022).More information on the collection and processing of wastewater samples can be found in Xiao et al (2022).Fitting to wastewater data, as opposed to incidence or mortality data, allows us to avoid underreporting issues related to clinical reporting.The B.1.1.7 (Alpha) variant was detected in Massachusetts in January 2021 (Massachusetts Department of Public Health, 2021), while the B.1.617(Delta) variant was found in the state in April 2021 (Markos, 2021).It should be noted that Massachusetts (population size 7 million) began vaccinating healthcare workers on December 15, 2020 during Phase 1 of the state's vaccination plan (Massachusetts Department of Public Health, 2022).For simplification purposes, we assume that individuals in the susceptible (S ) and recovered (R l ) compartments are vaccinated at a rate v and that the vaccine offers immediate protection from both strains.Individuals who have recovered from the emerging strain are not tracked or vaccinated for several reasons.In the presented model, these individuals are removed from the population and thus do not impact infection dynamics.It has been shown that two vaccine doses provided significant protection against the Alpha and Delta variants with respect to infection and hospitalization (Gram et al, 2022).Although protection against infection has been found to wane over time, Gram et al (2022) found that, after 120 days, vaccine efficacy against Delta decreased from 92.2% to 64.8% in those aged 12 to 59 years.Vaccine efficacy in individuals over 60 years of age saw decreases in efficacy from 90.7% to 73.2% and 82.3% to 50% for Alpha and Delta, respectively (Gram et al, 2022).Due to the limited time-scale of vaccination in the model and the scope of this study, we assume protection does not wane.
Based on data on fully-vaccinated individuals (defined as those who received all doses of the vaccine protocol) from the U.S. Centers for Disease Control and Prevention (U.S. CDC), compiled by Our World in Data (Mathieu, 2022; U.S. Centers for Disease Control and Prevention, 2022), we fix the per capita vaccination rate at v = 0.0038 per day with vaccination beginning on January 5, 2021 due to the three week time period between first and second doses (Massachusetts Department of Public Health, 2022).The calculation of v is shown in Figure 8.In order to fit system (1) to the wastewater data, we add a compartment C V (t) denoting the cumulative viral RNA copies in the wastewater following the formulations in Saththasivam et al (2021).Hence, the dynamics of the cumulative virus released into the wastewater is governed by where α denotes the fecal load per individual in grams per day, δ denotes the viral shedding rate per gram of stool, and (1−η) denotes the proportion of RNA that arrives to the wastewater treatment plant.Because the wastewater data is daily and C V (t) is The emergence of a virus variant cumulative, the objective function be minimized is given by , new viral RNA entering the sewershed on day t n ).Parameter estimation is carried out using fmincon and 1000 MultiStart runs in Matlab.For comparison purposes, both the ODE and DDE versions of the model were fit to the data.Initial values for I 1 and I 2 are estimated by using the initial viral RNA data and the estimated values of α, δ, and η; that is, the constraint and assuming that I 1 (0) ≥ I 2 (0).For the model with time delay, the same constraints are used for the initial histories.Values for estimated and fixed parameters are listed in Table 2. Figure 9 depicts model simulations without time delay using the best-fit parameters when compared to daily wastewater data (Figure 9a) and seven-day average case data (Figure 9b).The ODE version of the model predicts peak new infections on December 29, 2020, preceding the daily reported case data by 11 days.Due to the unreliability in the case data, however, this 1.5 week difference may be reasonable.Furthermore, the model projects approximately six times more new cases than the reported case data at their respective peaks.
Best-fit simulations with time delay are shown in Figure 10.Here, the model predicts daily incidence peaking on January 4, 2021, approximately five times higher than the reported cases on January 9, 2021, a difference of 5 days.Unlike the ODE version, the inclusion of time delay allows the model to capture the decline of the Alpha wave, but both the ODE and DDE versions of the model are unable to capture the Delta wave.

Sensitivity analysis
In this section, we carry out a local sensitivity analysis to explore which parameters are the most important to model dynamics.We use a normalized sensitivity analysis so that the sensitivity coefficients are not affected by parameter magnitude.Here, the normalized sensitivity coefficients are given by (Saltelli et al, 2000): where p and Y denote the parameter and response of interest, respectively, and ∆p is the perturbation size.Each parameter is varied by 1% individually from the values listed in Table 2 while all other parameters are fixed.Here, the response variable Y is cumulative cases evaluated at steady state.We ignore the parameters related to wastewater (α, δ, and η) since they do not impact disease dynamics in the analysis.Results are shown in Figure 11.The height of the bars indicates how sensitive the response variable is to the parameter; the direction of the bars (or sign of the sensitivity coefficient) indicates the direction of correlation.
The ODE and DDE versions of the model display significant sensitivity to the strain-specific contact rates (β 1 , β 2 ) and the strain-specific recovery rates (γ 1 , γ 2 ); the DDE version of the model has increased sensitivity to the initial number of those infected with strain 1 compared to the model without time delay.Furthermore, model dynamics, independent of time delay, are only slightly (if at all) impacted by changes in the strain-specific mortality rates (d 1 , d 2 ).

Discussion
In this paper, we have constructed a mathematical model describing two-strain virus dynamics with temporary cross-immunity.Although this general framework is applicable to many diseases, we put our model into the context of the COVID-19 pandemic and connected infectious individuals with wastewater data.The model produces rich long-term dynamics that include: (1) a state where the two strains are not infectious  2 enough and are cleared from the population, (2) two competitive exclusion states where one of the strains is more infectious than the other and ultimately forces the other to extinction, and (3) a coexistence state where the two strains coexist together.By using a quasi-steady state argument for S we reduced the four dimensional system (1) to a two dimensional system (12).This simpler system exhibited a competitive exclusion equilibrium where the first strain forces the second strain to extinction and a coexistence equilibrium.Results and open questions are summarized in Table 1.
The model presented in this study uses a time delay to account for cross-immunity between two strains and is shown to be a harmless delay since it doesn't influence the stability of the boundary equilibrium points (Gopalsamy, 1983(Gopalsamy, , 1984;;Driver, 1972).However, the time delay's influence on the stability of the coexistence equilibrium is an open question.This time delay acts as a definitive period for immunity as opposed to a continuous or distributed waning of protection (Pell et al, 2022).For comparison, we simulate the ODE version of the model (1) with the β 2 R l I 2 terms replaced by β 2 R l I 2 in order to study the effects of waning immunity, as shown in Figure 12.As → 0 (i.e. the waning period for cross-immunity increases) it is shown that the emergent strain requires more time to be established in the population if all parameters between the two strains are equal.Additionally, we can interpret the term εβ 2 I 2 R l as the number of new breakthrough infections that occur per time unit.As ε increases to 1, the more likely breakthrough infections will occur.Using a similar model that does not account for demography, Boyle et al., showed that the rapid turnover from one variant to another is influenced by two components: the increase in transmissibility and the breakthrough infections (Boyle et al, 2022).They deduce that emergent strains are the ones that are best at evading immunity (Boyle et al, 2022).Our simulations in Figure 12 further support this.
We fit the model (1) to wastewater data from the greater Boston area in order to show that the model can capture two-strain dynamics in the real world.Using wastewater data for fitting, as opposed to clinical data, allows us to avoid issues related to under-reporting or reporting lags of case data.We fit the model with and without time delay and found that incorporating time delay allowed the model to better follow the trend of the data for the first wave.However, regardless of the inclusion of time delay, the model did not qualitatively capture the second wave in the data (although the model with time delay performed slightly better).This is due to the models not accounting for the vaccination program that began around December.Although the issue of parameter identifiability is present (and beyond the scope of this study), we ultimately show that this four-dimensional model is able to capture complex two-strain dynamics.A local sensitivity analysis was carried out on the values obtained via curve-fitting and indicated that cumulative infections are sensitive to strain-specific contact and recovery rates.
This paper may also be viewed as an extension of the work done by Fudolig and Howard (2020).While Fudolig and Howard did consider cross-immunity because they focused on SARS-CoV-2 and influenza strains co-circulating, we incorporated a time delay to account for one SARS-CoV-2 strain providing temporary immunity to another.Although that study included a compartment for vaccinated individuals, more direct comparisons may be made by setting the vaccination rate of their model (p) and the time delay of the model presented here (τ) to zero.The authors derived the same local stability results for the disease-free and emergent strain (strain 2) equilibrium, and also found that the reproduction number of the emergent strain (strain 2) must be sufficiently small in order for the local stability of the established strain boundary equilibrium (Fudolig and Howard, 2020).Furthermore, we provide global stability results for the boundary equilibria.Our bifurcation plane for the full model, shown in Figure 3, mirrors that of Fudolig and Howard (2020).In addition, we provide an analogous bifurcation plane for the transient system (12) in Figure 6.
SARS-CoV-2-infected individuals always go through a latent period, where they are yet to be transmissible clinically.This duration is related to the number of infectious viruses (or the within-host viral load) and should not be confused with the sub-clinical symptomatic phase, which can follow the latent period (Ke et al, 2021;Heitzman-Breen and Ciupe, 2022).Existing models examining SARS-CoV-2 transmission often consider latency, which better integrates epidemic data (Phan et al, 2023;Patterson and Wang, 2022;Eikenberry et al, 2020).However, for our analytical purposes, the inclusion of latency can complicate the mathematical analysis but usually has a small effect on the basic reproduction number and often does not affect global stability (Patterson and Wang, 2022;Van den Driessche, 2017;Feng et al, 2001).A similar simplification to facilitate model analysis was also done by Boyle et al (2022).Thus, we made the simplifying assumption to not include latency in our current model.
In general, immunity against one strain may not confer protection for a different strain, if the two strains are sufficiently different from one another.However, while the initial infection may be due to a single strain, mutations occur during the course of infection and may allow for the development of antibodies to various mutations of the initial strain, which may include the particular second strain.Yet more paradoxically, antibodies obtained from one strain may enhance the infection of another, which is known as the antibody-dependent enhancement of infection phenomenon (Junqueira et al, 2022;Maemura et al, 2021;Wan et al, 2020;Nikin-Beers and Ciupe, 2015).The evolutionary dynamic of SARS-CoV-2 itself is interesting and quite complex and should vary from individual to individual.Instead, the motivation for our model comes from the scenario when a mutant strain begins to emerge while another strain is dominant, as is the case of Alpha and Delta variants.In particular, taking into account the timing (e.g., the beginning of Delta vs. the end of Alpha) and scale differences in the number of infected individuals from each variant, we assume individuals recovered from Alpha can lose immunity and get infected with Delta during this time period.On the other hand, we assume individuals who recovered from Delta may not get infected with Alpha.Due to this reason, we chose not to include vaccinations of individuals recovered from strain 2. In particular, if an individual is recovered from the emerging strain, regardless of the particular emerging variant, they may have some protection from the dominant strain.By the time the protection of this individual wanes, there should be much fewer individuals infected by the originally dominant strain to consider reinfection as a viable path of infection.Throughout the course of the SARS-CoV-2 pandemic, we have never observed a strain become dominant for multiple periods.Future work may consider extensions to these aspects of our model.
Ultimately, the model developed here, although simple in appearance, exhibits rich dynamics and, with the inclusion of wastewater-based epidemiology, is capable of capturing interactions of two strains circulating in the community.Future extensions of the model may include more than two strains and use standard incidence.For example, a model with N strains may include N infectious compartments but N or fewer recovered compartments, depending on how cross-immunity is modeled.It may also be desirable to include a mutation factor to study the emergence mechanisms of various strains.Another fruitful direction would be to more realistically model the temporary cross-immunity period using a distributed delay framework.

Fig. 1
Fig.1Schematic of the general model, system (1).Solid arrows correspond to disease-related transitions, dashed lines correspond to demographic transitions (birth and natural death) and the squiggle arrow corresponds to the time delay for an individual that has recovered from strain 1 to be susceptible to strain 2.

Fig. 3
Fig. 3 System (1) stability and existence regions of equilibria in the R 1 R 2 plane.(a) Parameterized by β 1 .(b) Parameterized by d 1 .(c) Parameterized by γ 1 .For example, starting in the E 0 stability region and increasing β 2 ultimately produces a bifurcation as E 0 loses stability and E 2 gains stability.This is also illustrated with panel (a) of Figure 2. Starting in the E 1 stability region and increasing β 2 produces a bifurcation as E 1 loses stability, E c appears and ultimately for higher β 2 values E 2 becomes stable.See panel (b) of Figure2for the bifurcation diagram.We note that the curve given by (7) can only be plotted as a function of β 1 , d 1 , γ 1 , a and d.We only show plots for the first three since the latter two have similar geometries.
(x * , y * ) = U c .The system is the exact same as system (10) and hence a Lyapunov function is V(x, y) = x − x * ln(x) + y − y * ln(y).

Fig. 6
Fig. 6 The R2 R1 -plane for the transient model (system 12).The model exhibits 3 different dynamics: 1) unbounded solutions, 2) coexistence of the two virus strains and 3) competitive exclusion of the 2nd strain.(a) Stability plane parameterized by β 1 ; (b) Stability plane parameterized by d 1 ; (c) Stability plane parameterized by γ 1 ; If starting in the competitive exclusive region where I 1 is the long-term winner and then increasing β 2 (thus increasing R2 ) we see that a bifurcation occurs when R2 = R1 γ 1 I 1 (0) d R1 + 1 and the second strain can coexist with the first strain.Increasing β 2 even more will ultimately lead to another bifurcation where I 2 becomes unbounded (by Proposition 13).

Fig. 7
Fig.7Long-term dynamics of system (12).(a) coexistence steady state where I 1 and I 2 coexist.(b) competitive exclusion of I 2 by I 1 .We do not plot S or I 1 since they are held constant at S M and I 1 (0), respectively.

Fig. 8
Fig. 8 Line of best-fit (blue) compared to Massachusetts vaccination data.The per capita vaccination rate v = 0.0038 is given by the slope of the best-fit line divided by the total state population of 7 million.

Fig. 9 a
Fig. 9 (a) Best-fit model without time delay compared to wastewater data (SSE = 8.6230).(b) Model output of total daily new cases of I 1 and I 2 compared to seven-day average of new reported cases.Dotted lines indicate date of maximum reported cases for the data (orange) and the model (blue).(c) Model output of strain 1 (solid line) and strain 2 (dashed) line over time

Fig. 11
Fig. 11 Local normalized sensitivity analysis with respect to cumulative steady state cases for the model (1) (a) without time delay, and (b) with time delay.Parameters are varied by 1% one at a time.Baseline values are listed in Table 2