Experimental and mathematical insights on the interactions between poliovirus and a defective interfering genome.

During replication, RNA viruses accumulate genome alterations, such as mutations and deletions. The interactions between individual variants can determine the fitness of the virus population and, thus, the outcome of infection. To investigate the effects of defective interfering genomes (DI) on wild-type (WT) poliovirus replication, we developed an ordinary differential equation model, which enables exploring the parameter space of the WT and DI competition. We also experimentally examined virus and DI replication kinetics during co-infection, and used these data to infer model parameters. Our model identifies, and our experimental measurements confirm, that the efficiencies of DI genome replication and encapsidation are two most critical parameters determining the outcome of WT replication. However, an equilibrium can be established which enables WT to replicate, albeit to reduced levels.

Introduction

Co-infections, the simultaneous infection of a host by multiple pathogen species, are frequently observed [1, 2]. The interactions between these microrganisms can determine the trajectories and outcomes of infection. Indeed, competition between pathogen species or strains is a major force driving the composition, dynamics and evolution of such populations [1, 3]. Three types of competition among free-living organisms have been defined from an ecological point of view: exploitation, apparent, and interference competition [1-3]. Exploitation competition is a passive process in which pathogens compete for access to host resources. Apparent competition occurs when the a population increases in the number of individuals. This, in turn, results in the increase in number of predators. In this scenario, other group of individuals can be indirectly affected by the the increase in the number of predators hunting both groups in the area. Thus, competition that is not due to using shared resources, but to having a predator in common [4]. In the context of the ecology of infection, apparent competition is associated with the stimulation of host immune responses, which acts as predator [3]. Interference competition represents a direct attack inhibiting the growth, reproduction or transmission of competitors, either chemically or mechanically [5].

In the current study we aim to understand better the interference competition process between two RNA virus genomes co-infecting a single cell. Thus, we examine the competition between full length wild-type (WT) poliovirus type 1 (PV1) and a replication competent replicon RNA derived from poliovirus type 1. This replicon was engineered by deleting the entire region encoding for capsid proteins. Poliovirus, the causative agent of poliomyelitis, is a positive-sense single-stranded RNA enterovirus belonging to the family Picornaviridae [6]. Upon cell infection, WT poliovirus initiates a series of processes that leads to the production of structural and nonstructural viral proteins and genome amplification. Structural proteins encapsidate the viral RNA, which leads to infectious viral particle production and cell-to-cell spread. Indeed, only encapsidated poliovirus genomes can survive outside the cells and can bind to new cells to initiate infections.

As an RNA virus, poliovirus is characterized by a high level of genome plasticity and evolution capacity, due to both high replication rate and error-prone nature of viral RNA polymerase [7, 8], which generates a large proportion of mutants in the viral population, called viral quasispecies [9-11]. In addition, defective genomes, carrying genomic RNA deletions, are thought to be produced also by RNA polymerase errors during RNA replication [12-14]. Natural poliovirus DI particles have been observed carrying in-frame deletions within the P1 region of the genome encoding structural capsid proteins [15, 16]. The DI genome used in this study features similar deletion of the P1 region, while expression of nonstructural proteins is not affected. When co-infecting cells with the WT, it can exploit capsid proteins produced by the WT to form DI particles, and in this way spread to new cells. However, it is not able to reproduce progeny in the absence WT helper virus [16].

From an ecological perspective, the WT can be viewed as a cooperator, producing capsid proteins as public goods. The DI particles are non-producing cheaters, bearing no production cost while exploiting the capsid protein products from the WT [1, 3]. Hence, co-infection should enable DI particles to replicate and spread, while hindering WT growth and propagation by interference competition.

DI particles have been identified in a number of viral species, such as vesicular stomatitis virus [17], poliovirus [18], ebola virus [19], dengue virus [20] or influenza virus [21, 22]. Given that DI particles compete and reduce WT production, they can modulate the outcome of infection [18, 22]. A recent study showed that defective viral genomes can contribute to attenuation in influenza virus infected patients [23], and they can protect from experimental challenges with a number of pathogenic respiratory viruses [24, 25].

Several studies have examined the interactions between defective particles and their parent viruses using a mathematical formalism. The majority of the prior work has focuses on the intercellular competition process in various settings (bioreactor production systems, passaging experiments, etc.) [26-28]. Much less attention has been devoted to the investigation of the intracellular dynamic of the WT-DI system. Stauffer Thompson and colleagues [29] investigate the impact of DIP dosing on the production of new DIP during a WT-DIP co-infection in Vesicular Stomatitis Virus. Differently from the model we present, their model does not explicitly take into account the competition for cellular resources. Laske and colleagues introduce a very detailed kinetic model of the dynamic of influenza with its DIs. [30]. While similar intracellular kinetic models exist for poliovirus [31], their application for the biological scope of this paper brings an unnecessary burden from the point of view of computation and parameter estimation.

In this study, we examined the competition between poliovirus WT and DI genomes within cells during one infection cycle, using a coarse grained dynamic model. DI particles are spontaneously generated during a significant number of virus infection, but the factors affecting DIs generation and propagation are not well understood. We then developed an ordinary differential equations (ODEs) mathematical model that captures the dynamics of DI and WT replication and encapsidation. We aimed to understand the processes within a single replication cycle that determine interference of the defective genome with WT within a cell. Our study indicates that DI and WT genomes compete for limiting cellular resources required for their genome amplification and for capsid proteins required for particle morphogenesis. The model was further used to evaluate the potential outcome of the interaction between DIs and WT viruses over a large range of parameter values and initial conditions (multiplicities of infection, temporal spacing and order of infection). Using DI variants carrying mutations that change their kinetic of replication, we further examine the predictive value of our model and its value to improving DI designs. Thus, the mathematical model described facilitate the understanding and interpretation of the biological impact of DI particles in the context of virus infection.

Results

Interference of WT poliovirus production by DI genomes

Initially, we evaluated whether DI genomes, carrying a deletion of the entire region encoding for capsid proteins, could affect WT virus replication (Fig 1A). The DI genome used in this study does not produce capsid proteins and is thus unable to encapsidate its genome and spread to other cells. However, it retains full capacity to produce nonstructural viral proteins and replicate its genomic RNA. WT poliovirus and DI genomic RNAs were transfected by electroporation to HeLaS3 cells and infectious titers of WT virus were determined over time by plaque assay (see Material and methods). As a control we used a replication-incompetent defective RNA lacking the capsid-encoding region, a part of 3D-polymerase encoding region, and the entire 3’ nontranslated region (NTR). HeLa S3 cells transfected by only WT genomes produced nearly 1 × 107 PFU/ml WT virus 9 hours after transfection, while co-transfection of WT genomes together with DI RNAs at a ratio WT:DI = 1:4 resulted in 100-fold decrease of WT titers (Fig 1B). The non-replicating defective RNA did not affect WT virus production, suggesting that replicating DI genomes are required for effective interference, as previously reported [32].

Fig 1

DI genomes inhibit WT virus production via genomic replication and encapsidation.

(A) Structure of the WT poliovirus, DI(Δ1175–2956), and DI(Δ1175–2956)(Rep-). The DI genome has an in-frame deletion in its P1-encoding region. The PvuII-cut DI genome is used for a non-replicating RNA control (Rep-). (B) Growth curves of WT poliovirus after transfection of WT genomes and/or DI genomes (sampling time-points: 0, 3, 6 and 9 hours after transfection). (C, D) Copy numbers of total (C) or encapsidated (D) genomes over time after transfection (sampling time-points: 0, 2, 3.5, 5, 7 and 9 hours post transfection). Blue and red squares indicate copy numbers of WT and DI genomes, respectively. Solid lines indicate single-transfected (WT or DI genome-transfected) samples, while dotted lines indicate double-transfected (both WT and DI genomes-transfected) samples. n = 3, mean ± standard deviations (SD) (B-D).

Quantification of the copy number of WT and DI genomes following co-transfection

Next, we examined the interaction between DI and WT genomes by varying the ratio of each RNA used to initiate transfection. As DI RNA demonstrated an increase in transfection efficiency, probably due to its shorter genome length, we adjusted DI and WT RNA concentrations to achieve equal number of transfected cells. Given that DI genomes are ∼2,000 nucleotide shorter (∼1/4 shorter) than WT genomes, it is expected to replicate faster than WT, and thus the copy number of DI genomes are higher than that of WT genomes. Transfection of shorter genomes is also more efficient than larger RNAs. We thus optimized our protocol to deliver equal copy number of DI and WT genomes into the transfected cells. We transfected 5μg of WT to 1.25μg of DI genomes, and we collected RNA samples at given timepoints (t = 0, 2, 3.5, 5, 7 and 9 hours after transfection). The average number of genomes in a single cell was estimated as the total number of genomes divided by the number of transfected cells (successfully transfected cells were 25.7%). Replication of WT and DI decreased ∼7 hours after co-transfection, but this effect was not observed in the cells transfected only with WT or DI (Fig 1C). Thus, replication of WT genomes was inhibited by DI, and the number of accumulated DI genomes also decreased in the presence of WT. This suggests that WT and DI genomes compete with each other for one or more limiting factors required for replication. Nonetheless, DI genomes have an small advantage and accumulated at slightly higher level than WT genomes (Fig 1C). To determine the numbers of encapsidated WT and DI genomes, we treated supernatants of infected cells with a mixture of RNase A and RNase T1. Viral RNAs encapsidated in virus particles are resistant to RNase treatment, while naked RNAs are degraded. Therefore, the number of RNase-resistant genomes is a measurement of particles containing WT and DI RNA. Of note, after the RNase treatment, DI genomes were not detected in samples obtained from cells transfected with DI RNA alone (Fig 1D), indicating, as expected, that without trans-encapsidation by capsid proteins provided by WT DI RNA is sussceptible to RNase treatment. The decrease of encapsidated WT genomes between singly and co-transfected cells conditions was two-fold larger than that of WT genomes without RNase-treatment, consistent with the idea that DI genomes hamper WT genome encapsidation (Fig 1C and 1D, compare the difference between plain and dashed blue lines at 9 hours after transfection in Fig 1D to the difference in Fig 1C).

These results support the idea that DI RNAs replicate faster than WT genomes, due to their shorter genome, consistent with previous observations [33, 34]. Interestingly, co-transfection results in a reduction in replication of both WT and DI genomes most likely due to competition for some host-cell limiting factor needed for genome amplification. In addition, capsid proteins produced by WT genomes limit DI and WT virus production as DI genomes compete for these proteins and thus further inhibit WT production. To further examine the mechanism of defective interference and quantitatively evaluate the effect of co-replicating DIs, we designed a simple mathematical model that describes the DI/WT genome interactions.

Mathematical description of intracellular competition

