Classification of In Vitro Phage–Host Population Growth Dynamics

The therapeutic use of bacteriophages (phage therapy) represents a promising alternative to antibiotics to control bacterial pathogens. However, the understanding of the phage–bacterium interactions and population dynamics seems essential for successful phage therapy implementation. Here, we investigated the effect of three factors: phage species (18 lytic E. coli-infecting phages); bacterial strain (10 APEC strains); and multiplicity of infection (MOI) (MOI 10, 1, and 0.1) on the bacterial growth dynamics. All factors had a significant effect, but the phage appeared to be the most important. The results showed seven distinct growth patterns. The first pattern corresponded to the normal bacterial growth pattern in the absence of a phage. The second pattern was complete bacterial killing. The remaining patterns were in-between, characterised by delayed growth and/or variable killing of the bacterial cells. In conclusion, this study demonstrates that the phage–host dynamics is an important factor in the capacity of a phage to eliminate bacteria. The classified patterns show that this is an essential factor to consider when developing a phage therapy. This methodology can be used to rapidly screen for novel phage candidates for phage therapy. Accordingly, the most promising candidates were phages found in Group 2, characterised by growth dynamics with high bacterial killing.


Introduction
Bacteriophages (phages) are viruses that specifically infect bacteria. They are estimated to be the most abundant organisms on Earth (~10 31 entities) and play a major role in shaping the microbial communities [1,2]. Phages are unable to replicate independently of a susceptible bacterial host, and their host range is determined by a combination of various factors, including host-binding protein specificity and bacterial phage-resistance mechanisms [3,4]. Bacteria can readily evolve resistance to phage infections though different mechanisms, which can result in distinct resistance phenotypes [5]. These can differ in whether the resistance is partial or complete, in the level of fitness cost associated with resistance, and in whether the mutation can be countered by the infecting phage. Consequently, these differences determine the effect of the phage infection on the bacterial population dynamics and the resulting community structure [6,7].
The therapeutic use of phages (phage therapy) represents an urgently needed alternative or supportive antibacterial agent to antibiotics to control bacterial pathogens [8][9][10][11]. and the Comprehensive Antibiotic Resistance Database (CARD) [28]. Virulence genes were identified using ABRicate with data from the Ecoli_VF database. Prophage regions were identified using the PHAge Search Tool Enhanced Release (PHASTER) tool [29]. The Rapid Annotation using Subsystem Technology (RAST) server and the SEED viewer, version 2.0 [30] were used for identification of the coding sequences (CDSs) and initial annotation of the APEC genomes.

Bacteriophage Isolation, Purification, and Enumeration
A total of 18 lytic coliphages were used in this study. These were selected from the collection based on their genomic diversity. Phages were isolated from poultry faecal material using E. coli laboratory strain K514, sequenced, and identified as described before [31]. Phage stocks were stored at titres of 1.2 × 10 7 to 4.5 × 10 10 plaque-forming units (PFU)/mL at 4 • C. Working stocks used for phage infectivity and phage-host growth experiments were kept at titres of~10 10 PFU/mL.

Phage Infectivity and Phage-Host Growth Dynamics
Bacterial overnight cultures were used, and the cell concentration was adjusted tõ 10 8 colony-forming units (CFU)/mL for every experiment. Bacterial solutions were inoculated with 20-µL phages yielding initial MOIs of 10, 1, or 0.1. All bacterial reduction curves were generated using 96-well plates with working volumes of 200 µL. Experiments were performed in triplicate. For the experiments of the susceptible combinations, the experiment was performed on a duplicate plate at another time. A well of phage-free bacterial culture and a well of bacteria-free phage culture were included on every plate as control experiments in addition to one media blank for reference. The optical density (OD) was measured for a wavelength of 600 nm (OD600) with the Thermo Fisher Scientific Multiskan GO Microplate Spectrophotometer, v1.01.12, and the data were recorded using SkanIt software, v6.0.2.3. The OD600 was measured with fast measurement mode and no pathlength correction or use of transmittance, and the measurements were taken immediately after inoculation and then at regular intervals of 30 min afterward for 22 h. The incubation temperature was 37 • C, and shaking was continuous at a medium speed.
Growth curves were obtained by plotting OD600 values after baseline adjustment against time. Phage infectivity was defined based on endpoint measurements. Successful phage infection was defined as OD600: <0.2, somewhat successful infection was defined as OD600: 0.2-0.5, and failed phage infection as OD600: >0.5. The phage-host growth dynamics were assessed based on measurements throughout the experiment.

