Co-circulation of SARS-CoV-2 and Influenza under vaccination scenarios

The interaction and possibly interference between viruses infecting a common host population is the problem addressed in this work. We model two viral diseases both of the SIRS type that have similar mechanism of transmission and for which a vaccine exists. The vaccine is characterized by its coverage, induced temporal immunity and efficacy. The population dynamics of both diseases considers infected individuals of each disease and hosts that are susceptible to one but have recovered from the other. We do not incorporate coinfection. We postulate two main transmission factors affecting the effective contact rates: i) that the virus with higher reproduction number can superinfect the one with lower reproduction number and ii) that there is some protection induced by vaccination against the weaker virus that reduces the probability of infection by the stronger virus. Our results indicate that coexistence of the viruses is possible in the long term. The time-dependent effective contact rate may induce either alternating outbreaks of each disease or synchronous outbreaks. We also found the existence of bi-stability triggered by a backward bifurcation, conducive to scenarios where, at the population level, vaccine application may promote persistence of both diseases provided the effective coverage and vaccine efficacy are low.


Introduction
The COVID-19 pandemic has entered worldwide, into a second phase of growth. In the Northern hemisphere, this phase coincides with the winter season which normally presents an increase of influenza and influenza-like illnesses (ILIs). However, this year has been atypical at least in two ways. One is the SARS-CoV-2 pandemic becoming the dominant and most prevalent respiratory viral infection from the beginning of 2020. The other is the characteristic absence of a significant number of influenza cases; influenza activity has almost null in the southern hemisphere and now in the northern, while the SARS-CoV-2 pandemic is active [1]. Figure 1 exemplifies this in the case of Mexico However, it is plausible that in the winter months (Northern hemisphere) with the many holidays that occur along this period, may finally produce influenza outbreaks that could seriously compromise the hospital capacity in Mexico and elsewhere [2,3]. The occurrence of a syndemic episode with cocirculating influenza and SARS-CoV-2 viruses is a potential and serious reality.
At the time of writing (early winter 2020), several SARS-CoV-2 vaccines are about to be deployed and vaccination against influenza and SARS-CoV-2 has already started. The COVID-19 vaccines will not be deployed to the general population in many parts of the world until later in 2021 [4,5] and, from now to that unknown date, humanity must rely in the use of non-pharmaceutical interventions (NPI) to prevent contagion and to control a likely rebound in the number of cases. In several parts of the world, we are not succeeding. Mask use, contact trace/isolation and social distancing have been relaxed. On the other hand, vaccines are not the silver bullet or panacea that will resolve the COVID-19 crises around the world. Vaccines may have a protective effect on COVID-19 via immune responses or indirect effects by reducing the burden of viral and bacterial respiratory diseases on individual patients among other things [6]. The impact of vaccination on the transmission of SARS-CoV-2 will start slowly and build up over several months to reach target coverage levels. The duration of immunity induced by a given COVID-19 vaccine is not yet clear [7]. Recently it has published that 4-10 weeks after SARS-CoV-2 infection, IgG titer decrease by 50% [8]. Because of the vaccine characteristics [5], we need to evaluate them taking into account, together with the transmission dynamics of the disease, vaccine properties such as efficacy, existing population immunity, coverage and temporal immunity, both natural and vaccine induced. Vaccines also confront problems arising from public opinion. There are, as reported in the news, population sectors that either reject taking the vaccine or do not trust the official programs that are being designed to deploy it.
Motivated by this situation, we present here a mathematical model geared to investigate the possible interference and competitive interactions between influenza and SARS-CoV-2 in the presence of a vaccine for both viruses.
The interaction between viral species follows known patterns of coexistence. [9] were the first to show that the inclusion of superinfection makes coexistence between competing Figure 1: Atypical influenza evolution. Panel A shows the number of weekly influenza cases reported from January 5, 2019, to December 12, 2020. It can be seen that the number has significantly dropped during 2020. On the other hand, panel B shows COVID-19 cases from February 22, 2020, to December 5, 2020. It is clear that the number of cases is rising, entering a second wave. In panel C it can be seen that the number of cases in 2020 (red curve) is atypically low compared to previous years since the 22nd epidemiological week (late May 2020). species possible. Superinfection is a process by which a competitive hierarchy exists among a group of species that compete for the same resources. The hierarchy facilitates coexistence and prevents competitive exclusion of the weaker species in a process that is mainly driven by the relative magnitude of the reproduction numbers of both species [10] that reflect the abilities to use a contested common resource. Seasonal influenza has a median basic reproduction number of 1.28 [11]; while for SARS-CoV-2, the median R 0 is 2.79 [12,11]. This difference in transmission potential supports the assumption that while competing for hosts, their common resource, the likelihood of long-term coexistence and cocirculation is high.
The process of superinfection has been identified as a promoter of coexistence of pathogen strains in a given host population. The structure of the models that have been used to theoretically demonstrate this property [9,13,14], has been also applied to explain the organization of community structure in ecology for species that inhabit a common landscape but are not necessarily closely related [15]. For infectious diseases the concept of superinfection has pervaded theoretical explanations for coexistence, in the same host population, of variants of a given pathogen [16,17]. The pioneering model in [13] versed on infectious diseases and was used to explore a plausible hypothesis for the length of the latent period of HIV before the onset of AIDS. Simultaneously, these same authors addressed the problem of community structure postulating a trade-off between colonization and extinction [18] in the presence of a hierarchy of competitive abilities for the exploitation of the resources in a common landscape [19]. The possibility of coexistence then, is driven essentially by the hierarchy of basic reproduction numbers where a species with higher reproduction number is capable of displacing from the host a pathogen of a different strain or species but with lower reproduction number. The opposite is not possible, that is, the displacement of a species with a higher R 0 by another species with a lower one. For acute respiratory infections, this competitive hierarchy that henceforth we will call superinfection, has been postulated to explain the alternating dynamics between influenza and RSV (respiratory sincytial virus). It is known that the reproduction number of influenza is higher than that of RSV. This fact, coupled with weather variability (seasonality), induces alternating patterns where influenza and RSV infections have a limited temporal overlap [16].
As winter settles down in the Northern hemisphere, the current SARS-CoV-2 pandemic has generated preoccupation since it, very likely, will cocirculate with other respiratory infections, particularly influenza and RSV [20,21,22,23,24], although interactions with other respiratory disease are certainly possible [25].
Despite that there are several reports on coinfections, that is the coincidental infection on the same patient between SARS-CoV-2 and other respiratory pathogen, particularly influenza [26,27,28], this is not happening yet in community context. We center in the cocirculation of both viruses and on the conditions that must hold, from the ecological perspective, in order to achieve coexistence or the extinction of one of the viruses.
In what follows we present a model for the interaction of influenza and SARS-CoV-2 in the presence of vaccination where the vaccine is characterized by coverage, efficacy and induced temporary immunity.

