Multiobjective learning and design of bacteriophage specificity

Machine learning-directed design of bacteriophages

A phage’s RBP is a significant contributor to its species targeting and strain selectivity, but the molecular details of these interactions are largely unknown. In our previous work, we generated large-scale mutagenesis data mapping variants in the terminal 82 residues of the T7 RBP across five E. coli strains^23,24. These hosts were selected due to their known differences in susceptibility to wild-type T7. Data from this work represent both single mutations in the RBP²³ and higher-order mutations selected with guidance from metagenomic databases²⁴. To augment the information density of the training data, we created an additional dataset for multi-host fitness of bi- and tri-mutant RBPs with combinations of single mutants that were beneficial on any host (Supplementary Fig. 1a, b). These combined data provide a valuable resource for understanding how the T7 RBP sequence influences infectivity and selectivity.

We apply supervised sequence-function modeling^30,31 to learn the mapping from T7 RBP sequence to strain-specific virulence. To evaluate which modeling approaches were best suited for these datasets, we examined the extrapolation performance of different approaches by training models on T7 RBP sequences with 1–3 mutations and measuring their predictive performance on data for held-out higher-order mutants with an average of ~6 mutations. We explored three model architectures: a simple linear regression (LR) model, a convolutional neural network (CNN) that learns non-linear features across an input sequence^31,32, and a transformer-based protein language model (ESM2)²⁵ fine-tuned for phage fitness prediction. Models were trained using a standard mean square error (MSE) regression loss function³³. These model architectures have all been applied for protein engineering, but we were curious if the loss function could be further optimized. We constructed a custom loss function, termed a thresholded MSE loss, which accounts for inaccurate fitness measurements caused by the dropout of non-functional variants during selection (Supplementary Fig. 2a). We found that more complex modeling approaches, such as neural networks (CNN and ESM2) and custom loss functions, improved the model’s predictive accuracy over linear regression or standard MSE regression (Supplementary Fig. 2b). Based on these results, we chose to move forward with the CNN and fine-tuned ESM2 architectures in combination with the thresholded MSE loss function.

To guide the optimization of synthetic phages, we compiled a single dataset consisting of all data (train and extrapolation sets from above) and trained ensembles of CNN and ESM2 models with random training and validation datasets. The final CNN and ESM2 ensembles showed high fitness prediction accuracies (weighted Spearman ~0.8) when evaluated on a shared held-out test set of 10% of all data, which included sequences not seen by any models during training and across all three training datasets (Supplementary Fig. 3a–c). We used these models to guide the mutant selection process by employing simulated annealing to search for optimal RBP sequence configurations^31,34,35(Fig. 1a). We also tested if incorporating structural information from ProteinMPNN could boost performance by increasing the likelihood that a given mutant will fold into the wild-type backbone structure²⁶. Altogether, we designed phage mutants using four approaches: CNN, CNN + MPNN, ESM, and ESM + MPNN.

Employing each of these four sequence generation methods, we optimized phages for 26 different objectives: (1) high infectivity for each of the five E. coli strains individually, (2) 20 different specificity combinations, and (3) generalist phages optimized toward all five training strains (Fig. 1b). These phage variants contained 3, 6, 9, and 12 mutations compared to the wild-type T7 RBP which enabled assessment of model performances at diverse mutational distances. In total, we designed, synthesized, and tested a total of 7845 mutant phages. Each design approach led to distinct mutational patterns in the resulting sequences, despite being trained on the same experimental data (Fig. 1c, d, Supplementary Fig. 4a). We highlight an example of the mutational preferences between models by showing the sequence motifs for divergent positions from subsampled clusters (Fig. 1d). Notably, the addition of structural constraints from ProteinMPNN reduced an effect observed in the structure-agnostic CNN and ESM models in which low levels of background mutations were observed throughout the protein (Fig. 1e). Inclusion of structural constraints also increased the preponderance of conservative mutations, measured by mean BLOSUM80 score³⁶ (Supplementary Fig. 4b, p<5e-60).

One concern when using machine learning to generate new protein variants is that models will regenerate sequences preexisting in the training data. Although each model generally favored mutations at similar positions in the T7 RBP, the overall distributions of mutations for the ML-generated phages were highly divergent from the training data (Fig. 1e). Even though we did not directly prohibit generated sequences from overlapping with the training data, only one bacteriophage was an exact match for a previously observed phage. We found a small relationship between the mutational distance and frequency of conservative mutations, but this effect size was minimal beyond comparing three mutants to higher-order mutants (Supplementary Fig. 4c, p<5e-30). This suggests that after including a few choice mutations, additional mutations were, on average, slightly more conservative.