A deterministic mathematical model describing the intracellular competition between DI and WT genomes was developed, adapted from an existing competition model for human immunodeficiency virus (HIV) [35]. In order to describe appropriately the intracellular dynamics of poliovirus, we explicitly account for limiting resources depleted by the virus during replication, slowing down the growth of the population over time. This slowdown was experimentally shown by [36], who reported an exponential growth of viral RNA up to the third hour of infection, followed by a linear increase and then a plateau. The effect of limiting resources on poliovirus replication has been investigated by [37] using a mathematical model to explain the observed saturation in viral replication dynamics. Resources depleted by the virus may include phospholipids and de novo synthesized membranes for the formation of replication organelles [38, 39], pro-virus host factors required for virus translation, replication or encapsidation [40], or even the supply of amino acids for protein synthesis or nucleotides for genome replication [37]. Our model considers a generic set of resources (R) at the virus disposal during replication. It describes the changes, over the course of infection of a cell, of the numbers of WT (G) and DI (G) RNA genome copies, of free capsids produced by the WT (C) and of limiting resource units (R) depleted by the genomes for their replication. The set of ODEs is the following:

Model parameters are summarized in Table 1 and can be described through three distinct stages of the viral cycle: replication, capsid synthesis and encapsidation. A flow diagram of the model is available in Fig 2. Limiting resources are produced at a constant linear rate λ and captured by WT (G) and DI (G) genomes at a rate ε per unit of resource (uor) per minute. One uor and one viral genome, by definition, form one replication complex (c = 1 uor ⋅ genome−1, [41]). Conditionally on the capture of a resource unit, a WT genome replicates and turns into θ genomes, before the replication complex disintegrates (we set the condition θ > 1 in order for virus genomes to replicate). We assume that this happens quickly compared to the other processes. The replication rate of WT genomes is given by the product θεR(t). The replication rate of DI genomes is faster than that of WT genomes by a fixed factor P (P > 1), which likely depends of their shorter genome size requiring less time for the polymerase to complete a round of template copying [33, 34]. Because DI genomes lack the genes responsible for capsid proteins synthesis, only WT genomes are capable of producing free capsids (C), with the capsid-to-genome accumulation ratio η. WT genomes are then encapsidated (i.e. packaged into free capsids) at rate κ. DI genomes are set to encapsidate faster by a fixed factor ω (ω > 1), based on experimental observations from Fig 1C and 1D. One viral genome encapsidates into one capsid to form a virion (c = 1 genome ⋅ capsid−1). Finally, α, β, and γ represent the decay rates of, respectively, viral genomes, free capsids and resources.

Table 1

Notations used in the model and model parameters.

Notation	Definition	Unit a
	Observed variables
gWTtot , gDItot	WT, DI total genome number	gen
vWT, vDI	WT, DI virion number	gen
	State variables
GWT, GDI	WT, DI genome number	gen
C	Free capsid number	caps
R	Resource units	uor
CWT, CDI	WT, DI virion number	vir
	Model parameters		Best value c	Confidence Interval d	Sensitivity analysis e
c g	Number of viral genomes per capsid	gen ⋅ caps−1	1b
c r	Number of uor per viral genome	uor ⋅ gen−1	1b
θ	Genome replication factor	-	1.121	[1.000; 26.983]	[1, 1.25, 1.5, 1.75, 2]
ε	Resource capture rate by viral genomes	(uor ⋅ min)−1	2.021 ⋅ 10−6	[1.355; 3.000]×10−6	[0.8, 1.3, 1.8, 2.3, 2.8]×10−6
κ	Encapsidation rate	(gen ⋅ min)−1	5.097 ⋅ 10−6	[4.900; 5.282]×10−6	[2, 3.5, 5, 6.5, 8]×10−6
α	Viral genome decay rate	min−1	0	-	0
η	Capsid to genome accumulation ratio	caps ⋅ gen−1	5.260 ⋅ 10−2	[5.058; 5.439]×10−2	[2, 3.5, 5, 6.5, 8]×10−2
β	Capsid decay rate	min−1	2.589 ⋅ 10−3	[1.659; 3.476]×10−3	[1, 1.75, 2.5, 3.25, 4]×10−3
P	DI-to-WT replication ratio	-	1.075	[1.062; 1.087]	[1, 1.15, 1.3, 1.45, 1.6]
ω	DI-to-WT encapsidation ratio	-	2.185	[1.943; 2.453]	[1, 1.75, 2.5, 3.25, 4]
λ	Resource production rate	uor ⋅ min−1	10.008	[0.497; 88.924]	[1, 1.25, 1.5, 1.75, 2]×10−3
γ	Resource decay rate	min−1	7.471 ⋅ 10−4	[7.426; 88.582]×10−4	[3.5, 5.25, 7, 8.75, 10.5]×10−8
L	Logistic’s maximum	min−1	3.078 ⋅ 10−2	[3.061; 3.099]×10−2	-
s	Logistic’s steepness	min−1	3.234 ⋅ 10−2	[3.234; 3.234]×10−2	-
t 0	Logistic’s midpoint	min	318.459	[318.459; 318.459]	-

a gen: genomes; caps: capsids; uor: units of limiting resource; vir: virions; ‘-’: dimensionless

b Fixed parameter values

c Best optimized value shared by the 85 first sets (identical until fourth significant digits) for reduced model optimization (P, ω, κ, η, α, β, L, s, t0), and arbitrarily chosen best set for full model optimization (ε, θ, λ, γ). These best values are used in Fig 3 and S2 Fig, and S1 Table (model ).

d Range of variation over 150 best sets of parameter values

e The set of best values estimated from the full model and taken for the sensitivity analysis and all of the theoretical modelling work (i.e. Figs 4, 5 and 6, S3 and S5 Figs) was: θ = 1.192, ε = 1.814 ⋅ 10−6(uor ⋅ min)−1, λ = 1.000 ⋅ 10−3uor ⋅ min−1 and γ = 7.157 ⋅ 10−8min−1.

Fig 2

Flow diagram of the model (Eqs (1)–(6)).

State variables: resource units number R; wild-type (WT) poliovirus naked genome number G; defective interfering (DI) naked genome number G; WT-produced capsid proteins number C; encapsidated WT genome number C; encapsidated DI genome number C. Model parameters are defined in Table 1. Color code is blue for WT, red for DI and brown for resources. Segments represent genomes and diamonds represent capids.

Fig 3

Experimental data and model fit for the naked and encapsidated wild-type (WT) and defective interfering (DI) genomes.

Evolution of the number of WT and DI (A-B) naked genome copies and (C-D) encapsidated genome copies with time, from 2 to 9 hours post transfection (hpt). (A & C) show data in dually infected cells whereas (B & D) show data in singly infected cells. WT and DI results are shown in blue and red color, respectively. Crosses indicate experimental data for 3 replicates per sampling time at 2, 3.5, 5, 7 and 9 hpt (data drawn from Fig 1C and 1D, see Material and methods for details). Dashed curves represent the fit of the reduced model with logistic equation and solid curves show the fit of the full model (Eqs (1)–(6)).

Fig 4

Sensitivity indices for the proportion of WT virions at cell lysis (Φ) and the impact of the most important parameters.

A: Main, total and most important interaction indices for Φ. Main indices correspond to single parameter effect (black parts), and total indices add the effect of the factor in interaction with all other factors (second-order interactions, full bars). Hatched areas represent the importance of the strongest pairwise interaction between defective interfering (DI) to WT encapsidation ratio ω and DI to WT replication ratio P. B: Heat map representing the impact of P and ω on Φ, all other parameters being fixed to their best estimated value. The cross and lines in (B) show estimated parameter values from fitting the model to the experimental data.

Fig 5

Investigating DI mutant genomes.

A: PV3 and PV1-DI particles were serially passaged 8 times. B: Titers of PV3 (PFU/mL) and PV1-DI particles (IU/mL) included in the samples were determined. C: Mutations in the DI genomes positively selected over passages. D: One-step replication cycle. Dynamics of the ratio of PV1-DI to (wild-type) WT PV1 naked genome copies (total RNA) over time (3, 4.33, 6.33 and 8.17 hours post transfection), for distinct DI mutants (Ori: original non-mutant DI genome). E: Estimation of the ratio of PV1-DI to WT PV1 replication rate P from a least square regression between experimental and model predicted PV1-DI to WT PV1 naked genome ratios. Confidence intervals represent values of P at 10% of the minimum least square value. Model predictions of P against experimental DI to WT naked genome ratios at the last sampling time-point.

Fig 6

Model-predicted impact of infection delay and initial multiplicity of infection (MOI) on WT and DI burst sizes.

A: WT (blue) and DI (red) predicted burst sizes (encapsidated genomes at 9 hours post infection) for various delays in DI infection of a cell. Negative delays correspond to cases where DI infects first and positive delays to cases where WT infects first. The dashed gray line marks the case of simultaneous infection by WT and DI. B-C: Heat maps representing the predicted relative burst sizes of (B) WT, , and (C) DI, , as a function of initial WT (x-axis) and DI (y-axis) MOIs (input). WT burst sizes were normalized by WT burst size resulting from WT:DI = 1:0 MOIs and DI burst sizes by DI burst size resulting from WT:DI = 1:1 MOIs (black squares).

The number of encapsidated WT genomes (i.e. WT virions, C) and DI genomes (i.e. DI virions, C) were measured experimentally (Fig 1D) and can easily be derived from Eq 2 as the loss of free capsids due to encapsidation:

Further, burst sizes, which are defined as the number of virions at cell lysis, i.e. 9 hours post transfection (hpt, [42, 43]), can be written as for WT virions and for DI virions.

The system of Eqs (1)–(6) encompasses 10 parameters and six variables, among which only four variables could be experimentally measured. Therefore, the mathematical model presents a classical problem of parameter identifiability, specifically regarding parameters ε and θ that appear as a product, with parameter ε figuring separately in Eq (4) corresponding to the unmeasured variable R. To solve this problem, we built a reduced model by assuming that the decrease in resources due to viral uptake follows a logistic decreasing function. We use a two-step optimization procedure to estimate the parameters. Briefly, we first estimate a set of parameters using the reduced model. Then, we use the parameters that are common between the reduced and the full models as “best guess” and run the estimation on the full model. Using this procedure, we are able to estimate the model parameters with reasonably high confidence (see Material and methods and Table 1).

We additionally conducted a model selection procedure in order to validate that the features introduced in our model were improving fitting statistics chosen as (i) the R-squared between experimental and fitted data, and (ii) the log-likelihood and (iii) Akaike information criterion (AIC) of a linear model explaining the experimental data with the fitted data. The different versions of the model that were tested are detailed in S1 Text and the results are presented in S1 Table. In what follows we will present results from both the reduced and full models for comparison. We report a full comparison of the models in the S1 Text.