Model set-up
We formulate a mathematical model considering the simultaneous presence of two viruses (influenza and SARS-CoV-2) and vaccination for each one. We do not distinguish between influenza subtypes and, therefore, we approximate its dynamics as if it were a single viral population. This allows us to take influenza as a disease with recovery, for a host can be yearly infected with a new strain that we call, generically, as influenza. SARS-CoV-2 is also assumed to induce temporal immunity. The model then, is a coupled system of two SIRS (susceptible-infected-recovered-susceptible) equations ( Figure 2). Figure 2: A compartmental mathematical model for two coupled SIRS. s, i, y, r i , r y , y r i , i ry and r represent the population of susceptible, infected by influenza, infected by SARS-CoV-2, recovered from influenza, recovered from SARS-CoV-2, recovered from influenza but infected with SARS-CoV-2, recovered from SARS-CoV-2 but infected with influenza, and recovered from both viruses. Here, r k includes population recovered after infection and successfully vaccinated by virus k, where k represents influenza (i) or SARS-CoV-2 (y). Dashed blue lines represent vaccination dynamics for both viruses and dashed red line denotes superinfection process.
We have normalized the total population so that N = 1. The model is the following: In eq. (1), s, i, y, r i , r y , y r i , i ry and r represent the population of susceptible, infected by influenza, infected by SARS-CoV-2, recovered from influenza, recovered from SARS-CoV-2, recovered from influenza but infected with SARS-CoV-2, recovered from SARS-CoV-2 but infected with influenza, and recovered from both viruses (see Table 1).
Note that the recovered population for each virus can be infected (vaccinated) by (against) the other virus. Likewise, it is assumed that SARS-CoV-2 may displace or competitively exclude influenza from its host, but with somewhat diminished efficacy in the 5 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted January 4, 2021. Recovered from SARS-CoV-2 y r i Recovered from influenza but infected with SARS-CoV-2 i ry Recovered from SARS-CoV-2 but infected with influenza r Recovered from both viruses  [29,30,31]. This situation is explored in the last section.