We constructed a pooled library of synthetic phages for experimental testing using ORACLE²³. ORACLE is a phage genome engineering tool we previously developed that uses Cre recombinase-guided insertion of a variant library at a targeted genomic locus (Fig. 2a). The resulting phage population represents the pre-selection library of RBP variants, which was then incubated independently with each strain of interest in triplicate. We sequenced pre- and post-selection libraries and calculated the relative change in abundance of each library member to determine phage fitness on each host. These calculations yield continuous measurements of relative virulence for each phage in the library against each host of interest.

We evaluated the predictive capacity of models using the average precision (AP/AUPRC) classification³⁷ metric based on a model’s ability to classify whether phages successfully met their given objective. We define success as the ability to infect the target host for infectivity, higher virulence on target hosts than anti-target (minimized) hosts for specificity, and ability of phage to infect all screened hosts for generalist phage. Each model demonstrated high accuracy at low mutational distances, and ESM + MPNN had the highest accuracy at all mutational distances (Fig. 2b). To assess where in the RBP models made mistakes, we calculated the Kullback-Leibler Divergence (KLD) between the distributions of mutation frequencies for non-functional phage compared to all engineered phage at each mutational distance³⁸. We found that each model made mistakes at different positions, due to increased KLD at different regions in the primary sequence (Supplementary Fig. 5a). Additionally, mistakes were generally made in similar positions regardless of the number of mutations (Supplementary Fig. 5b). In both cases, mutations that were predicted to be functional but were actually deleterious were primarily made at less frequently mutated positions (Supplementary Fig. 5a, b)

We mapped the measured fitness values for each host onto a dimension reduction of the mutant RBP sequences to visualize the fitness relationships between engineered bacteriophages (Fig. 2c). Through this analysis, distinct patterns are apparent in the fitness and specificity landscapes, such as clusters of sequences with multi-host specificity and high similarity in the sequence-fitness effects for susceptible strains (10G, BL21, BW25113) or resistant strains (BW25113ΔrfaD, BW25113ΔrfaG). There also appear to be several regions of sequence space where few variants are functional on any host, representing objectives and/or mutational distances that proved challenging for the models to predict.

Machine learning enables the optimization of diverse phages with high virulence

When utilized in therapeutic or biosanitation applications, bacteriophages are often administered as cocktails of diverse phages, which helps restrict target pathogen evolution^17,18,39,40. Though this work is not application focused, we considered what properties may be useful to optimize for different applications. One basic objective in cocktail preparation is to select phages with high virulence against the target host, which can help ensure rapid clearance of the pathogen. Thus, our first design objective, which we denote as “infectivity”, was to engineer phages with high virulence on individual target E. coli strains (10G, BL21, BW25113, ΔrfaD, and ΔrfaG).

We define success in optimizing phage infectivity as the ability of a phage to infect its intended host above the limit of detection. The majority (80–95%) of triple mutants successfully infected their target strain, while variants at high mutational distances showed reduced success rates (Fig. 3a). Nonetheless, over 30% of the 12-mutant phages generated by the CNN and ESM+MPNN models retained the ability to infect the target strain. We also consider overall success when pooling phage by different design criteria, such as models that use a CNN vs ESM or whether a model includes structural constraints. Sequences designed with the CNN were more likely to be functional than those designed with ESM at an edit distance of 9 (p<0.01), but not other distances (Supplementary Fig. 6a). Focusing instead on the impact of structural constraints, there was no significant difference in infectivity design outcomes at any mutational distance between models that included ProteinMPNN and those that did not (Supplementary Fig. 6b). No improvement in success when including ESM or MPNN may be a result of overemphasizing protein stability over dynamic function⁴¹ or could indicate that the sequence-function models had sufficient data to model single-host virulence and adding additional constraints hindered success.