Model predictions fit well the experimental measurements in Fig 1C and 1D, with best R2 = 0.974 for the reduced model and 0.964 for the full model (Fig 3). In dually transfected cells, both versions of the model reproduce well the number of naked and encapsidated genome copies. In singly transfected cells, the reduced model somewhat underestimates the number of naked genome copies while fitting satisfactorily the number of encapsidated genome copies. Conversely, the full model describes well the number of naked genome copies but overestimates the number of encapsidated genome copies. Overall, either model is able to describe the reduction on genome replication when WT and DI are co-transfected as compared to singly transfected cells (compare same color plain or dashed curves between Fig 3A and 3B). We hypothesized that this effect is the consequence of competing for limiting resources necessary for replication. The model also shows that the WT is more severally affected by competition for resources than the DI, as observed experimentally (Fig 3A, compare red and blue plain curves) thanks to the higher replication rate of DI genomes (P = 1.075, Table 1). Model predictions are also consistent with the experimental observation that DI genomes are more efficiently encapsidated than WT genomic RNA (Fig 3C, compare red and blue plain curves). This is due to the higher encapsidation rate of DI genomes (ω = 2.185, Table 1). Indeed, if we introduced in the model the same encapsidation rate for WT and DI genomes (i.e. ω = 1, model in S1 Text, eqs. (S10)–(S13)) we observed an underestimation of DI encapsidated genomes as compare with the experimental measurements. We conclude that DI genomes are encapsidated more efficiently than WT RNA. It is possible that the faster RNA replication of DI genomes provide an advantage on the competition for capsid proteins during encapsidation (S1 Fig). Thus, the model describes the experimental observation that WT genomes are encapsidated at a decreased rate in dual transfection compared to single WT transfections (compare Fig 3C and 3D, blue curves).

All parameter estimates are narrowly defined by the fitting procedure (Table 1) except for θ, λ and γ as they tend to correlate with each other, yielding non-uniqueness of best-fit parameter values (S2(A), S2(C) and S2(D) Fig, see also the distribution of ε values in S2B Fig). Also, a strong log-to-log relationship was found between the resource production to decay ratio λ/γ and the replication factor θ (S2(E) Fig).

Both the reduced and full models feature a time-dependent virus replication rate (S2(F) Fig). In the reduced model, this is given by the logistic function Λ(t) (Eq 10 in Material and Methods), and in the full model by the product θεR(t). In both models, the best fit yields approximately the same time-dependent replication rate, starting at 3.07 ⋅ 10−2 min−1 for the reduced model and at 3.04 ⋅ 10−2 min−1 (θελ/γ) for the full model, and decreasing with time towards 0. According to the reduced model, the time of half-decay is 318 minutes (t0 in Table 1), which corresponds to 5.3 hours post transfection.

The predictive power of the full model with best estimated parameter values was tested on independent experimental measurements of WT burst size corresponding to various initial DI-to-WT ratios for which the model had not been trained (S3 Fig). Relative experimental and predicted WT burst sizes were normalized by their respective value for WT-only transfection. The model was able to predict experimental outputs well, albeit underestimating WT output for some DI-to-WT input ratios. The largest underestimation was observed for the DI-to-WT input ratio of 0.5, and this discrepancy vanished as the input ratio increased.

Model predictions

The aim of our work is to understand the competition dynamics of WT and DI genomes during co-infection. To achieve this goal, we used the model described above to study how changes in parameter values around their experimental estimates can impact the outcome of the competition. Additionally, we also used the model to evaluate the effect of initial infection conditions, such as initial genome copy numbers of WT and DI and a time delay of cell infection on their respective burst sizes.

Sensitivity analysis

A sensitivity analysis was performed to identify parameters that have a significant effect on the output variable of interest, which we set as the proportion of WT virions at the time of cell lysis, Φ (Eq 12). Parameters were varied by ±50% of a set of best fit value based on experimental data, with five equally spaced values for each parameter (Table 1). The decay rate of genomes was not varied, because it was estimated to be negligible. The results indicate that the DI-to-WT replication ratio, P, and the DI-to-WT encapsidation ratio, ω, as well as their second-order interaction, were the most influential factors for the variation of Φ, explaining 76%, 17% and 7% of the variance, respectively (Fig 4A). All the remaining factors and their second-order interactions had a negligible effect (less than 1% of the variance). Hence, our model predicts that only parameters associated with the DI design have a strong impact on the degree of suppression of WT by DI.

To further examine the effect of P and ω, we varied both parameters from their best estimated value (Fig 4B). As expected from the global sensitivity analysis, P was found to be more important than ω for the production of WT virions, the gradient of Φ being steeper along P-axis than along ω-axis. Within the tested range of parameters P and ω, the value of Φ varied between 2% and 50%. The reference value of Φ corresponding to best-fit parameter estimates from experimental data was 23%. Therefore, we can predict that a DI particle with a lower replication rate, or, to a lesser extent, with a lower encapsidation rate than the DI particle used in the present work would weaken its competitivity with the WT virus, potentially leading to an increase in the proportion of WT virions at cell lysis of up to 27%. Conversely, a DI particle characterized by a higher replication rate or a higher encapsidation rate would strengthen its competitivity, potentially leading to a decrease in Φ of up to 21%.

Validation of model predictions using DI variants

To evaluate predictions of the model, we next investigated the competition between WT virus and DI variants. Our sensitivity analysis provides a quantitative measurement of the WT/DI competition as a function of the replication advantage of the DIs (DI-to-WT replication ratio P). As assessed in the previous Section, the most significant determinant for the outcome of infection is the ratio between DI and WT replication, P. Moreover, the simplest model that explains the difference between these mutants and the WT assumes no difference in encapsidation efficiency. The only nostructural protein involved in encapsidation is 2C, and none of the mutations we consider occur in that protein. To validate these predictions we isolated DI variants with defined replication advantage. In this experiment, we used WT poliovirus type 3 (PV3) and PV1-DI to discriminate their genomes by sequencing. Briefly, PV1-DI and WT PV3 were serially passaged for 8 times (Fig 5A and S1 Text). We determined PV3 and DIP titers at each passage to gain information on the dynamics of the co-infection experiments (Fig 5B). Both DI and PV3 compete with each other, resulting in a 10-folds reduction on PV3 titer, and the establishment of an equilibrium in which both PV3 and DIs are maintained for 8 passages (Fig 5B, red and green squares).

These experiments also serve to determine whether mutations accumulate during co-passaging in either DI or PV3 genomes. This information was used to construct DI particles with defined replication advantage with respect to the original DI (DI Ori). No change in encapsidation efficiency is expected. The mutation rates in RNA viruses can range from 10−4 to 10−6 per base, while conventional next-generation sequencing (NGS) can only detect variants at a frequency of about 1 in 100–500 (due to sequencing errors). Consequently, many variants in a viral population which exist at low frequencies cannot be distinguished from noise using conventional NGS. Therefore, we used Circular Sequencing (CirSeq), a technique developed to improve accuracy of NGS [44, 45]. We engineered several of the mutations identified into the DI cDNA, which facilitates identification of rare variants detection in the population. We identified eight mutations positively selected over passages (Fig 5C).

Next, we engineered DI variants carrying the identified mutations and determined the outcome of competition between WT virus and DI variants. HeLaS3 cells were co-transfected with WT poliovirus and DI variants and the ratio of DI-to-WT genome copies was estimated by RT-PCR over the course of infection (Fig 5D and S1 Text). In this way we were able to estimate genome replication P using the model, with all other parameters fixed to a set of estimated best values from the main experiment (Table 1 and S1 Text), and we established that the model prediction P and the experimental determined WT/DI ratios correlate with a Pearson coefficient of 0.96 (Fig 5E). Thus, these data confirms that, as predicted by our sensitivity analysis, the relative DI/WT replication ratios predict the outcome of infection. Thus our model serves as a useful tool to identify the most important determinants on the DI/WT interaction and enable to quantitatively estimate key parameters controlling modulation of the WT infection.

Effect of the multiplicity of infection and the timing of co-infection on WT burst size

We use our model to investigate the impact of varying initial conditions including (i) the time difference between WT and DI infection of the cell and (ii) the initial quantities of WT and DI on the WT and DI burst sizes, and , respectively. First, the time for DI infection compared to WT was varied from -7 to +7 hours post WT virus infection (Fig 6A). Negative delay values indicate that DI infects first, while positive values indicate that WT infects first. WT was allowed to produce at least one virion at cell lysis when DI was infecting the cell no more than 1.37 hours prior to the WT virus (delay = −1.37 hours). From this delay value, WT burst size increased as a steep logistic function, reaching a plateau at virions from a delay for DI infection of around 4 hours post WT infection. The curve of DI burst size was bell-shaped, reaching a maximum of virions at a delay for DI infection of 0.4 hours post WT infection. The delay window allowing DI genome to be encapsidated and produce at least 1 virion is narrow, from -3.43 to 4.37 hours. Most importantly, DI burst size superseded WT burst size only until the delay of 0.62 hours. The difference between DI and WT burst sizes was the largest when DI infected the cell 0.03 hours prior to WT ( virions).

We then investigated the impact of varying initial WT and DI multiplicities of infection (MOIs, corresponding to the number of viral genomes successfully entering a cell and initiating infection) from 0 to 1000 on WT and DI relative burst sizes (Fig 6B and 6C). WT burst sizes were normalized by WT burst size obtained for WT:DI = 1:0 MOIs, while DI burst sizes were normalized by DI burst size obtained for WT:DI = 1:1 MOIs. WT relative burst size varied from 0 to 1.11 depending on MOIs. DI relative burst size ranged from 0 to 1.26.

For WT virions to be completely suppressed, DI must infect the cell 2 hours prior to WT infection or at higher MOI than the WT, complete suppression taking place at co-infection of 100 DI and 1 WT. Conversely, the optimal initial conditions for maximizing the DI burst size arrive when WT is initially present in a slightly larger quantity than DI. The DI needs enough WT to exploit its capsids and produce virions. Because (i) the DI replicates and encapsidates faster than the WT, (ii) only the WT produces free capsids and (iii) replication and capsid production result in resource depletion, it is better for DI virion production to have the WT infect a cell in a slightly higher quantity than the DI. In that case, the WT has a slight initial advantage over the DI and can use resources to produce free capsids. In return, the DI can exploit those free capsids at its own advantage as it replicates and encapsidates more efficiently.

We also examined the cross-effect of the time delay and the variation of initial MOIs (S5(A) Fig). Globally, the shapes of WT and DI burst size curves as a function of delay are very similar between different WT and DI MOIs. The observed effect is a shift of the curves on the delay axis. Equal WT and DI MOIs generate very similar WT and DI burst size curves as a function of delay time. At these equal MOIs, DI competes more efficiently than WT upon co-infection, i.e. simultaneous infection of the cell by DI and WT (S5(A) and S5(B) Fig) and maximal DI burst size occurs when DI infects the cell around 0.3–0.5 hours after the WT (S5(A) Fig). As WT MOI gets larger than DI MOI, both the delay time that maximizes the DI to WT burst size difference and the peak of DI burst size shift towards delays where DI infects the cell before the WT (S5(A), S5(B) and S5(C) Fig). Giving an initial MOI advantage to the WT results in the DI having to infect the cell sooner in order to maximize production. On the opposite, as DI MOI gets larger than WT MOI, the shift is towards delays where WT infects the cell before the DI. Giving the DI an additional advantage in terms of MOI than it already has by being faster at replicating and encapsidating can be harmful for DI production, but this can be compensated by an earlier infection of the WT.

Discussion

