Co-circulation of two viral populations under vaccination

The interaction and possibly interference between viruses infecting a host population is addressed in this work. We model two viral diseases with a similar transmission mechanism and for which a vaccine exists. The vaccine is characterized by its coverage, induced temporary immunity, and efficacy. The population dynamics of both diseases consider infected individuals of each illness and hosts susceptible to one but recovered from the other. We do not incorporate co-infection. Two main transmission factors affecting the effective contact rates are postulated: i) the virus with a higher reproduction number can superinfect the one with a lower reproduction number, and ii) there exists some induced (indirect) 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, even considering the absence of superinfection. Influenza and SARS-CoV-2 are employed to exemplify this last point, observing that the time-dependent effective contact rate may induce either alternating outbreaks of each disease or synchronous outbreaks. Finally, for a particular parameter range, a backward bifurcation has been observed for dynamics without vaccination.


Introduction
The interaction between viral species follows known patterns of coexistence. [1] were the first to show that the inclusion of superinfection makes coexistence between competing species possible. Superinfection is a process where a competitive hierarchy exists among a group of species that compete for the same resources. In general (in the absence of vaccines or treatments), this 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 [2] that reflect the abilities to use a contested common resource.
The process of superinfection has been identified as a promoter of the coexistence of pathogen strains in a given host population. The structure of the models that have been used to theoretically demonstrate this property [1,3,4], has also been applied to explain the organization of community structure for species that inhabit a common landscape but are not necessarily closely related [5]. For infectious diseases, the concept of superinfection has pervaded theoretical explanations for coexistence, in the same host population, of variants of a given pathogen [6,7]. The pioneering model in [3] is about infectious diseases and has been used to explore a plausible hypothesis for the length of the latent period of HIV before the onset of AIDS. Simultaneously, these same authors [8] addressed the problem of community structure postulating a trade-off between colonization and extinction [9] in the presence of a hierarchy of competitive abilities for the exploitation of the resources in a common landscape. 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 syncytial 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 [6]. Seasonal influenza has a median basic reproduction number of 1.28 [10]; while for SARS-CoV-2, the median R 0 is 2.79 [11,10]. This difference in transmission potential supports the assumption that while competing for hosts, their common resource, the likelihood of long-term coexistence and co-circulation is high and that this balance may be associated with superinfection.
Vaccines have a protective effect via immune responses or indirect effects by reducing the burden of viral and bacterial respiratory diseases on individual patients, among others [12]. However, in general, the impact of vaccination on the transmission of a viral disease starts slowly and builds up over several months to reach target coverage levels. Vaccines are evaluated in terms of efficacy, existing population immunity, coverage and temporary immunity, both natural and vaccine-induced, reduction of mortality or infection risks [13]. Vaccines are applied in order to eradicate a disease. However, this aim depends on many factors and may not be fully achieved.
The present work addresses the general theoretical problem of the co-circulation and long-term persistence or eventual extinction of two viral respiratory infections subject to vaccination. In addition, we explore: i) the order relation between basic reproductive number and vaccine reproduction number, and ii) the existence of backward bifurcation when vaccination is not considered. Finally, we exemplify some of our results with the case of the developing ecological interaction between SARS-CoV-2 and influenza.
The paper is organized as follows. In Section 2 we formulate our mathematical model. In Section 3 we develop the local analysis, including cases of backward bifurcation. In Section 4, we address a particular case, co-infection dynamics between influenza and SARS-Cov-2. Finally, in Section 5 we draw some conclusions about this work.

