Quasineutral multistability in an epidemiological-like model for defective-helper betacoronavirus infection in cell cultures

Applied Mathematical Modelling 137 (2025) 115673 Available online 3 September 2024 0307-904X/© 2024 The Authors. Published by Elsevier Inc. This is an open access article under the CC BY-NC license (http://creativecommons.org/licenses/by-nc/4.0/). Contents lists available at ScienceDirect Applied Mathematical Modelling journal homepage: www.elsevier.com/locate/apm Quasineutral multistability in an epidemiological-like model for defective-helper betacoronavirus infection in cell cultures Juan C. Muñoz-Sáncheza,b,1, J. Tomás Lázaro c,d,e,f,1, Julia Hillunga, María J. Olmo-Ucedaa, Josep Sardanyése,f, Santiago F. Elena a,f,g,∗ aInstitute for Integrative Systems Biology (I2SysBio), CSIC-Universitat de València, Paterna, 46980 València, Spain bDepartament de Física Teòrica, Universitat de València, Burjassot, 46100 València, Spain cDepartament de Matemàtiques, Universitat Politècnica de Catalunya (UPC), 08028 Barcelona, Spain dInstitute of Mathematics, UPC-BarcelonaTech (IMTech), 08028 Barcelona, Spain eCentre de Recerca Matemàtica (CRM), Cerdanyola del Vallès, 08193 Barcelona, Spain fDynamical Systems and Computational Virology, CSIC Associated Unit CRM-I2SysBio, Spain gSanta Fe Institute, Santa Fe, NM 87501, USA A B S T R A C T It is well known that, during replication, RNA viruses spontaneously generate defective viral genomes (DVGs). DVGs are unable to complete an infectious cycle autonomously and depend on coinfection with a wild-type helper virus (HV) for their replication and/or transmission. The study of the dynamics arising from a HV and its DVGs has been a longstanding question in virology. It has been shown that DVGs can modulate HV replication and, depending on the strength of interference, result in HV extinctions or self-sustained persistent fluctuations. Extensive experimental work has provided mechanistic explanations for DVG generation and compelling evidences of HV-DVGs virus coevolution. Some of these observations have been captured by mathematical models. Here, we develop and investigate an epidemiological-like mathematical model specifically designed to study the dynamics of betacoronavirus in cell culture experiments. The dynamics of the model is governed by several degenerate normally hyperbolic invariant manifolds given by quasineutral planes -i.e., filled by equilibrium points. Three different quasineutral planes have been identified depending on parameters and involving: (i) persistence of HV and DVGs; (ii) persistence of non-infected cells and DVG-infected cells; and (iii) persistence of DVG-infected cells and DVGs. Key parameters involved in these scenarios are the maximum burst size (𝐵), the fraction of DVGs produced during HV replication (𝛽), and the replication advantage of DVGs (𝛿). More precisely, in the case 0 <𝐵<1 +𝛽the system displays tristability, where all three scenarios are present. In the case 1 +𝛽<𝐵<1 +𝛽+𝛿this tristability persists but attracting scenario (ii) is reduced to a well-defined half-plane. For 𝐵>1 +𝛽+𝛿, the scenario (i) becomes globally attractor. Scenarios (ii) and (iii) are compatible with the so-called self-curing since the HV is removed from the population. Sensitivity analyses indicate that model dynamics largely depend on DVGs production rate (𝛽) and their replicative advantage (𝛿), and on both the infection rates and virus-induced cell deaths. Finally, the model has been fitted to single-passage experimental data using an artificial intelligence methodology based on genetic algorithms and key virological parameters have been estimated. Contents 1. Introduction.............................................................................. 2 2. Summary of the experimental results .............................................................. 4 3. Mathematical model......................................................................... 4 * Corresponding author at: Institute for Integrative Systems Biology (I2SysBio), CSIC-Universitat de València, Paterna, 46980 València, Spain. E-mail address: [email protected] (S.F. Elena). 1Equal contribution. https://doi.org/10.1016/j.apm.2024.115673 Received 13 February 2024; Received in revised form 24 July 2024; Accepted 29 August 2024 Applied Mathematical Modelling 137 (2025) 115673 2 J.C. Muñoz-Sánchez, J.T. Lázaro, J. Hillung et al. 4. Results ................................................................................. 5 4.1. Planes of equilibria, stability and basins of attraction............................................... 6 4.2. Cell culture dynamics: estimation of the amount of infective particles .................................... 10 4.3. Case 𝑚 >1: all cells are simultaneously infected.................................................. 14 4.4. Sensitivity analysis .................................................................... 16 4.5. Model fitting to HCoV-OC43 data: parameters’ estimation............................................ 17 5. Discussion ............................................................................... 19 CRediT authorship contribution statement................................................................ 20 Declaration of competing interest ..................................................................... 20 Data availability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 Acknowledgements .............................................................................. 20 Appendix A. Time-series examples for the Π𝑉𝐷 , Π𝐶𝐷𝐷and Π𝐶𝐶 𝐷 ................................................. 21 Appendix B. First order variational equation with respect to initial conditions........................................ 21 Appendix C. Additional fittings results................................................................. 22 Appendix D. Model solutions for some 𝜄∕𝛼values.......................................................... 23 References.................................................................................... 23 1. Introduction RNA viruses can quickly adapt and trigger epidemics by crossing species barriers due to their high mutation rates, fast replication, and large population sizes [1–3]. However, a high mutation rate is a double-edge sword, as many mutations during infection result in defective viral genomes (DVGs). DVGs cannot complete the infectious cycle by themselves thus depending on the viral proteins synthesized by a wild-type helper virus (HV). The generic term DVG includes point mutations, hypermutated genomes, deletions, insertions, and genomic reorganizations [4]. Huang and Baltimore coined the term defective interfering particles (DIPs) for viral particles containing DVGs and normal structural proteins encoded by the HV [5]. DIPs were first identified in the late 1940s by Von Magnus and Gard [6]based on the negative impact they exerted on virus accumulation. DIPs rely on a HV for replication, disrupting HV accumulation and impacting viral pathogenesis [7,8]. Several studies have shown that viruses rich in DIPs reduce virulence [9], induce high interferon levels [10], aid viral persistence [11,12], and modulate infection capacities as shown for different SARS-CoV-2 strains [13]. Recent high-throughput sequencing techniques have uncovered the emergence of a plethora of DVGs within a single infected host [14–17]. Notably, these studies have shown that distinct subsets of prevalent DVGs recurrently appear along time, suggesting intricate dynamics within the viral population. These dynamics encompass competition, and possibly compensation or cooperation, among various DVGs. Positive selection favours the most competitive DVG variants, indicating their relative fitness concerning the HV and other DVGs. Despite the generation of hundreds or even thousands of distinct DVGs during infections, the majority are lost due, among other factors, to population bottlenecks occurring during in vivo transmissions among individuals [18]or during in vitro diluted serial passages [19]. However, DVGs might persist for long periods of time in immunosuppressed hosts or in those with comorbidities [20], or if a high ratio between viral particles and susceptible cells (a parameter known as the multiplicity of infection, MOI) is experimentally imposed [21–23]. During infections of in vitro cell cultures, as those that have motivated this modelling work, DVGs accumulate when viral populations are repeatedly passed at a high MOI, while they do not accumulate if MOI is low or strong bottlenecks are imposed at each transmission event [21,22,24]. The higher the MOI, the more likely DVGs and HV would coinfect the same cell and thus persist in the population [21,22,24]. While the molecular recombination mechanisms by which DVGs are generated are well understood, their role in modulating the outcome of viral infections is sometimes unclear. Understanding the characteristics and functions of DVGs is thus crucial for comprehending the complexity of viral infections and for developing strategies to control or mitigate their impact. DIPs act as true hyperparasites interfering with the HV replication [23,25], competing for resources and reducing its accumulation and transmission efficiency [5,22,26–28]. Experimental in vitro studies with cell cultures have shown that DIPs may engage in an arms race with the HV [29–32]. Shorter genomes earn an advantage in terms of replication speed compared to the HV [33]. Additionally, there is evidence of a stronger form of interference in which the DIPs compete more effectively for the viral replication machinery [30,34]. Since their discovery, and given their effect on the accumulation of the HV, DIPs have attracted the attention of researchers as potential antiviral candidates [28,35,36], known as therapeutic interfering particles (TIPs). Following this idea, Xiao et al. created a TIP by deleting the capsid-coding region of poliovirus [37]. Remarkably, the administration of this TIP to mice triggered a broad antiviral response against diverse respiratory viruses, including enteroviruses, influenza A virus, and SARS-CoV-2. The broad-spectrum antiviral effects of this synthetic TIP was attributed to local and systemic type I interferon responses. Notably, a single dose not only safeguarded animals from SARS-CoV-2 infection but also stimulated the production of SARS-CoV-2 neutralizing antibodies, protecting against reinfection [37]. Another successful application of TIPs was reported by Chaturvedi et al. [38]. Following a synthetic biology approach, these authors generated artificial particles that had significant anti SARS-CoV-2 effect in cell cultures and primary lung organoids, reducing viral accumulation 10to 100-fold. Furthermore, intranasal application of these TIPs in infected hamsters suppressed the virus by 100-fold in the lungs, reduced pro-inflammatory cytokine expression, and prevented severe pulmonary edema [39]. Interestingly, the prediction of SARS-CoV-2 inhibition by a single TIP administration was obtained from a within-host mathematical model based on differential equations [39]. Applied Mathematical Modelling 137 (2025) 115673 3 J.C. Muñoz-Sánchez, J.T. Lázaro, J. Hillung et al. Fig. 1. (a) Experimental setup for coronaviruses HCoV-OC43 and MHV in cell 6-wells plates [15]. (b) Passage dynamics for helper virus (HV) and defective viral genomes (DVGs), with large titer fluctuations for the experiments performed in [15]. DVG abundance estimated using [63]. (c) Schematic diagram of the dynamical system modelling withinand between-cell virus dynamics for the cell culture experiments performed in [64](see Section 2and [15]). The model considers infection of a susceptible population of host cells, 𝐶, by HV (red) and/or by DVGs (yellow), producing HV only-infected cells 𝐶𝑉, DVGs only-infected cells 𝐶𝐷, and double-infected cells 𝐶𝐷𝑉 . Within-cell replication involves the amplification of viral genomes, which produce DVGs at a rate 𝛽. Infected cells are lysed releasing HV and DVGs to the medium. (d) Time series for HCoV-OC43 viral particles infecting BHK-21 cells with different MOIs in cell cultures. (e) State variables of the model. Mathematical models describing trans interactions between viral genomes are found in the literature [40,41]. The dynamics of DIPs have been extensively studied as an extreme case of complementation [41–48]. For example, the early work by Szathmáry [43], presented structured deme models to provide a description of the coexistence of virus segments considering HV and DIPs, sensitive and resistant viruses together with DIPs, covirus pairs (i.e., virus that exist as two or more separated particles all of which must be present for the complete replication cycle to occur), and virus–covirus systems. A deeper analysis of the model with HV and DIPs was later developed in [44]by considering cell populations infected by particles differing in number. Different dynamics, such as stable fixed points, periodic orbits, and strange chaotic attractors were identified with this model. Later, Kirkwood and Bangham [49] developed a differential equations model to analyze a system consisting of well-mixed host cells, HV, and DIPs under serial passage dynamics. The model successfully explained various dynamic behaviours observed in cell culture experiments: fluctuations in virus accumulation during successive passages and self-curing involving the simultaneous extinction of the HV and DIPs as shown experimentally [46,50]. More recently, passage experiments of baculoviruses in moth larvae have also provided experimental evidence of chaotic dynamics between HV and DIPs containing long genomic deletions [32]. Among coronaviruses, deletions are the most common type of DVGs [13,15], formed through recombination due to conserved homology in specific regions and/or RNA structures [51–54]. Indeed, SARS-CoV-2 deletion DVGs have been pervasively found both in cell cultures and in patients [13,20]. In this regard, asymptomatic patients tended to have lower DVG loads than symptomatic ones, and highly diverse populations of DVGs were observed after long-term COVID-19 in an immunosuppressed patient, suggesting a relationship between interferon responses and DVGs [20]. Given the pervasiveness of DVGs in betacoronavirus populations [20]and the lack of evidences supporting their sustained interference activity (and hence their potential development as TIPs), Hillung et al. [15] performed long-term in vitro evolution experiments to explore the role of DVGs in betacoronaviruses dynamics of diversification and evolution. Two virus models were chosen for this study, the human coronavirus OC43 (HCoV-OC43) and the murine hepatitis virus (MHV). To better understand the population dynamics for this experimental system, here we develop and investigate a mathematical model gathering key interactions at the level of within-passage infection dynamics (Fig. 1). The model shows the presence of degenerate normally hyperbolic invariant manifolds (quasineutral invariant manifolds). The existence of such manifolds implies that the orbits Applied Mathematical Modelling 137 (2025) 115673 4 J.C. Muñoz-Sánchez, J.T. Lázaro, J. Hillung et al. in the phase space reach planes of equilibria, depending on parameters. That is, orbits are strongly attracted to a given equilibrium population which is not a single point attractor but a curve or plane filled with equilibria and different initial conditions reach different equilibrium population values. Research on quasineutral planes is rather limited. These neutral surfaces have been characterised in predator-prey models with Holling type III functional responses [55]and in socio-economical models [56]. Moreover, quasineutral states governed by lines or curves of equilibria have been identified in prey-predator models given by differential [57]and partial [58] differential equations, in Lotka-Volterra competition models [59], in strains’ competition models of disease dynamics [60], and in models of RNA genomes replication [61]. More recently, a quasineutral curve was found in a mathematical model for an autocatalytic replicator with an obligate parasite [62]. For this later case, a bistability mechanism determined whether a given initial condition achieved the quasineutral curve or co-extinction. The model investigated in here exhibits a variety of quasineutral objects (mainly planes) displaying, for some parameter combinations, tristability between different scenarios as a function of the initial conditions. 2. Summary of the experimental results To investigate the dynamics of defective viral genomes (DVGs) accumulation, we performed serial passages with two betacoronaviruses: HCoV-OC43 in either baby hamster kidney cells (BHK-21) or human large intestine carcinoma cells (HCT-8), and MHV in murine liver cells (CCL-9.1). Passages involved stochastically varying inocula size within two wide but disjoint intervals that can broadly be defined as low and high multiplicity of infection (MOI). At every passage, infectious viral titer (defined as concentration of plaque forming units, PFUs/mL) was measured as a proxy to helper virus (HV) accumulation by plaque assays, while the amount and type of the DVGs component of the evolving populations was evaluated at four equidistant passages by high-throughput RNA sequencing (RNA-seq). For illustrative purposes, here we will focus in the case of HCoV-OC43 in BHK-21 cells. In this particular combination, 49 serial passages were performed, and the DVG component evaluated every 12 passages. Fig. 1c shows the dynamics for the first 50 passages. For the high MOI treatment, the median MOI at the onset of each passage was 25 PFU/cell (IQR: 131.94), while for the low MOI treatment, it was 3.75 ×10 −4 PFU/cell (IQR: 0.025). See Hillung et al. [15]for more experimental details. Each passage involved the replication of the HV and DVGs in a cell culture dish with susceptible host cells. In order to characterize the infection dynamics in the cell culture and test the validity of the model here investigated, a second set of experiments was performed by Hillung et al. [64]. Specifically, the dynamics of virus accumulation was monitored within a single infectious passage. To do so, three independent confluent monolayers of BHK-21 cells were inoculated at two viral MOI in the order of units per cell with HCoV-OC43. Then, HV accumulation was evaluated as above at the time points indicated in Fig. 1b. 3. Mathematical model In this section we introduce the dynamical system modelling experiments by Hillung et al. [15,64] summarised in Section 2. The model describes the infection dynamics taking place within a single cell culture dish (Fig. 1a), considering a helper virus (HV), which infects and replicates within a population of susceptible cells, and produces defective viral genomes (DVGs). By assumption, a DVG can only infect a cell but cannot replicate or lysate cells on its own, since it needs the products from the HV. It is also assumed that DVGs have a replication advantage in cells coinfected with the HV. The presence of a DVG inside a cell, by contrast, hinders the replication of the virus, resulting in a fitness penalty. An important assumption made is that viral accumulation dominates against virus’ particles decay, i.e., the rate of genomic RNA synthesis is much greater than its degradation rate (previous research on RNA viruses supports this assumption [65]). Furthermore, the model does not consider spontaneous cells’ decay and presumes that cell death is driven by virus-induced lysis. All these assumptions are consistent with the experimental time scale (see below) and aim to focus on the dynamics that arise from the infection process. The state variables of the model are given by infecting particles, the HV (𝑉) and DVGs (𝐷), and four different types of cells, including susceptible (𝐶) and infected cells 𝐶𝑝, where the subscript indicates the particle (or particles) that has (have) infected the host cell, 𝑝 ∈{𝑉, 𝐷, 𝐷𝑉 }(see Fig. 1e). Without affecting any theoretical conclusion and to avoid the use of large numbers (for instance, the number of cells moves around 106) in the propagation of numerical errors, we have scaled all state variables by dividing them by 𝐶(0), the initial number of cells. The resulting system is dimensionless and it gives rise to cell variables ranging in the interval [0, 1], with 𝐶(0) =1, and viruses and DVGs taking values in an interval that depends on some parameters of the system. The infection process is modelled as follows. When a 𝐶encounters a 𝑉or a 𝐷it becomes an infected cell 𝐶𝑉or 𝐶𝐷, respectively. This process happens at an infection rate 𝜄whose inverse could be considered as the average time between infections. Each type of these infected cells can be coinfected or superinfected with the alternative particle type, resulting in a double-infected cell 𝐶𝐷𝑉 . 𝐶𝐷cells, as the DVG is not replicating itself, will remain in the same stage until a superinfection with 𝑉occurs, if so. The next process will be the replication of HV and DVGs and the lysate of the infected cells. 𝐶𝑉will be lysed resulting into 𝜂new HV particles and 𝛽𝜂 DVGs as a result of errors during the replication process. The total number of particles released from the lysis, or burst size, is 𝐵=𝜂(1 +𝛽). If replication occurs inside a 𝐶𝐷𝑉 cell, then the DVG will use the replication machinery of the HV to replicate, competing for common resources and resulting in a decrease of the HV accumulation. The viral production after the 𝐶𝐷𝑉 lysis will be reduced to 𝜂∕𝜅with 𝜅>1. On the other hand, as DVGs have a replication advantage denoted by 𝛿, if the virus offspring after lysate 𝐶𝐷𝑉 is 𝜂∕𝜅then the DVGs produced will be 𝜂(𝛿+𝛽)∕𝜅. This latter term accounts for the HV replication errors resulting into additional DVGs. Thus, the total offspring of the coinfected cells 𝐶𝐷𝑉 is 𝜂(1 +𝛿+𝛽)∕𝜅. Under the assumption that the total burst size of 𝐶𝑉and 𝐶𝐷𝑉 is the same, as both are the same type of cells, one could derive the penalty coefficient as 𝜅=1 +𝛿∕(1 + 𝛽). Applied Mathematical Modelling 137 (2025) 115673 5 J.C. Muñoz-Sánchez, J.T. Lázaro, J. Hillung et al. Table 1 Model parameters. All parameters are dimensionless except the infection and lysis rates with inverse time dimensions. HV: helper virus; DVG: defective viral genomes. Parameter Description Range 𝐵Number of particles released after cell lysis (burst-size) >0 𝜂Number of HVs produced per cell >0 𝛽Fraction of DVGs produced per HV due to erroneous replication [0,1] 𝜅Replication penalty for HV in cells coinfected with DVGs >1 𝛿Replication advantage of DVGs >1 𝜄Infection rate >0 𝛼Virus infection-induced cell death rate >0 𝑚Multiplicity of infection (MOI) >0 For the sake of clarity, the previous processes are represented stoichiometrically with the set of reactions 𝐶+𝑉𝜄 ⟶𝐶𝑉,(1) 𝐶+𝐷𝜄 ⟶𝐶𝐷,(2) 𝐶𝑉+𝐷𝜄 ⟶𝐶𝐷𝑉 ,(3) 𝐶𝐷+𝑉𝜄 ⟶𝐶𝐷𝑉 ,(4) 𝐶𝑉 𝛼 ⟶𝜂𝑉 +𝜂𝛽 𝐷, (5) 𝐶𝐷𝑉 𝛼 ⟶ 𝜂 𝜅𝑉+𝜂(𝛽+𝛿) 𝜅𝐷. (6) From these reactions, and using the law of mass action, the following set of autonomous ordinary differential equations (ODEs) is derived:  𝐶=−𝜄𝐶(𝑉+𝐷),(7)  𝐶𝑉=𝜄𝐶𝑉 −𝐶𝑉(𝜄𝐷+𝛼),(8)  𝐶𝐷=𝜄(𝐶𝐷−𝐶𝐷𝑉),(9)  𝐶𝐷𝑉 =𝜄(𝐶𝐷𝑉+𝐶𝑉𝐷)−𝛼𝐶 𝐷𝑉 ,(10)  𝑉=𝛼𝜂(𝐶𝑉+ 𝐶𝐷𝑉 𝜅)−𝜄𝑉(𝐶+𝐶𝐷),(11)  𝐷=𝛼𝛽𝜂(𝐶𝑉+ 𝐶𝐷𝑉 𝜅)+𝛼𝛿𝜂𝐶𝐷𝑉 𝜅−𝜄𝐷(𝐶+𝐶𝑉).(12) Here, 𝛼, 𝛽, 𝛿, and 𝜄are independent parameters, while 𝜂and 𝜅are derived from the stoichiometric relations through 𝜂=𝐵 1+𝛽and 𝜅=1+ 𝛿 1+𝛽.(13) Notice that the replication advantage of the DVG versus its HV, 𝛿, appears in the term 𝛼𝛿𝜂𝐶 𝐷𝑉 ∕𝜅. We assume 𝛽∈[0, 1], 𝛿>1, and 𝐵, 𝜄, 𝛼>0(see Table 1). The assumption 𝛽∈[0, 1] relies on the biological implausibility of 𝛽>1, since DVGs arise due to deletions of the HV genomes. Hence, 𝛽can be interpreted as the fraction of DVGs produced from HV during cell infection. As we mentioned above, the model assumes that the amplification of both HV and DVGs dominate over their degradation within the experimental time-scales. This assumption is grounded in the experimental measures of accumulation and degradation rates for HCoV-OC43 done in a single passage by Hillung et al. [64]. As it can be seen from this system, the scaling performed to the state variables formally affects only the expression of the parameter 𝜄. These data show that during the first 60 hours post-inoculation (hpi), the effect of degradation was negligible compared to that of production. Only after exhaustion of productive cells, degradation becomes relevant, though the rate of degradation was still 78.14% slower than the rate of production [64]. In any case, the possible consequences and impact on dynamics of relaxing this assumption are discussed below. 4. Results In the following sections, we compute the equilibria of system Eqs. (7)-(12)and discuss their stability. Phase diagrams for relevant biological parameters are numerically obtained.2Then, the basins of attraction of equilibria are numerically computed (Section 4.1). 2Numerical integrations of the ODEs have been performed using the Runge-Kutta-Fehlberg-Simó method of 7𝑡ℎ-8𝑡ℎ order with automatic step size control and local relative tolerance 10−13 . Applied Mathematical Modelling 137 (2025) 115673 6 J.C. Muñoz-Sánchez, J.T. Lázaro, J. Hillung et al. Section 4.2 provides further results, including those cases with the degradation of viral particles. Section 4.3 explores the particular case of multiplicity of infection (MOI) >1where most cells get infected simultaneously. Next, Section 4.4 contains a sensitivity analysis of the solutions with respect to parameters and to initial conditions. Last but not least, in Section 4.5 the experimental time series of the HV dynamics in cell cultures generated by Hillung et al. [64]have been fitted with the mathematical model using artificial intelligence. 4.1. Planes of equilibria, stability and basins of attraction As it is commonly done in Dynamical Systems Theory, the study of the dynamics of system (7)-(12)is first based on the computation of its equilibrium solutions and the analysis of their local stability. Lemma 1 (The origin and its local stability). The origin  =(0, 0, 0, 0, 0, 0), the full extinction of all the populations, is always an equilibrium point of the system (7)-(12). Moreover, for any value of the parameters, the eigenvalues of its jacobian matrix at the origin are 0(with multiplicity 4and semisimple) and −𝛼(double and semisimple as well). Besides the origin, other equilibria are found and discussed in the following proposition. Proposition 1 (Planes of equilibria). System (7)–(12)has three planes formed by equilibrium points (i.e., all the points forming the planes are fixed by the dynamics) that we label as Π𝐶𝐶 𝐷 , Π𝐶𝐷𝐷, and Π𝑉𝐷. The origin trivially belongs to all of these planes of equilibria, but it does not share their biological interpretation. Because of this, when we refer to these planes the origin will not be considered a part of them, but studied aside. These planes involve different biological equilibrium scenarios and are defined as follows: (a) Persistence only of non-infected cells and defective viral genomes (DVG)-infected cells: Π𝐶𝐶 𝐷 ={(𝐶,𝐶 𝑉,𝐶 𝐷,𝐶 𝐷𝑉 ,𝑉,𝐷)=(𝐶,0,𝐶 𝐷,0,0,0) |||𝐶,𝐶 𝐷≥0,(𝐶,𝐶 𝐷)≠(0,0)} The spectrum of the jacobian matrix at these points inside the planes is {0 (double, semisimple),−𝛼,−𝜄𝐶,𝜀±}(14) where 𝜀±=− 𝜄(𝐶+𝐶𝐷)+𝛼 2±1 2√4𝛼𝐶𝜂𝜄 +(𝜄(𝐶+𝐶𝐷)−𝛼)2 +4𝛼𝐶𝐷𝜂𝜄 𝜅.(15) For any value of the parameters, the discriminant of (15)is always non-negative, and so all the eigenvalues, for any point in Π𝐶𝐶 𝐷 , are real. This means that Π𝐶𝐶 𝐷 is a so-called normally hyperbolic invariant manifold (NHIM). This plane involves no population of HV (either free or inside cells) or double-infected cells, i.e., a situation of self-curing driven by the quick and efficient outcompetition of the HV by the DVGs. Non-infected cells still remain in the system. (b) Persistence only of DVGs and DVG-infected cells: Π𝐶𝐷𝐷={(𝐶,𝐶 𝑉,𝐶 𝐷,𝐶 𝐷𝑉 ,𝑉,𝐷)=(0,0,𝐶 𝐷,0,0,𝐷)|||𝐶𝐷,𝐷≥0,(𝐶𝐷,𝐷)≠(0,0)}. The spectrum of the jacobian at any of its points is {0 (double, semisimple),−𝜄𝐷, −(𝜄𝐷 +𝛼),𝜆±},(16) where 𝜆±=− 𝜄𝐶𝐷+𝛼 2±1 2√(𝜄𝐶 𝐷−𝛼)2+4𝛼𝐶𝐷𝜂𝜄 𝜅.(17) As in the case above, the discriminant is also non-negative for any choice of the parameters. Hence, all the eigenvalues are real and, thus, Π𝐶𝐷𝐷is also a NHIM. This plane also involves self-curing since, asymptotically, no HV are found as DVGs have steadily and efficiently displaced them from the system. (c) Persistence only of free HV and DVGs: Π𝑉𝐷 ={(𝐶,𝐶 𝑉,𝐶 𝐷,𝐶 𝐷𝑉 ,𝑉,𝐷)=(0,0,0,0,𝑉,𝐷)|||𝑉,𝐷≥0,(𝑉,𝐷)≠(0,0)},(18) involving no cells and only HVs and DVGs in the medium. The spectrum of the Jacobian at these points is given by {0,0,−𝜄𝑉 ,−𝛼,−𝜄(𝑉+𝐷),−(𝜄𝐷 +𝛼)} (19) Applied Mathematical Modelling 137 (2025) 115673 7 J.C. Muñoz-Sánchez, J.T. Lázaro, J. Hillung et al. The two (semisimple) 0-eigenvalues come from the fact that, for any (𝑉, 𝐷)the point (0, 0, 0, 0, 𝑉, 𝐷)is an equilibrium point. All the eigenvalues are non-positive (and semisimple) and so Π𝑉𝐷 is locally attracting for any (𝑉, 𝐷) ≠(0, 0). This plane represents the most common outcome in which HV ends up killing all cells and DVGs are unavoidably present as byproducts of HV replication. It is not difficult to check that there are no other equilibrium solutions out from the origin and these three planes. Let us discuss in more detail the local stability of the equilibrium points contained in the planes. Regarding those on Π𝐶𝐷𝐷, notice that their local stability depends on the sign of the eigenvalue 𝜆+(the rest are all negative, except those 0coming from the fact that any point on it with arbitrary 𝐶𝐷, 𝐷is also an equilibrium). Hence, the expression 𝜆+≥0⇔𝜄𝐶𝐷+𝛼≤√(𝜄𝐶𝐷−𝛼)2+4𝛼𝐶𝐷𝜂𝜄 𝜅 can be squared (no spurious solutions are introduced since the discriminant is always non-negative), and it leads to (𝜄𝐶𝐷+𝛼)2−(𝜄𝐶𝐷−𝛼)2≤4𝛼𝐶𝐷𝜂𝜄 𝜅⇔4𝜄𝐶𝐷𝛼≤4𝛼𝐶𝐷𝜂𝜄 𝜅⇔𝜅≤𝜂⇔1+ 𝛿 1+𝛽≤𝐵 1+𝛽. That is, 𝜆+≥0⇔𝐵≥1+𝛽+𝛿. (20) So, if 𝐵>1 +𝛽+𝛿all the points on Π𝐶𝐷𝐷are unstable and if 𝐵<1 +𝛽+𝛿, then they are all locally attracting. Concerning the equilibrium points forming Π𝐶𝐶 𝐷 we have, like in the previous case, that their stability depends only on the sign of the eigenvalue 𝜀+. Thus, squaring again, we obtain 𝜀+≥0⇔4𝛼𝐶𝜂𝜄 +(𝜄(𝐶+𝐶𝐷)−𝛼)2 +4𝛼𝐶𝐷𝜂𝜄 𝜅≥(𝜄(𝐶+𝐶𝐷)+𝛼)2 ⇔4𝛼𝐶𝜂𝜄 +4𝛼𝐶𝐷𝜂𝜄 𝜅≥(𝜄(𝐶+𝐶𝐷)+𝛼)2 −(𝜄(𝐶+𝐶𝐷)−𝛼)2 ⇔4𝜄𝛼(𝐶+𝐶𝐷)≤4𝛼𝐶𝜂𝜄 +4𝛼𝐶𝐷𝜂𝜄 𝜅⇔(𝜂−1)𝐶≥(1− 𝜂 𝜅)𝐶𝐷. Since 𝜂−1= 𝐵 1+𝛽−1,1− 𝜂 𝜅=1− 𝐵 1+𝛽+𝛿, it follows that 𝜀+≥0⇔(𝐵 1+𝛽−1 )𝐶≥(1− 𝐵 1+𝛽+𝛿)𝐶𝐷.(21) From conditions (20) and (21)three possible cases arise (summarized in Table 2): (i) 𝐵>1 +𝛽+𝛿; (ii) 1 +𝛽<𝐵<1 +𝛽+𝛿; and (iii) 0 <𝐵<1 +𝛽. It is worth remarking that 𝐵is one of the parameters that can be better inferred experimentally by counting PFUs as described in Section 4.3. The plane of equilibria Π𝑉𝐷 (taking out the origin) is locally attracting in all three cases, so the following discussion will concern only the planes of equilibria Π𝐶𝐶 𝐷 and Π𝐶𝐷𝐷. (i)Case 𝑩>𝟏+𝜷+𝜹. On one hand, this case involves 𝜆+>0, and so all the equilibrium points forming Π𝐶𝐷𝐷are, simultaneously, unstable. On the other, 𝐵>1 +𝛽+𝛿>1 +𝛽and so 𝐵 1+𝛽−1>0and 1− 𝐵 1+𝛽+𝛿<0. This implies condition (21)and so all the equilibria forming Π𝐶𝐶 𝐷 become unstable. Biologically, this is the most commonly expected situation: 𝐵contains both HV and DVGs; the greater 𝛽+𝛿, the more DVGs are produced at the expense of HV production. (ii)Case 𝟏+𝜷<𝑩<𝟏+𝜷+𝜹. From (20)it follows that 𝜆+<0and, consequently, all the points constituting Π𝐶𝐷𝐷are locally attracting. Moreover, one has that 𝐵 1+𝛽−1>0and 1− 𝐵 1+𝛽+𝛿>0. Thus, from expression (21), it derives that the points (𝐶, 0, 𝐶𝐷, 0, 0, 0) ∈Π 𝐶𝐶 𝐷 are: Applied Mathematical Modelling 137 (2025) 115673 8 J.C. Muñoz-Sánchez, J.T. Lázaro, J. Hillung et al. Table 2 Stability analyses: attracting planes in terms of the values of the parameters 𝐵, 𝛽, and 𝛿. Case Locally attracting equilibrium points 𝐵>1+𝛽+𝛿Π𝑉𝐷 1+𝛽<𝐵<1+𝛽+𝛿Π𝑉𝐷 ,pointsinΠ𝐶𝐶 𝐷 such that (22)holdsandΠ𝐶𝐷𝐷 0<𝐵<1+𝛽Π𝑉𝐷 ,Π𝐶𝐶 𝐷 and Π𝐶𝐷𝐷 (iia) locally attracting if they belong to the half-plane 𝐶 𝐶𝐷 < 1− 𝐵 1+𝛽+𝛿 𝐵 1+𝛽−1 ,(22) (iib)or unstable (precisely, a saddle) if they fall in 𝐶 𝐶𝐷 > 1− 𝐵 1+𝛽+𝛿 𝐵 1+𝛽−1 .(23) Notice that condition (22) can be equivalently written as 𝐶<1− 𝜂 𝜅 𝜂−1𝐶𝐷, with 1 <𝜂<𝜅. Since 𝐶, 𝐶𝐷are positive, any choice (𝜂0, 𝜅0)inside the open sector bounded by the lines 𝜂=1and 𝜂=𝜅(with 𝜂, 𝜅>1) gives rise to a region of (𝐶, 𝐶𝐷) ∈(0, +∞) ×(0, +∞) of attracting equilibrium points. Such regions vary from the void situation (for 𝜂=𝜅), with no point inside, to the one where all the points are attractors (for 𝜂→1+). In the extreme cases where 𝜂=𝜅, the production of virus per cell would be compensated by the penalty that the HV receives by infecting a 𝐶𝐷cell. (iii)Case 𝟎<𝑩<𝟏+𝜷. Clearly, 𝐵<1 +𝛽+𝛿and therefore all the equilibrium points forming Π𝐶𝐷𝐷are locally attracting. Besides, one has that 𝜀+in equation (15)is strictly negative and, consequently, all the points in Π𝐶𝐶 𝐷 are also locally attracting. Remark 1. From Proposition 1and its local stability analysis above, the following statements should be highlighted: 1. Recall that the planes Π𝑉𝐷 , Π𝐶𝐶 𝐷 , and Π𝐶𝐷𝐷are highly degenerate in the sense that they are formed by equilibrium points. This fact implies the existence of a couple of zero eigenvalues of the jacobian matrix at any of the equilibrium points. 2. The equilibrium points (0, 0, 0, 0, 𝑉, 𝐷)constituting the plane Π𝑉𝐷 are all locally attracting for any value of the parameters and for any value of 𝑉and 𝐷. Π𝑉𝐷 is a NHIM (attracting in this case). 3. The local stability of the equilibrium points (𝐶, 0, 𝐶𝐷, 0, 0, 0) ∈Π 𝐶𝐶 𝐷 depends on the sign of its eigenvalue 𝜀+=𝜀+(𝐶, 𝐶𝐷). In particular, as seen in (ii) above, there exists a line on Π𝐶𝐶 𝐷 , namely 𝐶 𝐶𝐷 = 1− 𝐵 1+𝛽+𝛿 𝐵 1+𝛽−1 , separating those points that are stable from those that are unstable. 4. The local stability of the equilibrium points (0, 0, 𝐶𝐷, 0, 0, 𝐷) ∈Π 𝐶𝐷𝐷depends only on the value of 𝐶𝐷. Precisely, it is attracting if and only if its eigenvalue 𝜆+=𝜆+(𝐶𝐷) <0 ⇔𝐵<1 +𝛽+𝛿, a condition that relies only on the parameters of the system. Moreover, one has that 𝜆+(𝐶𝐷) =𝜀+(0, 𝐶𝐷). 5. Since all the eigenvalues of all the jacobian matrices around the equilibrium points are real, the bifurcations (in stability) that these points undergo are always of transcritical type, that is, given by a change in the sign of a real eigenvalue. The results shown in Table 2lead to some multistability scenarios that can be numerically analysed in terms of the parameters. For example, let us consider the orbits of system (7)-(12)with initial conditions (𝐶(0),𝐶 𝑉(0),𝐶 𝐷(0),𝐶 𝐷𝑉 (0),𝑉(0),𝐷(0)) = (1,0,0,0,𝑚𝑞𝑉0,𝑚(1 − 𝑞𝑉0)),(24) and where 𝑚=𝑉(0) + 𝐷(0) 𝐶(0) =𝑉(0) + 𝐷(0) Applied Mathematical Modelling 137 (2025) 115673 9 J.C. Muñoz-Sánchez, J.T. Lázaro, J. Hillung et al. Fig. 2. Phase diagrams in the plane (𝑞𝑉0, 𝐵). The (approximate) 𝜔-limits are numerically obtained for an orbit with initial conditions (1, 0, 0, 0, 𝑚𝑞𝑉0, 𝑚(1 −𝑞𝑉0)) in terms of 𝐵and the ratio 𝑞𝑉0=𝑉(0)∕(𝑉(0) +𝐷(0). Different colours show the 𝜔-limit given by: Π𝑉𝐷 (blue), Π𝐶𝐷𝐷(olive), and Π𝐶𝐶 𝐷 (black). An equilibrium is assumed to be reached when ‖𝐹‖2<10−12 , where 𝐹is the vector field of system (7)-(12). Dashed lines correspond to 𝐵=1 +𝛽and 𝐵=1 +𝛽+𝛿. The results are displayed for different values of 𝑚: from left to right 𝑚 =0.01, 𝑚 =0.1, 𝑚 =1, and 𝑚 =10. In all the panels, we have used 𝛽=0.2, 𝛿=2, and 𝜄∕𝛼=10. Fig. 3. (a) Same as in Fig. 2but fixing 𝐵=1.5, 𝛽=0.2, 𝜄∕𝛼=10, and tuning 𝛿. (b) Phase diagrams in the plane (𝑞𝑉0, 𝛿)for fixed 𝑚 =0.5and different values of 𝐵: from left to right and top to bottom 𝐵=1.1, 𝐵=1.5, 𝐵=2and 𝐵=3. The other parameters are the same as in (a). (c) Examples of trajectories with variables larger than zero at the three planes of equilibria. Numerical values for the different graphs: (i) (𝑞𝑉0, 𝑚) =(0.85, 0.1); (ii) (𝑞𝑉0, 𝑚) =(0.85, 0.5); (iii) (𝑞𝑉0, 𝑚) =(0.45, 0.1); (iv) (𝑞𝑉0, 𝑚) =(0.45, 0.5); (v) (𝑞𝑉0, 𝑚) =(0.15, 0.1); (vi) (𝑞𝑉0, 𝑚) =(0.15, 0.5). The labels (i),...,(vi) of the time series (at the bottom) correspond to locations in the attraction basin diagrams (top). Other examples of trajectories for all the variables can be found at Fig. A1 in Appendix A. is a fixed MOI. Observe that, in this sense, the parameter 𝑞𝑉0provides the virus proportion with this MOI, that is, the ratio 𝑉(0)∕(𝑉(0) + 𝐷(0)). For any choice (𝑞𝑉0,𝐵), we compute the approximate 𝜔-limit3of its corresponding orbit and assign it a colour depending on the plane of equilibria (Π𝑉𝐷, Π𝐶𝐷𝐷, or Π𝐶𝐶 𝐷 ) it achieves. This study leads to different diagrams in terms of the values of 𝑚. Some of them have been depicted in Fig. 2. Diagrams computed under changes in 𝛽of around a 10% −30% provide qualitatively similar results to the ones in Fig. 2. On the contrary, an increase in the value of 𝛿(from 1.02 to 5in these examples) exhibits substantial growth of the region with final points on the plane Π𝐶𝐶 𝐷 (results not shown). Similar explorations can be carried out by setting the values of 𝛽and 𝐵and studying the behaviour when 𝛿≥1varies. As a sample, some of these plots and diagrams are depicted in Fig. 3for fixed values of the parameters and varying 𝛿. It is noteworthy to analyse the evolution of the basins of attraction of the planes Π𝐶𝐶 𝐷 , Π𝑉𝐷, and Π𝐶𝐷𝐷of Table 2and their biological interpretation in terms of 𝐵, 1 +𝛽, and 1 +𝛽+𝛿. For instance, the case 𝐵<1 +𝛽+𝛿exhibits tristability, in which self3Recall that, roughly speaking, the 𝜔-limit of an orbit 𝜓(𝑡)can be defined as 𝜔(𝜓) = lim𝑡→+∞ 𝜓(𝑡). Applied Mathematical Modelling 137 (2025) 115673 16 J.C. Muñoz-Sánchez, J.T. Lázaro, J. Hillung et al. Thus, adding both contributions, one obtains 𝑉f=𝑉(𝑖) f+𝑉(𝑟) f=𝜂(𝐶0−𝐷0)+ 𝜂 𝜅𝐷0+(𝑉0−𝐶0) 𝐷f=𝐷(𝑖) f+𝐷(𝑟) f=𝛽𝜂(𝐶0−𝐷0)+ 𝜂 𝜅(𝛽+𝛿)𝐷0, and its sum 𝑉f+𝐷f=(𝑉0−𝐶0)+(1+𝛽)𝜂(𝐶0−𝐷0)+ 1+𝛽+𝛿 𝜅𝜂𝐷0=𝑉0+𝐶0(𝐵−1 ), where it has been taken into account that 𝐵=(1 +𝛽)𝜂and that 1+𝛽+𝛿 𝜅𝜂=(1+𝛽)𝜂=𝐵. 4.4. Sensitivity analysis A crucial piece of information that mathematical models provide, beyond predictions and qualitative interpretations, is to get insights into the relevance of their parameters determining or affecting its dynamics (sensitivity). It is well known in Dynamical Systems Theory that this information can be obtained from the so-called variational equations. In its general form, they supply information about the dependence of any particular solution with respect to its initial conditions. But besides, variational equations with respect to parameters provide knowledge on the sensitivity of such solutions with regard to these parameters. These two complementary sources of information can be both computed using a similar procedure. Let us first introduce the study regarding parameters and show, afterwards, the one corresponding to the initial conditions. A practical way of getting variational equations in the first case is to differentiate Eqs. (7)-(12)with respect to the chosen parameter, to commute such derivative with the one with respect to the time variable 𝑡, and to solve the resulting differential equation. As an illustrative example, we detail the case of the variational with respect to the parameter 𝛽. Let us denote by (𝐶(𝑡; 𝛽0), 𝐶𝑉(𝑡; 𝛽0), … , 𝐷(𝑡; 𝛽0)) the solution of equations (7)-(12)with initial condition (𝐶0, … , 𝐷0). Let now (𝐶(𝑡; 𝛽), 𝐶𝑉(𝑡; 𝛽), … , 𝐷(𝑡; 𝛽)) be the solution of the same Cauchy problem but for a parameter 𝛽close to 𝛽0. We wonder how close the solution (𝐶(𝑡; 𝛽), … , 𝐷(𝑡; 𝛽)) will evolve (in time) with respect to (𝐶(𝑡; 𝛽0), … , 𝐷(𝑡; 𝛽0)). To do it, we compute the variational equation with respect to 𝛽along (𝐶(𝑡; 𝛽0), … , 𝐷(𝑡; 𝛽0)), given by the following system of ODEs: 𝑑 𝑑𝑡 (𝜕𝐶 𝜕𝛽 )=𝜕 𝜕𝛽  𝐶=−𝜄(𝑉+𝐷)𝜕𝐶 𝜕𝛽 −𝜄𝐶 𝜕𝑉 𝜕𝛽 −𝜄𝐶 𝜕𝐷 𝜕𝛽 , 𝑑 𝑑𝑡 (𝜕𝐶𝑉 𝜕𝛽 )=𝜕 𝜕𝛽  𝐶𝑉=𝜄𝑉 𝜕𝐶 𝜕𝛽 −(𝜄𝐷 +𝛼) 𝜕𝐶𝑉 𝜕𝛽 +𝜄𝐶 𝜕𝑉 𝜕𝛽 −𝜄𝐶𝑉 𝜕𝐷 𝜕𝛽 , 𝑑 𝑑𝑡 (𝜕𝐶𝐷 𝜕𝛽 )=𝜕 𝜕𝛽  𝐶𝐷=𝜄𝐷 𝜕𝐶 𝜕𝛽 −𝜄𝑉 𝜕𝐶𝐷 𝜕𝛽 −𝜄𝐶𝐷 𝜕𝑉 𝜕𝛽 +𝜄𝐶 𝜕𝐷 𝜕𝛽 , 𝑑 𝑑𝑡 (𝜕𝐶𝐷𝑉 𝜕𝛽 )=𝜕 𝜕𝛽  𝐶𝐷𝑉 =𝜄𝐷 𝜕𝐶𝑉 𝜕𝛽 +𝜄𝑉 𝜕𝐶𝐷 𝜕𝛽 −𝛼𝜕𝐶𝐷𝑉 𝜕𝛽 +𝜄𝐶𝐷 𝜕𝑉 𝜕𝛽 +𝜄𝐶𝑉 𝜕𝐷 𝜕𝛽 , 𝑑 𝑑𝑡 (𝜕𝑉 𝜕𝛽 )=𝜕 𝜕𝛽  𝑉=(𝛼𝜕𝜂 𝜕𝛽 (𝐶𝑉+ 𝐶𝐷𝑉 𝜅)−𝛼𝜂 𝜅2 𝜕𝜅 𝜕𝛽 𝐶𝐷𝑉 ) −𝜄𝑉 𝜕𝐶 𝜕𝛽 +𝛼𝜂𝜕𝐶𝑉 𝜕𝛽 −𝜄𝑉 𝜕𝐶𝐷 𝜕𝛽 +𝛼𝜂 𝜅 𝜕𝐶𝐷𝑉 𝜕𝛽 −𝜄(𝐶+𝐶𝐷)𝜕𝑉 𝜕𝛽 , 𝑑 𝑑𝑡 (𝜕𝐷 𝜕𝛽 )=𝜕 𝜕𝛽  𝐷=(𝛼(𝜂+𝛽𝜕𝜂 𝜕𝛽 )(𝐶𝑉+ 𝐶𝐷𝑉 𝜅)−𝛼𝛽𝜂 𝜅2𝐶𝐷𝑉 𝜕𝜅 𝜕𝛽 +𝛼𝛽 𝜕𝜂 𝜕𝛽 𝐶𝐷𝑉 𝜅−𝛼𝛽𝜂 𝜅2𝐶𝐷𝑉 𝜕𝜅 𝜕𝛽 ) −𝜄𝐷 𝜕𝐶 𝜕𝛽 +(𝛼𝛽𝜂 −𝜄𝐷) 𝜕𝐶𝑉 𝜕𝛽 +1 𝜅(𝛼𝛽𝜂 +𝛼𝛽𝜂) 𝜕𝐶𝐷𝑉 𝜕𝛽 −𝜄(𝐶+𝐶𝑉)𝜕𝐷 𝜕𝛽 , with 𝜕𝜂 𝜕𝛽 =− 𝐵 (1 + 𝛽)2,𝜕𝜅 𝜕𝛽 =− 𝛿 (1 + 𝛽)2, and where the variables 𝐶, 𝐶𝑉, … , 𝐷are evaluated at (𝑡; 𝛽0)and 𝛽=𝛽0. The initial conditions of this system are 𝜕𝐶 𝜕𝛽 (0) = 𝜕𝐶𝑉 𝜕𝛽 (0) = 𝜕𝐶𝐷 𝜕𝛽 (0) = 𝜕𝐶𝐷𝑉 𝜕𝛽 (0) = 𝜕𝑉 𝜕𝛽 (0) = 𝜕𝐷 𝜕𝛽 (0) = 0. Applied Mathematical Modelling 137 (2025) 115673 17 J.C. Muñoz-Sánchez, J.T. Lázaro, J. Hillung et al. Fig. 7. Time solutions of the variational equations with respect to parameters 𝛽, 𝛿, and 𝜄∕𝛼(top) and initial conditions 𝐶(0), 𝑉(0), 𝐷(0) (bottom) for high multiplicity of infection (MOI) (𝑚 = 100) around the solution plotted in Fig. 6with no degradation and values 𝑞𝑉0=0.75, 𝐵= 500, 𝛽=10 −6, 𝛿=10, and 𝜄∕𝛼=0.1. The solutions 𝜕𝐶(𝑡)∕𝜕𝛽, 𝜕𝐶𝐷(𝑡)∕𝜕𝛽, ... 𝜕𝐷(𝑡)∕𝜕𝛽 of this system provide the (first order) time evolution of the variation of the variables 𝐶, 𝐶𝑉, … , 𝐷for 𝛽close to 𝛽0. Thus, for instance, in the case of 𝑉(𝑡; 𝛽)we have 𝑉(𝑡;𝛽)=𝑉(𝑡;𝛽0)+ 𝜕𝑉 (𝑡;𝛽0) 𝜕𝛽 (𝛽−𝛽0)+((𝛽−𝛽0)2), and therefore, for small values of |𝛽−𝛽0|, we get that 𝑉(𝑡;𝛽)∼𝑉(𝑡;𝛽0)+ 𝜕𝑉 (𝑡;𝛽0) 𝜕𝛽 (𝛽−𝛽0), where 𝜕𝑉 𝜕𝛽 (𝑡; 𝛽0)has been obtained solving the variational equation above. Large values (positive or negative) of 𝜕𝑉 (𝑡; 𝛽0)∕𝜕𝛽 at a time, say, 𝑡 =𝑡∗would correspond to significant variation (growth or decay) in the value of 𝑉(𝑡; 𝛽)with respect to 𝑉(𝑡; 𝛽0)for 𝛽∼𝛽0. A similar sensitivity analysis can be performed for the rest of the state variables and with respect to the parameters 𝐵, 𝛿, and 𝜄∕𝛼. Furthermore, since the experimental data measurements are reduced to values of 𝐶, 𝑉, and 𝐷, we focus our attention only on the variation of these variables with respect to the selected parameters. An illustrative example of the solution of these variational equations (with respect to the parameters 𝐵, 𝛽, 𝛿, and 𝜄∕𝛼) can be found in Fig. 7. On one hand, the variationals for 𝑉are always in the negative domain, suggesting that increases in three parameters always translate in reductions of HV accumulation. On the other, the variationals for 𝐷are always in the positive domain, indicating that increases in the magnitude of any of the three parameters result in more accumulation of DVGs at the cost of the HV. Interestingly, the variational equations with respect to 𝛿show a small range of values, thus suggesting a weak dependence of the dynamics on the advantage DVGs have on HV in terms of replication. In sharp contrast, the dependence on 𝛽is strong and even stronger in the case of 𝜄∕𝛼. Together, these observations support the idea that DVGs mostly gain advantage by interfering with the virus (via 𝛿) and, interestingly, reducing the efficiency by which free virus results in virus-producing infected cells (see (34)and discussion therein). Variational analysis can also be extended to the initial conditions of the system. The derivation of the variational equations of a solution with respect to its initial conditions is deferred to a brief explanation in Section Bin the Appendix. As in here, our analysis only focused on the variables 𝐶, 𝑉, and 𝐷. Fig. 7illustrates the effect of 𝐶(0), 𝑉(0) and 𝐷(0) on the result of the process. Here, the situation is more complex than in the case of parameters. The number of available cells at the beginning of the infection has a strong positive effect both in the 𝑉and 𝐷, although the accumulation of DVGs is remarkable more sensitive to 𝐶(0) than the accumulation of HVs. In contrast, the effect of the starting number of HVs and of DVGs has minor effects on the outcome of the process. As expected, starting with more 𝑉particles slightly benefits the HV and weakly penalizes the DVGs. This has to be understood in a situation in which 𝑚 = 100 and 𝑞𝑉0=0.75. Likewise, increasing the number of 𝐷particles in the inoculum slightly penalizes the HV in the same small magnitude that favours DVGs accumulation. 4.5. Model fitting to HCoV-OC43 data: parameters’ estimation To assess the accuracy of the model in reproducing real time-series data and to estimate model parameters, we have fitted the model using artificial intelligence to the data reported in Hillung et al. [15], shown in Fig. 1b. The aim of such work was to compare Applied Mathematical Modelling 137 (2025) 115673 18 J.C. Muñoz-Sánchez, J.T. Lázaro, J. Hillung et al. Table 3 Ranges of the parameters of the model explored by the genetic algorithm. Parameter 𝐵∈[10 1,104] 𝛽∈[10 −8,100] 𝛿∈[1,200] 𝛼∈[10 −3,101]h−1 𝜄∈[10 −5,100]h−1 𝛾∈[0.01,0.10] h−1 𝑉0∈{[10 5,107]for 𝑚𝑉=1.8∥[10 6,108]for 𝑚𝑉=3.8} 𝐷0∈{[1,107]for 𝑚𝑉=1.8∥[1,108]for 𝑚𝑉=3.8} the decay of the HCoV-OC43 in cell culture lysates and in fresh media. For this reason the experiment was conducted for several hours after cellular extinction, including degradation processes in our model. A genetic algorithm (GA) has been used to estimate the vector of parameters (𝐵, 𝛽, 𝛿, 𝛼, 𝜄, 𝛾), and the initial conditions 𝑉0, 𝐷0better fitting the experimental data. The size of the population of vectors of parameters was fixed at 600, and the number of generations was set to 104. The exploration of the parameter space was refined through successive searches, and the parameter’s range constrained by biological and experimental considerations. Within each generation, the 5% of the parameters with the lowest values (see below) were designated as the elite population, remaining constant for the next generation. The next-best 80% underwent parameter crossover, while the remaining 15%, representing the least favourable parameters, experienced random mutations. The parameter space was subdivided into a logarithmic scale to ensure a more balanced exploration, given that the intervals spanned several orders of magnitude. Exceptions were made for 𝛿and 𝛾, where linear scale searches were performed due to the nature of their intervals. Concerning the degradation rate, 𝛾, a narrow interval was chosen around the reported values at [15]. The explored parameter ranges are listed in the Table 3. The cost function used to compute the differences between the experimental and the simulated data is defined as follows: =log(𝑁 ∑ 𝑖=1 (log(𝑉𝑖⋅ 𝑉−1 𝑖))2 +𝜎(log 4 ∑ 𝑗=1  𝐶𝑗(𝑡= 65);𝑘, 𝑥∗,𝑝))(39) with 𝜎(𝑥;𝑘,𝑥∗,𝑝)= 𝑝 1+𝑒−𝑘(𝑥−𝑥∗).(40) Here, {𝑡𝑖, 𝑉𝑖}𝑁 𝑖=1 represent the experimental data and  𝑉𝑖denotes the values for the viral population obtained with the mathematical model for a given set of parameters. Additionally, ∑4 𝑗 𝐶𝑗(𝑡 = 65) corresponds to the total number of cells at 𝑡 =65hpi, a time point estimated to have a low remaining cell count. The function is formulated as the logarithm of a 𝜒2-function for the logarithm of the experimental data, augmented by a penalty term. The penalty term takes the form of a sigmoidal function with adjustable parameters and is incorporated to account for the approximate time of cell death on the culture, enhancing alignment with the experimental dataset. The parameters of the penalty term have been empirically adjusted to guide the GA in selecting parameter families where the numerical integration of the model significantly reduces the cell count at 𝑡 =65hpi (with 𝑘 =5, 𝑥∗=−1, 𝑝 =10). External estimates of the cell count over time would enable the removal of this penalty term. To assess the robustness of the final result, it was decided to conduct five identical batches, thereby obtaining five distinct populations of optimized parameters. While performing this procedure, we observed that after 104generations, the parameters of the elite population were nearly identical. Consequently, we opted to extract the parameter vectors from each batch and compute the mean and standard deviation across batches. For the dataset obtained from inoculation at viral MOI = 1.8, the estimated parameters that best fitted the experimental data were found to be 𝐵=78, 𝛽=10 −8, 𝛿=1, 𝛼=0.1827 h−1, 𝜄 =0.0027 h−1 and 𝛾=0.0154 h−1, with corresponding initial conditions 𝑉0=1, 516, 364 and 𝐷0=1. Similarly, for the ones obtained at viral MOI = 3.8, the estimated parameters yielding the optimal reproduction of experimental data were identified as 𝐵= 523, 𝛽=10 −8, 𝛿= 200, 𝛼=0.1898 h−1, 𝜄 =0.0006 h−1 and 𝛾=0.0361 h−1, with corresponding initial conditions 𝑉0=3, 676, 350 and 𝐷0=1. The standard deviation of the optimal parameters in both cases was less than 0.005% in all cases. The fittings of the mathematical model to the experimental data using the parameters obtained with the GA are displayed in Fig. 8. The optimal combination for both datasets, as evident from the resulting optimal combination, shows that 𝛽and 𝛿do not find an optimal value within the search range. However, examining the value for 𝛿, it is noteworthy that under such similar experimental conditions, their results are remarkably different. Consequently, we decided to analyze how the function behaved while keeping the remaining parameters fixed and varying both 𝛽and 𝛿(Appendix Fig. C1). As a result, we found that the function hardly changed its value within the parameter search range for 𝛽and 𝛿. This observation arises because one of the parameter combinations optimizing is the absence DVGs (note that 𝐷0does not also find an optimal value within the search range). In the nearly DVGs-free scenario (small 𝛽and 𝐷0), the value of 𝛿becomes irrelevant. These results underscore the flexibility of the model and the necessity of acquiring additional data for the fitting process for a rigorous interpretation of the parameters. Measures that could prove helpful would include, for instance, counting total cell numbers over time using life imaging techniques and teasing apart infected from non-infected cells by using viruses tagged with fluorescent Applied Mathematical Modelling 137 (2025) 115673 19 J.C. Muñoz-Sánchez, J.T. Lázaro, J. Hillung et al. Fig. 8. Experimental data (circles) and the corresponding fitting of the mathematical model (solid lines) using the parameters values optimised by the genetic algorithm. We show the cases for inoculation conducted at viral multiplicity of infection (MOI) 1.8(blue) and 3.8(red). proteins and living from dead cells by using dead-specific dyes such as trypan blue or propidium iodide. Lastly, to precisely estimate the amount of DVGs, it would be essential to extract total RNA from the cultures at various time points, subjecting them to RNA-seq analysis followed by bioinformatic analyses using tools such as DVGfinder [67]to identify different DVGs and estimators as described in [63]to quantify their abundance. 5. Discussion Defective viral genomes (DVGs) have been often viewed as replication byproducts, stemming from cell culture passaging conditions. However, a renewed interest in this fraction of the total viral population suggests that their existence might have more functional or biological significance than previously thought [4]. Given the notably high recombination rates seen in some viruses, including betacoronaviruses, it is plausible that DVGs could influence the adaptive evolution of the viable virus population [13]. The question of whether DVGs confer or not a selective advantage to the viral population and if the conditions of high multiplicity of infection (MOI) or localized co-infection hold biological relevance beyond cell culture remains to be explored. Serial passage experiments serve for many different applications in Virology, not only studying in vitro evolution [68]. The two most common applications are (i) the maintenance and amplification of viral stocks in the laboratory and (ii) the classic way of producing attenuated life vaccines by adapting the target virus to cells from a host as different as possible from the one to be immunised. DVGs have been found in live-attenuated vaccines for polio, measles and influenza viruses (reviewed in [4]). However, their impact on the development of protective immunity and vaccine efficacy has not been formally evaluated. Given their potential to interfere with and stimulate the immune system, there is speculation that DVGs could improve vaccine efficiency while ensuring virus safety by limiting its replication and spread. If this hypothesis holds true, it becomes crucial to carefully control the amount of DVGs in vaccine preparations to prevent complete interference and a significant reduction in the virus’ effectiveness. Our results suggest that the number of DVGs generated from helper virus (HV) replication (𝛽) and the replicative advantage of DVGs (𝛿) are the two more relevant parameters to be manipulated in order to optimize the ratio HV:DVGs in terms of vaccine efficiency. Multitude of mathematical models for HV-DVGs have been investigated at different scales [32,38,40–48]. Interestingly, these models have reproduced experimentally observed dynamics such as deterministic chaos [32]and self-curing [46,49], which involves the simultaneous extinction of the HV and defective interfering particles (DIPs), result that was observed in vitro by Jacobson et al. [50]and Stauffer Thompson and Yin [46]. Our model also identifies combinations of parameters in which self-curing arises: large 𝛿 values would in turn allow for the relationship 1 +𝛽+𝛿>𝐵to be fulfilled, shifting from the scenario of a single attracting plane Π𝑉𝐷 to a more complex scenario with three attracting planes and tristability, including the planes Π𝐶𝐶 𝐷 and Π𝐶𝐷𝐷(being these planes degenerate NHIMs) representing virus-free self-curing solutions. Model predictions are as valid and realistic as the assumptions from which they build up. Our first strong assumption, already discussed at large, was that at the time scales of the evolution experiments, virus production over-weights virus degradation. This assumption is well supported by experimental data (see [64]and references therein as well as [65]). Our model was designed to shed light on the dynamics of betacoronavirus populations during serial passages and thus could not be applied to a much more complex in vivo situation. For example, we are collapsing all possible DVGs into a single category, while in reality, viral populations contain a large fraction of very diverse DVGs that are in a dynamic equilibrium, with some appearing and disappearing at every transmission event and others persisting for long periods of time [15–17,20]. Our approach aligns with the current mathematical models mentioned in the previous paragraph that typically focus on a single dominant DVG. Further developments to incorporate the potential cooperation and competition among multiple types of DVGs are needed. Our model is similar to those proposed by Frank [45]and by Liang et al. [48]in the sense that they incorporate 𝑉, 𝐷, 𝐶, 𝐶𝐷, 𝐶𝐷𝑉 , and 𝐶𝑉as state variables. However, there are two major differences. Firstly, their models consider a spatial distribution of cells and thus use reaction-diffusion systems. Secondly, they incorporate an additional category of cells, namely 𝐶∗ 𝑉, considered as cells that were infected by 𝑉some time ago, becoming 𝐶𝑉, and are already producing HV but cannot be superinfected by DVGs. In this case, the fate of 𝐶∗ 𝑉is to die out after some time producing only 𝑉. This makes a notable difference with our model, as we opted for Applied Mathematical Modelling 137 (2025) 115673 20 J.C. Muñoz-Sánchez, J.T. Lázaro, J. Hillung et al. simplicity by suppressing this category as the modulation of the quotient 𝜄∕𝛼embeds the scenario in which 𝐶𝑉lysates previous to superinfection with 𝐷(Appendix Fig. D2). Although we recognize that superinfection exclusion [69]might be a relevant process for many viruses, it is not a universal one, and no evidence of such process have been reported for coronaviruses. The model presented here provides two new results in theoretical Virology: (i) HV-DVGs replication and infection dynamics governed by quasineutral planes (degenerate 2D NHIMs formed by equilibrium points in this case); (ii) scenario of tristability formed by these degenerate NHIMs. Concerning point (i) and, as far as we know, quasineutral dynamics in RNA virus dynamics have been up to now found in degenerate one-dimensional objects [61]. Hence, our findings extend this result on degenerate manifolds to planes. Concerning point (ii), these mentioned one-dimensional manifolds are usually found in scenarios of monostability and bistability. We here provide an example where three different quasineutral planes can be achieved depending on initial conditions. The existence of these degenerate planes is subject to non spontaneous degradation of viral particles and cells (𝛾=0), scenario that may seem unlikely in a real experiment. However, experimental estimates have indicated that degradation remains very low and that the dynamics within the time scale used in the experiments is dominated by virus replication over degradation [64]. Hence, we conjecture that the long time delays arising with parameter values close to the ones at which these degenerate NHIMs are found [62]could be observed in real experiments. Moreover, the passage experiments showed extremely large fluctuations between passages (some of them of about 1-2 orders of magnitude), as we show in Fig. 1b (see [15]for further details). These extremely large fluctuations could be due to the combined effect of stochastic sampling between passages and the dynamics tied to degenerate NHIMs. Under this scenario, different initial conditions (stochastically varying at each passage) may give place to different stationary values of HV and DVGs due to influence of NHIMs or to their remnants for those cases with some (and small) contribution of degradation. In conclusion, the mathematical model studied in this manuscript fits well the dynamics of virus accumulation within singlepassages i.e., withinand between-cell dynamics in a cell culture, in cell cultures of betacoronaviruses. Our parameter sensitivity analysis also suggests that the most relevant parameters to explain the observed dynamical patterns are the production of DVGs from the HV (𝛽) and the infection to lysate rates (𝜄∕𝛼). The combination of experiments and biologically inspired modelling provides a powerful tool to identify parameters to be optimized for future developments of DVGs as therapeutic interfering particles or as adjuvants in life attenuated vaccines. CRediT authorship contribution statement Juan C. Muñoz-Sánchez: Conceptualization, Data curation, Formal analysis, Investigation, Software, Visualization, Writing – original draft, Writing – review & editing. J. Tomás Lázaro: Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Software, Visualization, Writing – original draft, Writing – review & editing. Julia Hillung: Data curation, Methodology, Resources, Writing – review & editing. María J. Olmo-Uceda: Conceptualization, Data curation, Software, Writing – review & editing. Josep Sardanyés: Conceptualization, Data curation, Formal analysis, Methodology, Software, Visualization, Writing – original draft, Writing – review & editing. Santiago F. Elena: Conceptualization, Data curation, Formal analysis, Funding acquisition, Methodology, Project administration, Resources, Supervision, Writing – original draft, Writing – review & editing. Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Data availability Data will be made available on request. Acknowledgements We thank José A. Oteo for valuable suggestions and critical reading of the manuscript. JCM has been funded by grant ACIF/2021/296 (Generalitat Valenciana). MJO was funded by contract FPU19/05246 by MCIU/AEI/ 10.13039/501100011033 and “ESF invests in your future”. JTL has been funded by the projects PGC2018-098676-B-100 and PID2021-122954NB-I00 funded by MCIU/AEI/10.13039/501100011033/ and “ERDF a way of making Europe”, and by the grant “Ayudas para la Recualificación del Sistema Universitario Español 2021-2023”. JTL also thanks the Laboratorio Subterráneo de Canfranc, the I2SysBio and the Institut de Mathématiques de Jussieu-Paris Rive Gauche (Sorbonne Université) for their hospitality as hosting institutions of this grant. We also thank the MCIU/AEI/10.13039/501100011033/, through the María de Maeztu Program for Units of Excellence in R&D (CEX2020001084-M) and CERCA Programme/Generalitat de Catalunya for institutional support. JS has been also supported by the Ramón y Cajal grant RYC-2017-22243 funded by MCIU/AEI/10.13039/501100011033 and “ESF invests in your future”. SFE was supported by CSIC PTI Salud Global grant 202020E153 and by grants SGL2021-03-009 and SGL2021-03-052 from European Union Next Generation EU/PRTR through the CSIC Global Health Platform established by EU Council Regulation 2020/2094. Many computations were performed on the HPC cluster Garnatxa at I2SysBio (CSIC-UV). JCMS, JTL, MJOU, and SFE acknowledge the support of the Santa Fe Institute, where part of this research was developed. Applied Mathematical Modelling 137 (2025) 115673 21 J.C. Muñoz-Sánchez, J.T. Lázaro, J. Hillung et al. Appendix A. Time-series examples for the 𝚷𝑽𝑫 , 𝚷𝑪 𝑫𝑫and 𝚷𝑪𝑪 𝑫 Fig. A1. Time series to extinctions for cells (end state at Π𝑉𝐷 ) and for viral infective particles (IP) (end state at Π𝐶𝐷𝐷) and double extinction (end state Π𝐶𝐶 𝐷 ) Parameters: 𝐵=10, 𝛽=0.5and 𝛿=20and 𝜄∕𝛼= 100. Initial conditions: (a) 𝑉0=10 −3 and 𝐷0=10 −3; (b) 𝑉0=10 −3 and 𝐷0=0.9; and (c) 𝑉0=10 −3 and 𝐷0=0.5. 𝐶0=1. Appendix B. First order variational equation with respect to initial conditions In this section we present a well-known derivation of the so-called first variational equation with respect to initial conditions. To this end, let us consider the following Cauchy problem: 𝑥 =𝑓(𝑥),𝑥(0) = 𝑥0,𝑥∈ℝ𝑛,=𝑑 𝑑𝑡,(B.1) which, without loss of generality, we can assume autonomous and satisfying all the required conditions of smoothness and derivability, ensuring the existence and uniqueness of solutions. We denote by 𝑥(𝑡; 𝑥0)its unique solution. We are interested in the computation of the first-order (in 𝜀) variation of the solution 𝑥(𝑡; 𝑥0+𝜀)of the Cauchy problem 𝑥 =𝑓(𝑥),𝑥(0) = 𝑥0+𝜀, (B.2) where, abusing of notation, 𝑥0+𝜀 =𝑥0+𝜀Id, Id being the 𝑛-dimensional identity matrix. If 𝜀is small enough, we can write, using Taylor, 𝑥(𝑡;𝑥0+𝜀)=𝑥(𝑡;𝑥0)+𝜀𝜕𝑥(𝑡;𝑥0) 𝜕𝑥0 +(𝜀2).(B.3) Differentiating this expression with respect to 𝑡we get 𝑥(𝑡;𝑥0+𝜀)= 𝑥(𝑡;𝑥0)+𝜀𝑑 𝑑𝑡 (𝜕𝑥(𝑡;𝑥0) 𝜕𝑥0)+(𝜀2).(B.4) On the other side, using that 𝑥(𝑡; 𝑥0+𝜀)is a solution of (B.2)and using Taylor again, we have that 𝑥(𝑡;𝑥0+𝜀)=𝑓(𝑥(𝑡;𝑥0+𝜀)) = 𝑓(𝑥(𝑡;𝑥0)) + 𝜀𝜕𝑓(𝑥(𝑡;𝑥0)) 𝜕𝑥0 +(𝜀2) =𝑓(𝑥(𝑡;𝑥0)) + 𝜀𝐷𝑓(𝑥(𝑡;𝑥0)) 𝜕𝑥(𝑡;𝑥0) 𝜕𝑥0 +(𝜀2), where 𝐷𝑓 denotes the differential matrix of 𝑓with respect to 𝑥 =(𝑥1, 𝑥2, … , 𝑥𝑛). Equating terms of order 𝜀1in the latter expression with those in (B.4)it turns out that 𝜕𝑥(𝑡; 𝑥0)∕𝜕𝑥0satisfies the ODE 𝑑 𝑑𝑡 (𝜕𝑥(𝑡;𝑥0) 𝜕𝑥0)=𝐷𝑓(𝑥(𝑡;𝑥0)) 𝜕𝑥(𝑡;𝑥0) 𝜕𝑥0 .(B.5) Applied Mathematical Modelling 137 (2025) 115673 22 J.C. Muñoz-Sánchez, J.T. Lázaro, J. Hillung et al. Regarding its initial condition, we know that (substituting 𝑡 =0in expression (B.3)) 𝑥(0;𝑥0+𝜀)=𝑥(0; 𝑥0)+𝜀𝜕𝑥(0; 𝑥0) 𝜕𝑥0 +(𝜀2)⇔𝑥0+𝜀=𝑥0+𝜀𝜕𝑥(0;𝑥0) 𝜕𝑥0 +(𝜀2), from where, equating again powers in 𝜀, it follows that 𝜕𝑥(0;𝑥0) 𝜕𝑥0 =Id,(B.6) where Id is the identity matrix. Equations (B.5) and (B.6)are usually called first variational equations around a solution 𝑥(𝑡; 𝑥0). They provide information about the local dynamics (tangential and normal) around a solution 𝑥(𝑡; 𝑥0)of (B.1). This first variational equation is linear and homogeneous. Recurrently, one can compute the variational equations of higher order (in 𝜀). All of them are also linear but non-homogeneous. Appendix C. Additional fittings results Fig. C1. (a) Mean (solid line) and minimum (doted line) value of at each generation for data from inoculation at viral multiplicity of infection (MOI) 1.8(blue) and 3.8(red). (b) Value of for different (𝛽, 𝛿) combinations being the rest of the parameters selected as optimum at the genetic algorithm. The same study has been performed with 𝐷0=0for both datasets obtaining the same results. Applied Mathematical Modelling 137 (2025) 115673 23 J.C. Muñoz-Sánchez, J.T. Lázaro, J. Hillung et al. Appendix D. Model solutions for some 𝜾∕𝜶values Fig. D2. Some numerical integrations of the systems for 𝐵= 100, 𝛽=0.01 and 𝛿=10. From left to right the 𝜄∕𝛼quotient is 0.01, 0.1, 1 and 10. This figure aims to illustrate how the result depends on these parameters, both in terms of cellular type abundance and the helper virus (HV) - defective viral genomes (DVGs) final ratio (IP: infective particles). References [1] S. Duffy, L.A. Shackelton, E.C. Holmes, Rates of evolutionary change in viruses: patterns and determinants, Nat. Rev. Genet. 9 (2008) 267–276, https://doi .org / 10 .1038 /nrg2323. [2] R. Sanjuán, M.R. Nebot, N. Chirico, L.M. Mansky, R. Belshaw, Viral mutation rates, J. Virol. 84 (2010) 9733–9748, https://doi .org /10 .1128 /JVI .00694 -10. [3] R. Belshaw, R. Sanjuán, O.G. Pybus, Viral mutation and substitution: units and levels, Curr. Opin. Virol. 1 (2011) 430–435, https://doi .org /10 .1016 /j .coviro . 2011 .08 .004. [4] M. Vignuzzi, C.B. López, Defective viral genomes are key drivers of the virus–host interaction, Nat. Microbiol. 4 (2019) 1075–1087, https://doi .org /10 .1038 / s41564 -019 -0465 -y. [5] A.S. Huang, D. Baltimore, Defective viral particles and viral disease processes, Nature 226 (1970) 325–327, https://doi .org /10 .1038 /226325a0. [6] P. Von Magnus, S. Gard, Studies on interference in experimental influenza, Ark. Kemi, Mineral. Geol. 24 (1947) 4. [7] J. Xu, Y. Sun, Y. Li, G. Ruthel, S.R. Weiss, et al., Replication defective viral genomes exploit a cellular pro-survival mechanism to establish paramyxovirus persistence, Nat. Commun. 8 (2017) 799, https://doi .org /10 .1038 /s41467 -017 -00909 -6. [8] E. Genoyer, C.B. Lopez, Defective viral genomes alter how Sendai virus interacts with cellular trafficking machinery leading to heterogeneity in the production of viral particles among infected cells, J. Virol. 93 (2018) e01579-18, https://doi .org /10 .1128 /JVI .01579 -18. [9] D.R. Cave, F.M. Hendrickson, A.S. Huang, Defective interfering virus particles modulate virulence, J. Virol. 55 (1985) 366–373, https://doi .org /10 .1128 /jvi .55 . 2 .366 -373 .1985. [10] F.J. Fuller, P.I. Marcus, Interferon induction by viruses. IV. Sindbis virus: early passage defective-interfering particles induce interferon, J. Gen. Virol. 48 (1980) 63–73, https://doi .org /10 .1099 /0022 -1317 -48 -1 -63. [11] B.K. De, D.P. Nayak, Defective interfering influenza viruses and host cells: establishment and maintenance of persistent influenza virus infection in MDBK and HeLa cells, J. Virol. 36 (1980) 847–859, https://doi .org /10 .1128 /jvi .36 .3 .847 -859 .1980. [12] J.C. Kennedy, R.D. Macdonald, Persistent infection with infectious pancreatic necrosis virus mediated by defective-interfering (DI) virus particles in a cell line showing strong interference but little DI replication, J. Gen. Virol. 58 (1982) 361–371, https://doi .org /10 .1099 /0022 -1317 -58 -2 -361. [13] C. Campos, S. Colomer-Castell, D. Garcia-Cehic, J. Gregori, C. Andrés, et al., The frequency of defective genomes in Omicron differs from that of the Alpha, Beta and Delta variants, Sci. Rep. 12 (2022) 22571, https://doi .org /10 .1038 /s41598 -022 -24918 -8. [14] J. Gribble, L.J. Stevens, M.L. Agostini, J. Anderson-Daniels, J.D. Chappeli, et al., The coronavirus proofreading exoribonuclease mediates extensive viral recombination, PLoS Pathog. 17 (2021) e1009226, https://doi .org /10 .1371 /journal .ppat .1009226. [15] J. Hillung, M.J. Olmo-Uceda, J.C. Muñoz-Sánchez, S.F. Elena, Accumulation dynamics of defective genomes during experimental evolution of two betacoronaviruses, Viruses 16 (4) (2024) 644, https://doi .org /10 .3390 /v16040644. [16] M.A. Rangel, P.T. Dolan, S. Taguwa, Y. Xiao, R. Andino, et al., High-resolution mapping reveals the mechanism and contribution of genome insertions and deletions to RNA virus evolution, Proc. Natl. Acad. Sci. USA 120 (2023) e2304667120, https://doi .org /10 .1073 /pnas .2304667120. [17] E. Jaworski, A. Routh, Parallel ClickSeq and Nanopore sequencing elucidates the rapid evolution of defective-interfering RNAs in Flock House virus, PLoS Pathog. 13 (2017) e1006365, https://doi .org /10 .1371 /journal .ppat .1006365. [18] J.T. McCrone, R.J. Woods, E.T. Martin, R.E. Malosh, A.S. Monto, et al., Stochastic processes constrain the within and between host evolution of influenza virus, eLife 7 (2018) e35962, https://doi .org /10 .7554 /eLife .35962. [19] I.S. Novella, S.F. Elena, A. Moya, E. Domingo, J.J. Holland, Repeated transfer of small RNA virus populations leading to balanced fitness with infrequent stochastic drift, Mol. Gen. Genet. 252 (1996) 733–738, https://doi .org /10 .1007 /BF02173980. [20] T. Zhou, N.J. Gilliam, S. Li, S. Spandau, R.M. Osborn, et al., Generation and functional analysis of defective viral genomes during SARS-CoV-2 infection, mBio 14 (2023) e0025023, https://doi .org /10 .1128 /mbio .00250 -23. [21] M. Stampfer, D. Baltimore, A.S. Huang, Absence of interference during high-multiplicity infection by clonally purified vesicular stomatitis virus, J. Virol. 7 (1971) 409–411, https://doi .org /10 .1128 /JVI .7 .3 .409 -411 .1971. Applied Mathematical Modelling 137 (2025) 115673 24 J.C. Muñoz-Sánchez, J.T. Lázaro, J. Hillung et al. [22] A.S. Huang, Defective interfering viruses, Annu. Rev. Microbiol. 27 (1973) 101–117, https://doi .org /10 .1146 /annurev .mi .27 .100173 .000533. [23] J.J. Holland, L.P. Villareal, M. Breindl, Factors involved in the generation and replication of rhabdovirus defective T particles, J. Virol. 17 (1976) 805–815, https://doi .org /10 .1128 /JVI .17 .3 .805 -815 .1976. [24] S.A. Felt, E. Achouri, S.R. Faber, R. Sydney, C.B. López, Accumulation of copy-back viral genomes during respiratory syncytial virus infection is preceded by diversification of the copy-back viral genome population followed by selection, Virus Evol. 8 (2022) veac091, https://doi .org /10 .1093 /ve /veac091. [25] T.A. Damayanti, H. Nagano, K. Mise, I. Furusawa, T. Okuno, Brome mosaic virus defective RNAs generated during infection of barley plants, J. Gen. Virol. 80 (1999) 2511–2518, https://doi .org /10 .1099 /0022 -1317 -80 -9 -2511. [26] A.D.T. Barrett, N.J. Dimmock, Modulation of Semliki Forest virus-induced infection of mice by defective-interfering virus, J. Infect. Dis. 150 (1984) 98–104, https://doi .org /10 .1093 /infdis /150 .1 .98. [27] L. Roux, A.E. Simon, J.J. Holland, Effects of defective interfering viruses on virus replication and pathogenesis in vitro and in vivo, Adv. Virus Res. 40 (1991) 181–211, https://doi .org /10 .1016 /s0065 -3527(08 )60279 -1. [28] C.M. Smith, P.D. Scott, C. O’Callaghan, A.J. Easton, N.J. Dimmock, A defective interfering influenza RNA inhibits infectious influenza virus replication in human respiratory tract cells: a potential new human antiviral, Viruses 8 (2016) 237, https://doi .org /10 .3390 /v8080237. [29] F.M. Horodyski, S.T. Nichol, K.R. Spindler, J.J. Holland, Properties of DI particles resistant mutants of vesicular stomatitis virus isolated from persistent infections and from undiluted passages, Cell 33 (1983) 801–810, https://doi .org /10 .1016 /0092 -8674(83 )90022 -3. [30] N.J. DePolo, J.J. Holland, Very rapid generation/amplification of defective interfering particles by vesicular stomatitis virus variants isolated from persistent infections, J. Gen. Virol. 67 (1986) 1195–1198, https://doi .org /10 .1099 /0022 -1317 -67 -6 -1195. [31] N.J. DePolo, C. Giachetti, J.J. Holland, Continuing coevolution of virus and defective interfering particles and of viral genome sequences during undiluted passages: virus mutants exhibiting nearly complete resistance to formerly dominant defective interfering particles, J. Virol. 61 (1987) 454–464, https://doi .org / 10 .1128 /jvi .61 .2 .454 -464 .1987. [32] M.P. Zwart, G.P. Piljman, J. Sardanyés, J. Duarte, C. Januário, et al., Complex dynamics of defective interfering baculoviruses during serial passage in insect cells, J. Biol. Phys. 39 (2013) 327–342, https://doi .org /10 .1007 /s10867 -013 -9317 -9. [33] J. García-Arriaza, S.C. Manrubia, M. Toja, E. Domingo, C. Escarmís, Evolutionary transition towards defective RNAs that are infectious by complementation, J. Virol. 78 (2004) 11678–11685, https://doi .org /10 .1128 /JVI .78 .21 .11678 -11685 .2004. [34] A.K. Pattnaik, G.W. Wertz, Cells that express all five proteins of vesicular stomatitis virus from cloned cDNAs support replication, assembly, and budding of defective interfering particles, Proc. Natl. Acad. Sci. USA 88 (1991) 1379–1383, https://doi .org /10 .1073 /pnas .88 .4 .1379. [35] A.C. Marriott, N.J. Dimmock, Defective interfering viruses and their potential as antiviral agents, Rev. Med. Virol. 20 (2010) 51–62, https://doi .org /10 .1002 / rmv .641. [36] T. Notton, J. Sardanyés, A.D. Weinberger, L.S. Weinberger, The case for transmissible antivirals to control population-wide infectious disease, Trends Biotechnol. 32 (2014) 400–405, https://doi .org /10 .1016 /j .tibtech .2014 .06 .006. [37] Y. Xiao, P.V. Lidsky, Y. Shirogane, R. Aviner, C.T. Wu, et al., A defective viral genome strategy elicitgs broad protective immunity against respiratory viruses, Cell 184 (2021) 6037–6051.e14, https://doi .org /10 .1016 /j .cell .2021 .11 .023. [38] S. Chaturvedi, G. Vasen, M. Pablo, X. Chen, N. Beutler, et al., Identification of a therapeutic interfering particle -a single-dose SARS-CoV-2 antiviral intervention with a high barrier to resistance, Cell 184 (2021) 6022–6036.e18, https://doi .org /10 .1016 /j .cell .2021 .11 .004. [39] S. Chaturvedi, N. Beutler, G. Vasen, M. Pablo, X. Chen, et al., A single-administration therapeutic interfering particle reduces SARS-CoV-2 viral shedding and pathogenesis in hamsters, Proc. Natl. Acad. Sci. USA 119 (2022) e2204624119, https://doi .org /10 .1073 /pnas .2204624119. [40] H. Gao, M.W. Feldman, Complementation and epistasis in viral coinfection dynamics, Genetics 182 (2009) 251–263, https://doi .org /10 .1534 /genetics .108 . 099796. [41] J. Sardanyés, S.F. Elena, Error threshold in RNA quasispecies models with complementation, J. Theor. Biol. 265 (2010) 278–286, https://doi .org /10 .1016 /j .jtbi . 2010 .05 .018. [42] C.R.M. Bangham, T.B.L. Kirkwood, Defective interfering particles: effects in modulating virus growth and persistence, Virology 179 (1990) 821–826, https:// doi .org /10 .1016 /0042 -6822(90 )90150 -p. [43] E. Szathmáry, Natural selection and dynamical coexistence of defective and complementing virus segments, J. Theor. Biol. 157 (1992) 383–406, https://doi .org / 10 .1016 /s0022 -5193(05 )80617 -4. [44] E. Szathmáry, Co-operation and defection: playing the field in virus dynamics, J. Theor. Biol. 165 (1993) 341–356, https://doi .org /10 .1006 /jtbi .1993 .1193. [45] S.A. Frank, Within-host spatial dynamics of viruses and defective interfering particles, J. Theor. Biol. 206 (2000) 279–290, https://doi .org /10 .1006 /jtbi .2000 . 2120. [46] K.A. Stauffer Thompson, J. Yin, Population dynamics of an RNA virus and its defective interfering particles in passage cultures, Virol. J. 7 (2010) 257, https:// doi .org /10 .1186 /1743 -422X -7 -257. [47] L. Chao, S.F. Elena, Nonlinear trade-offs allow the cooperation game to evolve from Prisoner’s Dilemma to Snowdrift, Proc. R. Soc. B 284 (2017) 20170228, https://doi .org /10 .1098 /rspb .2017 .0228. [48] Q. Liang, J. Yang, W.T.L. Fan, W.C. Lo, Patch formation driven by stochastic effects of interaction between viruses and defective interfering particles, PLoS Comput. Biol. 19 (2023) e1011513, https://doi .org /10 .1371 /journal .pcbi .1011513. [49] T.B.L. Kirkwood, C.R.M. Bangham, Cycles, chaos, and evolution in virus cultures: a model of defective interfering particles, Proc. Natl. Acad. Sci. USA 91 (1994) 8685–8689, https://doi .org /10 .1073 /pnas .91 .18 .8685. [50] S. Jacobson, F.J. Dutko, C.J. Pfau, Determinants of spontaneous recovery and persistence in MDCK cells infected with lymphocytic choriomeningitis virus, J. Gen. Virol. 44 (1979) 113–122, https://doi .org /10 .1099 /0022 -1317 -44 -1 -113. [51] W.Y. Liao, T.Y. Ke, H.Y. Wu, The 3’-terminal 55 nucleotides of bovine coronavirus defective interfering RNA harbor cis-acting elements required for both negativeand positive-strand RNA synthesis, PLoS ONE 9 (2014) e98422, https://doi .org /10 .1371 /journal .pone .0098422. [52] P.A. Jennings, J.T. Finch, G. Winter, J.S. Robertson, Does the higher order structure of the influenza virus ribonucleoprotein guide sequence rearrangements in influenza viral RNA?, Cell 34 (1983) 619–627, https://doi .org /10 .1016 /0092 -8674(83 )90394 -x. [53] K. Saira, X. Lin, J.V. DePasse, R. Halpin, A. Twaddle, et al., Sequence analysis of in vivo defective interfering-like RNA of influenza a H1N1 pandemic virus, J. Virol. 87 (2013) 8064–8074, https://doi .org /10 .1128 /JVI .00240 -13. [54] E.Z. Poirier, B.C. Mounce, K. Rozen-Gagnon, P.J. Hooikaas, K.A. Stapleford, et al., Low-fidelity polymerases of alphaviruses recombine at higher rates to overproduce defective interfering particles, J. Virol. 90 (2016) 2446–2454, https://doi .org /10 .1128 /JVI .02921 -15. [55] M. Farkas, E. Sáez, I. Sźantó, Velcro bifurcation in competition models with generalized Holling functional response, Miskolc Math. Notes 6 (2005) 185–195, https://doi .org /10 .18514 /MMN .2005 .115. [56] A. Bocsó, M. Farkas, Political and economic rationality leads to velcro bifurcation, Appl. Math. Comput. 140 (2003) 381–389, https://doi .org /10 .1016 /S0096 - 3003(02 )00235 -7. [57] M. Farkas, Zip bifurcation in a competition model, Nonlinear Anal., Theory Methods Appl. 8 (11) (1984) 1295–1309, https://doi .org /10 .1016 /0362 -546X(84 ) 90017 -8. [58] J.D. Ferreira, L.A.F. de Olivera, Zip bifurcation in a competitive system with diffusion, Differ. Equ. Dyn. Syst. 17 (2009) 37–53, https://doi .org /10 .1007 /s12591 - 009 -0003 -0. [59] Y.T. Lin, H. Kim, C.R. Doering, Features of fast living: on the weak selection for longevity in degenerate birth-death processes, J. Stat. Phys. 148 (2012) 646–662, https://doi .org /10 .1007 /s10955 -012 -0479 -9. Applied Mathematical Modelling 137 (2025) 115673 25 J.C. Muñoz-Sánchez, J.T. Lázaro, J. Hillung et al. [60] O. Kogan, M. Khasin, B. Meerson, D. Schneider, C.R. Myers, Two-strain competition in quasineutral stochastic disease dynamics, Phys. Rev. E 90 (2014) 042149, https://doi .org /10 .1103 /PhysRevE .90 .042149. [61] J. Sardanyés, A. Arderiu, S.F. Elena, T. Alarcón, Noise-induced bistability in the quasineutral coexistence of viral RNA under different replication modes, J. R. Soc. Interface 15 (2018) 20180129, https://doi .org /10 .1098 /rsif .2018 .0129. [62] E. Fontich, A. Guillamon, T. Lázaro, T. Alarcón, B. Vidiella, et al., Critical slowing down close to a global bifurcation of a curve of quasi-neutral equilibria, Commun. Nonlinear Sci. Numer. Simul. 104 (2022) 106032, https://doi .org /10 .1016 /j .cnsns .2021 .106032. [63] J.C. Muñoz-Sánchez, M.J. Olmo-Uceda, J.A. Oteo, S.F. Elena, Quantifying defective and wild-type viruses from high-throughput RNA sequencing, bioRxiv, https:// doi .org /10 .1101 /2024 .07 .23 .604773, 2024. [64] J. Hillung, T. Lázaro, J.C. Muñoz-Sánchez, M.J. Olmo-Uceda, J. Sardanyés, et al., Decay of HCoV-OC43 infectivity is lower in cell debris-containing media than in fresh culture media, microPubl. Biol. (2024), https://doi .org /10 .17912 /micropub .biology .001092. [65] F. Martínez, J. Sardanyés, S.F. Elena, J.A. Daròs, Dynamics of a plant RNA virus intracellular accumulation: stamping machine vs. geometric replication, Genetics 188 (2011) 637–646, https://doi .org /10 .1534 /genetics .111 .129114. [66] J.M. Cuevas, A. Moya, R. Sanjuán, Following the very initial growth of biological RNA viral clones, J. Gen. Virol. 86 (2005) 435–443, https://doi .org /10 .1099 / vir .0 .80359 -0. [67] M.J. Olmo-Uceda, J.C. Muñoz-Sánchez, W. Lasso-Giraldo, V. Arnau, W. Díaz-Villanueva, et al., DVGfinder: a metasearch tool for identifying defective viral genomes in RNA-Seq data, Viruses 14 (2022) 1114, https://doi .org /10 .3390 /v14051114. [68] S.F. Elena, R. Sanjuán, Virus evolution: insights from an experimental approach, Annu. Rev. Ecol. Evol. Syst. 38 (2007) 27–52, https://doi .org /10 .1146 / ANNUREV .ECOLSYS .38 .091206 .095637. [69] M. Hunter, D. Fusco, Superinfection exclusion: a viral strategy with short term benefits and long-term drawbacks, PLoS Comput. Biol. 18 (2022) e1010125, https://doi .org /10 .1371 /journal .pcbi .1010125.