We used a combination of mathematical modeling and empirical measurements to better understand of the interaction between a cooperator, the WT, producing capsid proteins as public goods, and a cheater, the DI, exploiting those capsid proteins from the WT. This type of direct competition affecting the growth of a competitor represents a case of interference competition. In addition, in co-infected cells, we also observed that WT and DI compete for shared resources necessary for replication, a phenomenon known as exploitation competition [1-3]. Hence our study identifies two different types of competition occurring during the co-infection of a cell by WT and DI.

Huang & Baltimore [46] argue that defective particles may influence the outcome viral infections. DI RNAs are produced by RNA replication errors, and are maintain in the virus population by high-multiplicity of infection, where DI and WT virus co-infected cells [47]. Internal sequences of the original vRNA of DI RNA segments are deleted, but cis-acting elements are retained, which enables DI replication. Naturally occurring immunostimulatory defective viral genomes (iDVGs) are generated during respiratory syncytial virus (RSV) replication and are strong inducers of the innate/natural antiviral immune response to RSV in mice and humans [25].

Mechanism of defective interference

Our study provides evidence for three main conclusions. First, the number of WT or DI genomes, taken individually, is lower in dually infected cells compared to singly infected cells (Fig 1C), indicating that they compete for a limiting resource for replication. Second, in dually infected cells, DI genomes replicate faster than WT genomes (Fig 1C), showing the advantage of their shorter genome size [33, 34]. Furthermore, we identified DI variants with increase replication fitness to confirm the model prediction indicating that the WT/DI replication ratio is a major determinant for the outcome of infection. Third, our study also demonstrated WT genomes encapsidation is inhibited by co-infecting DIs (Fig 1D), which further decrease WT virions production.

We have designed a minimal mathematical model able to capture key features of the DI/WT interaction during a single-cell replication cycle. We accounted explicitly for depletion of cellular resources and available capsid proteins produced by the WT virus. The model accurately describes the experimental data, and to predict new data on which the model had not been trained (S3 Fig). Fitting the available experimental data enables estimating parameters within biologically realistic ranges providing a better how the interaction of these processes affect WT/DI interaction during co-infection.

If DI’s replicate faster than WT simply because its shorter genome, the ratio of WT to DI genome replication should correspond to the genome lengths, that is 7515bp/5733bp = 1.311. However, the best optimized value of the corresponding parameter P was lower (1.075, Table 1). This underlines the complexity of the replication process and the value to develop reliable model to examine these type of processes. For example, poliovirus replication is step-wise process in which formation of replication complexes, follow for synthesis of negative- and positive-stranded viral RNAs may affect the overall kinetics of the process [48]. Thus, each of these steps might lower the average difference in replication rate between the WT and DI genomes over the course of infection of a cell. Indeed, when our model optimization procedure to the reduced model was performed keeping P fixed to the ratio in genome lengths, we observed a worse fit to the data (S6 Fig).

Our analysis suggests that some limiting resources are required and depleted during replication of viral genomes. Those resources might be lipids recruited by the viral machinery for the formation of replication complexes [39, 40] or host proteins involved in in different steps of viral replication [37]. When both WT and DI co-infect a cell, they are in competition for the exploitation of those shared and limiting resources. As DI replicates faster, it depletes resources faster than WT, affecting WT replication compared to WT-only infections [3]. Additionally, virus replication may represent a disruption for the cell physiology, which may lead depletion of resources at later stages of the replication cycle [37]. Consistent with this idea, fitting the models to the experimental data we estimated an initial replication rate of 3.04 ⋅ 10−2 − 3.07 ⋅ 10−2min−1, representing initial replication steps where the resources are not limiting. As resources are depleted, the model predicts that replication rate will tend to 0 following a logistic function (S2(F) Fig).

Insights on model fit to experimental data

The reduced model underestimates the number of viral genome copies in singly infected cells. The logistic equation considers the time-dependent genomic replication to be the same in both dually and singly infected cells. Under this assumption depletion of resources required for replication is identical in singly or dually-infected cells. As a result, genome copies are well predicted in dually-infected cells but underestimated in singly-infected cells, given that resource consumption is lower in single-infected cells. In contrast, the full model is able to recapitulate the impact of resource depletion on RNA genome production in both dually and singly infected cells, as a specific variable for resources and mass action terms are considered within the full model.

However, the full model overestimates the number of encapsidated WT genomes in singly infected cells. This result suggests that the WT virus encapsidates less efficiently when it is alone than in DI co-infected cells. Our model assumes that the encapsidation rate of WT is the same in singly- and dually-infected cells (κ). However, it is possible that DI could facilitate WT encapsidation. Indeed, during co-infection, two viruses can exploit a common pool of resources equally [49]. Given that the WT and DI genomes encode non-structural proteins that involved in capsid processing (3C/3CD) or encapasidation (2C), co-infection may generate an excess of these viral factors increasing the effiency of WT encapsidation [50, 51].

Our experimental observations indicate that the number of encapsidated genomes decrease on average from 7 to 9 hours post transfection (from 410 to 385 for WT in singly infected cells, from 115 to 81 for WT in dually infected cells, and from 355 to 336 for DI in dually infected cells). This effect is not observed in the model, because it assumes relatively slow decay of encapsidated genomes, which may take days [52], compared to the time scale of the experiment, which takes hours. A possible explanation is inactivation of virus particles overtime, an process that has not been incorporated in our model.

Nonetheless, the cross-validation experiment showed an overall good predictive power of the model, although it underestimated the relative WT output when the DI was transfected at low DI-to-WT input ratios, i.e. 0.46, S3 Fig). In addition, full model simulations with best parameter values overestimates the number of WT encapsidated genomes in singly infected cells while the estimation is more accurate in dually infected cells (Fig 3C and 3D). This discrepancy stems from the normalization process. The WT output was normalized using single WT infection in the cross-validation analysis, resulting in underestimations of WT relative output in dual infection condition. Despite these small deviations between model prediction and experimental measurements, even with differing measurements between experiments and simulations (see Material and methods), the model is able to recapitulate WT outputs for various inputs, which reinforces its robustness and predictive ability.

Identification of parameters affecting WT virion production

A sensitivity analysis showed that the WT/DI replication rate (P) and encapsidation rate (ω) are the most important parameters determinind the outcome of co-infection (Fig 4A and 4B). These parameters relate to the DI design and they can be engineered to modulate WT infection. Thus, the model provides estimations on the modifications to the DI genome enhancing its replication and/or encapsidation which lead to a significant 10-fold decrease in WT virion production. This result is consistent with the hypothesis that the mechanisms of DI interference with WT replication relies on exploitation competition [1-3] for resources required for genome replication and direct competition [5] for capsid proteins.

Impact of initial conditions

The model further predicts a crucial impact of temporal spacing and order of infection, as well as initial MOIs, on the proportion of WT virions at cell lysis (Fig 6 and S5 Fig). At equal MOIs, the DI particle needs to infect a cell within approximately a 2 hours window before or after the WT in order to be able to generate progeny, DI virions, and up to approximately 30 minutes after the WT in order to reduce WT yields (Fig 6A). Upon co-infection, the DI particles will maximize their virion production when WT and DI initial MOIs are approximately equivalent, and the WT particles will maximize their virion production when WT MOI is larger than DI MOI (Fig 6B and 6C). At equal MOIs, the difference between DI and WT virion production is maximized at approximately simultaneous co-infection of a cell. When WT MOI is larger than DI MOI, this difference maximization is obtained when the DI infects the cell before the WT. Conversely, when DI MOI is larger than WT MOI, it is obtained when the WT infects the cell before the DI (S5 Fig).

These results are in agreement with those of [53], who found that the extent of interference, assessed by the yield of WT poliovirus, is inversely proportional to the percentage of DI in the inoculum, and that it is also affected by varying the time interval between primary and secondary infection of a cell. In their viewpoint article, [2] highlight the importance of initial conditions, such as relative initial frequencies and the time of DI emergence in a population. An earlier review [54] also emphasized this aspect, while focusing on the exclusion of one virus strain by another. Experimental studies have shown inhibition of superinfection by a resident strain, in bacteria [55] and in viruses [56, 57]. Furthermore, the DI interference efficiency may also be affected by the the ability of picornaviruses rapidly induce resistance of the host cell to superinfection [58].

Importantly, preaccumulated replicon RNAs are not trans-encapsidated by capsids made from a coinfecting helper virus [59], consistent with the idea that encapsidation is coupled to RNA replication. Thus, if DI is the first to infect a cell, the genomes replicated before WT superinfection would not get trans-encapsidated by WT capsids.

Limited resource and co-evolution

The competition between WT poliovirus and DI particles within cell can be analysed in light of evolutionary game theory. For a game between WT cooperators and DI defectors, the pay-off matrix features a fitness of zero in the case of a population composed only of DIs, because they are unable to reproduce [60]. With such a feature, the evolution of a mixed WT and DI population is predicted to result in a polymorphic equilibrium, despite the greater pay-off that would result if the population was composed only of WT cooperators [60]. These strategies of cooperation and defection are common in viruses, as co-infection of the same host cell induces competition for shared intracellular products [60].

Resource availability can have important consequences on the dynamics and evolution of mixed pathogen populations [2, 61, 62]. For example, playing on resource availability could allow to slow the evolution of resistance to antimicrobial drugs [63]. In the case of DI particles, their presence within-cells infected by a WT virus is decreasing the number of resources available to the WT for replication and encapsidation, lowering WT virions production. Hence DI particles could be used to control WT infections, by lowering WT viral load, thereby facilitating further action of the immune system and/or drugs to clear the infection [46].

It has been suggested that the potential of DI particles to compete against a WT virus can be exploited to develop a new type of therapeutic antiviral strategy based on defective interfering particle competition [64]. Our model and the sensitivity analysis we have performed suggest that parameters P, the DI-to-WT replication ratio, and ω, the DI-to-WT encapsidation ratio, are the first and second most important parameters impacting the proportion of WT virions at cell lysis. Therefore, a rational strategy to strengthen interference activity of DI genomes and thus reduce the production of WT virions is to modify DI genomes towards higher replication rate and encapsidation efficiency. Such improvements may be realized by taking advantage of the evolvability of DI genomes. Serial co-passages of WT and DI particles followed by genetic analyses, as done in Fig 5A, 5B and 5C, would allow for the screening of other mutations providing higher replication or encapsidation of the DI. Also, the production of shorter DI genomes could lead to its faster replication.

Improving the interference at the intracellular level may cause less inhibition of WT viral load at the intercellular level, as there could be trade-offs. A reduced production of WT particles within-cells could result in a decreased MOI of WT viruses for the next infection cycle, and also to a decreased probability of a cell being co-infected. Furthermore, the narrow window of delay of co-infection for the DI to outcompete the WT as shown in Fig 6A also suggests the importance of simultaneous infection. Interestingly, recent studies show several possibilities for how co-infection is or can be favored [65-67]. Notably, the existence of vesicles containing multiple copies of virions as well as bacteria binding virions may increase the probability of simultaneous co-infection [65, 66]. Erickson et al. [67] reported that poliovirus binds lipopolysaccharide of bacteria, allowing co-infection of mammalian cells even at a low MOI. Finally, the potential of the combined DI-WT system to synergistically trigger an effective innate immune response (e.g. interferon) is also a potential avenue of investigation for the rational design of an antiviral therapy.

