Exploring the therapeutic potential of defective interfering particles in reducing the replication of SARS-CoV-2

SARS-CoV-2 still presents a global threat to human health due to the continued emergence of new strains and waning immunity amongst vaccinated populations. Therefore, it is still relevant to investigate potential therapeutics, such as therapeutic interfering particles (TIPs). Mathematical and computational modelling are valuable tools to study viral infection dynamics for predictive analysis. Here, we expand on the previous work by Grebennikov et al. (2021) on SARS-CoV-2 intra-cellular replication dynamics to include defective interfering particles (DIPs) as potential therapeutic agents. We formulate a deterministic model that describes the replication of wild-type (WT) SARS-CoV-2 virus in the presence of DIPs. Sensitivity analysis of parameters to several model outputs is employed to inform us on those parameters to be carefully calibrated from experimental data. We then study the effects of co-infection on WT replication and how DIP dose perturbs the release of WT viral particles. Furthermore, we provide a stochastic formulation of the model that is compared to the deterministic one. These models could be further developed into population-level models or used to guide the development and dose of TIPs. Author summary SARS-CoV-2 continues to evolve, with new strains or sub-strains being identified thanks to efforts to monitor the virus. Consequently, new strains threaten human health as current vaccinations may not adequately protect against future strains. It is therefore important to understand the roles that additional therapeutics could play in protecting against these future strains. Therapeutic interfering particles (TIPs), otherwise referred to as defective interfering particles (DIPs), could provide an additional treatment option against future strains. Previous models have examined the role of DIPs at the within-host level during co-infection with wild-type virus, but have paid little attention to intra-cellular dynamics. Here we extend the previous intra-cellular replication model of SARS-CoV-2 by Grebennikov et al. (2021) to include co-infection of WT virus with DIPs. We show that DIPs lead to a reduction in the WT virus in a dose-dependent manner, with higher doses leading to up to 10-fold reduction in total WT virus released from a cell depending on the multiplicity of infection (MOI). We find these results to be consistent for both deterministic and stochastic formulations of the model. Our approaches could be developed into a within-host model or population-level model, which could then be used to guide therapeutic DIP doses.


Introduction
In December 2019, a new infectious disease was reported to the World Health Organisation (WHO) that would later be identified as a novel coronavirus (SARS-CoV-2) [1].By 30 January 2020, the WHO declared SARS-CoV-2 a "public health emergency of international concern" [2], as it rapidly spread to 113 countries.By the 11th of March 2020, it had caused 118,319 infections and 4,292 deaths.
Consequently, the WHO declared SARS-CoV-2 a pandemic [3,4], and as of the 29th of July 2022, about 572 million infections and over 6 million deaths have been recorded worldwide.During the early stages of the pandemic, treatment options were limited to chloroquine and remdesivir [5,6].However, since then several effective vaccines have been developed that provide protection and reduce transmission, with many countries rolling out mass vaccination programs [7].Although vaccines for SARS-CoV-2 now exist, the emergence of new strains due to mutations has led to further concerns about vaccine effectiveness [8,9].This fact, together with that of waning immunity and the existence of individuals who are unable to be vaccinated or out-right refuse to do so, highlight the need for additional therapeutics and prophylactics [10,11].
One such potential therapy is viral interfering particles.During viral replication, mutants lacking essential parts of the viral genome arise [12,13], which are unable to replicate in the absence of wild-type (WT) virus.These are known as defective interfering particles (DIPs).DIPs can be exploited to make therapeutic interfering particles (TIPs), which inhibit the replication of WT virus by outcompeting WT gene segments for resources required during viral replication and assembly [14,15].TIPs/DIPs have been investigated for several viruses, including HIV, Ebola, influenza, and SARS-CoV-2 and have been found to cause a two-fold reduction in viral titres [14][15][16].However, caveats exist in their production; for instance, which sections of the viral genome are to be removed to allow for replication at a faster rate than WT, they are virus-specific, and little is known about how mutations change replication dynamics [13,17].
From a mathematical modelling perspective, a long-standing effort exists to describe transmission dynamics at the population and within-host levels (see Ref. [18] and Refs.[19,20].However, little effort has been devoted to investigating the intra-cellular replication kinetics of DIPs in the presence of WT virus.Grebennikov et al. [21] have recently proposed a SARS-CoV-2 intra-cellular replication dynamics model.This model allowed for the quantification of viral genomes and proteins during the replication cycle.We wish to exploit this model to explore co-infection with DIPs and the effect of DIPs on the replication dynamics of the WT virus.In particular, in this study, we formulate a mathematical model of SARS-CoV-2 replication in a cell co-infected with DIPs.As in Ref. [21], we will follow a deterministic approach to calibrate model parameters.We shall use sensitivity analysis to study the impact parameters have on the release of both WT and DIP viral particles.We also introduce a stochastic description of this model to compare to the deterministic one.We shall also investigate how initial doses of each virus affect viral particle production (WT and DIPs) to quantify DIP-related inhibition of WT replication and the reliance of DIPs on the WT replication machinery.

Materials and methods
The kinetics of the corresponding biochemical reactions are described in the deterministic mathematical model introduced in the mathematical model section of this paper.The system of ordinary differential equations (ODEs) is formulated under the assumption of mass action kinetics, Michaelis-Menten approximations, and on the biological scheme presented in Fig 1 .The model can, in principle, be defined as a stochastic process.

Sensitivity analysis
The mathematical model included parameters which encode the biological mechanisms under investigation.Since many parameters required calibration, it is important to identify which have the greatest effect on model outputs.Global sensitivity analysis allowed us to evaluate the results of simultaneous changes in parameter values [22].For implementing this approach, consider the vector of parameters θ = (θ 1 , θ 2 , . . ., θ n ) such that the model output is described as Y = g(θ).We use the Sobol approach to determine global sensitivities [23].Each parameter θ i can be considered a random variable with an associated range.Since Y is a function of these variables, it is also a random variable with variance V (Y ).We were interested in the conditional variance V (Y |θ i = θ * i ).However, since the value of θ * i is not known, we instead considered the average conditional variance, , where the expectation is with respect to θ i , and the variance is taken over all remaining parameters θ j , j ̸ = i.The law of total probability gives from which the first order Sobol index for parameter θ i is defined as We also investigated the result of multiple fixed parameter values.If we let ) be the expected reduction in the variance by fixing all parameters except θ j , then the total effect of parameter θ i can be defined as