For hosts with high susceptibility to wild-type T7 (10G, BL21, BW25113), engineered phages show minimal improvement compared to wild-type T7 (Fig. 3b). This is likely because T7 is already highly optimized for these hosts, leaving limited room for further enhancement. In contrast, wild-type phage infectivity on resistant hosts (ΔrfaD and ΔrfaG) is more than four orders of magnitude lower relative to susceptible hosts²³. Regardless of this natural deficiency, many phages optimized for these resistant hosts are several orders of magnitude more virulent than wild-type T7 on these resistant hosts. In general, there are few common mutational motifs within each objective (Fig. 3c). We found that only a few mutations are necessary to improve phage infectivity and increasing the mutational distance was unlikely to enhance virulence further (Fig. 3b, Supplementary Fig. 6c). Considering that we observe high success without the emergence of common motifs, it appears that there are many potential optima to modify the T7 RBP for enhanced virulence on single hosts, demonstrating high adaptability of T7 phage.

Directed design of high-order bacteriophage specificity

While it is possible to apply directed evolution to enhance single-host virulence, engineering phages with multi-host targeting specificities is challenging and not easily achieved through traditional directed evolution methods^14,42. Multi-task machine learning models offer a promising approach for engineering proteins with precisely defined combinations of properties. Applying the same models as above, we evaluated their ability for multi-task design by optimizing phages for 20 different specificity objectives, each targeting a unique combination of 2–5 target and anti-target strains (Fig. 1b). Success for this specificity task is defined by a phage showing higher relative virulence on all target strains than all anti-target strains (scores calculated from selection and sequencing assay). When pooling sequences by base model (CNN/ESM) or structural constraints, we found that the ESM-based models outperformed their CNN counterparts at edit distances of 6 and 9 (Supplementary Fig. 7a, p<0.01). Inclusion of structural constraints with ProteinMPNN improved success rates for specificity by 9, 3, and 11% at mutational distances of 6, 9, and 12 (Supplementary Fig. 7b, p<0.002). Overall, the best approach for designing specificity was ESM+MPNN, which had higher success rates compared to all models at 9 mutations and higher success rates than CNN and ESM base models at 12 mutations (Fig. 4a, p<0.02).

We denote each specificity objective as a sequence of five characters representing the five strains: 10G, BL21, BW25113, ΔrfaD, and ΔrfaG. In this notation, a “1” indicates that virulence toward that strain was maximized, a “0” indicates that virulence was minimized, and strain positions marked as “x” were not considered (Fig. 1b). For example, the task “x101x” specifies an objective that maximizes fitness for BL21 and ΔrfaD, minimizes fitness for BW25113, and does not consider 10G or ΔrfaG. For a more granular assessment of the specificity of engineered phages for each objective, we plotted the average fitness of the target strains against the average fitness of the anti-target strains. In this analysis, successful specificity designs are expected to occupy the upper left sector, and the complexity of the objectives increases from left to right across the subplots. Many of these phages demonstrated the intended specificity profile, with some showing high virulence on all target strains and infection below the limit of detection for anti-target strains (Fig. 4b, Supplementary Fig. 8a–c).

Based on results from our high-throughput phage selection assay, we successfully designed phages for all 20/20 specificity objectives. For 19 of these objectives, our designs achieved diverse specificities with only 3 mutations from wild-type T7, indicating high adaptability in the natural T7 RBP (Fig. 4b). The 00×1x objective only had successful designs at higher mutational distances (Supplementary Fig. 8a). In several instances (x0011, 00010, 0xx1x, 00×1x, xx011, 00011, x1011, x101, 1101x, 11011) mutants achieved the complete opposite specificity as wild-type T7 with only a few mutations. Overall, it was very uncommon for designs to have a specificity opposite to the intended task, and the more common mode of failure for design was for a phage to be completely non-functional.

We observed that some specificity objectives were much more challenging to design (00010, 0xx1x, 00×1x, 00011), and most of the phages for these objectives were simply nonfunctional. One contributing factor is likely that simultaneous minimization of 10G and maximization of ΔrfaD is a challenging task due to similarities in these strains. Low success for these objectives is also in part because, during optimization for these tasks, many sequences were predicted to have a large gap in scores between maximized and minimized objectives, but overall, low scores for both the target and anti-target strains. This over-minimization toward anti-target hosts is likely also a contributing factor to the progressive decline in success rates for higher complexity objectives (Supplementary Fig. 9). Over-minimization of negative function can be remedied through the inclusion of evolutionary or structural constraints (Supplementary Fig. 7a, b). Retrospectively, this effect could have also been mitigated by adjusting the scoring function to co-optimize for high fitness on target strains in addition to a large gap between maximized and minimized strains. Because the most common failure mode when optimizing for specificity was a complete loss of virulence, altering the optimization function or including structural constraints could enhance success rates, an important consideration for future specificity engineering projects.