Conclusion

While ecological studies for the control of pathogen populations mainly focus on preventing or slowing down the emergence of drug resistance [63, 68, 69] or on the evolution of virulence [70, 71], we take an original approach here by rather studying how to use cheater defective pathogens, competing more efficiently for shared resources, for the control of disease-inducing pathogens. Since we learned the mechanisms of intracellular interference, in a future work we would like to apply these findings for the study of the competition between WT and DI at the larger level of the tissue, embedding intracellular knowledge. It would allow us to draw guidelines to optimize DI particles at this level, based on WT viral load inhibition, and further confirm their efficiency in vivo.

Materials and methods

Competition experiment between defective genomes and wild-type genomes

Cells

HeLaS3 cells (ATCC CCL-2.2) provided by R. Geller and J. Frydman (Stanford University) were maintained in 50% Dulbecco’s modified Eagle medium and 50% F-12 medium (DMEM/F12) supplemented with 10% newborn calf serum (NCS), 100 U/ml penicillin, 100 U/ml streptomycin and 2 mM glutamine (Invitrogen).

Construction of viral cDNA plasmids

The cDNA plasmid prib(+)XpA, encoding the genome of poliovirus type 1 Mahoney strain under T7-promoter and hammerhead ribozyme sequences, was reported previously [72]. Plasmid prib(+)XpA was digested by NruI and SnaBI (New England Biolabs) and ligated to produce prib(+)XpA lacking the poliovirus capsid-encoding region from 1175 to 2956 (prib(+)XpA-Δ1175–2956).

In vitro RNA transcription

Plasmids prib(+)XpA and prib(+)XpA-Δ1175–2956 were digested by EcoRI or PvuII. Linearized plasmids were used as templates to obtain WT and DI genomic RNAs by in vitro transcription. In vitro transcribed RNAs were purified by phenol-chloroform extraction and the quality of purified RNAs was analyzed by electrophoresis on a 1% agarose gel in tris-acetate-EDTA Buffer (TAE).

Transfection of defective interfering and wild-type genomes

Monolayer of HeLaS3 cells was trypsinized and washed three times in D-PBS. Cells were resuspended in 1 ml D-PBS and the number of cells were counted on a hemacytometer, followed by adjusting the concentration to 1 × 107 cells/ml. 800 μl of cells and virus RNAs (5μg of WT genomes and/or different amounts of DI genomes described later) were combined in a chilled 4-mm electroporation cuvette and incubated 20 minutes on ice. Cells were electroporated (voltage = 250 V, capacitance = 1000 μF) using Gene Pulser II (Bio-Rad), washed two times, and recovered in 14 ml prewarmed (37°C) DMEM/F12 medium with 10% NCS. Samples were distributed on 24 well plates (250 μl/well).

Samples were collected at different time points (0, 3, 6, 9 hours for titration, and 0, 2, 3.5, 5, 7, 9 hours for RNA extraction) after electroporation. For titration and evaluation of encapsidated RNAs, samples were then frozen and thawed three times, followed by centrifugation at 2,500 g for 5 minutes, and supernatants were collected. Samples for evaluation of encapsidated RNAs were further treated with mixture of RNase A (20 μg/ml) and RNase T1 (50 U/ml) (Thermo Fisher Scientific) for three hours. Samples were stored at −80°C.

Titration of virus samples

Monolayers of HeLaS3 cells in 6-well plates were infected with 250 μl of serially diluted virus samples at 37°C for 1 hour and then overlaid with DMEM/F12 including 1% agarose. After 48 hours of infection, infected cells were fixed by 2% formaldehyde and stained by crystal violet solution. Titers were calculated by counting the number of plaques and multiplying their dilution rates.

RNA extraction

250 μl of samples was added to 750 μl of TRI-reagent LS (Sigma Aldrich), and RNAs were extracted following the kit protocol. Briefly, 200 μl of chloroform was added to each sample, shaken vigorously, and incubated at room temperature for 10 minutes. Then samples were centrifuged at 12,000 g for 15 minutes at 4°C. The upper aqueous phase was transferred to a fresh tube and 0.5 ml of isopropanol was added. After incubation at room temperature for 10 minutes, samples were centrifuged at 12,000 g for 8 minutes at 4°C to precipitate RNAs at the bottom of the tube. The supernatant was removed and the residue was washed by 1 ml of 75% ethanol. After centrifugation at 7,500 g for 5 minutes at 4°C, the RNA pellets were dried for 5–10 minutes. RNAs were resuspended in nuclease-free water.

Reverse transcription

2.5 μl of RNA samples was mixed with 0.5 μl of 2 μM primer (5’-CTGGTCCTTCAGTGGTACTTTG-3’), 0.5 μl of 10 mM dNTP mix, and 2.5 μl of nuclease-free water. Samples were incubated at 65°C for 5 minutes, and then placed on ice for 1 minute. After adding 10 μl of cDNA synthesis mix (1 μl of 10× RT buffer, 2 μl of 25 mM MgCl2, 1 μl of 0.1 M DTT, 1 μl of RNaseOUT and 1 μl of Superscript III RT enzyme), samples were incubated at 50°C for 50 minutes, and then at 85°C for 5 minutes to terminate reactions. 1 U of RNase H was added to each sample, followed by incubation for 20 minutes at 37°C. Then, 0.1 U of Exonuclease I was added to each sample, followed by incubation at 37°C for 30 minutes and at 80°C for 15 minutes to terminate reactions. cDNA samples were stored at −20°C.

Design of primers and Taqman probes

Primers and Taqman probes for droplet digital PCR assay were designed with PrimerQuest Tool (Integrated DNA Technologies). The primers and probe for WT genomes are 5’-CCACATACAGACGATCCCATAC-3’, 5’-CTGCCCAGTGTGTGTAGTAAT-3’, and 5’-6-FAM-TCTGCCTGTCACTCTCTCCAGCTT-3’-BHQ1. The primers and probe for DI genomes are 5’-GACAGCGAAGCCAATCCA-3’, 5’-CCATGTGTAGTCGTCCCATTT-3’, and 5’-HEX-ACGAAAGAG/ZEN/TCGGTACCACCAGGC-3’-IABkFQ.

Droplet digital PCR assay

2 μl of serially diluted cDNA samples was mixed with 10 μl of 2× ddPCR supermix for probes, 1 μl of 20× WT primers/probe, 1 μl of 20× DI primers/probe, and 6 μl of nuclease-free water. 20 μl reaction mix of each sample was dispensed into the droplet generator cartridge, followed by droplet production with QX100 droplet generator. Then PCR was performed on a thermal cycler using the following parameters: 1 cycle of 10 minutes at 95°C, 30 cycles of 30 sec at 94°C and 1 minute at 60°C, 1 cycle of 10 minutes at 98°C, and held at 12°C. Positive and negative droplets were detected by QX100 droplet reader.

Model reduction

As our mathematical model (Eqs 1–4) presents a classical problem of parameter identifiability, we built a lower dimensional model to solve this problem by assuming that the decrease in resources due to viral uptake for replication follows a logistic decreasing function. This assumption was verified by analyzing the curves of R(t)θε as a function of time on a first set of “blind” optimizations. Thus, we can recast the model using the following lower dimensional description:

The logistic function (Eq 10) is characterized by the curve’s maximum value L and steepness s, and the time of the sigmoid’s midpoint t0. While this reduced version only decreases the number of parameters to be estimated by one (L, s and t0 instead of θ, ε, λ and γ), it partially solves the identifiability problem by removing biologically interpretable parameters and just assuming a logistic function for resource uptake and replication. An analytical justification of this model choice based on time scale separation can be found in the S1 Text.

Fit to experimental data

The model was fitted to the experimental data of Fig 1C and 1D in order to estimate model parameters describing our biological system. De novo RNA replication can only be detected around 2 hours post transfection [31, 73–75]. We thus used experimental data from 2 hours post transfection for parameter estimation. Indeed, there are several steps of poliovirus infection cycle before replication can start, including translation of positive-sense genomes [48] and transition from a linear, translating RNA to a circular RNA competent for replication [31, 73–75]. As our model does not account for those first steps, we only used experimental data from 2 hours post transfection for parameter estimation.

Raw experimental data are WT and DI total RNA copy number ( and , Fig 1C) and RNAse treated genome copy number (v and v, Fig 1D). The former corresponds to the total number of genomes (naked and encapsidated) and the latter to the number of encapsidated genomes. The numbers of WT and DI naked (i.e. non-encapsidated) genomes are thus: and . Additionally, as raw data was obtained at the cell population level, it was normalized by the average number of successfully transfected cells in order to get the average number of naked and encapsidated WT and DI genomes per cell.

In all, experimental data comprises three replicates of independent populations sampled at 2, 3.5, 5, 7 or 9 hours post transfection. Three different conditions were tested: (i) cells dually transfected by WT and DI genomes, (ii) cells transfected by WT genomes only and (iii) cells transfected by DI genomes only. Transfected volumes of WT and DI genomes were calibrated to a ratio of WT:DI = 4:1 in order to approximately obtain a ratio of 1:1 after transfection.

Parameter estimation was achieved through nlminb optimization function in R software [76] embedded in an iterative process. Each optimization consisted in minimizing the sum of the squared errors between experimental and simulated normalized data points for all variables and conditions. The least square function is as: with r indicating the replicate number and t the sampling time (t = 3.5 to 9 hours post transfection). Initial conditions for all variables in all infection conditions were obtained from averaging experimental observations over the 3 replicates at t0 = 2 hours post transfection. In Eq (11), we define: with x, resp. y, being either experimental (g, resp. v) or numerical (G, resp. C) naked, resp. encapsidated, genomes data and i for WT or DI.

The iterative process was applied as follows (see also S4 Fig for a schematic representation). Boundaries on parameter values were defined based on a first set of “blind” optimizations, the intervals still remaining large and realistic. For each parameter p, let us denote p the lower boundary and p the upper boundary. For the first iteration, random values of parameters were drawn from uniform distributions, as p ∼ Unif(p,p), defining the starting point for optimization. Let us denote the optimized parameter value. In subsequent iterations, the starting point for each parameter was then randomly drawn from a uniform distribution, as . In all, 20 iterations were conducted, and this iterative procedure was implemented 250 times, each time with a different random starting point in the first iteration. Thus, 250 × 20 = 5000 optimizations were performed in total.

The goodness of fit was evaluated by ordinary least square (Eq (11)) and sum of residuals R2 between experimental and simulated normalized data.