Model development and calibration
In formulating the new mathematical model, we introduced several additional parameters that relate to the kinetics of DIPs and the loss of non-structural proteins due to DIPs using trans-elements from WT virus for their replication.These variables are summarised in Table 1.Grebennikov et al. provided parameter estimates for the WT virus [21].These values are summarised in Table 2.
The remaining parameters were estimated using approximate Bayesian computation (ABC) [51].The ABC algorithm allows a user to define a set of prior beliefs about parameter distributions, π(θ), and combine this with model simulations and data to arrive at a posterior distribution π(θ|D).Given a sample parameter set, θ * ∼ π(θ), a user can simulate data D * ∼ π(D|θ * ) and compare them to the experimental data, D. If the simulated data are within a given threshold distance, ε (with distance measure d(•, •)), from the experimental data, D, then the sample parameter set (θ * , D * ) is accepted.Otherwise, the parameter set is rejected and this process is continued until N parameter sets are accepted [51].We made use of an Euclidean distance measure, defined as where T is the set of time points within the experimental data set, D, and M is the March 5, 2024 5/34 Data sets for defective interfering particles are limited, with little investigation of the intra-cellular replication kinetics of WT virus in the presence of DIPs.Chaturvedi et al. [15] investigated two SARS-CoV-2 DIPs as TIPs.Both DIPs had shorter genomes, around 6%-10% than the WT virus.Chaturvedi et al. [15] performed a virus yield-reduction assay by transfecting Vero cells with TIP or control RNAs (one µg/million cells) 24 hours before infection with SARS-CoV-2 at an MOI=0.05, and harvesting culture supernatants for titration at various time-points (24,48, or 72 hours post-infection).They discovered that these particles lead to a 1.5 − 1.2 log fold reduction in virus produced compared to control samples.We compared the fold reduction generated by therapeutic interfering particle two (TIP2) [15], at 24 and 48 hours, to the fold reduction from our mathematical model of [V wt released ] against the original model proposed by Grebennikov et al. [21].These fold reductions are summarised in Table 3.For the ABC rejection method, given that the choice of a suitable ε is difficult, we sampled 10 6 parameter sets and kept the top 0.1% that minimise the distance measure d(•, •).We assumed uniform prior distributions for the parameters within the search ranges summarised in Table 4.

Stochastic simulation algorithm
To perform stochastic simulations of the model formulated as a Markov chain, the most popular exact approach is Gillespie's direct method [52].The Markov chain specifies the propensity a m for the m-th jump process (i.e., the respective elementary reaction rate), which changes the variables by a discrete amount when that process takes place.The propensity a m defines the probability p m = a m dt that the m-th process is triggered in the infinitesimal time interval [t, t + dt).At each step of the simulation, two random numbers r 1 , r 2 ∼ U (0, 1) are generated to sample the time of the next jump process, τ , and the index r m of the process to perform: where a 0 = M m=1 a m is the total sum of propensities.
In this work, we used the rejection stochastic simulation algorithm (RSSA) [53].
This method estimates the upper and lower bounds on the propensities a m , a m , instead of calculating the exact values, a m , and uses the third random number, r 3 ∼ U (0, 1), for a rejection test to check if the exact value is needed to be computed (see details in Ref. [53]).The propensity values are updated only when necessary, therefore, the algorithm is practical when the propensity computation is time-consuming, e.g., for non-linear process rates parameterised with Michaelis-Menten functions.Additionally, two dependency graph structures are defined to reduce the number of propensity updates and accelerate computations: the first one specifies for each process which variables are affected when the corresponding jumps occur, and the second one specifies for each variable the process indices with propensities dependent on the value of the variable.Note that certain search strategies for the candidate process can also be  rate of loss of N SP s by trans elements 8.61 × 10 −3 10 [−2.9,−1.17]from positive sense WT RNA, h −1 Median estimates of unknown model parameters values for the model presented in Eqs. ( 4)-( 27), with variables defined in Table 1, as well as relevant search ranges.
March 5, 2024 8/34 jump processes by their propensity bounds).Alternatively, one can use approximate methods to significantly speed up computations, such as the tau-leaping method [52], or the other hybrid methods that make use of SDE or ODE approximations [53,54].In this work, however, we used the exact RSSA method, as the performance of the parallelised code to compute the ensemble of stochastic trajectories was acceptable.