Parameter
Parameter description β k (t) time-dependent effective contact rate of virus k α reduction in susceptibility for superinfection φ k effective coverage rate of vaccine against virus of type k q k effectively vaccinated population proportion EVP T time to achieve a coverage rate φ k given q k p i protective effect against infection by the coronavirus η k recovery rates from virus k for primary infections ω k recovery rates from virus k for secondary infections θ k loss of immunity rate of recovered from infection with virus k θ loss of immunity rate once an individual had both infections Table 2: Model parameters with k = i (influenza) and k = y (SARS-CoV-2).
We define p i , 0 ≤ p i ≤ 1, a parameter that measures the protective direct or indirect effects that being influenza-vaccinated conveys against infection by the SARS-CoV-2 [6]. We consider that this protective effect can be associated with a reduction in the effective contact rate of SARS-CoV-2 when in contact with an influenza-vaccinated host r i . We assume that there are reinfections by both viruses due to the non-lasting immunity conferred by previous infections. Thus, θ k is the loss of immunity rate of the recovered individuals from infection with virus k, η k and ω k are the recovery rates from virus k for primary and secondary infections, respectively, and θ is the loss of immunity rate once an individual 6 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted January 4, 2021. ; had both infections (see Table 2).
We are not defining a compartment for vaccinated individuals apart from the one for natural recovered individuals. Instead, we collect immunized individuals both with vaccine induced protective immunity, and natural immunity. The compartments r i and r y contain, therefore, those individuals that have either been vaccinated against influenza or SARS-CoV-2, respectively or that, alternatively, are recovered from a natural infection for either of these two viruses. This choice of compartments diminishes the system's dimensionality. To finish this section we remark that the coverage rate, φ k where k denotes either influenza (i) or SARS-CoV-2 (y) is an effective coverage rate, meaning that the vaccinated individuals are effectively protected. In other words, φ k incorporates vaccine efficacy. We further discuss this point in the following section.

Model parametrization
We estimate the baseline parameters from bibliographical sources. According to [22], for seasonal influenza, the median reproduction number for influenza is 1.28 and that for SARS-CoV-2 is 2.79, but other sources report large variability [32]. Likewise, [11] citing several sources, reports an R 0 for influenza in the range 1.06-3.4 with a mean of 1.68. The same source gives an incubation period in the range of 1-6.3 days with a mean of 2.61 and an infectious period (η i ) range of 1-9 days with a mean of 4.58 days. For SARS-CoV-2, [33] reports an incubation period of 3 to 4 days and an infectious period (η y ) of 4-5 days but [11] gives an incubation period in the range 1.9 to 14.7 days with a mean of 5 days, and an infectious period in the range 7-35 days with a mean of 15.2 days.
With respect to vaccination, susceptible individuals are constantly leaving this compartment at a rate −φ k S. Thus the probability of having been vaccinated at time t is 1 − exp(−φ k t). Therefore, if we wish that a proportion q k of the susceptible population is vaccinated at time T , then we set a vaccination rate such that The parameter q k is the effectively vaccinated proportion or EVP in the following sense: if the vaccine has a σ k % efficacy and we apply the vaccine toq k % of the population, then only a fraction q k = σ kqk is effectively protected. Influenza vaccine efficacy varies every year.
In 2019-2020 the [34] reports σ i = 0.29 but in 2010-2011 σ i = 0.6. [35,36] have argued that efficacy declines because of waning immunity that may last 6 months for influenza A(H1N1) and influenza B and at least 5 months for influenza A(H3N2) viruses. Also, we have postulated a parameter p i that mimics a protective role conveyed by influenza vaccination against SARS-CoV-2 infection. This hypothesis is based on the work of [37] that reports that protective influenza vaccination does not negatively affect the risk of contracting coronaviruses. We explore the possibility that the effect is positive, diminishing the risk of SARS-CoV-2 infection.