This optimization procedure was applied in two steps. In the first step, the optimization was performed on the reduced version of the model (Eqs 5–10), thus estimating nine parameters: the six parameters corresponding to (i) the DI (P and ω), (ii) the production of capsids (η), (iii) the encapsidation process (κ) and (iv) the decay rates of genomes and capsids (α and β); and the three parameters of the logistic function representative of the time-dependent replication rate (L, s and t0). In the second step, the six redundant parameters between the reduced and full version of the model (P, ω, η, κ, α and β) were fixed to their best estimated value obtained during the first step. The remaining four parameters (θ, ϵ, λ and γ) were estimated by optimizing the full version of the model (Eqs 1–6). The first 85 best sets of parameters that provide the most accurate fits to the experimental data are used for sensitivity analysis of the full model (Table 1 and Figs 4, 5 and 6, S3 and S5 Figs)) and to confirm the theoretical validity of the reduced model (Fig 3 and S2 Fig, and S1 Table).

DI variant identification, CirSeq and Analysis of allele frequencies

Monolayer of HeLaS3 cells was trypsinized and washed three times in D-PBS. Cells were resuspended in 1 ml D-PBS and the number of cells were counted on a hemacytometer, followed by adjusting the concentration to 1 × 107 cells/ml. 800 μl of cells and virus RNAs (5μg of poliovirus type 3 WT genomes and 5μg of DI genomes) were combined in a chilled 4-mm electroporation cuvette and incubated 20 minutes on ice. Cells were electroporated as described above. Supernatants containing PV3 WT virus and trans-encapsidated DI genomes were collected 24 hr post transfection and used to infect a fresh plate of HeLa cells. The infection procedure was repeated for 8 passages.