To summarize design outcomes, we display the target and measured fitness for the best-designed phages across all models and edit distances (Fig. 4c). Scores for strains that were not considered during mutant selection are concealed for clarity. Each modeling approach produced the best mutant phage for at least one specificity objective, and the best mutants were found across every mutational distance. While we do observe phages that are selective for each task, objectives minimizing 10G and maximizing rfaD (0xx1x, 00×1x, 00010, and 00011) had the weakest specificities.

To investigate the mutational drivers underlying diverse phage specificities, we analyzed a group of objectives that share a common design goal: minimizing BW25113 while maximizing ΔrfaD. This set includes x0011, xx01x, x101x, and 11011. Among these, xx01x, x101x, and 11011 consistently favor charged or polar substitutions at positions 479, 540, and 544—three of the most frequently mutated sites in phages optimized for these objectives (Fig. 4d, e). In contrast, x0011 tends to preserve wild-type residues at these positions, instead favoring mutations at positions 480, 518, and 520. Since x0011 differs from the others principally in its additional requirement to minimize BL21 virulence, these patterns suggest the existence of at least two distinct optima within the core motif: one that tolerates or enhances activity against BL21, and one that actively suppresses it. More broadly, the emergence of divergent mutational preferences among objectives sharing a common design subtask illustrates the ability of multi-task sequence-fitness models to navigate these trade-offs when engineering multi-part specificities, a clear benefit of multitask modeling.

Structurally, some objectives show a higher frequency of mutations at the exterior loops of the RBP, which is expected to be the primary binding interface to the LPS receptor^23,43,44. However, many mutations were observed across the entire RBP domain for all objectives (Supplementary Fig. 10a–e). Of note, for objective 11100, the top three most frequently mutated positions were located on the interior loops, which are distal to the exterior loops that are expected to primarily interact with the cell surface (Fig. 4e, Supplementary Fig. 10e). The presence of key mutations in non-intuitive places in the protein structure highlights the importance of mutational scanning beyond key domains and the benefit of applying data driven modeling.

Machine learning-directed design of generalist phages capable of infecting new hosts

A major constraint of natural phages is their individually narrow host range, which limits their utility against non-predisposed populations. This limitation is particularly problematic in clinical settings, where timely and precise strain identification is often challenging, resulting in the frequent use of broad-spectrum antibiotics. Rapidly engineering generalist bacteriophages with broad host range is therefore of high interest. Ideally, a researcher could start with a model phage system and identify mutations that render novel or enhanced virulence toward new targets. In line with this goal, the final design objective was generality.

These phages were optimized to exhibit high fitness on each of the five strains for which we had data. Characterization of engineered generalists revealed that these phages were frequently capable of infecting all training strains with approximately 80% success rates at an edit distance of 3 (Supplementary Fig. 11a,c). At edit distances of 6 and 9 mutations, the CNN and ESM + MPNN design approaches produced generalist phages more frequently than the other two approaches (p≤0.025), each achieving >50% success for infecting the five target strains at 9 mutations (Supplementary Fig. 11c).

While high success rates were observed in designing generalists towards the training strains, the ideal advantage of designing generalist function lies in the potential for phages optimized toward known strains to infect a broader range of strains with unknown susceptibility profiles, accelerating screening for novel fitness. To evaluate the true promiscuity of these phages, we challenged the engineered phages against eight additional E. coli strains with alternative modifications to the surface receptors, including deletions to genes involved in LPS modification (BW25113ΔrfaC, BW25113ΔrfaF, BW25113ΔrfaJ, BW25113ΔrfaI, BW25113ΔlpcA), unrelated iron transporter genes (BW25113ΔtonA and BW25113ΔtonB), and full deletion of the LPS primary receptor (ClearColi).