Mathematical model set-up
We formulate a mathematical model considering the simultaneous presence of two viruses and vaccination for each one. We do not distinguish between viral subtypes and, therefore, we approximate viral dynamics as if each of them were a single viral population. Both diseases are assumed to present temporary immunity and, thus, the possibility of reinfections. The model is a coupled system of two SIRS (susceptible-infected-recovered-susceptible) equations ( Figure 1). Figure 1: 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 virus i, infected by virus y, immune from virus i, immune from virus y, immune from virus i but infected with virus y, immune from virus y but infected with virus i, and immune from both viruses. Here, r k includes population recovered after infection and successfully vaccinated by virus k, where k represents i or y viruses. Dashed blue lines represent vaccination dynamics for both viruses and dashed red line denotes superinfection process.
We assume a constant total population, normalized to the total population N . In the presence of vaccination against both viruses, the equations stand as: In eq. (1), s, i, y, r i , r y , y ri , i ry and r are defined as described in Figure 1. Note that the 3 . 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) immune population for each virus can be infected (vaccinated) by (against) the other virus. With the aim of approximate the seasonal component that drives ILIs [14,15,16], for numerical analysis, we assume that the effective contact rates for both viruses, β i (t) and β y (t), are time-dependent. It is also assumed that virus y may potentially competitively exclude virus i from its host. The effective contact rate of this transmission route is modulated by α, with 0 ≤ α ≤ 1 that we define as the superinfection coefficient.
On the other hand, we define p i , 0 ≤ p i ≤ 1, to measure the indirect protective effects of ivaccinated individual against infection by virus y (e.g., as described for influenza and SARS-CoV-2 in [12], additionally [17], found a significant reduction in the odds of testing positive for COVID-19 in patients who received an influenza vaccine compared to those who did not receive the vaccine. In pediatric populations, seasonal influenza and pneumococcal vaccination may have protective value in symptomatic COVID-19 diseases [18]). This protective effect is modelled by a reduction of the effective contact rate of virus y when in contact with an immune individual from virus i (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 related to virus k. To simplify the model, we assume that θ, the loss of immunity rate once an individual becomes immune from both viruses, is a constant and does not depend on the order of the infections. Other parameter descriptions are shown in Table 1  We do not distinguish between vaccinated and recovered individuals. Instead, we collect immunized individuals into a single compartment; r i and r y contain, therefore, those individuals that have either been vaccinated against virus i or virus y, respectively or that, alternatively, are recovered from a natural infection for either of these two viruses. This modeling choice reduces the system's dimensionality.
Effective target coverage definition We define the vaccination rate, φ k where k denotes either virus i or virus y as an effective vaccination rate. In other words, φ k incorporates vaccine efficacy. Under 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 E k of the susceptible population is vaccinated at time T , then we set a vaccination rate such that q E k is the effective target coverage or ETC with time horizon T k : if the vaccine has σ k % efficacy and we apply the vaccine to q k % of the population, then only a fraction q E k = σ k q k is effectively protected where k = i denotes virus i and k = y denotes virus y.

Local 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 [19].

Reproduction number
The disease-free equilibrium of eq. (1) always exists and it is given by E 0 = (s * , i * , y * , r * i , r * y , y * ri , i * ry , r * ) where E0 = s * , 0, 0, φis * µ + θi + φy , φys * µ + θy + φi , 0, 0, s * θ + µ φiφy µ + θi + φy with where H = (µ + θ) (µ + θ i + φ y ) (µ + θ y + φ i ). Note that the susceptible population at the diseasefree equilibrium depends on the parameters for coverage and immunity for both viruses . Given the interaction of both viruses, their reproduction numbers give information on conditions for coexistence, competitive exclusion, or extinction. In what follows, 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 compute the vaccine reproduction number. We proceed as in [20] to obtain: Here, s * , r * i and r * y are defined in eq. (3). Then, the vaccine reproduction number is given by the spectral radius of matrix F V −1 : 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 May 26, 2021. ; https://doi.org/10.1101/2020.12.29.20248953 doi: medRxiv preprint where the vaccine reproduction numbers for virus i (R vi ) and virus y (R vy ) are: Remark: In the absence of vaccination for either one of the viruses, and for constant effective contact rates, the basic reproductive numbers are the classical expressions of an SIR epidemic, that is From eqs. (5) and (6), the order relation kept by R vk and R 0k , with k = i, y, is not obvious. For example, Figure 2 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, ω y = 1/16, p i = 0.05. Varying both coverages q E k from 1% to 99% (eq. (2)), it is possible that R vk > R 0k with k = i or y, i.e, vaccine application may enhance rather than reduce the occurrence of an outbreak. It is therefore, important to find the conditions that guarantee that only a reduction occurs, i.e., that R vk < R 0k . In Figure 2, the time horizon to reach vaccination coverages against viruses i and y is fixed on four and three months, respectively. We further characterize the relationship between these reproductive numbers.

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 May 26, 2021. ; Lemma 3.2 Let R vk and R 0k be defined as in eqs. (5) and (6), respectively; and M i = s * + (ηi+µ) (ωi+µ) r * y and M y = s * + M i and M y can be interpreted as the effective proportion of susceptible individuals available to infection by virus i and y, respectively. This lemma underlines the dependence of the order relation between the reproductive numbers R vk and R 0k , with k = i, y, on the pool of susceptible individuals for each virus.
The proof is given in Appendix A.

When R 0k < R vk
In lemma 3.2, we have shown that there is a range of parameter combinations such that R 0k < R vk , with k = i, y, that is, where vaccination is not effective in reducing transmission. In this section we explore this case. We fix µ = 0.000039139, α = 0.5, 001603, φ y = 0.002341 and θ y = 1/180. Thus, M i = 1.131837 and M y = 1.076793. Figure 3 shows the asymptotic equilibria, when t = 36500, of i and y as function of R vi and R vy . Both effective contact rates values are constant and defined such that 0.4 ≤ R vi ≤ 1.6 and 0.5 ≤ R vy ≤ 2, with the initial condition (s(0), i(0), y(0), r i (0), r y (0), i ry (0), y ri (0), r(0)) = (0.8498, 0.0001, 0.0001, 0.1, 0.05, 0, 0, 0). Figure 3 illustrates that when R vk > 1, with k = i, y, both disease coexist. Moreover, there exist combinations of parameter values such that R vi < 1 < R vy which also implies coexistence ( Figure 3A). This phenomenon can be interpreted as a kind of rescue effect of one virus by the other. Figure 3B also shows this phenomenon. Finally, when R vk < 1, with k = i, y, both disease go extinct. Figure 4 shows a viral coexistence. Here, β i and β y are 0.19 and 0.07, respectively. This case shows that even when both basic reproductive numbers R 0k are less than one, both viruses persist in the absence of vaccination. This pattern strongly indicates the existence of bi-stability and of a backward bifurcation [21]. Figure 5 confirms the existence of bi-stability. Here, we fix initial conditions (s(0), i(0), y(0), r i (0), r y (0), i ry (0), y ri (0), r(0)) = (0.8499 − y 0 , 0.0001, y 0 , 0.1, 0.05, 0, 0, 0). When A backward bifurcation means that the reproduction number being less than unity becomes only a necessary, but not sufficient condition, for disease elimination. To further explore the existence of the backward bifurcation shown in Figures 4 and 5, we performed the numerical continuation of the equilibrium points in appropriate parameter regions.

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 May 26, 2021.   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 May 26, 2021. ; Figure 5: Bi-stability when R 0k < 1, with k = i, y. Red line represents coexistence of both viruses. Blue line shows extinction. The inset is a zoom of the extinction dynamics.
of β k so that R 0k is below and above unity. All numerical results were implemented in Matcom (Matlab tool) [22]. Figure 6 illustrates the numerical continuation of equilibrium points for viral coexistence and extinction of both viruses (in the absence of vaccination), projected on the plane β k and virus k for k = i, y. Parameter values are β i = 0.19, α = 0.5, θ i = 1/365, θ y = 1/180, θ = 1/365, η i = 1/5, η y = 1/14, ω i = 1/21, ω y = 1/56, p i = 0. Figure 6A shows that the continuation extinction curve of both viruses given by s is locally stable (solid line) up to the value β y = 0.0714676, corresponding to a reproduction number of R 0y = 1.0054. In curve s as the virus y infection rate increases, the disease-free equilibrium curve becomes unstable (dashed line). The coexistence curve of both viruses denoted by iy, is locally stable (solid line) for β y ∈ [0.05814, 0.2697] or, equivalently for 0.81 < R 0y < 3.77.
In Figure 6B, we show the numerical continuation of equilibrium points: for viral coexistence and extinction of both virus projected on the plane i − β y . The disease-free equilibrium branch s, is locally stable (solid line) until R 0y = 1.0054. The stability interval of the branch with the two virus present, iy (solid line) is 0.81 < R 0y < 3.79. Figs. 6C-D illustrate the numerical continuation of equilibrium points (coexistence and extinction of both virus) projected on the planes y − β i and i − β i , respectively. In both, the disease-free equilibrium branch s is locally stable (solid line) until R 0i = 1.09, the stability of the coexistence curve starts at R 0i = 0.502.
In all cases, we show that a stable coexistence equilibrium exists together with a disease free equilibrium when R 0k < 1. The usual causes of backward bifurcation in some standard deterministic models are imperfect vaccination [23,24], the existence of exogenous re-infections [25] or vaccinederived immunity waning at a slower rate than natural immunity [26], the role of re-infection [27], among others [28]. Since the recovery rate from the second infection for both viruses is smaller 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 May 26, 2021. ; than the recovery rate from the first infection, the influx to the pool of susceptibles occurs in a time window larger than expected if only one infection were present. This continuous and extended influx generates a backward bifurcation. Therefore the vaccination strategy must be efficient (large coverage in as short a time as possible), to prevent the increase in the pool of susceptibles that may lead to undesirable outcomes. In general, as superinfection increases in strength (α increases), we observe a corresponding slight decrease in the number of superinfected ( with virus i) hosts and an increase in the superinfector (with virus y). This is consistent with standard results, [1,3,4]. However, there exist parameter values where an increase of α produces an unexpected change in disease dynamics. Figure 7 shows the effect of increasing superinfection (α) for a particular initial conditions and parameter values. As an example we take (s(0), i(0), y(0), r i (0), r y (0), i ry (0), y ri (0), r(0)) = (0.8496, 0.0001, 0.0003, 0.1, 0.05, 0, 0, 0). Other parameters are as in Figure 5 to insure that R 0k < 1, with k = i, y. Apparently, there is a threshold value for α that switches the dynamics from coexistence to extinction in both viruses.

A particular case: Influenza and SARS-CoV-2
The year 2020 and early 2021 have 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 been almost null in the southern hemisphere and now in the northern, while the SARS-CoV-2 pandemic is active [29]. Figure 8 exemplifies this in the case of Mexico. The winter months (Northern hemisphere) did not produce influenza outbreaks concurrent with COVID-19 resurgence and, thus, hospital capacity in Mexico and elsewhere was not compromised [30,31]. A possible explanation for the absence of influenza is that the measures are taken to prevent COVID-19 (social distancing, mask use, etc.) also prevent influenza transmission. These measures have limited ability to stop the COVID-19 (aerosols) but may be quite effective at preventing influenza transmission. Nevertheless, the occurrence of a syndemic episode with co-circulating influenza and SARS-CoV-2 viruses is still a potential reality.
In this section, Eq. (1) is used to explore co-circulation dynamics between SARS-CoV-2 and influenza. Given that, to date, there exist no evidence to suggest that a COVID-19 infection has been observed to displace an influenza infection, we consider α = 0. Baseline parameters are given in Appendix B.1.

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 May 26, 2021. ;

Reproduction numbers, vaccine efficacy, effective coverage and temporary immunity
The dependence of the vaccine reproduction number (eq. (4)) on vaccination coverage and temporary immunity is explored now. For the rest of this section, we fix µ = 0.000039139, η i = 1/5, η y = 1/14, θ = 1/365, θ i = 1/365, ω i = 1/5 and ω y = 1/14. The ETC for vaccine k (q E k ) determines the corresponding φ k through eq. (2), with k = i, y. We fix the effective transmission rates to β i = 0.32 and β y = 0.15. Figure 9 shows R V as function of ETC for SARS-CoV-2 (q Ey ) achieved in three months (see eq. (2)) and the average duration of immunity by SARS-CoV-2 (θ −1 y ). The protective effect of the influenza vaccine (p i ) are also left free to vary. As expected, the maximum value of R v is achieved when there is not vaccination, that is, q Ey = 0. Likewise, we also observe that as increases p i the maximum value of R v is lower.
The effect of the vaccine efficacy and time horizon over vaccine reproduction numbers are presented in Appendix B.2.

Transmission reduction by vaccination (R vk < R 0k )
Numerical explorations of eq. (1) allows us to postulate the diagram in Figure 10. We observe that 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. Finally, if both vaccine reproductive numbers are greater than one, coexistence of both diseases ensues. In Appendix B.3 we present numerical simulations in support of our results of this section.

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 May 26, 2021. ; Figure 9: Vaccine reproduction number as a function of coverage and average duration of immunity for SARS-CoV-2. A) No protection from influenza vaccination (p i = 0). B) 50% protection (p i = 0.5). C) Full protection (p i = 1). Target coverage (TC) for influenza is 35% in 4 months with a vaccine efficacy of 50%. Note that the greater p i is, the higher the reduction in R V .

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 [6,7]. We incorporate seasonality using a periodic effective contact rate β k (t) = β k (1 + cos ωt) where ω = 2π 365 for 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). Due to the higher reproduction number, the lack of previous immunity, the absence of antivirals or scarce vaccines, the equal and homogeneous effect of NPIs on reducing the effective contact rate, henceforth we consider SARS-CoV-2 as the competitively dominant virus in this interaction with influenza. This behavior is similar to the RSV over influenza. Figure 11 illustrates alternation patterns between influenza and SARS-CoV-2. Here, parameter values are β i = 0.45, β y = 0.22, α = 0, θ i = 1/365, θ y = 1/180, θ = 1/365, η i = 1/5, η y = 1/14, ω i = 1/5, ω y = 1/14, p i = 0.5, = 0.8, q i = 0.2, q y = 0.3. Figure 11A shows that influenza and SARS-CoV-2 have stronger outbreaks every two years in an alternating sequence. This behavior is related to a reduction in susceptibility and an increase in the strength of the periodic forcing under 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 May 26, 2021. ; https://doi.org/10.1101/2020.12.29.20248953 doi: medRxiv preprint Figure 10: Summary of the asymptotic dynamics when varying both vaccine reproduction numbers. Scenarios for SARS-CoV-2 and influenza when R vk < R 0k , with k = i, y. a vaccination scheme. Note also that SARS-CoV-2 peaks are weaker every two years when they coincide with stronger influenza outbreaks and vice versa, suggesting competition for hosts.
Figs. 11B-C show that this alternating sequence can be preserved even under vaccination for both viruses. Thus, in this scenario, every two years, the amplitude of the primary infections of both viruses is alternatively greater and, also, that this behavior may indeed correspond to the competitive alternation of both viruses.
Other patterns, between both viruses, are shown in Section B.4.

