The dual nature of bacteriophage: growth-dependent predation and generalised transduction of antimicrobial resistance

Bacteriophage (“phage”) are both predators and evolutionary drivers for bacteria, notably contributing to the spread of antimicrobial resistance (AMR) genes by generalised transduction. Our current understanding of the dual nature of this relationship is limited. We used an interdisciplinary approach to quantify how these interacting dynamics can lead to the evolution of multi-drug resistant bacteria. We co-cultured two strains of Methicillin-resistant Staphylococcus aureus, each harbouring a different antibiotic resistance gene, with 80α generalized transducing phage. After a growth phase of 8h, bacteria and phage surprisingly coexisted at a stable equilibrium in our culture, the level of which was dependent on the starting concentration of phage. We detected double-resistant bacteria as early as 7h, indicating that transduction of AMR genes had occurred. We developed multiple mathematical models of the bacteria and phage relationship, and found that phage-bacteria dynamics were best captured by a model in which the phage burst size decreases as the bacteria population reaches stationary phase, and where phage predation is frequency-dependent. We estimated that one in every 108 new phage generated was a transducing phage carrying an AMR gene, and that double-resistant bacteria were always predominantly generated by transduction rather than by growth. Our results suggest a shift in how we understand and model phage-bacteria dynamics. Although rates of generalised transduction could be interpreted as too rare to be significant, they are sufficient to consistently lead to the evolution of multi-drug resistant bacteria. Currently, the potential of phage to contribute to the growing burden of AMR is likely underestimated.