The receptor deletion strains are expected to have variable degrees of truncations to the LPS, with ΔrfaC, ΔlpcA, and ΔrfaD resulting in expected truncations directly above the KDO sugar, followed by ΔrfaF, ΔrfaG, ΔrfaI, and ΔrfaJ, each truncating one additional level of sugars above this KDO base relative to each other^24,45. ΔtonA and ΔtonB should show similar variant-host susceptibility as BW25113, as these deletions have the same primary LPS receptor as BW25113⁴⁶. Many generalist phages exhibited high virulence on these non-training strains even without being directly designed for these tasks (Supplementary Fig. 11b, c). Approximately 35% of phages designed by the CNN or ESM + MPNN models were functional on all 13 E. coli receptor deletion strains at an edit distance of 9 (Fig. 5a). Phages that did not infect all 13 strains were typically only unable to infect 1–2 non-training strains at low mutational distances (Fig. 5b). At higher mutational distances, the most common failure mode was an inability to infect any of the 13 strains, suggesting a complete loss of phage viability (Supplementary Fig. 11d). The frequency of generalist behavior was ~2x higher for phages directly optimized as generalists compared to phages optimized for infectivity, and ~30x higher compared to specialist phages (Fig. 5c). This improvement in promiscuity emphasizes the advantage of machine learning-directed generalist design to optimize for new tasks.

The distribution of fitness scores for generalist phages was shifted below wild-type for susceptible hosts and above wild-type for resistant hosts (Fig. 5d). This indicates that, as with specificity, there is a fitness tradeoff for multi-host targeting compared to optimization for single strains. The most common mutation for all generalists was G480K, which occurred in ~50% of variants and was also identified to be highly beneficial for 3/5 of the host strains in the single mutant dataset⁸ (Supplementary Fig. 10f). Other than this one common mutation, the successful generalist phages were highly diverse, and mutations were not exclusively localized to the primary binding interface of the RBP but rather spread throughout the protein. Similar to what was observed with the infectivity tasks, this result suggests high plasticity of T7 phage even when optimizing for many hosts.

Secondary validation of designed phage functionality

Due to the nature of selection-based sequencing assays, the fitness scores for a given phage are a function of all library members and the length of selection. This effect leads to high scale variability in scores measured with high-throughput selection assays. As a result, the scores do not represent an absolute measurement of virulence, but rather a relative measure. To assess the absolute specificity of the engineered phage, we selected six optimized phages for clonal preparation and screening with a plate-based plaque assay, in which the number of plaques was counted in duplicate on relevant hosts and compared to a reference host to compute the efficiency of plating (EOP).

For several variants, we observe significant improvements in virulence relative to wild-type T7 RBP (Supplementary Fig. 12). For example, a clonal phage for objective xxxx1 has a 1.7×10⁴ fold improvement in EOP compared to wild-type T7 (Supplementary Fig. 12a, p<0.0005). For a nine mutant generalist phage (11111), we observed a 5.0×10³ fold increase in the minimum EOP on the five E. coli strains that the variant was designed toward (Supplementary Fig. 12f, p<0.05). This improvement was achieved by maintaining virulence on the susceptible 10G, BL21, and BW25113 strains at levels similar to wild type, while introducing mutations that improved virulence toward the resistant ΔrfaD and ΔrfaG strains 5.0×10³ and 110-fold, respectively.

In optimizing for specificity, we selected mutations predicted to maximize the virulence gap between target and anti-target strains. Notably, this strategy permits decreases in virulence on target strains as long as these effects produce greater depletions in virulence toward the anti-target hosts. This effect is apparent for a six mutant from the objective x101x, which had a 1.7×10⁴ fold improvement in specificity over wild-type (Supplementary Fig. 12c, p&lt;0.0005) but also incurred a ~five-fold decrease in virulence toward one of the target strains. Similarly, a variant from objective 11011 showed a decrease of 12x toward the 10G target strain and a decrease of 15x toward the BL21 target strain, yet these depletions allowed a net 2.7×10⁶ fold increase in specificity relative to T7 (Supplementary Fig. 12e, p<0.05).

Our validation also highlights the constraints imposed by relative fitness measurements in interpreting outcomes from selection-based sequencing assays. For challenging design objectives, such as xx01x and x0011, the selected variants did not display absolute specificity in plaque assays, with EOPs on anti-target strains exceeding those on target strains. Nonetheless, our models produced variants that exhibited substantial relative gains in specificity toward target over anti-target strains relative to wild-type: 9.5×10⁴-fold (Supplementary Fig. 12b, p<0.0005) and 590-fold (Supplementary Fig. 12d, p<0.005) for xx01x and x0011, respectively. Despite assay constraints, our models consistently produced variants that improved selectivity by several orders of magnitude, validating the effectiveness of leveraging insights from multifunctional sequence data for multiobjective design.