7
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted January 4, 2021. ; https://doi.org/10.1101/2020.12.29.20248953 doi: medRxiv preprint For SARS-CoV-2, one vaccine (Pfizer) is being deployed and others are about to be deployed. [38] assumes that efficacy will be in the range of 10-90% with more likely scenarios of 70%. The Pfizer, Moderna and likely the Astra-Zeneca vaccines have efficacies at the upper end of this interval. On the other hand, coverage has three basic scenarios: low 20%, medium 50%, and high 80%. Given the form in which we are modeling coverage, we define that the above percentages can be reached after T = 90, 180 and 365 days. For SARS-CoV-2, the duration of immunity ranges from half a year to lifelong, with a more likely scenario of one year. These estimates are largely based on data on immunity to other coronaviruses [39]. Currently, whether past infections will prevent severe COVID-19 on re-exposure to SARS-CoV-2 is not known. Finally, we postulate here that, given that the basic reproduction number of the coronavirus is higher than the corresponding one to influenza, there exists the possibility that the SARS-CoV-2 competitively displaces influenza from an infected host rendering it with a SARS-CoV-2 infection.

Results
In this section, we present the main results of our analysis.

Equilibrium analysis
First we briefly characterize some basic properties of the solutions of model eq. (1).
The proof is immediate and follows from Proposition A.17 in Appendix A of [40].

Disease-free equilibrium and local stability
With vaccination, the disease-free equilibrium of eq. (1) Note that the susceptible population at the disease-free equilibrium depends on the parameters for coverage and immunity for both viruses .
8 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted January 4, 2021. ; https://doi.org/10.1101/2020. 12.29.20248953 doi: medRxiv preprint For completeness, in the absence of vaccination for either one of the viruses, and for constant effective contact rates, the basic reproductive numbers are the classical expression for an SIR epidemic. Thus, assuming a constant effective contact rates The disease-free equilibrium E 0 (eq. (3)) is locally asymptotically stable if and only if T 0i < 1 and T 0y < 1.
To prove the claim we compute the jacobian matrix of eq. (1) by substituting s = where s * is the value of s at the disease-free equilibrium. The eigenvalues of J E 0 are given by those of A and C. A is a diagonal matrix and therefore its eigenvalues are negative if and only if T 0i < 1 and T 0y < 1. Matrix C does not depend on the transmission rates. Two of its eigenvalues are −µ − ω i and −θ − µ. The remaining ones are the roots of a polynomial of the form with strictly positive coefficients Applying the Descartes rule of signs, this polynomial has no complex and no real positive roots. Therefore the DFE is locally asymptotically stable if T 0i < 1 and T 0y < 1.
Remark: This result indicates that the disease-free equilibrium E 0 , loses its local stability whenever either T 0i or T 0y , are greater than one. Note that these threshold numbers 9 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted January 4, 2021. ; https://doi.org/10.1101/2020.12.29.20248953 doi: medRxiv preprint depend on s * . Vaccination plays a role in the local stability of E 0 through the equilibrium value of s * . The other coordinates of E 0 namely, r * i , r * y and r * depend also on the vaccine parameters.

Reproduction number
Given the interaction of both viruses, their reproduction numbers may present broader spectra of interrelation regarding the possibilities of coexistence, competitive exclusion or extinction. In this section, we give a first characterization for these properties. When the diseases have not yet invaded the host population, but the host is vaccinated against both viruses, we can compute the corresponding vaccine reproduction number. We proceed as in [41] to obtain: where s * , r * i and r * y are defined in (3). Then, the vaccine reproduction number is given by the spectral radius of matrix F V −1 . Thus, where the specific vaccine reproduction numbers for influenza (R vi ) and SARS-CoV-2 (R vy ) are: Remark: Since 0 < s * ≤ 1, from eq. (6), it is not clear the order relation between R vk and R 0k , with k = i, y. Figure 3 shows that taking µ = 0.000039139, β i = 0.3, β y = 0.2, η i = 1/5, η y = 1/14, θ = 1/365, θ i = 1/365, θ y = 1/180, ω i = 1/7 and ω y = 1/16, and vary EVP q from 1% to 99%, there exist values where possibly R vk > R 0k with k = i or y, that is to say, vaccine application may enhance the occurrence of an outbreak. It is therefore, important to give conditions that guarantee R vk < R 0k . This will be explored in section 5. For the moment, we have: Let T 0k be defined as in lemma 3.2. The vaccine reproduction numbers for influenza and SARS-CoV-2 R vi , R vy respectively, have the following characterization.
Remark: This case corresponds to the inequality T 0k > 1 and R vk < 1, when T 0k > R vk holds. Under vaccination, the disease-free equilibrium point is unstable.