33
To counter the rapidly increasing global public health threat of antimicrobial resistance (AMR), we 34 must urgently develop new solutions 1 . "Phage therapy" is one such tool which has recently seen a 35 renewed interest 2 . This relies on using bacteriophage (or "phage"), major bacteria predators and the 36 most abundant biological entities on the planet 3 , as antimicrobial agents. However, phage are also 37 natural drivers of bacterial evolution through horizontal gene transfer by "transduction" 4,5 . AMR 38 genes can be transferred by transduction at high rates, both in vitro and in vivo [6][7][8] . The dual nature of 39 phage (predation and transduction) makes them a double-edged sword in the fight against AMR, as 40 they are themselves capable of contributing to the spread of the problem they aim to solve, yet our 41 understanding of these dynamics and how to best represent them is limited. 42 There are two types of transduction; here, we focus on "generalised transduction", which occurs 43 during the phage lytic cycle, when non-phage genome DNA is mistakenly packaged in a new phage 44 particle ( Figure 1). The resulting transducing phage released upon lysis can then inject this genetic 45 material into another bacterium. Current guidelines for phage therapy recommend that exclusively 46 lytic phage should be used, removing the risk of the second type of transduction which relies on 47 bacteria not exposed or previously exposed to phage. The starting concentration of both single-158 resistant S. aureus parent strains (BE to erythromycin & BT to tetracycline) was 10 4 colony-forming 159 units (cfu) per mL, and the starting concentration of exogenous phage 80 (PL) was 10 4 plaque-forming 160 units (pfu) per mL. Double-resistant progeny (BET) are generated by transduction. The initial co-culture 161 was diluted in fresh media after 24h, to form a new co-culture with bacteria previously exposed to 162 phage. Phage were added in the new co-culture to reach a concentration of 10 4 pfu/mL. Error bars 163 indicate mean +/-standard error, from 3 experimental replicates. Where S1(t) and S2(t) represent the number of bacteria (in cfu/mL) from the chosen strains 1 and 2, 177 at times t = 0 or 24 hours. 178 BE did not show a significant fitness cost relative to BT over 24h of growth (mean relative fitness 1.02, 179 sd 0.03). The double-resistant progeny BET did not show a significant fitness cost relative to either 180 single-resistant parent strain (mean relative fitness to BE: 0.96, sd 0.06; mean relative fitness to BT: We obtained maximum growth rate estimates by fitting a logistic growth model to the in vitro data. 183 This is shown in Equation 2, where max is the maximum bacterial growth rate, Bmax is the carrying 184 capacity, and B the concentration of bacteria (in cfu/mL). 185 The median estimated maximum growth rates were 1.61 for BE (95% credible interval 1.59-1.63), 1.51 187 for BT (1.49-1.53) and 1.44 for BET (1.42-1.47), with a total carrying capacity of 2.76 x 10 9 cfu/mL (2.61 188 x 10 9 -2.98 x 10 9 ). 189 Investigation of possible phage-bacteria interactions using a flexible 190 modelling framework 191 Model structure 192 We designed a mathematical model to reproduce the in vitro phage-bacteria dynamics, including 193 generalised transduction of resistance genes. During our experiment, our co-culture contained up to 194 three strains of bacteria: the two single-resistant parents (BE, BT) and the double-resistant progeny 195 (BET). Although we were only able to count lytic phage (PL), based on the biology of generalised 196 transduction ( Figure 1) we know that there were also transducing phage carrying either the 197 erythromycin resistance gene (PE), or the tetracycline resistance gene (PT). Since we did not detect any 198 evidence of 80 lysogeny in our co-culture after 24h, we did not include this feature in the model. The 199 corresponding model diagram is shown in Figure 4a. The complete model equations can be found in 200 Methods. 201 Using this modelling framework, we explored a combination of different phage-bacteria interactions, 202 described below (Figure 4b-c). By fitting the models to our experimental data, we could rule out 203 certain interactions and suggest the best model to reproduce the phage-bacteria dynamics seen in 204 vitro. 205 Equilibrium analyses for the density and frequency-dependent models 250 Despite these being common methods to represent phage-bacteria interactions in mathematical 251 models, previous analyses have suggested that the density and frequency-dependent interactions 252 alone cannot capture the equilibrium levels we and others have seen 18,35 . We explore this in the 253 context of our own in vitro data using equilibrium analyses. 254 Assuming that transduction and the phage latent period are negligible, a simplified model 255 representing phage predation as a density-dependent process is shown in equations (6) and (7). 256 Where max is the maximum bacterial growth rate, Bmax is the carrying capacity,  is the phage 259 adsorption rate,  is the phage decay rate, and  is the phage burst size, with ' equal to  -1. To solve 260 for equilibrium (i.e. B = P = 0), equations (6) and (7) can be rewritten as equations (8)  Substituting this into equation (15) leads to a value for phage decay rate () (equation (17)). 283 = 10 9 * 10 −5 * ′ ≈ 10 4 * ′ (17) 284 Giving rise to the condition that the phage decay rate  must be approximately 10 4 times greater than 285 the burst size '. As at least one phage must be released upon bursting, the burst size ' is greater than 286 1. However, the phage decay rate , which represents the proportion of phage inactivated during one 287 time-step, must be less than 1 and hence this is impossible. 288 As for the frequency-dependent interaction, a simplified model of this process is shown in equations 289 (18) and (19). 290 With the condition that B = 0 equation (18) can be rewritten as equation (20). This condition must 293 hold in order for there to be a non-zero equilibrium. 294 Using our equilibrium values (B = 10 9 , P = 10 5 , Bmax = 2.8 x 10 9 , max = 1.52), this is equivalent to equation 296 Therefore,  must be approximately equal to 6733 to obtain a non-zero equilibrium as seen in our in 313 vitro data with our model fitted values (biologically plausible ones). However, the definition of  Even though these analyses rely on a simplified set of equations, using realistic parameter values we 316 have shown that a non-zero equilibrium, as we have seen in vitro, cannot be reproduced using models 317 with only a density or frequency-dependent interaction. Instead, phage-bacteria co-existence may be 318 explained by variations in phage predation parameters depending on bacterial resources availability, 319 or bacterial growth rate 14,17-22 . However, to the best of our knowledge a simple mathematical 320 expression linking phage predation to bacterial growth has not yet been developed. 321 322 Second phage-bacteria interaction: dependence of phage predation on bacterial growth 323 Here, we consider that a decrease in bacterial growth as bacteria reach stationary phase could firstly 324 affect the phage adsorption rate , due to changes in receptors on bacterial surfaces, which affect 325 opportunities for phage to bind ( Figure 4c). Secondly, this could affect phage production, and thus 326 burst size , as phage replication relies on bacterial processes and may decrease when bacterial 327 growth slows down (Figure 4c). Using a single phage predation multiplier, with the same principle of 328 logistic growth applied to bacteria, we allow either or both  and  to decrease as bacterial growth 329 decreases in our model (equations (32) and (33)). 330 335 reproduce the in vitro dynamics 336 Overall, we considered 6 different models, either density-or frequency-dependent, and with either or 337 both the phage adsorption rate and burst size linked to bacterial growth. Note that we did not include 338 a phage decay rate in these models, as this did not affect the dynamics of the system over 24h, for a 339 wide range of decay rates ( All models successfully reproduced the trends seen in vitro when the phage were started at either 10 3 346 and 10 4 pfu/mL (Figure 5a-b). However, only the two models where only phage burst size decreases 347 as the bacteria population approaches carrying capacity were able to reproduce the increase in phage 348 numbers seen in the later hours of the 10 5 pfu/mL dataset, despite all models having been fitted to 349 this dataset (Figure 5a-b). This was confirmed by calculating the average Deviance Information Criteria 350 (DIC) value for the models, which favours best-fitting models while penalising more complex models 351 (i.e. with more parameters) 36 . The two models where only phage burst size decreases as the bacteria 352 population approaches carrying capacity had the lowest DIC values, indicating that they were the 353 better-fitting models (Table 1). 354

Identification of the best-fitting phage-bacteria interactions to
Our initial experiments considered the dynamics over 24h for varying phage starting concentrations. 355 To test the ability of our model to recreate the dynamics under changing bacterial levels, we replicated 356 our transduction co-culture experiments with starting concentrations of 10 6 cfu/mL bacteria instead 357 of 10 4 cfu/mL, varying the starting phage concentration between 10 4 and 10 6 pfu/mL, and measuring bacteria and phage numbers after 24h of co-culture. We then used the estimated parameter values 359 (Table 1) to try to reproduce these 24h numbers of bacteria and phage. 360 Increasing the starting phage concentration led to an increase in the number of phage after 24h 361 (Figure 5c). For a starting phage concentration between 10 4 and 10 6 pfu/mL, increasing starting phage 362 numbers did not affect single-resistant parents BE and BT numbers after 24h, but led to a progressive 363 increase in double-resistant progeny BET numbers. Increasing starting phage numbers above 10 6 364 pfu/mL caused bacteria numbers after 24h to decrease. 365 Using the estimated parameter values (Table 1) with the model where only burst size is linked to 366 bacterial growth, we see that the density model cannot reproduce these dynamics as it predicts that 367 bacteria become extinct rapidly (Figure 5c). The frequency-dependent model is able to reproduce 368 these trends, but fails to recreate the exact same numbers of phage and bacteria, predicting a decline 369 in bacterial levels when the starting phage concentration increases above 10 5 pfu/mL, a lower 370 threshold than seen in the data (Figure 5c). The same overall trends are seen for the models where 371 only the adsorption rate is linked to bacterial growth, or both adsorption rate and burst size (

411
Parameter estimates for our best-fitting model (with a frequency-dependent interaction and a link 412 between phage burst size and bacterial growth only) suggest that the adsorption rate is 2.3 x 10 -10 413 (95% credible interval: 2.1 x 10 -10 -2.4 x 10 -10 ) which is the smallest estimate from the models (Table  414 1). On the other hand, the estimated burst size is relatively large at 76 (70 -83) phage, and is higher 415 than a previous in vitro estimate for 80 of 40 37 . However, due to the decrease in burst size when 416 bacteria are in stationary phase, we expect that this number would change depending on the 417 conditions under which it is measured (Figure 6a). Finally, the estimated latent period of 0.72h (0.69 -418 0.77) is slightly longer than a previous in vitro estimate of 0.67h 37 . Regarding the other models, we 419 note some biologically unlikely parameter estimates which further suggest that these models are 420 inappropriate, such as the low burst size for the models with only the adsorption rate linked to 421 bacterial growth (12 (10 -14) and 10 (8 -12)), or the high latent period for the models with both 422 adsorption rate and burst size linked to bacterial growth (0.93 (0.86 -0.99) and 0.88 (0.79 -0.96)) 423 (Table 1). rather than by growth (Figure 6d). This is because DRP only appear after 2 to 4h, while after 4h 436 bacterial growth rate starts decreasing as the total bacteria population approaches carrying capacity 437  (Table 1). (a) Phage burst size over time, by starting phage 443 concentration. As bacteria reach stationary phase after 8h, phage burst size decreases. In the 10 5 444 dataset, we see that burst size is predicted to increase again after 20h. This is due to bacterial numbers 445 decreasing as bacteria are being lysed by phage. Density-dependent models have been compared to data at less fine time scales (e.g. daily time points) 466 or over smaller time periods (e.g. less than 8h), where they were able to reproduce in vitro values 467 from experiments in chemostats, and have been helpful to improve our basic understanding of phage-468 bacteria dynamics 14-16 . However, here we show that this type of interaction is not able to reproduce 469 finer hourly dynamics, and does not perform consistently when varying concentrations of starting 470 phage and bacteria. Using this, alongside a critique of the mathematical implications of this process, 471 we argue that density-dependence is not a biologically accurate representation of phage predation, 472 as it fails to reproduce these dynamics at high or low numbers of phage and bacteria, which would 473 correspond to scenarios potentially seen during phage therapy. 474 Our work adds to the growing body of evidence that phage predation depends on bacterial growth 475 14,[17][18][19][20][21][22][23] . This has implications for antibiotic-phage combination therapy, as it suggests that 476 bacteriostatic antibiotics, which prevent bacterial growth, could reduce phage predation. This effect 477 has been previously seen in S. aureus 38 . In the environment, including in persistent infections, bacteria 478 spend most of their time in stationary phase 39 . This suggests that bacteria and phage may be able to 479 co-exist for prolonged periods of time in a broad range of settings, without the phage systematically 480 eradicating the bacteria. Under such conditions, phage may mediate horizontal gene transfer by particularly relevant for S. aureus, since approximately 20% of humans are colonised asymptomatically 483 by this bacterium at any given time 40 , and at least 50% of these carriers may also carry phage capable 484 of generalised transduction 41 , suggesting a constant background evolution rate for S. aureus in the 485 human population. Combined with environmental exposure to antibiotics which acts as a selective 486 pressure, this may contribute to the risk of multidrug resistant bacteria evolution. 487 Our experimental design is both a strength and a limitation of our study. Since we jointly designed the 488 experiments and models, we are confident that we have included in our mathematical model all the 489 organisms and interactions present in vitro. We are therefore confident in the conclusions on model 490 structure, however, the usage of such a specific experimental system with two bacterial strains of the 491 same genetic background and one phage limits the generalisability of our parameter values, as these 492 will likely vary for different bacteria and phage. Growth conditions will likely also differ between the 493 in vitro environment studied here, and in vivo conditions. Here, our model assumes that phage do not 494 decay, that bacteria do not become resistant to phage, and that they can grow indefinitely as they are 495 observed in a rich medium for 24h only, but over longer periods of time it may be necessary to revisit 496 these assumptions 42 . Finally, we assumed that the proportion of transducing phage created was 497 independent of the gene being transduced (ermB, on the bacterial chromosome, or tetK, on a 498 plasmid). This was supported by preliminary work (not shown), but should be further investigated to 499 improve our understanding of the factors that can facilitate or prevent transduction of different genes. 500 To answer all of these questions, future work should investigate both phage predation and 501 transduction dynamics over longer time periods, with different strains of bacteria and phage. 502 All our models captured certain aspects of the trends seen in vitro, but also underestimated phage 503 numbers between 5-7h by up to 20 times. This is likely a consequence of our experimental design. To 504 count lytic phage, we centrifuged and filtered the co-culture to remove bacteria. This we were able to fit our model to the dynamics for two MOI (0.1 and 10), and replicate those of a third 517 (1). However, when trying to use the same model for these same three MOI, but with a starting 518 bacterial concentration to 10 6 , we found differences between our model and values seen after 24h. 519 This indicates that MOI is not appropriate to summarise all the complexity of the underlying phage-520 bacteria dynamics. Future experimental studies should express their results as a function of their 521 starting concentration of phage and bacteria, not just MOI. 522 In any case, the failure of our model to replicate 24h values with a different starting bacteria 523 concentrations shows that, whilst we have reduced the model structure uncertainty, we are still not 524 fully capturing the phage-bacteria interaction. Currently, our model predicts that, for a starting 525 concentration of 10 6 bacteria, a starting concentration of 10 5 phage or more will be enough to cause 526 a decrease in bacterial numbers after 24h, while our data shows that the starting concentration of 527 phage must be higher than 10 6 for this to happen. In vitro, it is likely that slower bacterial growth 528 simultaneously affects the phage adsorption rate, latent period and burst size, each to varying extents 529 14,[17][18][19][20][21][22][23] . This would explain why we would need a higher starting concentration of phage for a higher 530 starting concentration of bacteria, to exert a strong enough predation pressure before bacteria reach 531 stationary phase, causing a reduction in phage predation. However, here we have only made the first 532 step in this process, having linked the burst size linearly to the bacterial growth rate, instead of trying 533 to link different phage predation parameters to bacterial growth using different functions. These 534 complexities need to be explored further, supported by in vitro work measuring phage predation 535 parameters at various time points. In S. aureus, wall teichoic acid (WTA) is the phage receptor 44,45 . 536 Lack of WTA glycosylation has been shown to induce phage resistance 46 , and changes in WTA 537 structure at different growth phases may be possible, since one of the genes involved in its synthesis 538 is repressed by a quorum sensing system 47 . However, to the best of our knowledge this has not yet 539 been investigated. 540 Despite being recognised as a major mechanism of horizontal gene transfer, thus far there have been 541 limited mathematical modelling studies on the dynamics of transduction of AMR 13 . Using our model, 542 we are able to estimate numbers of transducing phage which we cannot count in vitro, and see that 543 approximately 1 generalised transducing phage is generated per 10 8 lytic phage, consistent with 544 previous estimates 48, 49 . Here, we show that this number, which may seem insignificant, is enough to 545 consistently lead to the successful horizontal gene transfer of AMR, resulting in DRP after only 7h, 546 substantially less than the usual duration of antibiotic treatment. We also show that transduction is 547 the dominant mechanism to create new DRP throughout the entire experiment, rather than growth 548 of existing DRP. This echoes the conclusions of previously published work on the importance of 549 transduction, including in vivo experiments and with other Staphylococcus species 6,7,29,50 . 550 Our findings suggest that transduction is currently under-emphasised in the exploration of phage-551 bacteria dynamics. Future studies on this topic should not assume that transduction can be dismissed 552 by default, but instead investigate whether it is relevant in their system. This requires further in vitro 553 and in vivo monitoring to identify scenarios where transduction plays a significant role in the transfer 554 of AMR genes, likely depending on the environment, and characteristics of the bacteria and phage 555 present. This will require new experimental designs, since counting phage numbers can be difficult, 556 notably with clinical strains of bacteria. This should also be investigated in the presence of antibiotics, 557 where the importance of selection enters, increasing the fitness of the small numbers of DRP 558 generated by transduction. 559 In conclusion, the dual nature of phage (predation and transduction) leads to complex interactions 560 with bacteria. These dynamics must be clarified, to correctly evaluate the extent to which phage 561 contribute to the global spread of AMR. We must also understand this dual nature to guarantee a safe 562 design of phage therapy. Otherwise, ignoring transduction may lead to worse health outcomes in 563 patients if phage contribute to spreading AMR instead of overcoming it. Interdisciplinary work will be 564 essential to answer these urgent public health questions in the near future. 565

566
All analyses were conducted using the statistical software R 51 . The underlying code and data are 567 available in a GitHub repository: https://github.com/qleclerc/mrsa_phage_dynamics. 568

569
Strains and phage used 570 The Staphylococcus aureus parent strains used for our transduction experiment were obtained from 571 the Nebraska Transposon Mutant Library 52 . These were strain NE327, carrying the ermB gene 572 encoding erythromycin resistance and knocking out the 3 integrase gene, and strain NE201KT7, a 573 modified NE201 strain with a kanamycin resistance cassette instead of the ermB gene knocking out 574 the 2 integrase gene, and a pT181 plasmid carrying the tetK gene encoding tetracycline resistance 53 . 575 Growing these strains together in identical conditions as our co-culture below, but without the 576 addition of exogenous phage, does not lead to detectable horizontal gene transfer (HGT; data not 577 shown). To enable HGT, exogenous 80 phage was used, a well-characterised temperate phage of S. 578 aureus capable of generalised transduction 33 . To count lytic phage, S. aureus strain RN4220 was used, 579 a restriction deficient strain highly susceptible to phage infection 54 . 580 581 Transduction co-culture protocol 582 Pre-cultures of NE327 and NE201KT7 were prepared separately in 50mL conical tubes with 10mL of 583 Brain Heart Infusion Broth (BHIB, Sigma, UK), and incubated overnight in a shaking water bath (37°C, 584 90rpm). The optical densities of the pre-cultures were checked at 625nm the next day to confirm 585 growth. The pre-cultures were diluted in phosphate-buffered saline (PBS), and added to a glass bottle 586 of fresh BHIB to reach the desired starting concentration in colony forming units per mL (cfu/mL) for each strain, forming a master mix for the co-culture. CaCl2 was added at a concentration of 10mM to 588 the master mix. Phage 80 stock was diluted in phage buffer (50 mM Tris-HCl pH 7.8, 1 mM MgSO4, 589 4 mM CaCl2 and 1 g/L gelatin; Sigma-Aldrich), and added to the master mix to reach the desired 590 starting concentration in plaque forming units per mL (pfu/mL). Ten 50mL conical tubes were prepared 591 (one co-culture tube for each timepoint, from 0 to 8h and 16 to 24h), each with 10mL from the master 592 mix. Each co-culture tube was then incubated in a shaking water bath (37°C, 90rpm) for the 593 corresponding duration. 594 Bacteria counts for each timepoint were obtained by diluting the co-culture in PBS before plating 50μL 595 on selective agar, either plain Brain Heart Infusion Agar (BHIA, Sigma, UK), BHIA with erythromycin 596 (Sigma, UK) at 10mg/L, BHIA with tetracycline (Sigma, UK) at 5mg/L, or BHIA with both erythromycin 597 and tetracycline (10mg/L and 5mg/L respectfully). Note we plated 500μL instead of 50 on the plates 598 with both antibiotics, to increase the sensitivity of the assay. This allowed distinction between each 599 parent strain, resistant to either erythromycin or tetracycline, and the double resistant progeny (DRP) 600 generated by transduction. Plates were then incubated at 37°C for 24h, or 48h for plates containing 601 both antibiotics. Colonies were counted on the plates to derive the cfu/mL in the co-culture for that 602 timepoint. 603 Lytic phage counts for each timepoint were obtained using the agar overlay technique 55 . Briefly, the Polymerase chain reaction protocols 613 To confirm that DRP contained both the ermB and tetK genes, primers ermBF (5′-614 CGTAACTGCCATTGAAATAGACC-3′), ermBR (5′-AGCAAACTCGTATTCCACGA-3′), tetKF 615 (5′-ATCTGCTGCATTCCCTTCAC-3′), and tetKR (5′-GCAAACTCATTCCAGAAGCA-3′) were used. Strains 616 NE327 (only containing ermB) and NE201KT7 (only containing tetK) were used as positive and negative 617 controls. 618 To confirm that 80 lysogeny did not occur in our co-culture, we applied a previously published 619 method 33 and used a combination of four primers: SaRpmF (5′-GACTGAATGCCCAAACTGTG-3′) in the 620 S. aureus rpmF gene, SMT178 (5'-GGCTGGGAATTAATGGAAGATG-3′) in the 80 integrase, SaSirH (5′-621 TTAAGTAGCATCGTTGCATTCG-3′) in the S. aureus sirH gene, and SMT179 (5′-622 GAGTCCTGTTTGCGAATTAGG-3′) in the 80 ORF73 region. SaRpmF and SMT178 were used to amplify 623 the left prophage junction (attL), SaSirH and SMT179 to amplify the right junction (attR), and SaRpmF 624 and SaSirH to amplify the bacterial insertion site (attB). RN4220 was used as a negative control for 625 lysogeny, and JP8488, an RN4220 strain lysogenic for 80, was used as a positive control (obtained 626 from José Penadés and Nuria Quiles, Imperial College London). 627 All PCRs were conducted using OneTaq Hot Start Quick-Load 2X Master-mix, following the 628 manufacturer's protocol. Tested samples were homogenised in 20µl nuclease-free water, and 1.5µl of 629 each suspension was used as template for a total reaction volume of 25µl. To have a decline in the bacterial population, we must therefore satisfy the condition stated in 696 equation (48).
Finally, by subtracting 1 and multiplying by -2.8 x 10 9 , we obtain equation (50). 701 B > −( 1 1.5 − 1) * 2.8 * 10 9 (50) 702 Solving equation (50) leads to the condition that the number of bacteria must be greater than 9.3 x 703 10 8 in order for bacterial growth to be low enough that bacteria numbers decrease. In other words, it 704 is impossible for the effect of phage predation to decrease bacteria numbers below 9.3 x 10 8 . We 705 overcome this by applying a correction term to equation S1, leading to equation (51). This is equivalent to saying that bacteria infected by phage cannot replicate and hence is more 708 biologically realistic. 709 710 Link between bacterial growth and phage predation 711 We consider that reduced bacterial growth can lead to decreased phage predation, through reduced 712 adsorption () and/or burst size (). Equations (52) and (53) allow these parameters to decrease as 713 bacterial growth decreases, using the same principle of logistic growth as seen in equation (34). algorithm. For every iteration, this algorithm slightly changes the parameter values, runs the model, 723 assesses the resulting model fit to the data, and accepts or rejects these new parameter values based 724 on whether the model fit is better or worse than with the previous set of values. We run the algorithm 725 with two chains, and once convergence has been reached (determined using the Gelman-Rubin 726 diagnostic, once the multivariate potential scale reduction factor is less than 1.2 56 ), we generate 727 50,000 samples from the posterior distributions for each parameter. 728 In a first instance, we used our growth co-culture data, where phage are absent, to calibrate the 729 bacterial growth rate parameters max for each bacteria strain  ( ∈ {E, T, ET}), as well as the carrying 730 capacity Nmax using a simple logistic growth model (equation (9)). All other parameters related to 731