Lead Contact

Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Srivatsan Raman (sraman4@wisc.edu).

Materials Availability

This study did not generate new materials.

Data and Code Availability

• Processed sequencing data reported in this paper are deposited at GitHub and are publicly available as of the date of publication: https://github.com/naia-novy/multiobjective-protein-modeling/. DOI is listed in the Key Resources Table.

• All original code reported in this paper is publicly available as of the date of publication: https://github.com/naia-novy/multiobjective-protein-modeling/. DOI is listed in the Key Resources Table.

Microbes and culture conditions

All culturing experiments were conducted at 37°C with shaking in standard LB media. The E. coli bacterial strains used in this study are as follows: 10G, BL21, BW25113, BW25113ΔrfaD, BW25113ΔrfaG, BW25113ΔrfaC, BW25113ΔrfaF, BW25113ΔrfaI, BW25113ΔrfaJ, BW25113ΔlpcA, BW25113ΔtonA, BW25113ΔtonB, ClearColi BL21(DE3). Susceptible strains included E. coli BL21, a B strain derivative, BW25113, a K-12 derivative, and 10G, a DH10B derivative. All receptor deletion strains are resistant to phage infection and contain truncations at different points in the lipid polysaccharide receptor ^58–60. These deletions were made in the BW25113 strain and include BW25113ΔrfaD, BW25113ΔrfaG, BW25113ΔrfaC, BW25113ΔrfaF, BW25113ΔrfaI, BW25113ΔrfaJ, BW25113ΔlpcA. Deletion of these genes can reduce the ability of T7 to infect E. coli by several orders of magnitude²³. ClearColi BL21(DE3) contains a full deletion of the lipid polysaccharide in the BL21 background strain. BW25113ΔtonA and BW25113ΔtonB represent deletions in genes encoding an iron transporter and are not expected to affect phage infection. Deletion strains are referred to throughout simply as the gene that was deleted for brevity. No antibiotics were used during the selection experiments, but during library preparation, Kanamycin50 or Kanamycin50+Spectinomycin115 were used.

Machine learning model training and implementation

CNNs were implemented in PyTorch and trained from random initializations. The 35 million parameter model from the ESM2 suite was fine-tuned by appending a two-layer fully connected network after the embedding block and allowing training of the final embedding layer and these appended neurons. Model inputs for the LR and CNN models are one-hot encoded sequences, and all models predict per-strain fitness scores. An ensemble of 100 CNNs and an ensemble of 10 fine-tuned ESM2 models were trained for downstream phage engineering. Each model in the ensembles was trained with different random initializations. We used 100 or 10 models in the CNN and ESM2 ensembles, respectively. 80–10-10 training, validation, and test splits were used for each model, with the test dataset remaining constant for all dataset splits. A weighted MSE loss function was implemented based on the fitness scores and the inverse of the replicate standard error. This approach biases the model the fit closer to higher-quality data and reduces overfitting. For the thresholded regression, a threshold was set manually for each strain and dataset based on the distribution of fitness scores at the score below which phages were expected to be fully nonfunctional. Weights for scores below this threshold were set to the mean weight for each output to ensure the models did not overfit to functional variant data. Strains and each dataset were organized as separate outputs due to scaling differences across datasets. To determine the final predicted scores for each strain, the relevant scores for each dataset were averaged in a weighted fashion, where weights for each dataset were optimized during hyperparameter optimization. All hyperparameters were optimized using Optuna⁶¹. Model training was done using resources from the Center for High-Throughput Computing (CHTC) at UW-Madison.

Processing and generation of training data

This work utilizes three different datasets describing the impact of mutations of phage virulence on a panel of E. coli hosts. Two of these datasets were previously published. The first is a site saturation mutagenesis library, which contains all single mutations in the terminal 82 residues of the T7 RBP. This phage library was screened against 10G, BL21, BW25113, ΔrfaD, and ΔrfaG. A second dataset was collected as part of this work to augment the amount of training data available (Supplementary Fig. 1). We refer to this library as a combinatorial library, which was screened against 10G, BL21, BW25113, and ΔrfaG. The combinatorial library contains 12,830 combinations of double and triple substitutions of 122 individual substitutions throughout the terminal 82 residues of the T7 RBP. Substitutions were selected based on activity profiles in the original five hosts tested during site saturation mutagenesis of the tip domain. The top 25 candidates enriched on all five strains and a series of 20 substitutions that demonstrated enrichment on between one to four strains and depletion on the remainder of strains were used to generate the list of 122 candidate substitutions. The third library includes mutations spanning the same region of the RBP, but containing up to 10 mutations identified from metagenomic motifs. This library was screened against a large panel of E. coli strains, which includes data for 10G, BL21, BW25113, ΔrfaD, and ΔrfaG. To help normalize distribution variance across strains, data for each strain were z-score normalized using SciPy⁶².