10
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted January 4, 2021. ; Remark: To satisfy the condition T 0k > R vk two cases were obtained either both 11 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted January 4, 2021. ; basic reproductive numbers are greater than one T 0k > 1, R vk > 1 or both are less than one T 0k < 1, R vk < 1. This means that two different cases can occur: an scenario where both diseases are eradicated or the opposite case.
Remark: We have R vk > T 0k implying that the vaccine reproduction number can be higher than the basic reproduction number that leads to an unwanted scenario where vaccination may be counterproductive in the following ways: i) T 0k < 1, R vk < 1, ii) T 0k > 1, R vk > 1 and iii) T 0k < 1, R vk > 1. This last scenario shows that there might be cases where the disease-free equilibrium point is only a local attractor but the disease is still invading the host population.
In eq. (2), q k can be rewritten as σ kqk , where σ k is vaccine efficacy against virus type k andq k is the vaccine coverage for virus type k. R V can, therefore, be seen as a function of SARS-CoV-2 and influenza vaccine efficacy σ k , with k = i, y (figure not shown). In this case, the maximum of R V is reached, when p i is zero, for low SARS-CoV-2 vaccine 12 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted January 4, 2021. ; https://doi.org/10.1101/2020.12.29.20248953 doi: medRxiv preprint efficacy and medium influenza vaccine efficacy. As the protective effect of p i increases, this situation changes. The maximum of R V occurs when the vaccine efficacy is zero for both viruses (recall that natural temporal immunity is always present which makes R V still valid). For high vaccine efficacies, the R V is always less than one, as expected.

Asymptotic dynamics
In Section 3.1.2 we have shown that there is a wide range of parameter combinations such that R vk < R 0k , with k = i, y, that is, where vaccination is effective in reducing transmission. The other possibility, namely, R vk > R 0k , with k = i, y, may also occur. This case is explored in Section 3.2.1. So, assuming R vk < R 0k , with k = i, y, and fixing the parameters µ = 0.000039139, α = 0.5, η i = 1/5, η y = 1/14, θ = 1/365, θ i = 1/365, ω i = 1/7, ω y = 1/16, β i = 0.3, β y = 0.2, p i = 0.5 and θ y = 1/180 with the effective contact rates β i and β y constant, our numerical exploration of eq. (1) allows us to postulate the diagram shown in Figure 6A. In summary, when both vaccine reproductive numbers are less than one, both epidemics die out, when only one of the vaccine reproductive numbers is greater than one, then the virus associated with that reproductive number persists and the other goes extinct and, finally, if both vaccine reproductive numbers are greater than one, coexistence of both diseases ensues. Figure 6: Summary of the asymptotic dynamics when varying both vaccine reproduction numbers. Left panel shows the case for SARS-CoV-2 and Influenza when R vk < R 0k , with k = i, y. Right panel illustrates the case for two hypothetical viruses A and B when R 0k < R vk , k = A, B as in Section 3.2.1 Figure 7 shows the case where only SARS-CoV-2 persists. In this case, β i and β y are 0.3 and 0.2, respectively. In Figure 7A red lines represent dynamics with vaccination. Blue lines represent a counterfactual scenario without vaccination. We observe that the vaccine reduces the levels of both diseases, with the extinction of influenza ( Figure 7A). Figure 8 illustrates the case for coexistence of both viruses. In this case β i and β y are 0.55 and 0.2,

13
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted January 4, 2021. ;

Figure 7: Extinction of Influenza and persistence of SARS-CoV-2 with simultaneous vaccination. In the inset dots represent SARS-CoV-2 reproduction numbers, and crosses influenza reproduction numbers. Color codes: Red vaccination is applied, blue, no vaccination. A) Proportion of individuals infected with influenza. B) Proportion of individuals infected with SARS-CoV-2. Solid red lines represent disease dynamics when vaccination for both viruses is simultaneous. Dotted blue lines show dynamics without vaccination.
to make the reproductive numbers greater than one. Endemic equilibrium levels in the presence of vaccination are reduced compared to the no vaccine case (dotted blue lines). 14 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