Conclusions
Influenza and SARS-CoV-2 will likely be co-circulating in the near future in many countries. However, vaccination campaigns have not begun at the same time. For example, in the Northern hemisphere the influenza vaccination campaign started in the Fall 2020, but that for SARS-CoV-2 has begun in early 2021 in many part of the World. Therefore, it is important to carefully plan vaccination campaigns and to define a realistic and sufficient coverage to avoid the situation described in the first paragraphs of this section.
Our results show that the parameters related to SARS-CoV-2, such as the average time of loss immunity, effective target coverage, and protective effect against infection by the coronavirus, all are relevant in reducing the vaccine reproduction number (Figure 9). To date, there are some 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.

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

Figure 11: Alternation patterns. A) Annual prevalence cycles for primary influenza infections (blue line) and primary SARS-CoV-2 infections (dashed blue line). B) Annual cycles for primary (orange dashed line) and secondary (orange line) influenza infections. C) Annual cycles for primary (black dashed line) and secondary (black line) SARS-CoV-2 infections.
parameter estimates for SARS-CoV-2 that remain unknown like temporary immunity and the large uncertainity of vaccine availability for developing or poor countries.
Our model also shows, as expected, that asymptotic behavior is closely associated with the vaccine reproduction number for each type of virus. For example, when considering realistic parameters for influenza and SARS-CoV-2, Figure 10 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 reproductive numbers are above one. This may happen as a consequence, for example, of low coverage or the time horizon for achieving it is large.
We have also numerically explored the behavior of our model with time-dependent effective contact rates. We consider it relevant because COVID-19 is a new disease with a transmission route similar to other viral infections, and in consequence, seasonal variability can explain future behaviors when considering co-circulation dynamics.
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 the simulated scenarios, we have observed that the prevalence of secondary cases (hosts that are susceptible to one but have recovered from the other virus) of both viruses decreases. The simulations assume 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 the higher reproduction number) is much greater since the infection rate depends on the number of 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 May 26, 2021. ; contacts per unit time, but also on the infection probability per contact and the infectious period which are different for both viruses. Also, we have found that the transmission patterns lead to alternation of patterns where influenza and SARS-CoV-2 have stronger outbreaks every two years in an alternating sequence, suggesting competition for hosts (see Figure 11).
In a more theoretical and general case, our results show that the vaccine reproduction number for some parameter combinations can be higher than the basic reproduction number opening up the possibility of the undesired outcome where vaccination may have a negative public health impact than otherwise at the population level. This outcome underlines the importance of the design of the vaccination strategy and of the availability of vaccines. Figure 2 shows that a low ETC may push the vaccine reproduction number above the basic reproduction number, resulting in a higher prevalence than the case without vaccination. This low ETC case is unrealistic in most contexts 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.
Finally, in the absence of vaccination, we have shown that there are conditions under which the basic reproductive numbers do not need to be greater than one for both diseases to coexist ( Figure 4). The above confirms the existence of bistability and of a backward bifurcation where for this special case, the recovery rate from the second infection for both viruses is slower than the recovery rate of the first infection. For this same situation, when R 0k < 1 k = i, y, there are initial conditions where superinfection can switch the stability of equilibria: low α gives coexistence and higher α, extinction (see Figure 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 May 26, 2021. ; https://doi.org/10.1101/2020.12.29.20248953 doi: medRxiv preprint Then, the characteristic polinomial of the J E0 is We observe that a 1 , a 2 , a 3 , a 4 > 0 and a 1 a 2 a 3 − a 2 3 − a 2 1 a 4 > 0. Hence, p 1 (λ) satisfying the RouthHurwitz criterion. In consequence, all roots of p 1 (λ) have negative real parts. Likewise, it is clear that all roots of p k (λ) are negative if and only if R vk < 1, with k = i, y. Therefore, all eigenvalues of J E0 have negative real parts if and only if R v < 1.

Appendix B Influenza and SARS-CoV-2 B.1 Model parametrization
To explore numerical scenarios from Influenza and SARS-CoV-2, we estimate the baseline parameters from bibliographical sources. According to [32], for seasonal influenza the median reproduction 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 May 26, 2021. ; number is 1.28 and that for SARS-CoV-2 is 2.79, but sources report large variability [33]. Likewise, [10] 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, [34] reports an incubation period of 3 to 4 days and an infectious period (η y ) of 4-5 days but [10] 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.
Influenza vaccine efficacy varies every year. For 2019-2020 the [35] reports σ i = 0.29 but in 2010-2011 σ i = 0.6. [36,37] 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). Besides this reported efficacies we are postulation 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 [38] that reports that protective influenza vaccination does not negatively affect the risk of contracting coronaviruses. We explore the possibility that the effect is positive, thus conferring a reduction in the risk of SARS-CoV-2 infection.
For SARS-CoV-2, several vaccines have been deployed with efficacies in the range of 50-95% with more likely scenarios of 70%. The Pfizer, Moderna vaccines have efficacies at the upper end of this interval. Astra-Zeneca vaccine efficacy is around 75%. 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 set up scenarios where the above percentages are reached after T = 90, 180 and 365 days. We assume that SARS-CoV-2 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-infection to SARS-CoV-2 is not known.
B.2 Reproduction numbers, vaccine efficacy, effective coverage and temporary immunity

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 May 26, 2021. ;

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 May 26, 2021. ;  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 [40]. To mimic this situation in a simple, yet reasonable way, we aggregate all the influenza strains as a single influenza epidemic and allow for reinfections given the assumption that temporary immunity against influenza lasts one year. With this simplification, our model produces an annual pattern driven by the yearly weather variability.