Variant library generation

Variants were generated over 25,000 iterations of simulated annealing based on model predictions and iterative sampling of random mutations. A desired edit distance was specified before the search, and a decision to replace an existing mutation was made at each step. Beneficial mutations were always incorporated, and deleterious mutations were accepted with a probability that decreased over the annealing cycle. Median predictions from the CNN and ESM2 models were used for all optimization tasks. The mean ProteinMPNN predicted probability of each residue fitting the wild-type backbone was used to guide optimization for MPNN-constrained models (Protein_MPNN_v_48_010 weights). When ProteinMPNN was used to guide the optimization of sequences, the product of the CNN or ESM2 infectivity/specificity score and ProteinMPNN scores was used as the objective function. For infectivity objectives, the predicted score for the strain of interest was maximized directly. For specificity objectives, the gap difference between the lowest score to be maximized and the highest score to be minimized was used to guide mutation sampling via a utopia logistic gap function. Generalists were optimized based on the minimum predicted score of all 5 strains on which the model was trained. For our specificity objectives, we attempted to select a representative subset of the 211 meaningful combinations of functional profiles for up to 5 strains.

Ten times the number of desired phage variants (~75,000) were generated. Following the 25,000 iterations of generation, the resulting design pool was filtered to remove duplicate sequences and restricted to endpoint sequences that completed the full annealing cycle. To select the final variants for experimental validation, sequences were ranked using a lower confidence bound, defined as the ensemble objective score minus its standard deviation, to prioritize high-performing designs with low model uncertainty. To ensure the validated library captured diverse sequence space rather than a single local optimum, these top-ranked variants were one-hot encoded and grouped into 20 distinct clusters using K-means clustering. Fifteen phages were selected for each model and distance for infectivity and specificity. One hundred phages were selected for each generalist model and distance.

*After sampling several different scoring functions for specificity, we selected an optimization function based on a utopia logistic gap function. This function balances the trade-off between maximization and minimization while theoretically preventing over-optimization once a reasonable score is achieved in either direction, using a logistic (sigmoid) transformation. Before each operation, the sigmoid transformation is applied to normalize each score, ensuring that the true maximum and minimum scores are 1 and 0, respectively. Then, the following operation is used to compute the specificity score.

Phage library preparation

Genes encoding the mutant 82 amino acid terminal region of the RBP were ordered as an oligo pool from Twist and amplified using PCR. The resulting DNA fragments were cloned into a bacterial plasmid (pHT7Rec1 Template) containing a type IIS landing pad for the gene insert. The resulting library was used in the ORACLE phage library preparation method⁸. The first step of ORACLE is optimized recombination, in which the library is recombined into an acceptor phage with lox71 landing sites. Next, in the accumulation step, the resulting phage library is enriched for phages that have undergone this recombination step by allowing infection of the phage into a strain that contains Cas9 and a guide targeting a site in the genome of non-recombined phage (SC101_pCas9_Target3). During these two steps, the bacterial host expresses the wild-type T7 RBP at high levels to reduce library bias (pHT7Helper1 Pbr mCherry+gp17). The inclusion of this helper plasmid yields a genetically modified phage that presents the wild-type RBP. Finally, the enriched library is expressed by having the accumulated phage library infect a bacterial strain lacking this wild-type RBP helper plasmid for a single replication cycle.

Growth-based selection of phage

Bacterial cultures were grown to an OD600 of ~0.4 before adding the phage library at a multiplicity of infection (MOI) of 0.01 and shaking at 37 °C for two hours (or until culture lysis for ΔrfaD, ΔrfaG, ΔlpcA). After this selection, the cultures were placed on ice to halt further infections, and a small amount of chloroform was added to the culture to kill any remaining bacteria. Post-selection phage populations were harvested by filtering the lysate through a 0.22μm filter.