A special case
In this section we explore a marginal case that our model produces. We underline that this is a very specific case that appears for a limited range of parameter values and initial conditions, yet it seems worth exploring a bit further. To be precise, when R 0k < R vk with k = i, y, parameters are outside the range of interest for SARS-CoV-2 and influenza but the dynamics that appears is curious. Taking µ = 0.000039139, η i = 1/5, η y = 1/14, θ = 1/365, θ i = 1/365, θ y = 1/180, β i = 0.19, β y = 0.071, p i = 0, α = 0.5 and ω y = 1/35. EVP against SARS-CoV-2 and influenza are 1.4% and 1%, respectively; ω i varies. In this section we do not refer to any specific viral diseases but postulate two arbitrary viral strains A and B. Figure 6B now shows only two possibilities: either extinction or coexistence. If both vaccine reproduction numbers are less than one, then both diseases become extinct, whereas if any of the vaccine reproductive numbers is greater than one, then coexistence of both viruses occurs. Figure 9 shows the coexistence of both viruses when ω i = 1/22. The inset in Figure 9A shows the reproduction numbers for both viruses. In the absence of vaccination both diseases become extinct (dotted blue lines). A similar pattern is observed when R vy < 1 < R vi . When both vaccine reproduction numbers are greater than one and ω i = 1/25 we obtain the case shown in Figure 10. In contrast to the case shown in Figure 9, vaccination reduces both disease levels. However, without vaccination, both basic reproductive numbers R 0k are less than one but coexistence is the outcome. This strongly indicates the existence of bi-stability and the appearance of a backward bifurcation [42]. Figure 11 confirms the existence of bi-stability. In general, as α increases, we observe a slight decrease in the number of superinfected (virus A) hosts and an increase in the superinfector (virus B). This is consistent with 15 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted January 4, 2021. ;  standard results, [9,13,14]. However, there are conditions where an increase of α implies a unexpected change in the disease dynamics. Figure 12 shows the effect of increasing the superinfection index (α) in a very particular regions of initial conditions. We set, as example, (s(0), i(0), y(0), r i (0), r y (0), i ry (0), y ri (0), r(0)) = (0.84814, 0.0001, 0.00176, 0.1, 0.05, 0, 0, 0). Other parameters are as in Figure 11 implying R 0k < 1, with k = A, B. We can see that an increase in α changes the stability of the equilibrium points, causing that 16 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted January 4, 2021. ; for small α, coexistence occurs. For larger α the stability switches to the extinction of both diseases. We point out that this behavior occurs within the region of a backward bifurcation for which the parameter ranges are unrealistic for SARS-CoV-2 and influenza. However, it is an interesting theoretical result. In this very particular situation, superinfection is not a promoter of coexistence but of extinction. This counterintuitive results will be explored in a forthcoming work.

Numerical continuation in a neighborhood of R 0k = 1
To further verify the existence of the backward bifurcation found in the previous section, we performed a numerical continuation [43] of the equilibrium points in the parameter regions of interest. In Figure 13 we show the case for R 0k slightly above 1 and in Figure 14 we show the case for R 0k slightly below one. Both numerical results were implemented in Matcom (Matlab tool) and confirm the existence of a backward bifurcation.

Seasonal contact rates
Finally, we address the issue of the long-term dynamics of the interactions between the two viruses. Seasonal variability is important to explain intra-annual fluctuations of viral populations [16,17]. We incorporate seasonality using a periodic effective contact rate β k (t) = β k (1 + cos ωt) where ω = 2π 365 corresponding to an annual period, β k is the baseline constant effective contact rate for virus k and is the amplitude of the seasonal variation (strength of the periodic forcing, 0 < < 1).