Assessing the Effect of Phage Species, APEC Strain, and MOI
The two different response variables, PhageScore [32] and local virulence score [16], were derived to assess the effect of the phage species, APEC strain, and MOI on the growth dynamics. For the PhageScore method, the area under the growth curve (AUC) was determined for the complete study period for each treatment combination and for the corresponding control, i.e., without phage, and the ratio of the difference between the control and the treatment over the control (multiplied with 100) was calculated. For the local virulence score, the same was done, but only measurements until the stationary phase in the control group were taken into consideration, i.e., until the timepoint before the maximum OD value. A fixed effects model with a normally distributed error term was used, and the phage, bacterium, and MOI, as well as all two-way interactions, were included in the model. F-tests were performed at the 5% significance level to assess the effects of the different factors. Finally, the Pearson correlation coefficient between the PhageScore and the virulence score was calculated.

Classification of Phage-Bacterium Growth Dynamics
To classify the growth dynamics of phage-host interactions based on the OD measurements, we applied two statistical data mining techniques: nonmetric multidimensional scaling (NMDS) ordination and principal component analysis (PCA), with hierarchical agglomerative clustering. A NMDS ordination plot (Euclidean distance) using the vegan package in R (http://www.R-project.org/ (accessed on 18 April 2021)) [33] was applied to quantify and visualise the pairwise dissimilarity between samples of each timepoint. The stress value was determined to access how well the data were transformed. A stress value between 0.02 and 0.01 was considered an acceptable fit, and <0.01 was considered a good fit. Metadata was included using the envfit function to determine the effect of each factor (phage species, bacterial strain, and MOI). Using the factoextra package the right number of groups (of the growth dynamics pattern) was determined using bootstrap values = 100 and hierarchical agglomerative clustering ("ward.D" method). The results were visualised using the ggplot2 package. A PCA was performed as verification (sanity check) of how much of the variability in the data is explained by each factor and how much of the total variability is captured. Subsequent subclustering of the designated groups was performed as described above for the complete dataset.

Repeatability of Group Assignments
The dplyr package in R was used to determine if replicates of the same phagebacterium-MOI combination were placed in the same group and subgroup.

APEC Strains
Whole-genome sequencing (WGS) of the bacterial genomes yielded a total of 2,673,858-4,825,608 paired-end reads for each of the 10 isolates, with an average coverage of 77-142-fold. The characteristics based on the WGS analysis are summarised in Table 1.

Coliphages
The coliphages used in this study were previously characterised [31]. The characteristics are summarised in

Phage Infectivity
The infectivity of the 18 coliphages against each of the 10 APEC strains at MOI 10, 1, and 0.1 is shown in Figure 1. The levels of infectivity of the tested coliphages were highly variable, varying between 0% and 100% of the tested APEC strains. Variations in the degree of infection (successful, somewhat successful, or failed) were detected, with the most successful infections at MOI 10 and least at MOI 0.1. Phage P1 was shown to be the most infective, as it was the only one able to infect all APEC strains at all MOI tested, except for B7-MOI 0.1. Phages P2, P3, and P9 were able to infect nine of the APEC strains, excluding B4, B4, and B1, respectively. Phage P5 was able to infect eight of the APEC strains, excluding B1 and B4. Phage P8 was able to infect six of the APEC strains, not including B1, B4, B6, and B7. Phages P4, P6, and P7 were able to infect five of the APEC strains, excluding B1, B4, B6, B7, and B8. The nine phages P10-18 were not able to infect any of the 10 APEC strains.

Classification of Phage-Bacterium Growth Dynamics
The phage-host growth dynamics were classified based on the OD measurements from a total of 2869 phage-host combination experiments. The NMDS analysis grouped the data into three different groups ( Figure 4). Group 1 comprised resistant phage-host combinations characterised by bacterial growth, Group 2 comprised fully susceptible combinations characterised by bacterial killing, and Group 3 comprised in-between combinations. A two-dimensional plot was considered appropriate, as the generated stress value was below 0.01 (0.007). Growth dynamics curves associated with each combination experiment and phage-free control culture are shown in Supplementary Figure S1. The principal component analysis (PCA) captured a total of 97% of the variance and showed similar groupings (Supplementary Figure S2). The Pearson's correlation coefficient between the local virulence score and the PhageScore was equal to 0.94.

Classification of Phage-Bacterium Growth Dynamics
The phage-host growth dynamics were classified based on the OD measurements from a total of 2869 phage-host combination experiments. The NMDS analysis grouped the data into three different groups ( Figure 4). Group 1 comprised resistant phage-host combinations characterised by bacterial growth, Group 2 comprised fully susceptible combinations characterised by bacterial killing, and Group 3 comprised in-between combinations. A two-dimensional plot was considered appropriate, as the generated stress value was below 0.01 (0.007). Growth dynamics curves associated with each combination experiment and phage-free control culture are shown in Supplementary Figure S1. The principal component analysis (PCA) captured a total of 97% of the variance and showed similar groupings (Supplementary Figure S2).
NMDS on a two-dimensional graph was applied to investigate how the factors (bacterial strain, phage species, and MOI) affect the groupings ( Figure 5). The 10 phages P4 and P10-18 were associated with bacterial growth (Group 1), of which P11 had the greatest association (similar to the phage-free controls), followed by P17, P4, P16, P13, P18, P10, P15, P12, and P14 ( Figure 5A). Accordingly, the greatest bacterial growth was associated with the Demerecviridae phage belonging to the Tequintavirus genus. Myoviridae phages belonging to the Felixounavirus genus and Drexlerviridae phages belonging to the Guelphvirus or Warwickvirus genera were also associated with bacterial growth. Phages associated with more bacterial killing (Groups 2 and 3) included the eight phages P1-3 and P5-9, of which P1 had the strongest association, followed by P5, P9, P2, P3, P8, P6, and P7. A few P4 combinations were also found in Group 3. Accordingly, the greatest bacterial killing was associated with Myoviridae phages belonging to the Tequatrovirus genus. Myoviridae phages belonging to the Mosigvirus genus were also associated with bacterial killing, as well as two out of three Hanrivervirus phages (Drexlerviridae family). The five bacterial strains B1, B4, B6, B7, and B8 were associated with bacterial growth (Group 1). This included all three O1 serotype strains (B1, B4, and B8), as well as two O78 serotype strains ( Figure 5B). The other five strains: B2, B3, B5, B9, and B10, were associated with more bacterial killing and included two O2 serotype strains (B5 and B10) and three O78 serotype strains. The overall effect of the bacterial factor on the plot was less compared to the effect of the phage factor. A high MOI (MOI 10) was associated with more bacterial killing, and a low MOI (MOI 0.1) was associated with less killing ( Figure 5C).  NMDS on a two-dimensional graph was applied to investigate how the factors (bacterial strain, phage species, and MOI) affect the groupings ( Figure 5). The 10 phages P4 and P10-18 were associated with bacterial growth (Group 1), of which P11 had the greatest association (similar to the phage-free controls), followed by P17, P4, P16, P13, P18, P10, P15, P12, and P14 ( Figure 5A). Accordingly, the greatest bacterial growth was associated with the Demerecviridae phage belonging to the Tequintavirus genus. Myoviridae phages belonging to the Felixounavirus genus and Drexlerviridae phages belonging to the Guelphvirus or Warwickvirus genera were also associated with bacterial growth. Phages associated with more bacterial killing (Groups 2 and 3) included the eight phages P1-3 and P5-9, of which P1 had the strongest association, followed by P5, P9, P2, P3, P8, P6, and P7. A few P4 combinations were also found in Group 3. Accordingly, the greatest bacterial killing was associated with Myoviridae phages belonging to the Tequatrovirus genus. Myoviridae phages belonging to the Mosigvirus genus were also associated with bacterial killing, as well as two out of three Hanrivervirus phages (Drexlerviridae family). The five bacterial strains B1, B4, B6, B7, and B8 were associated with bacterial growth (Group 1). This included all three O1 serotype strains (B1, B4, and B8), as well as two O78 serotype strains ( Figure 5B). The other five strains: B2, B3, B5, B9, and B10, were associated with more bacterial killing and included two O2 serotype strains (B5 and B10) and three O78 serotype strains. The overall effect of the bacterial factor on the plot was less compared to the effect of the phage factor. A high MOI (MOI 10) was associated with more bacterial killing, and a low MOI (MOI 0.1) was associated with less killing ( Figure 5C).  Based on the NMDS analysis, Group 3, including 419 entries, was divided into five subgroups ( Figure 6). Only combinations with phages P1-P9 were found. The growth curves associated with each subgroup are shown in Supplementary Figure S3.
NMDS on a two-dimensional graph was applied to investigate how the factors (bacterial strain, phage species, and MOI) affect the Group 3 groupings (Figure 7). The association of phage P3 was in the direction of subgroup 3.1; however, the association was very weak. Phages P1 and P5 had a strong association with subgroup 3.2, while P2 and P3 had a weak association. A P5 association was in the direction towards subgroup 3.3. The remaining five phages, including P4, P6, P7, P8, and P8, were associated with subgroup 3.4. P6 was the only phage with a strong association ( Figure 7A).
Based on the NMDS analysis, Group 3, including 419 entries, was divided into five subgroups ( Figure 6). Only combinations with phages P1-P9 were found. The growth curves associated with each subgroup are shown in Supplementary Figure S3. represent combinations with the initial bacterial killing followed by exponential bacterial growth and the stationary phase. Subgroups 3.4 (n = 66) and 3.5 (n = 48) represent combinations with the initial bacterial growth followed by the stationary phase. Ellipses indicate a 95% confidence level based on a multivariate t-distribution. Stress < 0.02 = an acceptable fit.
NMDS on a two-dimensional graph was applied to investigate how the factors (bacterial strain, phage species, and MOI) affect the Group 3 groupings (Figure 7). The association of phage P3 was in the direction of subgroup 3.1; however, the association was very weak. Phages P1 and P5 had a strong association with subgroup 3.2, while P2 and P3 had a weak association. A P5 association was in the direction towards subgroup 3.3. The remaining five phages, including P4, P6, P7, P8, and P8, were associated with subgroup 3.4. P6 was the only phage with a strong association ( Figure 7A).
Bacterial strains B5 and B10 (O2 serotype) had a strong association with subgroup 3.1, and B2 (O78) had a weak association. B1 (O1), as well as O78 strains B6, B7, and B9, had an association with subgroup 3.2. B6 was the only stain with a strong association, as well as the only strain driving the plot towards subgroup 3.3. No strains were driving the plot towards subgroup 3.4. Strain B3 was found in this subgroup but showed a very weak/no association. O1 strains B4 and B8 had strong associations with subgroup 3.5 (Figure 7B). The MOI did not have a strong association with any of the subgroups ( Figure 7C). represent combinations with the initial bacterial killing followed by exponential bacterial growth and the stationary phase. Subgroups 3.4 (n = 66) and 3.5 (n = 48) represent combinations with the initial bacterial growth followed by the stationary phase. Ellipses indicate a 95% confidence level based on a multivariate t-distribution. Stress < 0.02 = an acceptable fit.
Bacterial strains B5 and B10 (O2 serotype) had a strong association with subgroup 3.1, and B2 (O78) had a weak association. B1 (O1), as well as O78 strains B6, B7, and B9, had an association with subgroup 3.2. B6 was the only stain with a strong association, as well as the only strain driving the plot towards subgroup 3.3. No strains were driving the plot towards subgroup 3.4. Strain B3 was found in this subgroup but showed a very weak/no association. O1 strains B4 and B8 had strong associations with subgroup 3.5 ( Figure 7B). The MOI did not have a strong association with any of the subgroups ( Figure 7C).

Description of Defined Growth Dynamics Patterns
Based on the bacterial growth (OD) detected in the 2869 phage-host combination experiments, three different growth dynamics patterns groups, including five subgroups, were defined ( Figure 8). The groups included: (1) combinations with a fully resistant

Description of Defined Growth Dynamics Patterns
Based on the bacterial growth (OD) detected in the 2869 phage-host combination experiments, three different growth dynamics patterns groups, including five subgroups, were defined ( Figure 8). The groups included: (1) combinations with a fully resistant bacterial growth pattern; (2) phage-host combinations with a fully susceptible pattern, showing minimal or no bacterial growth; and (3) combinations with one of five in-between patterns characterised by delayed growth, lower killing, or variable killing of the bacterial cells. For all dynamics patterns except Group 1, the phage effect on the bacterial growth kinetics was observed within only a few hours of incubation. When the stationary phase was reached, the bacterial density remained stable throughout all the coculturing experiments.
Group 1 (bacterial growth): The bacterial growth continued to increase during the first 7 h of incubation until the cultures reached the stationary phase. The final OD600 was 0.7. The pattern observed showed logistic growth as under standard conditions without phage present and could be explained by the presence of naturally phage-resistant strains, where the phage is unable to infect and has no effect on the growth. Group 2 (high level of bacterial killing): The bacteria were lysed and never recovered. A single small peak of bacterial growth followed by bacterial killing was observed in some cases but was not reflected in the average OD growth curve. Group 3.1 (initial bacterial killing followed by bacterial growth): Prolongation of the lag phase with no or low bacterial growth for~9 h was observed, followed by a slow increase in bacterial growth until the stationary phase was reached. The final OD600 was~0. 45. The greatest variations were seen in this subgroup (Supplementary Figure S3). Group 3.2 (initial bacterial killing followed by bacterial growth): Prolongation of the lag phase with no or low bacterial growth was observed, followed by a sharp increase in the bacterial density after~7 h of incubation before reaching the stationary phase. The final OD600 was~0.5. Group 3.3 was characterised in a similar way as Group 3.2, except for a higher final OD600 of~7.5, the highest observed for the seven different patterns. Group 3.4 (impaired bacterial growth): Increase of the bacterial growth was observed during the first~9.5 h of incubation until the cultures reached the stationary phase. Impaired growth was observed compared to the Group 1 and Group 3.5 patterns, with a final OD600 of~0.26. Group 3.5 was similarly characterised by impaired bacterial growth. Increase of the bacterial growth was observed during the first~12 h of incubation until the cultures reached the stationary phase. Impaired growth was observed compared to the Group 1 patterns, with a final OD600 of~0.45 (similarly to Group 3.1).
Microorganisms 2021, 9, x FOR PEER REVIEW 13 of 18 bacterial growth pattern; (2) phage-host combinations with a fully susceptible pattern, showing minimal or no bacterial growth; and (3) combinations with one of five in-between patterns characterised by delayed growth, lower killing, or variable killing of the bacterial cells. For all dynamics patterns except Group 1, the phage effect on the bacterial growth kinetics was observed within only a few hours of incubation. When the stationary phase was reached, the bacterial density remained stable throughout all the coculturing experiments. Group 1 (bacterial growth): The bacterial growth continued to increase during the first 7 h of incubation until the cultures reached the stationary phase. The final OD600 was ~0.7. The pattern observed showed logistic growth as under standard conditions without phage present and could be explained by the presence of naturally phage-resistant strains, where the phage is unable to infect and has no effect on the growth. Group 2 (high level of bacterial killing): The bacteria were lysed and never recovered. A single small peak of bacterial growth followed by bacterial killing was observed in some cases but was not reflected in the average OD growth curve. Group 3.1 (initial bacterial killing followed by bacterial growth): Prolongation of the lag phase with no or low bacterial growth for ~9 h was observed, followed by a slow increase in bacterial growth until the stationary phase was reached. The final OD600 was ~0.45. The greatest variations were seen in this subgroup (Supplementary Figure S3). Group 3.2 (initial bacterial killing followed by bacterial growth): Prolongation of the lag phase with no or low bacterial growth was observed, followed by a sharp increase in the bacterial density after ~7 h of incubation before reaching the stationary phase. The final OD600 was ~0.5. Group 3.3 was characterised in a similar way as Group 3.2, except for a higher final OD600 of ~7.5, the highest observed for the seven different patterns. Group 3.4 (impaired bacterial growth): Increase of the bacterial growth was observed during the first ~9.5 h of incubation until the cultures reached the stationary phase. Impaired growth was observed compared to the Group 1 and Group 3.5 patterns, with a final OD600 of ~0.26. Group 3.5 was similarly characterised by impaired bacterial growth. Increase of the bacterial growth was observed during the first ~12 h of incubation until the cultures reached the stationary phase. Impaired growth was observed compared to the Group 1 patterns, with a final OD600 of ~0.45 (similarly to Group 3.1).

Repeatability of Group Assignments
All replicates, including the phage-bacteria-MOI combinations with phage P10-P18, were found in same NMDS group (Group 1). Some discrepancies in the group

Repeatability of Group Assignments
All replicates, including the phage-bacteria-MOI combinations with phage P10-P18, were found in same NMDS group (Group 1). Some discrepancies in the group assignments were seen for replicates of specific combinations, including eight phage P1 combinations, 11 P2 combinations, three P3 combinations, nine P4 combinations, 12 P5 combinations, seven P6 combinations, three P7 combinations, seven P8 combinations, and 10 P9 combinations ( Figure 9). Most discrepancies included combinations with replicates grouped in Group 2 (bacterial killing) and Group 3.1 (in-between and bacterial killing), replicates grouped in Group 2 and Group 3.2/3.3 (initial bacterial killing followed by bacterial growth), or replicates grouped in Group 1 (bacterial growth) and Group 3.4/3.5 (in between dynamics and bacterial impaired growth) (Supplementary Tables S2 and S3).
Microorganisms 2021, 9, x FOR PEER REVIEW 14 of 18 assignments were seen for replicates of specific combinations, including eight phage P1 combinations, 11 P2 combinations, three P3 combinations, nine P4 combinations, 12 P5 combinations, seven P6 combinations, three P7 combinations, seven P8 combinations, and 10 P9 combinations ( Figure 9). Most discrepancies included combinations with replicates grouped in Group 2 (bacterial killing) and Group 3.1 (in-between and bacterial killing), replicates grouped in Group 2 and Group 3.2/3.3 (initial bacterial killing followed by bacterial growth), or replicates grouped in Group 1 (bacterial growth) and Group 3.4/3.5 (in between dynamics and bacterial impaired growth) (Supplementary Tables S2 and S3).

Discussion
It has become clear that successful phage therapy development and application, among others, depend on an understanding of the phage-host interactions and population dynamics [34]. This study presents the classification of bacterial growth dynamics in the presence of lytic phages using two statistical data mining techniques: NMDS and PCA. OD as the input is a fast and data-rich screening method for in vitro phage-host growth dynamics [35]. This approach captures the ongoing dynamics and produces a quantitative high-throughput data to determine the phage-host range, phage virulence/infectivity, Figure 9. Overview of the discrepancies of the groupings of the phage-APEC-MOI combination replicates. P = phage. B = bacterium. MOI = multiplicity of infection. White = no discrepancies. Black = one or more discrepancies. Sixty-nine technical replicates were included for each combination. Combinations with P10-18 were not included, as no discrepancies were detected.

Discussion
It has become clear that successful phage therapy development and application, among others, depend on an understanding of the phage-host interactions and population dynamics [34]. This study presents the classification of bacterial growth dynamics in the presence of lytic phages using two statistical data mining techniques: NMDS and PCA. OD as the input is a fast and data-rich screening method for in vitro phage-host growth dynamics [35]. This approach captures the ongoing dynamics and produces a quantitative high-throughput data to determine the phage-host range, phage virulence/infectivity, and bacterial phage resistance development [16,32,36]. These parameters can be important as pharmacodynamic (PD) parameters and also include part of the pharmacokinetics (PK), as it assesses the potential increase of the treatment dose. This would not be the case when relying on a single endpoint measurement. However, it is understood that ODs do not differentiate between viable and dead cells, and as such, there was no exact link between the OD values and bacteria viability but, nevertheless, a good proxy. The repeatability of the NMDS grouping was shown to be acceptable. Fully natural phage-resistant combinations were clearly identified. Most discrepancies observed between groupings of the replicates from the same phage-APEC-MOI combinations can be explained by grouping cut-off values. Replicates grouped in Group 1 and Group 3.4/3.5 included dynamics characterised by higher or lower bacterial growths. Replicates grouped in Group 2 and Group 3.1 included dynamics characterised by the initial bacterial killing, with no or low subsequent growth.
The optimal number of clusters/groups varied depending on the method used. In this study, we chose five clusters. However, subgroups 3.2 and 3.3 and subgroups 3.4 and 3.5 were relatively similar and could potentially be combined, resulting in only three Group 3 subclusters, and this would also create fewer discrepancies in the group allocations. The discrepancies due to biological variations can be expected due to the spontaneous emergence of phage-resistant variants after varying the incubation time (Group 2 vs. Groups 3.1, 3.2, and 3.3). It is also possible for a very small (partially) resistant sub-population of bacteria to be naturally present in the culture at the start of the experiment [17].
Various factors affecting the phage PK/PD have been described using mathematical and experimental models [15,17,18,37,38]. In this study, the influence of the factors (phage type, bacterial strain, and MOI) on the observed growth patterns was determined. Previous studies have highlighted the MOI influence on phage therapy, and recently, a fast microtiter plate assay for determination of the optimum MOI for a coliphage was further described [39]. However, in this study, we found the MOI to have a less significant effect on the phage-host growth dynamics outcome compared to the phage species. Furthermore, a recent study suggested that the description of MOI alone is not sufficient, as the concentration, particular to the bacteria, can significantly affect the results [40]. Therefore, in this study, the MOI at all the tested values was based on a fixed bacterial concentration.
A quantitative assessment of the phage lytic activity using the virulence score and PhageScore across a large dataset allowed direct comparisons of individual phages. In contrast to a single OD endpoint measurement and the well-established plaque assay, these methods (virulence score and PhageScore) captured the dynamics of phage infection, including bacterial (re)growth after prolonged growth inhibition or lysis. Additionally, compared to the overlay-based efficiency of plating assays and direct spot testing, these methods represent an accurate and less cumbersome and time-consuming approach and do not depend on the subjectivity and/or experience of the observer [41]. However, upscaling of this approach depends on the availability of a high-throughput plate reader. In accordance with previous findings, we found the two methods highly comparable, showing similar properties/values for the studied phages [32,42]. In this study, whenever the bacteria were able to grow and reach the stationary phase, the growth remained stable throughout all the experiments. If a second peak appeared or the growth started to decrease (after reaching the stationary phase), the comparability of the local virulence score and the PhageScore would be reduced. In this study, we only analysed the growth dynamics of cocultures of a single phage type and bacterial strain. However, both the virulence score and PhageScore have previously been used to compare phage combinations for use in phage therapy cocktails [32,41]. In future studies, the inclusion of mixed phage cultures (phage cocktails), preferably targeting different host receptors, may provide further indications of their potential as therapeutics against pathogenic target bacteria [17]. Accordingly, for future applications, the inclusion of phage, as well as bacterial, traits may be required for classification.
Given their great abundance and diversity, multiple candidate phages might be available to infect a target host; yet, we still lack a better understanding of which phage would perform best [43]. One approach to identify the cause(s) of treatment success is to compare the characteristics of phages with high success rates with those of phages with low success rates. The characteristics differing between these two groups of phages become candidates for causation. In this study, Tequintavirus phages (associated with bacterial growth) and Tequatrovirus phages (associated with bacterial killing) would be great candidates for comparisons. Accordingly, the inclusion of phage characteristics, such as the phage receptor, adsorption rate, latency period, burst size, and virion size, may provide further explanation of the phage-host dynamics and may help predict the phage therapy efficacy [44][45][46].
Phage therapy is, by its nature, a strongly selective treatment [17]. Accordingly, when selecting phages for therapeutical application, the emergence of phage-resistant bacteria should be taken into consideration [7]. Bacteria can develop resistance against phages through various mechanisms, including the modification of phage receptor-encoding genes, innate immune systems (such as CRISPR-Cas), and the presence of prophages in the bacterial genome [47,48]. Here, O1 serotype strains B1, B4, and B10 were associated with natural phage resistance/high levels of bacterial growth, whereas the strains with serotype O2 and serogroup O78:H4 were associated with phage susceptibility/bacterial killing. Moreover, all phage-host combinations including P10-18 were found in Group 1 (fully resistant combinations). These phages would be excluded as candidates for phage therapy targeting the selected APEC strains. P1 was the only phage not included in any Group 1 combinations and showed the greatest bacterial killing potential. Accordingly, phages only included in Group 2, characterised by high bacterial killing, are considered the most promising candidates for phage therapy. Multiple phage-host-MOI combinations were characterised by initial bacterial killing followed by bacterial growth (subgroups 3.1, 3.2, and 3.3), suggesting the emergence of phage-resistant bacterial variants. Whether phages found in these subgroups should be considered suitable for phage therapy depends on their specific applications, and further studies are needed to determine if the initial inhibition of bacterial growth for~7 h is sufficient to clear out the infection.
Although in vitro experiments do not capture many in vivo realities, such experiments can give significant insights into the phage-host dynamics and lead to interesting predictions, which could be useful in phage therapy and exploited in appropriately designed in vivo models [13]. Recently, a framework (Clinical Phage Microbiology) with recommendations for in vitro identification and the evaluation of phages intended for treatment was published [42]. One step of the framework pipeline included determination of the growth kinetics of liquid cultures and highlighted the need for a standardised quantitative assessment with reproducible scoring. The methodology applied here constituted such an assessment and may help to improve the standardisation of the quantitative evaluation of the phage candidates.

Conclusions
Our methodology assessing the host-phage interaction in vitro provided a highthroughput method for classifying the bacterial growth dynamics in the presence of virulent coliphages using measurements of bacterial growth by OD as the inputs. The established in vitro model was not only used to gain a better understanding of the phage PK/PD but can also be applied as a screening method for selecting new suitable phage candidates for therapeutic applications against pathogenic target bacteria. However, to fully understand the complexity of these phage-host dynamics, the underlying mechanisms behind these different interactions need to be deciphered.