B.5 Sensitivity analysis
For completeness, a variance based sensitivity analysis known as Sobol method [41,42] was conducted to evaluate the parameters' influence on the model's state variables. A global sensitivity analysis assumes that the output of a system is a function of a set of inputs (parameters). By 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.

(which was not certified by peer review)
The copyright holder for this preprint this version posted May 26, 2021.  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 B.5.1, while parameters η i = 1/5, η y = 1/14, ω i = 1/5 and ω y = 1/14 are held constant.   most important parameter, with its individual influence decreasing over time. φ i is the second most influential parameter followed by β y . This implies that changes in vaccination schemes are indeed important to control influenza. In the case of the SARS-CoV-2 (y) virus, the most dominant parameter is β y while the influence of the others is very small. Figure B.5.2: Sobol indices for secondary infections. It can be seen that the influence of the parameters on primary and secondary infections of influenza are very similar.In the case of SARS-CoV-2 secondary infections, the protective effect gained from influenza (p i ) is only below β y in terms of importance.

Parameter Lower bound Upper bound
We also perform a sensitivity analysis for the vaccine reproduction numbers R vi and R vy , Figure B.5.3. For each disease, the parameter that changes the most the reproduction number is the corresponding contact rate. The other parameters that show an important effect are the vaccination coverage φ i and φ y which, of course, depend on the corresponding ETC q Ei and q Ey .

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 May 26, 2021. ; https://doi.org/10.1101/2020.12.29.20248953 doi: medRxiv preprint

26
. 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 May 26, 2021. ; https://doi.org/10.1101/2020.12.29.20248953 doi: medRxiv preprint

27
. 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 May 26, 2021. ; https://doi.org/10.1101/2020.12.29.20248953 doi: medRxiv preprint . 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 May 26, 2021. ; https://doi.org/10.1101/2020.12.29.20248953 doi: medRxiv preprint . 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.