Software
The

Mathematical model of WT and DIP infection
The variables of the mathematical model characterising the life cycle of SARS-CoV-2 according to Figure 1 are listed in Table 1.

Cell entry and RNA release
The binding of infectious WT virion to the cellular trans-membrane protein ACE2 allows entry and release of the viral RNA into the host cell.We describe this process by equations specifying the rates of change of free-, receptor-bound, and fused virions, as well as the viral RNA genome in the cytoplasm: Here [V wt f ree ] is the number of extra-cellular free infectious virions, [V wt bound ], the number of virions bound to ACE2 and activated by TMPRSS2, [V wt endosome ], the number of virions in endosomes, and [gRN A wt (+) ], the number of ss-positive sense genomic RNA.A similar set of equations is used to describe the cell entry and RNA release of March 5, 2024 9/34 non-infectious viral defective interfering particles: Here [V

RNA transcription and DIP parasitism
The where Eq. ( 12) reflects the fact that non-structural proteins are translated only from the viral genomic RNA of infectious WT virions.Transcription of negative sense viral genomic and subgenomic RNAs described by Eq. ( 13) and Eq. ( 14) is regulated by the positive Translation and competition for nucleocapsid protein and other structural proteins where and Assembly and release of WT SARS-CoV-2 and DIPs New virions are assembled at the endoplasmic reticulum-Golgi compartment, where N-RNA complexes become encapsulated.These assembled virions can then exit an infected cell by exocytosis via the lysosomal pathway, budding, or cell death [56,57].
There is no competition associated with the release of new infectious and DIP virions, but the viral assembly rates, θ wt assemb and θ dip assemb , depend on the availability of structural proteins, since DIPs will likely require fewer of them than WT virions.The rates of change of the ribonucleocapsid, assembled and released infectious SARS-CoV-2 March 5, 2024 11/34 and DIPs are described below: and In this study, we wish to explore model behaviour for different initial conditions, [V wt f ree ](0) and [V dip f ree ](0), and thus, understand the replication dynamics of WT viral particles in the presence of DIPs, and how the initial dose of WT or DIP particles regulates infection and production kinetics of WT virions.We are also interested in investigating the sensitivities to model parameters of different outputs.

Stochastic Markov chain model
The deterministic model defined by Eqs. ( 4)-( 27) can be generalised to a stochastic one formulated as a discrete-state continuous-time Markov chain (DSCT MC).The stochastic model allows one to account for integer-valued variables, to obtain probability distributions rather than mean field estimates for the variables of interest, and to compute the probabilities of productive cell infection at low MOI [58].It is convenient to estimate model parameters for the system of ODEs and then, with a calibrated deterministic system, and a defined Markov chain model, perform stochastic simulations making use of Monte Carlo methods.We follow our previous effort on the stochastic modelling of SARS-CoV-2 [58] and HIV-1 [54] life cycles to formulate and simulate the Markov chain.The Markov chain corresponding to Eqs. ( 4)-( 27) is presented in Table 5.
It includes the state transition events and the propensities, a m , for the m-th jump process.The propensity a m defines the probability p m = a m dt that the m-th process takes place in the infinitesimal time interval [t, t + dt).This definition yields exponential distributions for the time between jumps and various Monte Carlo methods can be used to simulate the stochastic trajectories from these distributions [52,53].We note that the processes of ribonucleocapsid formation (m = 33, 34) and virion assembly (m = 37, 38) are formulated as single events, yet involve the simultaneous change of three different variables.In these processes, protein copy numbers are decreased by the corresponding number of protein molecules, n p , needed to form a complex or assemble a pre-virion particle (i.e., by n wt N , n dip N , n wt SP , or n dip SP , respectively).Alternatively, one can formulate the Markov chain (MC) with three separate processes for each assembly event, in which the protein molecules are decreased by only one molecule with the propensity multiplied March 5, 2024 12/34 by n p (see Ref. [58] for an example of the extended MC formulation).We have verified that the extended and reduced MCs produce similar statistics.This reduction can be viewed as a weighted sampling strategy used in probability-weighted dynamic Monte Carlo method (PW-DMC) to accelerate computations [53].

Sensitivity analysis
We was identified as an important parameter to minimise variation in the release of both WT [V wt released ] and DIPs [V dip released ].Consequently, transcription of negative sense WT genomic RNAs is an essential first step in producing positive stranded gRNA, which is then translated to form structural proteins S, M , and E ([SP ]), as well as nucleocapsid proteins ([N ]) which are required to form new viral particles.Parameters associated with WT virion or DIP assembly are also important to monitor to reduce variation in model outputs.Several of the parameters identified by the Sobol sensitivity analysis have been previously estimated in Ref. [21] and are summarised in Table 2. Other parameters required estimation and these are listed in Table 4.

Parameter calibration
In our extension of the model proposed by Grebennikov et al. in Ref. [21], we introduced several parameters which have not been previously quantified.To estimate their values, we performed Bayesian parameter calibration.Since experimental data sets on co-infection with DIPs is limited, we aimed to achieve the fold reduction experimentally quantified by Chaturvedi et al. in Ref. [15].We made use of an ABC rejection algorithm with 10 6 sample sets.As previously mentioned, since a choice of ε is March 5, 2024 13/34 ] Entry and RNA release (DIPs): ] ORF1 translation and competitive viral RNA replication: Translation and ribonucleocapsid formation: Assembly and release:   hard to determine, we instead took the 0.1% of parameter sets which minimise the Euclidean distance.We sampled the exponent of the search ranges shown in Table 4.
As a result, our sample size provided a large coverage of parameter space.We compared the fold log reduction between the reference solution of a model without DIPs and the one with DIPS to the data in Table 3, and Figure 4 illustrates the model output where we used the median values from the accepted 0.1% sample sets.From these median values we obtained a fold change of 1.08 (two d.p.) at 24 hours post-infection and 1.14 (two d.p.) at 48 hours post-infection, compared to the reference solution [21].Posterior histograms in Figure S1 showed that with the data set and the mathematical model, Bayesian inference has led to poor learning for all but one of the newly introduced parameters.Posterior distributions are still extremely wide, with k wt trans(+) being the only parameter with a narrow posterior distribution.This was due to lack of longitudinal data to compare modelled DIP replication dynamics with.
Figure 5 illustrates the time evolution for each variable in Table 1 given the median values found via ABC rejection.From the upper panels of Figure 5 we examined that the entry kinetics of WT virus into the cell are similar to those of the reference solution.DIPs, however, enter the cell at a faster rate than WT virions.It is important to remember that we assumed there are sufficient ACE2 receptors mediating viral entry; thus, there is no competition between WT and DIP for receptor binding.The number of non-structural proteins is greatly reduced (Figure 5   Assembly and release parametric (mean values) and non-parametric (medians, inter-quartile ranges) statistics 306 computed on an ensemble of 10 6 trajectories.Additionally, the histograms of the 307 simulated variable values at particular time points can be produced from the ensemble 308 for analysis (Figure S2).The means and medians follow approximately the deterministic 309 model outputs, while the inter-quartile ranges estimate the uncertainty of the 310 simulations due to stochastic effects caused by the discrete nature of the model variables.311 These stochastic effects are more prominent for variables which are present in a cell in 312 small numbers, and the assumption that their mean values can be approximated by the 313 deterministic model may not be satisfied.In particular, the deterministic model predicts 314 that an infection is productive for every positive initial dose of [V wt f ree ](0) = 10, while the 315 stochastic trajectories can become extinct due to stochasticity. Figure 7 illustrates the probability of productive infection as a function of the initial WT virion (MOI) and DIP doses.As can be seen in Figure 7 (left panel), the probability of a productive infection tends to one as the initial dose of the WT virus hits 20 viral particles.
However, the probability is affected by the initial dose of DIP particles (Figure 7, right panel), with this probability being reduced linearly as the dose increases.
The mean values of WT virion and DIP production 24 hours post-infection closely follow the outputs predicted by the deterministic model.However, a part of the trajectories simulated with the stochastic model become extinct.The probability of productive infection as a function of WT M OI and DIP 0 is shown in Figure 7.When this probability is close to one, a change in DIP 0 does not significantly modify it.For every WT M OI, an increase in DIP dose reduces this probability linearly.Figure S5 shows the dependence of this linear decay in probability, β wt , on WT M OI.

Dose response analysis
We examined the release kinetics, i.e., the abundance of WT virions compared to DIPs, as a function of the initial doses [V wt f ree ](0) = [V dip f ree ](0) = 10.However, one can expect that initial infection doses might vary from cell to cell.Therefore, we now examined the release kinetics of WT virions under different initial conditions.Figure 8  from the stochastic model presented in Figure S3 (upper panel), while the mean estimates (Figure S3, lower panel) showed marginally higher release in WT virus and lower release of DIPs.Additionally, for high doses of WT virus, a productive infection is almost guaranteed (Figure 7), but as shown in Figure S3, even if an infection is guaranteed the overall number of WT, and hence infectious particles, is reduced.Figure 9 shows viral particle release kinetics predicted by the deterministic model with fixed initial conditions for [V wt f ree ](0) ranging from 3 to 20 virions and varying initial conditions for DIPs [V dip f ree ](0) = 1-100.DIP release peaks at a MOI = 6 and then begins to decrease as the dose increases.An increase in dose continues to have an effect on the release of WT virions, so that for a MOI = 40 total WT virion production is < 30 virions released over the 24 hour time period considered.This highlights the ability of DIPs to compete (with an advantage) for replication resources with WT virions.Consequently, if the dose is high enough, DIPs sequester so many intra-cellular resources that WT production is significantly reduced.Finally, the non-linear effects of DIP MOI on WT virion and DIP production per cell suggest that there might be optimal dosing of DIPs when used as a therapeutic agent.The maximum effect can potentially be achieved at around 5 to 10 DIPs per cell as this would maximise the number of new DIPs produced by the infected cells, and these, in turn, will reduce the WT virion production in other infected bystander cells.

DIP dose effect on WT virion production
Given the predicted three-dimensional curves of model outputs as function of initial doses presented as heatmaps in Figure 8, we asked if the production of WT virions 24 hours post-infection W T 24 = [V wt released ](24) as function of initial doses M OI = [V wt f ree ](0) and DIP 0 = [V dip f ree ](0) can be approximated by a compact analytic expression.Figure 9 shows that W T 24 as a function of DIP 0 , for a fixed M OI, exhibits a decay that is slower than exponential (which would be displayed as a straight line on a logarithmic scale).Therefore, we used several analytical expressions with slower than exponential decay to fit the deterministic model predictions for WT virion production for fixed M OI = 10.These include: σ(M OI).The overall parameterisation is the following: where W T 24 (DIP 0 = 0, M OI) is the number of released WT virions 24 hours post-infection with zero DIP initial dose for a given WT M OI.
Figure S4 shows the fit of a generalised Pareto function (28) for M OI = 10 and the dependence of parameters ξ and σ on M OI.One can see that the fit follows the data closely for suitable numbers of produced WT virions (W T 24 > 1) and has some small discrepancies for W T 24 < 1 at large DIP doses DIP 0 > 40.The dependences of parameters ξ(M OI) and 1/σ(M OI) exhibit non-linear patterns.They can be approximated with a Hill function and a Dagum distribution p.d.f., respectively.However, when these analytic approximations are substituted in (28), the overall fit of ( 28) behaves approximately as an exponential decay function (data not shown).
Therefore, one should use the computed estimates of the parameters ξ and σ for every M OI, or approximate them with a higher degree polynomial that would follow the estimates closely, e.g., with a 30-degree Chebyshev polynomial as shown in Figure S4.
Overall, the relative error of the fit (weighted residual sum of squares) of the closed-form expression (28) reaches its peak for M OI ≈ 6, in the same region where the parameters ξ and σ shown a non-linear dependence on M OI.The root-mean-square deviation normalised to W T 24 (DIP 0 = 0, M OI), the maximum value of produced WT virions for each MOI, shows a similar increase near M OI ≈ 10, as well as a later increase for large MOIs.This can be explained since the discrepancy in the tail of a generalised Pareto distribution corresponds to larger numbers of W T 24 with an increase of M OI.In summary, we have provided a closed-form expression, (28), as a prediction of the effect of DIPs on productive cell infection, i.e., the expected mean number of WT virions produced in a productive infection scenario for a range of relevant MOIs.

Discussion
SARS-CoV-2 still presents a real threat to human health as a result of several compounding factors: emergence of new strains due to mutation, waning immunity amongst the vaccinated, and un-vaccinated individuals (for medical reasons or personal choice).Therefore, it is still important to investigate new treatment options, especially those that could be implemented early after infection, to alleviate pressure on healthcare systems.One such potential therapy is defective interfering particles.DIPs are virus-like particles with shorter genomes that require a wild-type (WT) virus to replicate.In this paper, we investigated the intra-cellular replication kinetics of WT virus in the presence of DIPs, making use of a mathematical model.To this end, we extended the model proposed by Grebennikov et al. in Ref. [21], which focused on the intra-cellular replication kinetics of SARS-CoV-2, to include co-infection with defective interfering particles, given their therapeutic potential [59,60].In particular, we investigated the ability of DIPs to reduce WT viral load by competing for resources required to replicate or encapsulate the viral genome to form new virions.Since DIP genomes lack key fragments, they need a "helper" virus, which encodes non-structural and structural proteins, for their replication.There is evidence of DIPs leading to cause a reduction in viral titres for several viruses including: influenza A, dengue fever and SARS-CoV-2 [14,15,61].With the emergence of new SARS-CoV-2 strains, the effectiveness of a DIP particle (derived from a particular viral strain) against novel ones remains to be investigated.

Mathematical models of WT virus and DIP co-infection have been investigated at
March 5, 2024 23/34 the within host-level and either consider: a standard infection model with target, eclipse phase and infected cells or include different localised areas of infection, such as the upper and lower respiratory tract [15,19].There are, however, no models (to the best of our knowledge) that examine the intra-cellular replication kinetics of SARS-CoV-2 in the presence of DIPs.Our aim was to assess the hypothesis that DIPs lead to a reduction not only in the number of released WT virions but also, negatively impact the transcription of positive sense genomic RNAs.Additionally, we investigated the effects of initial infection dose (WT and DIP) in the release of both new WT virions and DIPs.Since experimental data sets are extremely limited, it is important to note that the parameter values obtained in this manuscript are based on the data set used [15], and may not be globally identifiable.By globally identifiable we mean the identification of a unique parameter value from a data set.
The extension of the model presented by wt/dip tr(+) , the transcription rates for positive sense genomic RNA.The rates associated with cell entry, k f use and k uncoat , also lead to some variation in model outputs.Finally, if we examine as output WT and DIP release, we find their associated assembly rates, k wt assembl and k dip assembl , as the most sensitive parameters.
DIPs have potential as therapeutics, thus, it is important to explore how initial infection doses of WT and DIP alter the release of WT virus, to inform a treatment plan.We show that even a low MOI= 1 of DIPs can cause a reduction of approximately 50% in released WT virus compared to an infection in the absence of DIPs, with further reduction in released WT up to 10-fold for increasing MOI wt and MOI dip .Figure 8 illustrates how increasing the dose of DIPs leads to a reduction in the fraction of released WT virions, in relation to the initial WT infection dose.These trends are consistent with the results from the stochastic model also developed in this paper (Figure S3).The doses of both WT virus and DIPs also had an effect on the probability of a productive infection, which decreased with increased doses of DIPs, but is almost certain for high doses of WT virus.We also investigated the effect of initial MOI of DIPs given a fixed dose of WT virus (MOI=10) on viral particle release.Our results show that while DIP release peaks at an initial DIP dose of MOI=5, the release of WT virions decreases in a dose-dependent manner.Furthermore, by an initial DIP dose of MOI=40, WT virion release is effectively inhibited.
The deterministic and stochastic models we presented are a good first approximation to the kinetics of WT and DIP co-infection.Yet, there are a number of biological processes which have not been considered.First and foremost, we omitted the anti-viral response of the cell.While we need not consider the adaptive immune response since our time interval is 48 hours, the innate immune response would play a pivotal role [62,63].A family of cytosolic receptors, known as pattern recognition receptors (PRR), exists that detect viral RNAs to induce the production of type I interferons.
Type I interferons (or viral IFNs), which are secreted by infected cells, include IFN-α, IFN-β, IFN-ω and IFN-τ .These molecules are associated with activation of anti-viral cell states, which in turn lead to inhibition of viral replication and eventual viral clearance [64].Furthermore, innate immune responses have been shown to be induced by DIP binding to PRRs, providing additional stimuli and magnifying the anti-viral cellular response [59].As a consequence, it would, therefore be ideal to extend the March 5, 2024 24/34 proposed model to consider the role of an innate immune response.Another limitation of our model is that for WT virions, we do not distinguish between infectious and non-infectious particles.This would be important to understand the potential infectivity of the viral particles released.We also fail to characterize the natural generation of DIPs during the WT replication cycle (which is inherently characterised by mutations).This process would contribute to the release of other defective interfering particles, and would potentially reduce the number of WT virions released.However, a complete calibration of such a model would require a data set not currently at hand.
To conclude, we believe the model we have proposed shows the potential benefits of DIPs as a therapeutic tool to reduce WT virus production.We also have shown that even low doses of these particles can have a positive effect on limiting WT virus production and reducing the probability of a successful infection.This reduction continues, in a dose dependent manner, to greatly reduce WT virus production.Future work will focus on incorporating immune responses and the natural production of DIPs into the mathematical model presented here but will require further carefully curated data to assist in parameter estimation.
Figure S5 The estimation of the effect of DIP dose on the probability of productive infection as a function of MOI.Left: the effect of DIP initial dose on the probability of productive infection.For each MOI, the probability decreases linearly as DIP dose, DIP 0 , increases.Right: the dependence of the linear decay rate, β wt , on MOI.The decay rate, and therefore the effect of DIP dose, decrease with increase in MOI.

Fig 1 .
Fig 1.  Biological scheme of the competitive replication of the infectious SARS-CoV-2 and defective interfering viral particles.There are three types of arrows shown on the scheme: (i) arrows without a T-end (which are not green) indicate the synthesis processes, i.e., translation or transcription, (ii) arrows with T-shaped beginning indicate that the variable at which the arrow points with the arrow-end is increased while the variable at which the arrow points with a T-end is decreased (e.g., transitions of the entities from one state to another, or a decrease of non-structural proteins during the transcription activation), (iii) green arrows indicate the transport of the entities from one place to another (to or from double-membrane vesicles or to the cell membrane in vesicles), which is not modelled explicitly.All entities are subject to degradation, however, these processes are not shown to avoid cluttering the figure.
[30,49] d dip N −gRN A degradation rate of DIP ribonucleoprotein, h −1 0.268 10 [−1.16,0] [32, 35] k dip release rate of DIP virion release via exocytosis, h −1 105 10 [0.9,3.15][49, 50] d dip assembled assembled DIP virion degradation rate, h −1 4.89 × 10 −3 10 [−4,−0.62][26] k wt trans (−) rate of loss of N SP s by trans elements 5.39 × 10 −5 10 [−5,−3.7]from negative sense WT RNA, h −1 k wt trans (+) rate of loss of N SP s by trans elements 6.17 × 10 −3 10 [−2.22,−1.86]from positive sense WT RNA, h −1 k dip trans (−) rate of loss of N SP s by trans elements 4.72 × 10 −5 10 [−5.69,−3] from negative sense WT RNA, h −1 k dip trans (+) released WT viral genomic RNA undergoes translation into non-structural viral polyproteins, [N SP ], which operate to form the viral replication and transcription complex, i.e., the RNA-dependent RNA polymerase (RdRp).The main function of the RdRp replication complex is to generate a negative sense full-length genomic and subgenomic RNAs.As DIPs lack the ability of self-replication, the conditional transcription of DIP RNAs results in competition with WT SARS-CoV-2 for replication proteins [55].The use of WT virus trans elements by DIPs reduces [N SP ] availability for the transcription of WT viral RNA.The respective sets of equations have different structures, as detailed below.The abundance of non-structural proteins, [N SP ], the negative sense genomic and subgenomic RNAs, [gRN A wt (−) ], and positive sense genomic and subgenomic RNAs, [gRN A wt ], associated with the replication of the infectious virions are described by the following equations:

Fig 4 .
Fig 4. Time evolution of viral particle release for the parameter values estimated with the ABC method.Viral particle release kinetics predicted by the model with initial conditions [V wt f ree ](0) = 10 and [V dip f ree ](0) = 10 48 hours post-infection after model calibration using ABC rejection and data from Chaturvedi et al. [15].Median parameter values summarised in Table 4 were used for previously unknown parameter values.(Yellow line:) shows the reference solution to a model where DIPs are not considered in the replication dynamics.(Red line:) illustrates the production of WT virions [V wt released ] with DIPs (blue line) [V dip released ].
middle left panel), peaking at 7 hours with ≈ 20 molecules as opposed to the reference solution, which peaks at roughly 13 hours with ≈ 40 molecules.The production of [gRN A wt (−) ] halves and peaks earlier in the time course, with a greater number of DIP negative sense genomic RNA than WT.Consequently, we then saw an approximate fold reduction of positive sense genomic RNA, ribonucleocapsid proteins, assembled and released WT viral particles.March 5, 2024