Next-generation sequencing and data processing

Post-selection phage populations were used as a template for a PCR (KAPA HiFi, KK2102) directly off the phage genome using primers that targeted the mutated region in the RBP. These primers also added indexes corresponding to the biological replicate. A second PCR was conducted using the first PCR as a template, which added i5 and i7 indexes for Illumina sequencing. Sequencing was done at the UW-Madison Biotechnology Center on a NovaSeq X Plus shared 2×150 run. The resulting data was processed using custom Python scripts to extract the sequences of input and count the number of observations for each variant. Once pre- and post-selection counts for each variant were available, variants were filtered to exclude sequences with less than 25 pre-selection counts or less than 5 post-selection counts for all of the strains screened. Then, scores were calculated using the log ratio of post and pre-selection counts and maximum likelihood estimation to combine biological triplicates⁶³. All analyses utilize or show the combined triplicate scores. Count correlations for nearly all replicates frequently exceeded Pearson correlations of >0.95, except in instances of post-selection measurements where variant dropout was frequent. For the data shown in this study, scores were z-score normalized⁶² and additionally normalized between zero and one using the manually set threshold for phage function and the 98^th percentile score for that strain. This normalization approach was chosen so that zero would approximately represent non-functional phage and 1 would represent the maximum observed fitness on that strain (excluding outliers).

Computational analysis and statistical tests

Plotting and analysis were done using seaborn and matplotlib^64,65. Weighted Spearman metrics were calculated using a Python implementation of the weighted Spearman R package, where the weights were the same weights used during training⁶⁶. Dimension reduction was done on one-hot encoded RBP variants using PaCMAP⁶⁷ with the default parameters, except for n_neighbors and MN_ratio, which were set as 100 and 2, respectively. The coordinates of the dimension reduction were saved and reused in all dimension reduction figures. In Fig. 1c and 1d, circled regions were placed during interactive plotting, and the sequences were determined based on the coordinates of points within the circle. These circled regions were selected for their high concentration of sequences designed by each model and high abundance of generalist designs (to reduce motif divergences caused by differences in objective). A Mann–Whitney U test statistical test was used to compare distributions of BLOSUM scores. A Fisher's Exact Test was used to compare the success summed across sub-objectives for all pairs of models and mutational distances for infectivity, specificity, and generalist tasks. A student’s t-test was used to determine the significance between conditions in clonally validated phage data. During some model comparisons, we use the average precision (AP) metric to compare classification accuracy due to unbalanced datasets. This metric is calculated by taking the area under the precision-recall curve, where precision is defined as (true positives) / (true positives + false positives), and recall is defined as (true positives)/(true positives + false negatives)³⁷. Kullback-Leibler Divergence (KLD) was calculated using SciPy⁶².

Plaque assays

Plaque assays were performed to estimate phage titers by mixing phage within a top agar solution of 0.5% agar with 250μL of an overnight 10G+Helper bacterial culture. Phage was added at a dilution to enable plaque counting. These plates were incubated for 16 hours before counting plaques. The titers determined through plaque assays were used to guide the amount of phage required to achieve the target MOI during selection experiments and for clonal phage validation of experimental results. Efficiency of plaquing (EOP) was calculated for each phage individually by taking the ratio of plaquing efficiency on each relevant host and the plaquing efficiency on 10G+helper, which constitutively expresses the wildtype T7 RBP.

Construction of clonal phage variants:

To make clonal phage variants for validation, Gibson Assembly was performed using the NEBuilder HiFi DNA Assembly Master Mix (NEB E2621L). The T7 phage genome was split into five ~10kb fragments and mixed at a 1:1 molar ratio along with a synthetic gBlock carrying the gp17 variant sequence. Following the Gibson Assembly⁶⁸, drop dialysis was performed using a 0.025μM mixed cellulose ester membrane (VSWP01300), and E. coli 10G cells were then transformed using electroporation. The transformation mix was then recovered for 1 hour at 37 °C before 200μl of the transformation mix was directly plated in 0.5% LB+10G top agar⁶⁹. Following overnight incubation at 37 °C, individual plaques were picked and sequenced. Codon usage for mutations was determined by selecting codons optimized based on the T7 genome. Mutations used for clonal validation are noted in the table below. Only the terminal 82 residues of the RBP were mutated. Mutations are numbered according to the position in the full T7 RBP.