17
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted January 4, 2021. The disease-free equilibrium point curve s * is locally stable (green horizontal line) until β y = 0.071667 corresponding to R 0y = 1.0033. For higher SARS-CoV-2 infection rates, the disease free equilibrium becomes unstable. The stability interval of the coexistence curve denoted by x * iy , as a function of β y is y [0.0700, 0.2797] (green downward line) or 0.98 < R 0y < 3.91. This shows that without vaccination there are stable coexistence points despite the fact that R 0y < 1.
The higher reproduction number, but also the lack of previous immunity and absence antivirals/vaccines against SARS-CoV-2 makes it the competitively dominant virus in this interaction, like RSV over influenza. Figure 15A shows the time trajectory of both viruses in the absence of vaccination. Outbreaks alternate in time with influenza oscillations presenting narrower amplitudes. In contrast, SARS-CoV-2 outbreaks are broader. Both present inter-epidemic periods with very low prevalence. Epidemic outbreaks do not overlap. For secondary infections (those of recovered individuals from one disease but newly infected with the other), the alternation of patterns persists but prevalence of both viruses, particularly for influenza, decreases significantly (see Figure 15B-C). For primary and secondary infections, we show the biologically feasible case where SARS-CoV-2 is competitively superior to influenza, a decrease in α, meaning a decrease in the superinfection contact rate, provides a "more irregular" pattern among viruses, both in the periodicity and the ampli-18 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted January 4, 2021. ; https://doi.org/10.1101/2020.12.29.20248953 doi: medRxiv preprint iy and x * y , we find stability even though R 0y < 1. We can conclude that without vaccination the basic reproductive numbers associated with each virus need not be greater than one for the diseases to persist. tude of the epidemic outbreaks.
Every year influenza arrives in waves with each new influenza epidemic produced, in general, by a different viral strain in a process of lineage or variant replacement [44]. To mimic this situation in a simple, yet reasonable way, we aggregate all the influenza strains as a single influenza epidemic but allow for reinfections given the assumption that temporal immunity against influenza lasts one year. With this simplification, our model captures two main cyclic epidemic patterns: an annual pattern driven by yearly weather variability which is exploited by the dominant virus and a second biennially one, exploited by influenza. From alternating outbreaks, as the susceptibility to superinfection increases, various scenarios arise where influenza outbreaks are biennial and SARS-CoV-2 are annual,

19
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted January 4, 2021. ; once the superinfection continues to increase, both are synchronized being SARS-CoV-2 the dominant virus. These two cyclic behaviors are in the root of the superinfection mechanism for coexistence of both populations [16]. For influenza, it has been documented that incidence periodicity is generally biennial with different strains dominating alternate years [45]. This alternating pattern is caused by the immunity generated by the strains of influenza that circulate the previous year, by the grade of match and effectiveness of the influenza vaccine, but most important by the circulation of RSV subtype A or B. A predominantly higher circulation of RSV-A plus a predominance of influenza A/H3N2 will generates a flatter or lower impact influenza epidemic.
We introduce vaccination, (see eq. (2)), choosing φ k such that the times T to reach the desired EVP coverage q k , are equal to 4 and 3 months for influenza and SARS-CoV-2 respectively. In Figure 16A, we show the trajectories SARS-CoV-2 and influenza. Vaccination induces a synchronized oscillatory pattern in both viruses now strongly associated to the strength of the periodic forcing. This periodic behavior is also reflected in the secondary infections (see Figure 16B-C). Figure 16D shows that influenza rebounds every two years while SARS-CoV-2 does it every year. The effect of superinfection on the behavior of patterns can be observed by reducing the susceptibility to superinfection under a vaccination scheme. Note also that SARS-CoV-2 peaks are weaker every two years when they coincide with the influenza outbreaks ( Figure 16F) suggesting competition for hosts. Figure 16E shows that biennial behavior for influenza in both primary and secondary infections can be preserved under vaccination by increasing the strength of the periodic forcing and by 20 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted January 4, 2021. ; A plausible explanation of the Figure 16D-F behavior, is that, a greater seasonal forcing implies a greater amplitude in the outbreaks, which can be observed in the prevalence values for each type of virus. Reducing the effect of superinfection implies that more 21 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted January 4, 2021. ; influenza-infected individuals are not superinfected with SARS-CoV-2, which has an effect on reducing the prevalence of individuals recovered from SARS-CoV-2 but infected with influenza.