Fig 5 .
Fig 5. Deterministic model outputs.Time-dependent state variables of the mathematical model for the life cycle of SARS-CoV-2, including wild-type virions and defective interfering particles with initial conditions [V wt f ree ](0) = 10 and [V dip f ree ](0) = 10 over a 48 hour time course.(Upper left:) free WT or DIP virions bind and fuse to the cell ACE2 receptors, and (upper right:) virions entering endosomes and uncoating of viral positive sense genomic RNA.(Middle left:) transcription and translation to form a negative sense genome and ORF1 to form non-structural proteins (N SP s), which is then followed by (middle right:) the production of new positive sense genomic RNAs and translation of N protein.(Bottom left:) translation of structural proteins and formation of ribonucleocapsid molecules, which lead to (bottom right:) the assembly and release of new virions, both WT and DIP.(Dashed lines:) represent the reference model solution proposed in Ref. [21].

Fig 6 .Figure 6
Fig 6.Stochastic model outputs.Statistics of time-dependent state variables of the stochastic model with initial conditions [V wt f ree ](0) = 10 and [V dip f ree ](0) = 10 over a 24 hour time course based on an ensemble of 10 6 simulated trajectories.Solid lines: medians, dashed lines: mean values, filled areas: inter-quartile ranges.

34 Fig 8 .
Fig 8. Effects of varying initial dose on viral particle release.Top: Total WT virions released over the 24 hours post-infection for varying initial conditions of free WT virions [V wt f ree ](0) = 0-20 and free DIPs [V dip f ree ](0) = 0-20 from the deterministic model.Bottom: Total DIP particles released for varying initial doses.The isolines shown on the heatmaps as white lines coincide with the corresponding ticks in the colorbars.