For preparing CirSeq libraries, each passaged virus (6 × 106 PFU) was further expanded in parental cells seeded in four 150 mm dishes. The culture medium was harvested before the appearance of severe CPE, and the cell debris was removed by centrifugation at 3,000 rpm for 5 min. The virion in the supernatant was spun down by ultracentrifugation at 27,000 r.p.m, 2 hours, 4ºC and viral RNA was extracted by using Trizol reagent. Each 1μg RNA was subjected to CirSeq libraries preparation as described previously [44, 45]. The CirSeq pipeline allows error detection in RNAseq through consensus generation and quality filtering to overcome the intrinsic error rate associated with reverse transcription. The experimental and computational protocols are described in detail previously [44, 45]. Briefly, purified polyA+ RNA from infected cells is fragmented to yield 80–100 bp fragments, circularized, and subject to rolling-circle reverse transcription. This procedure yields tandem reverse transcripts that are used to correct reverse transcription errors. Variant base-calls and allele frequencies were then determined using the CirSeq v2 package (https://andino.ucsf.edu/CirSeq). Circularized repeats are oriented to the reference genome and variants are called from raw reads based on phred33 scores of 20. These tandem variant-called reads are then aligned to each other to generate consensus sequences with a theoretical error of 10−6. Technical replicates of passaged libraries, and individual sequencing lanes, were compared after CirSeq mapping and pooled for analysis of fitness. Raw reads are deposited at Bioproject PRJNA669406. All consensus, and mapped reads from CirSeq are available at https://purl.stanford.edu/gv159td5450. Positive selected mutations within the DI genome, which frequency increased over passages, were selected and analyzed for replication fitness and competition with WT virus (Fig 5D and 5E).

Variants passaging experiment

A PV1-DI construct encoding the Venus (a green fluorescent protein) gene in place of the P1 gene was used for the following experiments. We transfected 5μg of PV3 and 1.25μg of PV1-DI genomes, and collected viruses at 24 hours after transfection as the passage 1 (P1) sample. Then, 1.0 × 106 HeLaS3 cells were infected with the different amounts of the P1 sample (500, 100, and 5 μl), and collected 24 hours after infection. The virus samples were passaged 8 times in the same manner. Viral RNAs for each sample were then analyzed by CirSeq [44] to identify accumulated mutations in DI.

Competition experiment and analysis

An additional experiment was conducted on eight DI mutants and the parental PV1-DI construct encoding the Venus gene. To obtain DI particles for the experiments, a packaging cell line was established. HeLaS3 cells were transfected with a pcDNA4 plasmid encoding the PV1 P1 (capsid) gene, followed by selection with Zeocin. We used a clone stably expressing P1 proteins (HeLaS3/P1) to generate DI particles. HeLaS3/P1 cells were transfected with DI RNAs by electroporation, and DI particles were collected 24 hours after transfection. Each of the DI variant was put in competition with the WT virus, and the number of naked genome copies was measured at 3, 4.33, 6.33 and 8.17 hours post transfection, with three replicates per time-point. Only the ratio of DI-to-WT naked genome copies is kept for further analysis.

We assumed that the mutations only affect the replication of each DI variant, hence we set all model parameters to their best estimated values (see Table 1) except for the DI-to-WT replication ratio P that we varied between 1 and 1.8. For simplicity, we assumed that the ratio of DI to WT encapsidation rate (ω) remained the same as estimated from the main experiment, as we showed that P is a more sensitive parameter than ω (Fig 4). The full model was used for this analysis (Eqs 1–6). We set the initial conditions to 10 copies of DI and WT (each) naked genomes at 2 hours post transfection. We assumed that there was initially no capsids or encapsidated genomes, and the number of resource units was set as in the main experiment to λ/γ. We recorded the simulated numbers of DI and WT naked genome copies at the experimental time-points for varying values of P (step of 0.01). The best P value for each DI variant was found as the one minimizing the sum of squared differences between experimental and simulated ratio of DI-to-WT naked genome copies.

Cross-validation

We cross validate the results of our optimization procedure by assessing how well the model is able to predict the relative WT virus burst size for various WT to DI initial ratios (after transfection) for which it has not been trained. We first obtain an optimal set of model parameters on our time series experimental data (g and v) featuring initial WT:DI = 1:1. We then test five additional DI-to-WT initial ratios, ranging from 0 to 3.6. Initial conditions for model simulations were set as the average of experimental values for each of the five initial ratios.

In the cross-validation experiment, evaluation of the relative WT virus burst size was based on the count of plaque-forming units (PFUs). In the time-series experiment that was used for parameter estimation, the number of WT virions (v) was estimated by digital droplet PCR. Assuming that the ratio of WT infectious to non-infectious particles and the multiplicity of infection (MOI) of WT virus are both constant independently of initial conditions, the relative PFU of WT virus for each initial condition should be a good proxy of the relative WT burst size.

The experiment was conducted on a cell population, and then the measurements were normalized by the number of successfully transfected cells. In some cases, the average experimental MOIs were small, potentially leading to not all cells being co-infected by WT and DI genomes. We integrated this aspect in our simulated burst size calculations, with the probabilities that a cell would be infected by both DI and WT genomes or only by WT genomes. We assume that the number of DI genomes infecting a cell X results from a Poisson distribution of parameter the average DI MOI n, as X ∼ Pois(n). The probability that no DI genome enters a cell is thus . Conversely, the probability that at least one DI genome enters a cell is . The expressions are equivalent for the WT virus. Let us denote the WT burst size in WT-DI infection as and the WT burst size in WT-only infection as . We weighted WT simulated burst sizes as follows: .

WT PFU experimental values and WT burst size model predictions () were normalized for all initial ratios by their respective values in the absence of DI genome (i.e. DI-to-WT initial ratio of 0). For the experimental data the average over the three replicates was taken for normalization. The performance of the model to predict relative WT burst size was evaluated by R2 and p-value of a Pearson correlation test between experimental and simulated datapoints.

A sensitivity analysis was performed to assess the relative importance of each parameter on the proportion of WT virions at cell lysis, defined as:

Based on parameter estimation, each parameter was approximately varied by ±50% of its best estimated value. Based on these boundaries, each parameter was allocated a vector of five equidistant values (except for α that was not varied because it was estimated at 0, see Table 1). Then, all distinct combinations of parameter values were tested according to a full factorial design. In all, 59 = 1,953,125 simulations of the full model were performed. All simulations started at time 0 hours post transfection with 10 copies of WT and DI genomes (G(0) = G(0) = 10), no capsids nor encapsidated genomes (C(0) = C(0) = C(0) = 0), and R(0) = λ/γ. Simulations were conducted until 9 hours post transfection. An analysis of variance (ANOVA function in R software) was then conducted to assess the importance of each parameter and their second-order interactions on the variance of Φ.

Impact of delay and multiplicities of infection on WT and DI burst sizes

In the experiment, WT and DI genomes were co-transfected to cells and in quantities yielding an MOI ratio of approximately WT:DI = 1:1. We conducted two sets of simulations to study the impact of varying either (i) the time between cell infection by WT and DI or (ii) the MOIs of WT and DI on their burst sizes. All the simulations were conducted on the full version of the model (Eqs 1–6). In the first set of simulations, WT and DI burst sizes were recorded for various delays between primary and secondary infection of a cell, ranging from -7 to +7 hours for the time of DI infection compared to the WT. The MOI upon infection of the cell was set to 10 for both WT and DI (i.e. G (0) = 10 genomes and G (t) = 10 genomes, with t the delay for DI infection), the number of capsids and encapsidated genomes to 0 and R(0) to λ/γ. In the second set of simulations, WT and DI burst sizes were recorded for various WT and DI initial MOIs, ranging from 0 to 1000 virions. All the other variables were set as described for the study of delays. Then, WT burst sizes were normalized by the WT burst size corresponding to WT:DI = 1:0 initial MOIs (i.e. infection by the WT virus only at low MOI), and DI burst sizes by the DI burst size corresponding to WT:DI = 1:1 initial MOIs (i.e. infection by both WT and DIs at low MOI).

Experimental data and model fit for the naked and encapsidated wild-type (WT) and defective interfering (DI) genomes, with same encapsidation rate.

Evolution of the number of WT and DI (A-B) naked genome copies and (C-D) encapsidated genome copies with time, from 2 to 9 hours post transfection (hpt). (A & C) show data in dually transfected cells whereas (B & D) show data in singly transfected cells. WT and DI results are shown in blue and red color, respectively. Crosses indicate experimental data for 3 replicates per sampling time at 2, 3.5, 5, 7 and 9 hpt. Solid curves show the fit of the full model with same encapsidation rate for WT and DI naked genomes (Eqs. (S10)–(S13) in S1 Text).

(PDF)

Click here for additional data file.

Histograms and correlations of the genome replication factor (θ), the resource linear production (λ), the decay rate (γ), and the resource capture rate (ε).

The best 123 estimated values for each parameter are represented. A-D: Histograms of best estimated values for θ (A), ε (B), λ (C) and γ (D). E: Correlation between resource production to decay rate and genome replication factor. The best fit curve is shown in red and its equation is provided (Pearson p-value <2.2 ⋅ 10−16 and R2 = 0.897). F: Time-dependent replication rate given by the reduced model (Λ (t), dashed line) and the full model (θεR(t), plain line).

(PDF)

Click here for additional data file.

Cross-validation of the model (Eqs (1)–(6)).

Three experimental replicate values (black dots) of relative WT virus output are represented for various DI to WT input (proxy of multiplicities of infection) ratios. Red dots indicate predicted relative WT output starting with the same experimental input ratios. Experimental WT output corresponds to PFU while simulated WT output corresponds to burst size (number of encapsidated genomes at 9 hours post infection). All outputs were normalized by the output value (or the mean for experimental data) of WT:DI = 1:0 input ratio. R-squared and p-value of a Pearson correlation test between experimental and predicted WT outputs are given in the graphic.

(PDF)

Click here for additional data file.

Diagram of the parameter fitting procedure.

(PDF)

Click here for additional data file.

Model-predicted delay at various multiplicities of infection (MOIs).

A: Predicted impact of delay for DI particle infection of a cell (x-axis, in hours) on WT (blue) and DI (red) burst sizes (y-axis). One line represents one DI MOI and one column one WT MOI. The grey vertical line indicates no-delay (simultaneous infection). The red vertical line indicates the peak of DI burst size and the green vertical line the maximum difference of DI to WT burst size. B: Heat map of predicted delay for the maximum difference of DI to WT burst size (green lines in A). C: Heat map of predicted delay for WT and DI burst size curves intersection.

(PDF)

Click here for additional data file.

Experimental data and reduced model fit for the naked and encapsidated wild-type (WT) and defective interfering (DI) genomes, with P fixed to the ratio of WT to DI genome lengths.

Evolution of the number of WT and DI (A-B) naked genome copies and (C-D) encapsidated genome copies with time, from 2 to 9 hours post transfection (hpt). (A & C) show data in dually transfected cells whereas (B & D) show data in singly transfected cells. WT and DI results are shown in blue and red color, respectively. Crosses indicate experimental data for 3 replicates per sampling time at 2, 3.5, 5, 7 and 9 hpt. Solid curves show the fit of the reduced model with P fixed to the ratio of WT to DI genome lengths (P = 7515bp/5733bp, Eqs (7)–(10)).

(PDF)

Click here for additional data file.

Model selection.

For each model, described in S1 Text, the goodness of fit is evaluated by the squared Pearson correlation coefficient between experimental and fitted data (R2), and the quality of the models is evaluated by the log-likelihood (−2 ⋅ log(L)) and Akaike information criterion (AIC) from a linear model between experimental and fitted data.

(ZIP)

Click here for additional data file.

Model selection procedure and analytical study of the model.

(PDF)

Click here for additional data file.

20 Feb 2021

Dear Raul,

Thanks very much for submitting your manuscript "Experimental and mathematical insights on the interactions between poliovirus and a defective interfering genome" for consideration at PLoS Pathogens. As with all papers reviewed by the journal, your manuscript was reviewed by members of the editorial board and by three independent reviewers. In light of the reviews (below this email), we would like to invite the re-submission of a significantly-revised version that takes into account the reviewers' comments.

We cannot make any decision about publication until we have seen the revised manuscript and your response to the reviewers' comments. Your revised manuscript is also likely to be sent to reviewers for further evaluation.

When you are ready to resubmit, please upload the following:

[1] A letter containing a detailed list of your responses to the review comments and a description of the changes you have made in the manuscript. Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out.

[2] Two versions of the revised manuscript: one with either highlights or tracked changes denoting where the text has been changed; the other a clean version (uploaded as the manuscript file).

Important additional instructions are given below your reviewer comments.

Please prepare and submit your revised manuscript within 60 days. If you anticipate any delay, please let us know the expected re-submission date by replying to this email. Please note that revised manuscripts received after the 60-day due date may require evaluation and peer review similar to newly submitted manuscripts.

Thank you again for your submission. We hope that our editorial process has been constructive so far, and we welcome your feedback at any time. Please don't hesitate to contact us if you have any questions or comments.

Best regards,

Bert Semler

Guest Editor

PLOS Pathogens

Adolfo García-Sastre

Section Editor

PLOS Pathogens

Kasturi Haldar

Editor-in-Chief

PLOS Pathogens

​orcid.org/0000-0001-5065-158X

Michael Malim

Editor-in-Chief

PLOS Pathogens

orcid.org/0000-0002-7699-2064

***********************

Reviewer's Responses to Questions

Part I - Summary

Please use this section to discuss strengths/weaknesses of study, novelty/significance, general execution and scholarship.

Reviewer #1: There is a renewed interest in defective interfering (DI) particles (and other types of defective genomes) because they are produced and maintained in populations of many viruses (particularly RNA viruses), and because they may play regulatory roles, particularly in the interaction with the host immune system. In the present study, experimental data on WT and DI poliovirus are presented, and confronted with a relatively simple model (a reduced model to determine some parameters and a full one) of differential equations, based on one previously developed for HIV-1. The model considers as central parameters rates of replication and encapsidation to explain the outcome of WT-DI competitions. The predictive value of the model is tested with some experimental results, and reasons for agreements and disagreements are discussed.

Although the study is of potential interest, in its present form it fails to relate the findings to previous work on DI particles, and models on interference mechanisms. The authors fail also to underline the major conclusions, novelty, and implications of the experimental results and model. In addition, some experimental protocols should be described in greater detail.

1. Other previously described models and their predictions should be compared with the differential equation-based model presented in the study (Frank, J. Theor. Biol. 2000; Stauffer Thompson et al. J. Gen. Virol. 2009; Laske et al. Virus Res. 2016; Kirkwood and Bangham, PNAS, 1994; among others). For example, the deterministic model of Kirkwood and Bangham emphasizes unpredictable effects of DI activity, a point not discussed in the present study. The merits of the new model relative to previous ones, and its advantages for future work should be explained.

2. According to Figure 5 legend the poliovirus used was type 3 (PV3), while according to the Supplementary Methods the WT was PV3 while the DI was constructed on PV1. The poliovirus used as WT and for DI construction should be explained and, if different viral types were used, a justification is needed. Also, how the present DI structure compares with others previously characterized for poliovirus should be indicated.

3. A critical experimental value is the proportion of encapsidated RNA measured by resistance to a RNase mixture. Controls (or literature references of previous work by the authors where the controls are described) consisting of mixtures of encapsidated (protected) purified particles and free viral RNA (unprotected) to evaluate the reliability of the quantifications should be described. Indicate also the proportion of cells that were transfected.

4. There is some confusion about what the authors mean by “replication rate”. Having a shorter RNA, DIs are expected to complete a round of template copying earlier than WT RNA. This does not imply that the replicase (or replicative complex or organelle) travels faster along the RNA template when copying DI RNA than WT RNA. Unless a different rate of polymerase progression has been documented, it would be better to refer to replication being completed in a shorter time for DI than WT. An additional ambiguity becomes apparent when the authors state that “replication rate” in their model includes also translation of the positive sense RNA genome. To complicate matters, in section 4.1 the authors suggest that the DIs will replicate faster than the WT by a factor given by the ratio of WT to DI genome lengths, but the observed discrepancy may be attributed to the various processes linked to replication (which processes in addition to translation mentioned earlier?). In this section they also refer to “replication speed”.

5. The data presented in Fig. 5 require additional information. In panel A, the passages produced a number of mutant DIs. The reader assumes that no mutations occurred in WT but this should be stated. Fitness of variant DIs is mentioned but no information on their values and how was fitness defined and determined is provided. DI mutant 1952 in panel D is named 19522 in panel C. The meaning of the volumes given in panel B PV3 should de explained. The differences between panels E and F and their significance are not clear.

6. While the Results and figure legends will benefit from expansion, the Discussion is very long and with speculative arguments. At times, it is considerably naïve, for example when referring to evidence of dengue coinfections in India, when the literature includes many examples of coinfections. It is unlikely that the non-specialist reader can capture the main conclusions and novelty of the present study among a quite disorganized amount of information, mainly in the Discussion.

7. At times the distinction between what is predicted by the model without experimental verification and what has been experimentally studied is not clear. While the data in Figure S2 are an experimental verification of model predictions, those in Figures 6 and S4 do not seem to have an experimental verification with the poliovirus system. If this is correct, it should be indicated.

Reviewer #2: Summary:

The presented manuscript `Experimental and mathematical insights on the interactions between poliovirus and a defective interfering genome’ aims to investigate the intracellular time course dynamics of poliovirus and its defective viral genomes (DVGs) that interfere with the full-length virus (FLV) via two mechanisms. First, DVGs are replicons which need FLV’s capsid to form particles, thus hampering encapsidation of the FLV. Second, DVGs replicate more rapidly, thus depleting cellular resources at the expanse of the FLV. Briefly, the authors used a mathematical model to fit the time course genome concertation data of FLV and DVGs to estimate the model parameters. They conclude that the most important parameters that determine the outcome of FLV and DVG competition are DVG-to-FLV replication ratio, P, and DVG-to-FLV encapsidation ratio, omega. The authors then use additional experimental data to validate their model predictions.

Strengths/Weaknesses:

I personally find the idea of combining the experimental data on FLV and DVGs with mathematical modelling interesting and sufficiently novel. The study addresses a rather unexplored problem of numerical quantification of DVG parameters during the FLV-DVG competition as these could help to better design treatment protocols should DVG research progress to a phase of therapy development in the future, particularly the optimal dose and the time of DVG administration. On the other hand, successful DVG infection would have to be secured. Cell transfection is a limitation and it is unclear to me whether FLV-DVG multiple-cycle infection dynamics would be possible to track by the authors. The authors also present their assumption about the model based on the knowledge of their DVGs as a result claiming that the model predicts P and omega to be greater than 1 (e.g. Discussion, lines 386-387). If P=omega=1, there’s no competitive advantage and the outcome will depend on the initial conditions/time of delivery. Computational methods need to be described in more detail.

Reviewer #3: PPATHOGENS-D-21-00004_reviewer

Summary

(Given how covid has changed my working life, I'm reviewing this as I would review a preprint I need. As in, does it make sense, is it basically sound, and do I understand what they're saying? If yes to all three, then I don't think it needs much revision. So I read carefully, but I'm not looking to nitpick in this review.)

This is a really nice example of model and experiment working together. The question of interfering particles is important to understanding viral fitness, the experimental methods are clean, and the modeling is insightful and simple but not-too simple and directly relevant to the experiments. And the model makes successful, non-trivial out-of-sample predictions. Any modeling paper that actually predicts something it didn't fit to is good work (and much more rare than it should be!)

I'm familiar with the both the virology and model, but much better equipped to critique model methods than lab. The model is reasonable, as is the optimization procedure to fit it, and the results are clearly communicated.

In general, the paper is clearly written. I felt like I knew exactly what the paper was about from the intro and why you were interested in it on the first read (with one minor caveat below). This impression held up through the rest of it.

**********

Part II – Major Issues: Key Experiments Required for Acceptance

Please use this section to detail the key new experiments or modifications of existing experiments that should be absolutely required to validate study conclusions.

Generally, there should be no more than 3 such required experiments or major modifications for a "Major Revision" recommendation. If more than 3 experiments are necessary to validate the study conclusions, then you are encouraged to recommend "Reject".

Reviewer #1: See Part I

Reviewer #2: The presented mathematical model is a modification of a model from Rouzine and Weinberger (2013) in which a new variable was introduced to simulate cellular resource depletion by FLV and DVGs. I have the following questions regarding modelling, and I list them point by point:

1) Have the authors tried to fit the model with a fixed P? There is a note in the manuscript that the estimated ratio could correspond to 1.311 based on the FLV and DVG genome lengths. This would lower the dimensionality of the fitting problem. How would the predictions change?

2) What is the reason to assume that DVGs encapsidate faster than FLV? How do the predictions change if no such competitive advantage is assumed for DVGs?

3) When fitting the model, both naked and encapsidated genome abundancies define the objective function to be minimized by least squares. I am not an experimentalist and don’t understand the protocol much, but how do the authors differentiate between encapsidated genomes and empty capsids (e.g. Figure 3)? Are the data in Fig. 3 truly just encapsidated genomes? I think this is important to also write in the text as modellers could struggle to find this information in there.

4) What is the authors’ opinion on the possible assumption that DVG uses less cellular resources per genome than FLV, possibly due to its smaller size? Is that something they would consider in their model?

5) Line 222: Predictions also demonstrate that DVGs are more efficiently encapsidated than FLV (omega=2.185). I disagree with this sentence because omega>1 is the assumption the authors make, not a prediction. Can the model reproduce the data when omega=1? See my comment above.

6) Line 226-227: Again, not a prediction, interference is an assumption in the model which allows to fit the data.

7) In the equation (10) for Lambda(t), the factor s should have units of min^-1 and not be unitless as (t-t0) is in mins and L is unitless. Or is it L that should have units?

8) To be honest, I don’t follow the model reduction procedure. When I look at the models (1)-(4) and (7)-(10), they should at least yield the same units, that is genomes/time or capsid/time. When I look at eq. (2) for capsid dynamics, then e.g. the rate (eta x Lambda(t) x G_wt) has units (capsid x genome^-1 x nothing x genome) which results in just capsid as units, given Table 1. For this reason, the rates are not comparable. Also, lines 236-237, these time-dependent rates 3.07x10^-2 and 3.02x10^-2 are unitless?

9) Figure S2: Is this a reasonable validation? This is related to my question about free capsids – if the quantification of encapsidated genomes accounts for free capsids, PFUs might not correspond and normalization will mask this.

10) How exactly are the data from Figure 1 used in modelling? This part is a bit unclear to me.

11) Figure 5: What exactly is this figure saying? If I understand correctly, 8 DVGs were selected and transfected with the FLV and the DVG/FLV ratio was determined over the time. The DVG/FLV ratio at 8h post transfection was plotted against the model (which one?) prediction of P for each variant to demonstrate that higher DVG/FLV ratio implies higher P. Can you say the same for omega?

12) What were the criteria for those 8 DVGs to be selected? Have the authors used any bioinformatics methods or tools for variant calling? If so, they should be detailed in the Methods section, too.

13) Since I am not an experimentalist, it is difficult for me to assess which experiments can and cannot be technically done. Figure 6A and S4 show the impact of different time of delivery of DVGs on the FLV growth. If repeated transfection is possible, experiments to validate these predictions are desirable.

14) Lines 460-461: The authors write “This decay might also be real, but we chose not to include it in our model because of the too large number of parameters it has.” It is not acceptable to kill a parameter because it’s inconvenient.

Reviewer #3: Major comments

The idea to use the reduced model to pin down an unidentifiability before going back to the full model is clever, but I suspect the authors could go further with this line of reasoning. I interpret the fact that the reduced model works well (in fact slightly better than the full model) says G_wt(t) + G_di(t) is to an excellent approximation a linear function of R(t), such that dR/dt reduces to the logistic differential equation. Which is to say that the sum of Gs is very fast on the timescale R evolves, or that d(G_wt+G_di)/dt >> dR/dt . I don't quite have it in me tonight to formally work through this slow/fast separation of scales argument, but the authors might do so to better understand why the reduced dynamics should in fact be a good description of the system (http://www.bio-physics.at/wiki/index.php?title=Separation_of_Timescales; also http://www.scholarpedia.org/article/Multiple_scale_analysis). From the figures, we can see that the fast dynamics are the competition between G_wt and G_di, such that the total is roughly only a function of R while they play in between. The difference between the singly and dually-infected situations may simply come down to the initial condition (P) and may be derivable by hand (this is supported by the full model results for singly vs dually-infected cells). If this argument holds up, the reduced model may arguably be the better model and the main text can be simplified to only have one model in the figures and with a simpler description of model fitting. (THAT HAVING BEEN SAID, if the authors can't make this argument work, don't find it as enlightening as I do, or in fact I'm just wrong, the current model is publishable.)

Now Fig 5F is a nice prediction!

**********

Part III – Minor Issues: Editorial and Data Presentation Modifications

Please use this section for editorial suggestions as well as relatively minor modifications of existing data that would enhance clarity.

Reviewer #1: MINOR POINTS:

a. It is not obvious to this reviewer why it is expected that the mutations in DI will not affect the encapsidation efficiency.

b. The last sentence of Results is hard to understand.

c. End of Section 4.1: give all units for the range of values for the initial replication rate (is it 3.02-3.07.10E-2 or 3.02.10E-2-3.07.10E-2?)

Reviewer #2: Line 166: limiting resource are produced at a constant rate lambda, not linear, I suppose.

Page 12: Burst sizes are relative or absolute numbers? If the latter, the values are without units (e.g. line 320).

Page 12, line 334: It may be useful to interpret the results with respect to FLV minimization and not FLV maximization or DVG maximization since reduction of FLV should be the goal with DVG interference/treatments.

Line 682: The formula for g_wt = g_wt^tot – v_di should probably be g_wt = g_wt^tot – v_wt?

Line 693: It should probably be “…Each optimization consisted in minimizing the sum of squared errors between experimental and simulated normalized data…”. Least squares is a method.

Figure 4, line 6: omega is a ratio not a rate.

Supplementary Methods/Description of the models: I’m not sure how the models starting with eq. (13) up to (21) are relevant. The original model is enough as it describes all main determinants of the FLV-DVG dynamics that the authors want to address. And the reduced model is already in the Methods.

Line 1054: instead of “…minimizing the sum of the least square difference between…” should be “…minimizing the sum of squared differences between…”

Reviewer #3: Minor comments

Line 39: I really want an Oxford comma after "apparent"

Line 91: I was hoping for an em-dash saying which two steps... (I had to go back to the abstract to find replication and encapsidation, but it would be nice to be re-reminded and cued up for the results after the intro).

Figure 3 caption: Is "Experimental data and model predictions" the right phrase? The model was fit to this data and it's not an out-of-sample prediction. I would prefer simply "Experimental data and model" although I acknowledge "prediction" is (over)used for this type of figure by others.

Figure S1E: "-1log" is strange notation. Maybe just "-log"

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

Reviewer #3: No

Figure Files:

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, . PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at .

Data Requirements:

Please note that, as a condition of publication, PLOS' data policy requires that you make available all data used to draw the conclusions outlined in your manuscript. Data must be deposited in an appropriate repository, included within the body of the manuscript, or uploaded as supporting information. This includes all numerical values that were used to generate graphs, histograms etc.. For an example see here on PLOS Biology: http://www.plosbiology.org/article/info%3Adoi%2F10.1371%2Fjournal.pbio.1001908#s5.

Reproducibility:

To enhance the reproducibility of your results, PLOS recommends that you deposit laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. For instructions see

23 Jun 2021

Submitted filename: Poin-by point responses PLoS Pathogens DIP.docx

Click here for additional data file.

28 Jul 2021

Dear Raul,

We are pleased to inform you that your manuscript 'Experimental and mathematical insights on the interactions between poliovirus and a defective interfering genome' has been provisionally accepted for publication in PLOS Pathogens.  Good show, Doc!

Before your manuscript can be formally accepted you will need to complete some formatting changes, which you will receive in a follow up email. A member of our team will be in touch with a set of requests.

Please note that your manuscript will not be scheduled for publication until you have made the required changes, so a swift response is appreciated.

IMPORTANT: The editorial review process is now complete. PLOS will only permit corrections to spelling, formatting or significant scientific errors from this point onwards. Requests for major changes, or any which affect the scientific understanding of your work, will cause delays to the publication date of your manuscript.

Should you, your institution's press office or the journal office choose to press release your paper, you will automatically be opted out of early publication. We ask that you notify us now if you or your institution is planning to press release the article. All press must be co-ordinated with PLOS.

Thank you again for supporting Open Access publishing; we are looking forward to publishing your work in PLOS Pathogens.

Best regards,

Bert L. Semler

Guest Editor

PLOS Pathogens

Adolfo García-Sastre

Section Editor

PLOS Pathogens

Kasturi Haldar

Editor-in-Chief

PLOS Pathogens

​orcid.org/0000-0001-5065-158X

Michael Malim

Editor-in-Chief

PLOS Pathogens

orcid.org/0000-0002-7699-2064

***********************************************************

Reviewer Comments (if any, and for reference):

Reviewer's Responses to Questions

Part I - Summary

Please use this section to discuss strengths/weaknesses of study, novelty/significance, general execution and scholarship.

Reviewer #1: (No Response)

Reviewer #2: I have read the revised version of the manuscript entitled "Experimental and mathematical insights on the interactions between poliovirus and a defective interfering genome" and recommend it for publication as the authors have carefully addressed all my and other reviewers' comments. I have no further suggestions or questions.

**********

Part II – Major Issues: Key Experiments Required for Acceptance

Please use this section to detail the key new experiments or modifications of existing experiments that should be absolutely required to validate study conclusions.

Generally, there should be no more than 3 such required experiments or major modifications for a "Major Revision" recommendation. If more than 3 experiments are necessary to validate the study conclusions, then you are encouraged to recommend "Reject".

Reviewer #1: (No Response)

Reviewer #2: (No Response)

**********

Part III – Minor Issues: Editorial and Data Presentation Modifications

Please use this section for editorial suggestions as well as relatively minor modifications of existing data that would enhance clarity.

Reviewer #1: (No Response)

Reviewer #2: (No Response)

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

22 Sep 2021

Dear Dr. Andino,

We are delighted to inform you that your manuscript, "Experimental and mathematical insights on the interactions between poliovirus and a defective interfering genome," has been formally accepted for publication in PLOS Pathogens.

We have now passed your article onto the PLOS Production Department who will complete the rest of the pre-publication process. All authors will receive a confirmation email upon publication.

The corresponding author will soon be receiving a typeset proof for review, to ensure errors have not been introduced during production. Please review the PDF proof of your manuscript carefully, as this is the last chance to correct any scientific or type-setting errors. Please note that major changes, or those which affect the scientific understanding of the work, will likely cause delays to the publication date of your manuscript. Note: Proofs for Front Matter articles (Pearls, Reviews, Opinions, etc...) are generated on a different schedule and may not be made available as quickly.

Soon after your final files are uploaded, the early version of your manuscript, if you opted to have an early version of your article, will be published online. The date of the early version will be your article's publication date. The final article will be published to the same URL, and all versions of the paper will be accessible to readers.

Thank you again for supporting open-access publishing; we are looking forward to publishing your work in PLOS Pathogens.

Best regards,

Kasturi Haldar

Editor-in-Chief

PLOS Pathogens

​orcid.org/0000-0001-5065-158X

Michael Malim

Editor-in-Chief

PLOS Pathogens

orcid.org/0000-0002-7699-2064