Sensitivity analysis
For completeness, a variance based sensitivity analysis known as Sobol method [46,47] was conducted to understand the parameters' influence on the model's state variables. This is a global sensitivity analysis that starts by assuming that the output of a system is a function of a set of inputs (parameters). By assuming that the vector of parameters is a random variable, the output is also a random variable. The total variability of the output, induced by the variability of the inputs, is decomposed in proportions associated with individual or sets of parameters. The higher the proportion of variability caused by changes in a specific parameter, then the higher the sensibility of the model to that parameter. θ, θ i , θ y , φ i , φ y , p i , α, β i , β y vary uniformly in the ranges presented in Table 3, while parameters η i = 1/5, η y = 1/14, ω i = 1/7 and ω y = 1/16 are held constant. The results are presented in Figure 17 and Figure 18.  We also perform a sensitivity analysis for the vaccine reproduction numbers. The sensitivity analysis for R vi and R vy is presented in Figure 19. It is clear that, for each disease, the parameter that changes the most the reproduction number is the corresponding contact rate. The only other parameter that shows an important effect is the vaccination coverage φ k which, of course, depends on the corresponding EVP q k .

22
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted January 4, 2021. ; Figure 17: Sensitivity analysis. Left panels show the first order Sobol indices for the individual effects on influenza prevalence i (upper left panel) and SARS-CoV-2 prevalence y (lower left panel) at each point in time. Right panels show total Sobol indices, for the effect of each parameter and its interactions with all other parameters. For influenza, β i is the most important parameter, with its influence slightly decreasing over time. β y is the second most influential parameter at the beginning of the epidemic, but it is eventually replaced by φ i . This implies that changes in vaccination schemes are indeed important to control influenza. In the case of SARS-CoV-2 , it can be seen that the most important parameter is β y , followed by θ y and φ y indicating that, over time, changes in the the vaccination schemes slowly gain importance.

23
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

24
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted January 4, 2021. ; https://doi.org/10.1101/2020. 12.29.20248953 doi: medRxiv preprint Our results show that the vaccine reproduction number for some parameter combinations can be higher than the basic reproduction number opening the possibility of the undesired outcome where vaccination may be worst than otherwise at the population level. This outcome underlies the importance of an adequate vaccination strategy when two viruses co-circulate. Figure 3 indicates that a low target EVP may push the vaccine reproduction number above the basic reproduction number, resulting in a higher prevalence when vaccinated compared to the case without vaccination. This low EVP scenario is highly unrealistic in most context but could be of consequence in critical settings such as war, social disturbance and natural disasters where coverage may fall short of the desired target. On the other hand, influenza and SARS-CoV-2 are co-circulating simultaneously in many countries. However, vaccination campaigns do not begin at the same time. As example, in Mexico the influenza vaccination campaign started on October 2020, but that for SARS-CoV-2 hast started on late December 2020. It is important, therefore, to carefully weight the consequences of not starting the vaccination campaigns and their target EVP since this decision can affect the final endemic levels of both diseases. Our model shows, as expected, that asymptotic behavior is closely related to the vaccine reproduction number for each type of virus. For example, when considering realistic parameters for influenza and SARS-CoV-2, Figure 6 shows that if the vaccine reproduction number for influenza is greater than one and the vaccine reproduction number for SARS-CoV-2 is less than one, then influenza persists, and SARS-CoV-2 is eradicated. Coexistence sets in when both vaccine reproduction numbers are greater than one. When the vaccination rate is low we have shown, albeit in an special case, that the basic reproductive numbers do not need to be greater than one for both diseases to coexist ( Figure 11). For this same situation, when R 0k < 1 k = A, B, there are initial conditions where the superinfection index can switch the stability of equilibria: low α gives coexistence and higher α, extinction (see Figure 12).
We have also numerically explored the behavior of our model with time-dependent effective contact rates. In general we observe that influenza epidemics have less amplitude and show inter-epidemic periods with very low prevalence whereas SARS-CoV-2 epidemics are broader in amplitude and show a clear endemic phase between outbreaks. In all the scenarios that we simulated, we have observed that the prevalence of secondary cases (hosts that are susceptible to one but have recovered from the other) of both viruses, decreases. The simulations assumed an effective contact rate for influenza higher than that of SARS-CoV-2, however the force of infection of this last virus (the one with higher reproduction number) is much greater, since the infection rate depends of the number of contacts per unit time, but also on the infection probability per contact and the infectious period which are different for both viruses.
Our model has not specifically explored the role that non-pharmaceutical interventions will continue to have while vaccine coverage is still insufficient to grant an effective control of the epidemic. We presently are studying this factor.

25
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted January 4, 2021. ; https://doi.org/10.1101/2020.12.29.20248953 doi: medRxiv preprint