Fig 9 .
Fig 9. Total WT virions and DIPs released for increased initial doses of DIPs.Viral particle release kinetics predicted by the model with fixed initial conditions for [V wt f ree ](0) ranging from 3 to 20 virions, as a function of varying DIP doses [V dip f ree ](0) = 1-100.
(a) a Gompertz curve, and the probability density functions (p.d.f.) of (b) power-law, (c) Weibull, (d) Cauchy, (e) Burr, (f) Lomax, and (g) generalised Pareto heavy-tailed distributions.The error that was minimized is the sum of squares between the WT virion production, W T 24 , predicted by an analytic expression and predicted by the deterministic model for each DIP 0 ranging from 1 to 100.The generalised Pareto distribution (with a location parameter equal to zero) was chosen as the optimal analytic expression making use of the Akaike information criterion.The parameters of the generalised Pareto distribution ξ (shape) and σ (scale) can be fitted for different values of M OI, thus giving the functions ξ(M OI)

Table 1 .
Dynamical variables of the mathematical model for the life cycle of SARS-CoV-2, with defective interfering particles.

Table 2 .
Estimates of previously calibrated model parameter values.
dip f ree ] is the number of extra-cellular free DIPs, [V dip bound ], the number of DIPs bound to ACE2 and activated by TMPRSS2, [V dip endosome ], the number of DIPs in endosomes, and [gRN A dip (+) ], the number of ss-positive sense genomic RNA.DIPs for SARS-CoV-2 would require a functional spike (S) protein to successfully bind to ACE2 receptors and mediate cell entry.Consequently, we assume that the rates for k bind , k diss , k f use , and k uncoat are the same for both WT virus and DIPs.However, degradation rates related to cell entry will differ between WT and DIPs, since the shorter genome of DIPs might imply a different degradation rate.
DIPs compete with WT virions for packaging proteins, e.g., nucleocapsid N proteins ([N ]) [55].Structural S, envelope E, and membrane M proteins are translated from positive sense subgenomic RNAs in the endoplasmic reticulum (ER) and are considered in the mathematical model as a single population, [SP ].Nucleocapsid proteins, on the other hand, are translated in cytosolic ribosomes from positive sense RNAs.Both SP and N proteins are required for the formation of virus like-particles, WT or DIPs.It can be assumed that n dip SP ≤ n wt SP and n dip N ≤ n wt N , since the shorter DIP genome will require fewer N proteins for the formation of the ribonucleocapsid and construction of a viral particle.Translation of N and SP proteins are described by the following two equations: now evaluate how model outputs change with parameter values.To that end a Sobol global sensitivity analysis was performed on four different model outputs.We first considered the variability of WT genomic RNA, [gRN A wt ], and DIP genomic RNA, [gRN A dip ], as a result of modifying parameter values within a set range summarised in Table 2 and Table 4. Secondly, we investigated how parameter variability affects the release kinetics of both WT [V wt released ] and DIP [V dip released ] particles 48 hours post-infection.Understanding which parameters cause the most variability in our model will allow us to calibrate it with careful consideration to minimise output perturbations.Figure 2 illustrates the first and total order sensitivities for WT genomic RNA, [gRN A wt ], and DIP genomic RNA, [gRN A dip ], as outputs of the proposed model.For [gRN A dip ], the parameter k dip tr(−) was identified as generating the largest variation.is associated with the transcription of negative sense RNAs for DIPs, and thus, is essential in the formation of new positive sense genomic and subgenomic RNAs.The rate k dip tr(+) was also identified as a high sensitivity parameter, since it is associated with the transcription of positive sense RNAs.Consequently, k wt tr(−) was the second most important parameter in minimising variation in model output for [gRN A wt ], following the same reasoning as for DIP positive sense genomic RNA.A parameter that was of great importance, and not only caused large variation in model outputs of [gRN A] for WT or DIPs, but also [V released ], was the threshold parameter of non-structural proteins, K N SP .K N SP causes the most variation for [V dip released ] and [gRN A wt ] compared to any other parameter, and for [V wt released ] and [gRN A dip ] it is the second most important parameter.K N SP is associated with the transcription of both negative and positive sense genomic RNAs, and changes in the value of this parameter will modify the number of WT virions and DIPs released.k wt tr(−)

Table 5 .
The Markov chain models: individual transitions and their propensities.
[21]ennikov et al. inRef.[21]requirednew parameters to account for the kinetics of DIPs.Therefore, it was necessary to investigate the sensitivities of all model parameters.In particular, we made use of the Sobol sensitivity analysis to understand how variation in parameter values affects four different model outputs: [gRN A wt ], [gRN A dip ], [V wt released ], and [V dip released ].We found several parameters that have an effect on all four model outputs: K N SP , the threshold number of non-structural proteins, k wt/dip tr(−) , transcription rates of negative sense genomic RNA for WT virus and DIPs, respectively, and k