Title: Theoretical validation of growth curves for quantifying phage-bacteria interactions

Bacteria-infecting viruses, bacteriophages, are the most abundant biological entities on the planet, frequently serving as model systems in basic research and increasingly relevant for medical applications such as phage therapy. A common need is to quantify the infectivity of a phage to a given bacterial host (or the resistance of a host to a phage). However, current methods to quantify infectivity suffer from low-throughput or low-precision. One method that has the potential for high-throughput and high-precision quantification of phage-bacteria interactions is growth curves, where bacterial density is measured over time in the presence and absence of phages. Recent work has proposed several approaches to quantify these curves into a metric of phage infectivity. However, little is known about how these metrics relate to one another or to underlying phage and bacterial traits. To address this gap, we apply ecological modeling of phage and bacterial populations to simulate growth curves across a wide range of trait values. Our findings show that many growth curve metrics provide parallel measures of phage infectivity. Informative metrics include the peak and decline portions of bacterial growth curves, are driven by the interactions between underlying phage and bacterial traits, and correlate with conventional measures of phage fitness. Moreover, we show how intrapopulation trait variation can alter growth curve dynamics. Finally, we test the sensitivity of growth curve metrics to inoculum densities, and assess techniques to compare growth curves across different bacterial hosts. In all, our findings support the use of growth curves


Introduction
Bacteriophages are the most abundant and diverse biological entities on the planet (1-4). Not only do phages play key roles in global ecosystem functioning, they have long been used as laboratory models in molecular biology, ecology, and evolutionary biology, and are increasingly relevant for medical application as phage therapy (5). Phages depend on bacterial host cells, hijacking cellular machinery to replicate their genetic material and produce new virion particles. However, phages can vary immensely in their success at infecting and replicating in a given host cell. This can be driven both by variation in phage adaptation to infect and replicate in a given host, as well as variation in bacterial adaptation to resist phage infection and replication via numerous defense mechanisms (6,7).
Because of this variation, a common need is to quantify the infectivity of a phage on a given bacterial host (reciprocally, the resistance of a bacterial host to a given phage). Current approaches to quantify phage infectivity (and bacterial resistance) exhibit tradeoffs between throughput and precision (8,9). Specifically, some approaches (e.g. the efficiency-of-plaquing assay (10)(11)(12)) provide low-throughput but high-precision measures of phage infectivity, while other approaches (e.g. cross-streaks and spot-tests (13)) provide high-throughput but low-precision measures of infectivity, often only providing a binary measure of infectivity. Moreover, none of these approaches can be easily automated, for instance with assays conducted by liquid-handling robots.
One method that has the potential for scalable high-throughput quantification is growth curves of bacteria in the presence of phages. Growth curves are collected by growing bacteria in an automated spectrophotometer, which obtains measures that are proxies for bacterial cell densities (e.g., lightscattering/optical density) at defined time intervals. By collecting growth curves of bacteria cultured in the presence of phages, once can observe the effect of phages on the population density of bacteria over time. Indeed, growth curves are a widely-used assay of phage-bacteria interactions, often to measure the infectivity of phages on bacterial strains .
Recent work has proposed a number of metrics that can be calculated from growth curve data as quantitative measures of phage infectivity. Some approaches focus on specific portions of the growth curve data (Fig 1), including maximum density, area under the curve (AUC), or extinction time (42)(43)(44)(45)(46)(47)(48)(49)(50). Other proposals use multivariate analyses like principal component analysis to produce summary metrics (23,51). However, little is known about how these metrics relate to one another, how they vary depending on underlying phage and bacterial traits, and how they relate to conventional measures of phage infectivity.
To address this knowledge gap, we leverage an approach from another body of literature: the ecological modeling of phage-bacteria communities. There has been extensive theoretical work on phage-bacteria interactions (52,53), but little theoretical work has applied ecological modeling to understand how growth curves can be used to quantify phage infectivity (49). Moreover, direct implications of existing theoretical work for growth curve applications are limited, since theoretical work often focuses on longterm equilibria while experimental growth curve analyses use short-term dynamics.
Here, we bridge this gap between the two fields, using ecological modeling of population dynamics to inform applications of growth curves to quantify phage-bacteria interactions. We used common modeling approaches of susceptible bacteria, infected bacteria, phages, and nutrients to simulate growth curves of a single bacteria population and a single phage population with varying traits (Fig 1).
For instance, we manipulated traits like infection rate, phage burst size (number phages produced per infected cell), phage lysis time (the duration of time between host infection and lysis, when viral progeny are released), bacterial growth rate, and bacterial carrying capacity. We then applied proposed growth curve analysis metrics to these simulated data. Our findings show that many metrics provide parallel measures of phage infectivity, and that metrics often quantify generalized phage fitness. We also highlight the limitations of growth curve analyses in cases where bacterial resistance to phages arises and when comparing phages across different bacterial hosts. Simulating phage-bacteria growth curves. A. A modified SI-type system of delay differential equations was used to model the growth of bacterial and phage populations during growth curves. Susceptible bacteria (S) grow logistically (at rate u) up to a carrying capacity (k). Phage virions (P) infect susceptible hosts (at a rate a) and after a discrete lysis time (τ) release progeny virus particles (b per infection). Phages can also superinfect cells (at a rate z, relative to a). Lysis of cells releases nutrients back into the environment (at rate d). B. Simulation of this model shows initial growth of phage population and susceptible and infected bacterial populations, before bacteria go extinct and phage growth ceases. C. Experimentally, the typical available data output is the total bacterial population (S + I), plotted here. D. Metrics of the total bacterial population can then be extracted for comparison of different phages infecting bacteria, including peak density, peak time, time when density drops below a threshold ('extinction' time), and area under the curve (AUC).

Model implementation
We built a delay differential model of bacterial and phage growth (54) (Eqs. 1 -4). Our model explicitly tracks the density of four populations through time: susceptible bacteria (S), infected bacteria (I), free phage particles (P), and nutrients (N). Note that we defined N to be in the same units as S and I, such that consumption of one unit of nutrients produces one susceptible cell. Moreover, when d = 1, the model of bacterial growth reduces to logistic growth. For convenience, we let k be the carrying capacity of cells, and define = + + . All populations without a subscript are assumed to be the current value at time t. Description of all parameters and the default value used in our simulations are provided in Table 1.
This model produces growth curves under the assumption of no plasticity in any traits. However, susceptibility to phage infection is known to decline as bacterial growth slows and bacteria enter stationary phase (52,(55)(56)(57)(58)(59). To identify simulated growth curves where this assumption is most likely to be violated, we identified growth curves which approximately reached their carrying capacity if + ≥ 0.9 * at any point in time. In some plots in the main text these simulations are excluded (as noted in the figure legends), but plots with all simulations are included in the supplement.

Modeling variation in susceptibility to infection
To explicitly model variation in the susceptibility of bacteria to phage infection, we implemented two previously-published approaches (52). In the first, we modeled decreasing phage growth as nutrient density and bacterial growth rate decline via changes in the infection rate (60), burst size (57,60,61), and lysis time (26,52,62) (see Supplemental Results). In the second, we modeled transitions to a resistant nongrowing sub-population of cells (63) (Eqs. 5 -6) (See Supplemental Results for models where the resistant sub-population of cells has non-zero growth, and where transitions from S to R are non-linear). Description of all parameters and the default value used in our simulations are provided in Table 1.

Results
We first sought to test how the various different growth curve metrics relate to one another, since it remains unclear whether, and to what degree, different metrics are correlated. To do so, we simulated growth curves using phages with a wide array of trait values and calculated metrics for each growth curve. Across all curves, bacterial density increases, peaks, then declines to extinction. Perhaps unsurprisingly, we find that growth curve metrics are tightly correlated with one another (Figs 2, S1). For instance, there is a close relationship between the peak density of the bacterial population (peak density) and the time when that peak occurs (peak time) (Figs 2A, S1A). Similarly, peak time and area under the curve (AUC) are closely related (Figs 2B, S1B). In addition, peak time is approximately linearly related to the time when bacterial density drops below some threshold density ('extinction' time) (Figs 2C, S1C, S2). These patterns are robust to changes in bacterial traits (Fig S3). These findings suggest that changes in phage traits largely alter growth curve metrics along a single axis of variation, and that different growth curve metrics provide parallel measures of phage infectivity. Given the observed relationships between growth curve metrics, we sought to understand whether different portions of growth curve data are more or less informative for quantifying phage infectivity. When examining growth curves simulated with phages varying in their infection rates, we observed that the incline portion of the density dynamics is unaffected by infection rate (Fig 3A), or burst size or lysis time ( Fig S4). This suggests that phage infectivity generally does not alter bacterial population dynamics during the incline portion of the data. In contrast, the slope of the decline portion of the density dynamics is strongly correlated with the peak bacterial density (Figs 3B, S5). Thus, growth curve analyses should largely disregard the incline portion of the data, focusing on the more-informative peak and decline portions of the data for quantifying phage infectivity. We next sought to understand how growth curve metrics are related to the underlying phage and bacterial traits included in our model. In particular, it is unclear which traits have effects on growth curve metrics, and whether and how contributing traits interact with one another in determining the resulting metrics. To answer those questions, we simulated growth curves with phages with infection rates and lysis times along a grid of values, then calculated the peak time (of bacterial density) for each simulation. When plotted as a surface, we see that peak time is affected by both infection rate and lysis time (Fig 4). Moreover, infection rate and lysis time interact to produce a concave landscape, jointly determining peak time. Similar patterns are found for other phage ( Fig S6) and bacterial (Fig S7) traits in the model, as well as when considering other growth curve metrics (Figs S8, S9, S10). In all, these data show that growth curve metrics are determined by the interactions between all phage and bacterial traits, suggesting that growth curve metrics can act as a generalized measure of phage fitness. Given that growth curve metrics are tightly related to one another and driven by the interaction of bacterial and phage traits, we sought to understand how growth curve metrics relate to traditional measures of phage fitness. For instance, one common measure of phage fitness is the growth of the phage population over time. Indeed, past work has shown that the average phage growth rate during a growth curve has a negative relationship with bacterial extinction time both in vitro and in silico (49). We simulated growth curves using phages with a wide array of trait values, observing that our model reproduces this pattern (Figs 5A, S11A). Moreover, we now show how this relationship arises (Fig 5B,  S11B). Specifically, the final phage density is largely equal to the product of the phage's burst size and the peak bacterial density: once all the bacteria present at the peak are converted into phages, they comprise the bulk of the phages present at the end of the growth curve. At small burst sizes, this pattern is noisy; as burst size increases, it increasingly dominates the effect on final phage density, and the pattern becomes stronger. These findings suggest that the extinction time metric, and by extension all the other closely-related metrics (Fig 2), are correlated with conventional measures of phage fitness. We next sought to understand how relaxation of a key assumption of our model would alter growth curve dynamics. Thus far, our models were constructed under the assumption that bacterial susceptibility to phages is constant, both over time and within the bacterial population. Under such an assumption, the inevitable outcome of a phage-bacteria growth curve is bacterial extinction (Table S1). However, this assumption is unlikely to be realistic. In vitro bacterial susceptibility often declines as bacterial growth slows (52,(55)(56)(57)(58)(59), and experimental growth curves often exhibit incomplete lysis of bacterial populations that cannot be generated by our simple model ( (14), Fig S12). To assess the effects of relaxing this assumption, we used the two frequently-published approaches to incorporate variation in host susceptibility into our model (52). First, we modeled phage growth weakening as nutrients become scarce and bacterial growth slows (26,52,57,(60)(61)(62). Under these conditions, we find that bacterial populations can avoid extinction if they deplete nutrients sufficiently that phage infectivity becomes zero (Table S1). However, under none of these conditions are bacterial populations observed to experience partial lysis (Figs S14, S15, S16). Second, we introduced subpopulations of bacteria that vary in their susceptibility to phage infection to approximate intrapopulation variation (63). Under conditions where one sub-population is completely resistant to phage infection, the bacterial population avoids extinction by persisting in the resistant sub-population (Table S1). Moreover, transitions to a resistant sub-population can reproduce empirically-observed patterns of partial lysis of bacterial populations (Figs 6, S17, S18). These results suggest that slowing phage growth as bacterial growth slows can alter growth curve dynamics, but that partial lysis patterns are likely explained by the presence of sub-populations of bacteria with varying susceptibility to phage infection. Next, we set out to assess the sensitivities of growth curve metrics to the inoculum densities of bacteria or phages. We simulated growth curves with varying initial densities of bacteria or phages. We then compared how strongly growth curve metrics were affected by initial densities versus phage infectivity. We found that changes in the initial density of bacteria or phages can strongly alter growth curve metrics on the same order of magnitude as changes in infectivity (Figs 7, S19, S20, S21, S22). Thus, normalization of initial bacterial and phage densities is likely to be necessary when preparing a growth curve assay. In practice, this normalization can be subject to stochastic noise, where the true initial densities differ from the intended initial densities. Fortunately, the extent of this noise is minimal, and is only likely to obscure small differences in phage infectivity (Fig S23, S24, S25). Finally, we set out to test how to use growth curves to compare phage infectivity across different bacterial hosts. For instance, consider when different bacterial hosts vary in their susceptibility to a phage but also vary in their growth rate or carrying capacity. How can we eliminate the effects of differences in bacterial growth to infer differences in phage infectivity? One proposed approach is to use multivariate ordination methods like principal component analysis (PCA) to quantify phage infectivity (23,51). We simulated growth curves of varying bacteria with a phage and applied PCA, but found that the PCA axes did not correlate with phage infectivity (Fig S26). An alternate proposed approach is to use area under the curve (AUC) normalized by dividing by the AUC of a phage-free control (44,48,50). We simulated growth curves of varying bacteria with a phage and found that AUC generally correlates well with phage infectivity. However, under some conditions normalizing can counterintuitively worsen the signal of phage infectivity (Figs 8, S27). Specifically, when bacteria vary in their carrying capacities, the AUC of the phage-free control varies more than the AUC of the phage-bacteria curves, so normalizing worsens the signal of infectivity (Fig 8A). In contrast, when bacteria vary in their growth rates, the AUC of the phage-free control varies less than the AUC of the phage-bacteria curves, so normalizing improves the signal of infectivity (Fig 8B). These results suggest that normalizing AUC should only be done when phage-bacteria AUC values are expected to vary more than the control bacteria-only AUC values.

Discussion
Here, we set out to test how metrics of bacterial growth curves can be used to quantify phage infectivity and bacterial resistance. We found that many metrics of growth curves provide correlated information with each other on infectivity and resistance (Fig 2). Specifically, this information is provided in the peak and decline phases of the curves, whereas the initial growth phase is largely uninformative (Fig 3). We then showed that growth curve metrics act as an overall measure of phage fitness, incorporating effects from multiple phage and bacterial traits (Fig 4). Moreover, these metrics correlate with traditional measures of phage fitness (Fig 5). We tested how variation in infectivity can alter growth curve dynamics (Fig 6) and showed that growth curve metrics are largely robust to noise in initial conditions (Fig 7). Finally, we highlighted the limitations of proposed methods to compare infectivity and resistance across distinct bacterial hosts (Fig 8).
This work shows that growth curve metrics can be used to quantify general phage infectivity, but are unlikely to be useful to quantify specific phage life history traits like adsorption (cell-attachment) rate, burst size, or lysis time. We found no clear way to infer adsorption rate or lysis time from growth curve metrics. In principle, it could be possible to measure phage density after a growth curve and combine that information with peak bacterial density to estimate burst size (Fig 5B). However, in practice, single step phage growth curves remain likely to give more precise and accurate measures of specific phage traits.
Some patterns from these new findings also align with the results of past work in other study systems. For instance, our simulations show that the initial phase of growth curve data provides little information about the infectivity of a phage (Fig 3A). Similar findings have been reported in the epidemiology literature, where the initial growth of an epidemic provides little indication of the eventual severity of the epidemic (66,67). Additionally, our simulations recreated previously-reported patterns between phage growth rate and bacterial extinction time (Fig 5A). Indeed, our work expanded on these findings by showing how this pattern arises from an interaction between phage burst size and the peak of bacterial density (Fig 5B, S11B).
Using simulations, we were able to directly manipulate inoculation conditions and observe how they altered growth curve dynamics. However, in laboratory studies, normalizing phage inoculum densities can be difficult. Phages are commonly quantified by plaquing them on a bacterial host, but since phages vary in their plaquing efficiency across different hosts the true phage particle concentration may be far from the inferred concentration. Further work is necessary to improve on this shortcoming of plaquebased density estimation, perhaps exploring how alternate methods like qPCR or flow cytometry could be implemented for common-practice enumeration of phage density.
We showed that observed patterns of partial population lysis during growth curves (Fig S12) cannot be explained by declining phage infectivity as bacterial growth slows (Figs S14, S15, S16), but instead can be explained by subpopulations of bacteria that vary in their susceptibility to phage infection (Fig 6). Indeed, such mechanisms have previously been discovered and shown to underlie partial population lysis in some experimental growth curves (14). However, it remains unclear how to interpret partial lysis patterns with respect to phage susceptibility. Bacterial density after partial population lysis is driven not only by the phage infection rate, but also by the transition rate from susceptibility to resistance (Fig 6). Further work is needed to understand how these two traits are related to one another before a clear interpretation can be drawn from the density after partial population lysis.
Our work has established a broad foundation for using phage-bacteria growth curves to quantify phage infectivity and bacterial resistance, opening many avenues for future work. Given the limited theoretical work on phage-bacteria growth curves (this paper and (49)), future theoretical work could explore the effects of numerous biological phenomena, including stochasticity, coinfection exclusion, cooperation, failed infections, continuous intrapopulation trait variation, or lysogeny (phage integration into the host bacterial chromosome). Additionally, theoretical work could explore alternate uses of growth curves than quantifying phage infectivity, like quantifying the evolution of bacterial resistance to phage infection. At the same time, further empirical is needed. For instance, to understand the relationship between cell density and the measured proxies of cell density (e.g., optical density) (68)(69)(70), and to test theoretical predictions with phages whose life history traits have been experimentally quantified a priori.
In all, we have shown that measuring growth curves can be a powerful approach for quantifying phagebacteria interactions. Growth curves complement existing methods for measuring phage infectivity. They occur in liquid culture, which may better reflect the interactions of phages and bacteria in liquid (planktonic) environments than agar surface-based assays of infectivity. Moreover, growth curves are highly scalable, able to be collected in parallel for many replicates at once and easily implemented by automated platforms. Given that growth curves are already widely used heuristically in the phagebacteria literature, our findings suggest that they may be ripe for quantitative application from basic to applied questions.