EP4630966A2 - Intelligent design and engineering of proteins - Google Patents

Abstract

A pipeline for designing and engineering a protein is composed of directed evolution, sequencing, and machine learning analysis. This pipeline can explore the directed evolution sequences that are not present in libraries and can be used in facilitating the discovery of therapeutic proteins and imagine probes, and/or enhancing efficiency of CRISPR and enzyme activities. scFvs mutants with tight binding to and/or slow off-rate towards PD-L1 are discovered by this pipeline. These scFv could be used as potent antibodies against PD-L1 in cancer immunotherapy and/or in CAR-T cell therapy

Description

INTELLIGENT DESIGN AND ENGINEERING OF PROTEINS

CROSS-REFERENCE TO RELATED APPLICATION

This application claims the benefit of U.S. Provisional Application No. 63/431 ,556, filed on December 9, 2022, and U.S. Provisional Application No.: 63/432,836, filed on December 15, 2023, respectively, the entire content of each of the above referenced applications are incorporated herein by references.

STATEMENT REGARDING FEDERALLY SPONSORED RESEARCH OR DEVELOPMENT

This invention was made with government support under R01 GM1 11941 and R21AI1581 14, each of which is awarded by the National Institute of Health (NIH). The government has certain rights in the invention.

BACKGROUND

Identification of proteins and small molecules with desired properties is a task of crucial importance for the pharmaceutical and chemical industries. In recent years, machine learning assisted approaches have become more popular. Neural network/deep learning models are widely used owing to their flexibility, their ability to learn complex relationships from large datasets, and - owing to the development of libraries like PyTorch - ease of implementation.

Despite deep learning’s success, it suffers from several limitations. First, most deep learning architectures do not quantify their uncertainty on predictions, tend to be “over-confident” when extrapolating and are “black-box”, so that it is very hard to determine why the model makes a particular prediction for a specific input. Second, surprisingly, deep learning models have repeatedly been shown to lack robustness to so- called “adversarial attacks”, small perturbations to their input. Adding subtle noise not perceptible to the human eye, for example, could reliably cause deep learning models to misclassify a photo of a cat as guacamole, or a lionfish as eggnog. In biology, protein structure prediction models have been shown to exhibit the same vulnerability to adversarial attacks. Finally, deep learning often entails enormous computational cost. Many state of the art models use tens or hundreds of billions of learned parameters. Often a large fraction of their parameters can be removed without damaging performance, suggesting a more efficient approach is possible. For example, it was reported that 97% of the parameters of the BERT model were pruned in natural language processing and equivalent performance was yet achieved.

An alternative model architecture is Gaussian process (GP) regression. A GP is a Bayesian model which defines a multivariate normal distribution over possible functions mapping the input x variable to an output y. In regions where there is little or no training data, the model expresses high uncertainty. In regions where there is substantial training data, the model predicts outcomes with more confidence. The type of functions in the distribution is determined by a kernel function that needs to be selected and its hyperparameters “tuned” (as with any other model) for a specific problem. The kernel function measures the similarity of any two x-inputs.

GP models have at least four compelling advantages. First, like deep learning, if equipped with an appropriate kernel, a GP is able to approximate any relationship. Second, GP can calculate the marginal likelihood - the probability of the data averaged over all possible function values - in closed form. Consequently, kernel hyperparameters can be “learned” by maximizing marginal likelihood rather than the likelihood, reducing the risk of overfitting and thus making the model more robust. Third, the model’s predictions are generated using the similarity of new data points to those in the training set, where similarity is quantified by the kernel function. Unlike a deep learning model, a GP is not a black box, because one can determine how the model measures similarity between datapoints and generates predictions. The same kernel function that the GP model uses to make predictions can also be used to cluster the training data and retrieve the most “similar” datapoints in the training set, where “similarity” is determined by the kernel function we selected.

Fourth, a GP reliably quantifies its uncertainty and assigns high uncertainty to datapoints very dissimilar from its training set. Uncertainty is important for protein engineering and drug design where experimentally evaluating model predictions is expensive. In these scenarios, one would rather use only model predictions that are more likely to be reliable. GP would detect “distribution shifts”, where new data points are very different from the training set by assigning high uncertainty to its predictions, while a deep learning model may “fail silently”, generating low-accuracy predictions without providing any obvious sign of failure. Uncertainty can also be used for Bayesian optimization / active learning, where the model assists the practitioner in selecting which datapoints to experimentally evaluate next. Several techniques for estimating uncertainty in deep learning have been introduced in literature. However, the uncertainty estimates they provide are much more poorly calibrated than those provided by a GP.

GPs suffer from several limitations. The ability to choose the kernel function can be both a drawback and a benefit, since there are many problems to design an appropriate kernel (e.g., image classification) without enough prior knowledge. It may be possible in some cases to overcome this drawback by combining deep learning with Gaussian processes, and indeed this possibility has been explored by some practitioners in literature. For example, a deep learning model can be trained in an unsupervised fashion to learn a feature representation which can be used as input to a GP. In this strategy, the GP essentially serves as the last layer of the deep learning model, thereby hopefully garnering some of the benefits of both approaches.

The main challenge for GP is computational expense. A straightforward implementation scales with dataset size as O(W³). Indeed, the poor scaling of GPs is often cited as a reason to prefer deep learning. For example, when building a model for traffic prediction, it was argued that it was not possible to use a GP for their problem because “GPs are hard to scale to a large dataset”. A second major challenge is the lack of efficient kernels for sequence and graph data. Many kernels for graphs have been described, but typically exhibit quadratic or worse scaling in the size of the graph.

One approach to approximating a GP is the random Fourier feature (RFF) approach of Rahimi and Recht. Briefly, RFF methods approximate the kernel function via a random map corresponding to the kernel of interest, so that each input datapoint is converted into a “random feature” representation; one is then able to approximate the GP using Bayesian linear regression in the random feature space, and thereby achieve linear scaling with training set size.

Monoclonal antibodies are among the most successful biological therapeutics. Despite their impressive versatility and specificity, development of therapeutic antibodies continues to pose a variety of complex challenges. Typically, initial hits have insufficient affinity, and their binding must first be improved through in vitro affinity maturation, whereby repeated rounds of mutagenesis and selection for antibodies with improved affinity are performed. This process is frequently time-intensive and may take months to complete and cannot simultaneously optimize for other desirable properties like good solubility and low immunogenicity. Computational techniques that could aid in the faster identification of high-affinity antibodies with desirable properties would likely accelerate this process.

Traditional computational approaches to antibody binding affinity rely on estimation of free energy. These are often limited by high computational cost, low throughput and the limited reliability of the free energy estimates generated by these methods. Alternatively, machine learning techniques have been applied to both protein engineering and a variety of tasks in antibody design.

Machine learning-based approaches face at least two major challenges. The first one is the lack of estimated uncertainty in the predictions of the binding affinities or other properties. As the training data can only ever cover a small fraction of the sequence space, machine learning models typically perform poorly when asked to extrapolate far beyond the bounds of their training set. Deep learning models, while flexible and powerful, typically (outside of specific architectures) do not provide confidence intervals or estimates of uncertainty in their predictions. Gaussian process (GP) models, as discussed above, have been suggested as an alternative that does provide well-calibrated confidence intervals, but they scale poorly to large datasets without the use of approximations and are often infeasible for datasets larger than 5000 sequences.

The second key challenge is the selection of an appropriate representation for the input. Many different encoding schemes for proteins have been described in the literature, including most recently ones adopting language models; some of these are antibody- specific and some are general to protein sequences. There is, however, little consensus on which of these is most appropriate for a given problem. The classic one-hot encoding scheme is simple to implement but unnecessarily high-dimensional and uninformative since every amino acid is treated as being completely different from each other.

Thus, high-affinity antibodies are often identified through directed evolution, which may require many iterations of mutagenesis and selection to find an optimal candidate. Deep learning techniques hold the potential to accelerate this process, but the existing methods cannot provide the confidence interval or uncertainty needed to assess the reliability of the predictions.

Programmed death-ligand 1 (PD-L1 ) is a protein that binds to another protein called PD-1 as part of a pathway that in humans suppresses the immune response to a cell recognized as foreign. Many cancers overexpress PD-L1 to downregulate the immune response and prevent their own destruction by the immune system. Interfering with PD-L1 mediated inhibition of the immune response has been shown to be a highly effective approach for treating some cancers, and multiple monoclonal antibodies (mAbs) that bind PD-L1 have been approved by the FDA for a variety of indications, including Atezolizumab, Durvalumab and Avelumab. However, there is a need to develop more potent and high affinity antibodies against PD-L1 .

SUMMARY

The present disclosure provides a novel pipeline that overcomes the limitations in the state of art as stated above. In certain embodiments, the present disclosure provides a new pipeline that includes a novel machine-learning/AI method for designing and engineering proteins with desired properties. More specifically, in an application to identify high affinity antibodies, a new autoencoder model trained on over 3 million human B-cell receptor sequences to effectively encode antibody sequences is used. A variational Bayesian neural network is then developed to perform ordinal regression on the directed evolution sequences and quantify their likelihood to be tight binders against the antigen. Importantly, this model can also assess sequences that are not present in the directed evolution libraries and thus greatly expand the search space to uncover the best sequences for further evaluations.

In certain embodiments, the power of this pipeline was demonstrated by achieving a 17-fold improvement in the KD of the PD-L1 antibody Atezolizumab and this success illustrates the potential of the pipeline in facilitating general antibody development and other protein design/engineering tasks. In certain embodiments, the present disclosure provides a set of twenty-one (21 ) scFv sequences (a fusion of the heavy and light chain variable regions) that tightly binds to PD-L1 . These scFvs are all modified variants of the parent Atezolizumab mAb and are discovered through the machine learning-assisted pipeline disclosed herein for discovering antibodies that are tight binders to targets of interest. The measured binding affinity for the top candidate from this set of the twenty- one (21 ) is about 17-fold tighter than Atezolizumab, which means it has one of the tightest affinities of any current anti-PD-L1 mAbs. These scFvs appear to have a slower dissociation rate (off-rate) than Durvalumab and Avelumab, as well. Given its exceptionally tight binding, these scFvs may have potential either as a therapeutic or as a tool for in vitro oncology work.

Since the scFvs provided herein work by forming a tight, reversible, noncovalent interaction with the PD-L1 protein expressed on the surface of human cells, this may enable them to treat cancer by interfering with the PD-L1/PD-1 interaction that some cancers use to suppress the immune response, or alternatively to function as a tool for in vitro assays involving cancer cells that overexpress PD-L1 .

In addition to possible use as an antibody therapeutic or an in vitro tool, the scFvs provided herein could be useful for CAR-T - based oncology therapies. CAR-T cells are immune cells modified to recognize and attack cancer cells using proteins overexpressed on the surface of specific cancers. The field of CAR-T based therapeutics for cancer has recently seen significant growth, with the first FDA approval in 2017. CAR-T cells could potentially be engineered to express the exceptionally tight binding anti-PD-L1 antibody scFv of the present disclosure, and thereby recognize cancer cells that overexpress PD- L1. Indeed, CAR-T cells expressing an anti-PD-L1 nanobody have been reported to exhibit antitumor activity in vitro and in mouse models, suggesting this may be a viable strategy. Given the extremely tight binding exhibited by the scFvs of the present disclosure, these scFvs may be exceptionally well-suited to this task.

Thus, the present disclosure provides design and engineering proteins for therapeutics (such as antibodies), imaging probes (such as protein domains), enhancing efficiency of CRISPR and enzyme activities. Other systems, methods, features, and advantages of the present disclosure can be or become apparent to one with skill in the art upon examination of the following drawings and detailed description. It is intended that all such additional systems, methods, features, and advantages be included within this description, be within the scope of the present disclosure, and be protected by the accompanying claims. In addition, all optional and preferred features and modifications of the described embodiments are usable in all aspects of the disclosure taught herein. Furthermore, the individual features of the dependent claims, as well as all optional and preferred features and modifications of the described embodiments are combinable and interchangeable with one another.

BRIEF DESCRIPTION OF THE DRAWING

Many aspects of the present disclosure can be better understood with reference to the following drawings (also “Figures” or “FIGs”). The components in the drawings are not necessarily to scale, emphasis instead being placed upon clearly illustrating the principles of the present disclosure. Moreover, in the drawings, like reference numerals designate corresponding parts throughout the several views.

FIGs. 1 A-1 B. Number of epochs (passes over the full dataset) required to fit. In all cases, a relative tolerance of 1 e-6 is used as the threshold for convergence.

FIGs. 2A-2D. Relative error vs number of epochs: convergence rate of stochastic gradient descent and CG with and without preconditioning. AMSGrad was used here since in initial experiments, the closely related Adam stochastic gradient descent method popular for fitting neural networks was even worse (failed to reduce error significantly in less than 100 epochs).

FIG. 3. Schematic of the Fast-Conv-1d kernel as compared with a standard 3-layer convolutional neural network.

FIG. 4. A comparison of the acceleration achieved for conjugate gradients using a preconditioner constructed using SRHT or using the unmodified Algorithm S2.

FIG. 5. The number of iterations convergence across various benchmark datasets using various settings for preconditioner rank and three different variants on the preconditioner construction algorithm.

FIGs. 6A-6B. Number of epochs I iterations required for convergence with CG on different datasets as a function of beta I lambda^A2.

FIG. 7. Distribution of absolute percent error for marginal likelihood estimation across 112 evaluations using either 25 or 40 probe vectors.

FIG. 8. Kernel PGA on the QM9 test set, using the GraphConv-1d kernel.

FIGs. 9A-9B. Kernel k-means clustering of the QM9 training set, using the GraphRBF kernel.

FIG.10. Test query molecule 1 , together with the most similar training molecules.

FIG. 11 . Test query molecule 2, together with the most similar training molecules.

FIG. 12. Test query molecule 3, together with the most similar training molecules.

FIG. 13. Test query molecule 4, together with the most similar training molecules.

FIG. 14. Test query molecule 5, together with the most similar training molecules.

FIG. 15. Schematic illustration of the RESP computational and experimental pipeline, U-PDL1 and B-PDL1 are unlabelled and biotin-labeled PDL1 , respectively.

FIGs. 16A1-2 & 16B1-4. Sorting of the naive library by FACS and quality assessment of each enriched group by flow cytometry. FIGs. 16A1-2. Sorting of the naive Atezolizumab scFv library on the yeast surface compared to the WT scFv (after 29 hours of off-rate competition). The right (highest) gate was set to collect mutants with slower off-rates while the middle (medium ate was set for WT-level off-rates and the left (lowest gate set for faster off-rates. The y-axis represents scFv expression (measured by the level of V5 epitope/AF647) while the x-axis is PD-L1 binding/PE intensity. % values represent each group’s percentage of the total number of sorted cells. I6B1-4. Testing each sorted population for binding intensity relative to the WT after 8 hours of off-rate competition. LB are low-binders (faster off-rates), MB are medium binders (WT off-rates), HB are high binders (slower off-rates). Analyzed 10⁵ cells per group, mean fluorescent intensity (MFI) on the x-axis is given per plot, V5 expression (y-axis) from AF647 intensity while binding (x-axis) in PE intensity. FIG. 17. In silico directed evolution using approximate Bayesian inference and learned embeddings for efficient candidate sequence evaluation. Figure discloses ATSW as SEQ ID NO: 69, AESW as SEQ ID NO: 70, CTSW as SEQ ID NO: 71 , and CESW as SEQ ID NO: 72.

FIG. 18. The structure of the task-adapted autoencoder. The encoder module generates a learned representation for each input sequence; the decoder module tries to reconstruct the input sequence, while the classifier generates a binary prediction for decoy vs. human B-cell receptor in its training set.

FIG. 19A-19C. Scoring the human-ness of decoys and original sequences from the autoencoder training set. 47,772 sequences were randomly sampled from the autoencoder training set (half decoy, half human). These were then scored for human- ness using the AbLSTM model (FIG. 19A), the ANARCI tool (FIG. 19B), and the BioPhi model (FIG. 19C) from the literature. In all three cases, the model's score for decoys is significantly different from that for non-decoys, and the decoys are less human than the original sequences. In all three cases, using the two-sided Mann-Whitney U test as implemented in Python’s Scipy library version 1.5.4, the calculated p-value is 0.0 (meaning that it is approximately 0 given floating point error). The following conventions apply for each boxplot. The upper and lower bounds of the box are the 25th and 75th percentile of the data, and the whiskers are drawn at 1 .5x the interquartile range (the distance from the 25th percentile to the 75th percentile). The center is drawn at the median of the data, and the “notch” represents the 95% confidence interval on the median (as determined by nonparametric bootstrap). The diamonds represent “flier” points which lie outside 1 .5x the interquartile range. Four asterisks indicate the p-value is < 0.0001 . Source data are provided as a source data file.

FIGs. 20AI-2-20C. Analysis of pipeline results. FIGs. 2OA1-2. The per-position marginal distributions in the two main clusters. Each main cluster has several sub-clusters as apparent from the marginal distribution. Only positions mutated in either cluster are displayed. Each position is numbered first using numbering from the wild-type sequence and then in parenthesis using Chothia numbering. FIG. 20B illustrates the distribution of scores assigned to sequences in each binding category by the variational Bayesian model, which is designed to perform ordinal regression (classification on ranked categories). The experimental data does not directly measure off-rates but rather classifies sequences into three bins: RH01 (high off rate), RH02 (moderate off rate), RH03 (low off rate). The higher the score assigned by the ordinal regression model, the more confident the model that the sequence has a low off rate, while lower scores indicate greater confidence in a high off rate. FIG. 20C. The accepted scores vs iteration for a typical simulated annealing chain. The algorithm initially explores sequence space impartially and as the temperature is reduced is gradually forced to focus on the most promising regions it has found thus far. Source data are provided as a source data file for all panels.

FIG. 21. Location of mutations in the 21 mutants in the structure of the Atezolizumab heavy chain. Mutated residues are labeled in the heavy chain. Structure from RCSB Protein Data Bank (5XXY) in Mol* Viewer^{46 71}.

FIGs. 22A-22C. Experimental validation of the Kott and KD on the yeast surface. FIG. 22A. WT vs. Mutant 4 scFv dissociation after 92h at RT on the yeast surface, T1/2 is half-life. FIG. 22B. Comparison of WT Atezolizumab, Mutant 4, Durvalumab, and Avelumab scFv dissociation over 28h at RT on the yeast surface. FIG. 22C. Binding affinity (KD) measurements determined on the yeast surface between scFv and PD-L1 (3 independent measurements, also see FIGs. 29A-29F).

FIGs. 23AI-2 &23B. The variational Bayesian neural network architecture. FIG. 23AI-2. A comparison of a variational Bayesian neural network with a fully connected architecture. Both networks map a vector of input values x1 , x2,...xn to a hidden layer vector H1_1 , H1_2,...H1_n by a matrix multiplication followed by a nonlinear activation. In the fully connected network, however, each element of the hidden layer weight matrix is a learned value and once the network is trained it is a fixed value. In the Bayesian network, by contrast, each element of the hidden layer weight matrix is a Gaussian distribution specified by a learned mean and a learned standard deviation. To generate predictions, they can be sampled from the weight distributions, which provides an estimate of the uncertainty on the predictions. FIG. 23B. The structure of the ordinal regression scoring model used in the pipeline. FIGs. 24A-24B. WT Atezolizumab scFv binds human PD-L1 on the yeast surface. The y-axis of each panel represents expression of the scFv (V5 epitope/AF647), while the x-axis the binding intensity to human PD-L1/PE intensity (Left/24A: 0 nM PD-L1 vs Right/24B: 6.1 nM PD-L1 ). Numbers of double positive cells out of the total analyzed are given below each panel.

FIGs. 25A-25D. First screen of the 21 -member library after 8 hours off-rate competition (relative to WT, top plots 25A & 25B) and final screen of the library after 39 hours competition (bottom plots 25C & 25D). The right triangle represents the sort gate with associated % of cells collected per sort. The y-axis of each panel represents expression of the scFv (V5 epitope/AF647), while the x-axis the binding intensity to human PD-L1/PE intensity.

FIGs. 26A-I-2-26B. BLI binding analysis of soluble WT and Mutant 4 scFv against immobilized biotin PD-L1. FIGs. 26A1-2. Binding curve (association and dissociation) profiles for different scFv concentrations. WT was titrated over a range from 0.65 to 475 nM while Mutant 4 over a range of 0.219 to 160 nM. FIG. 26B. On/Off-rates and KD for WT and Mutant 4 scFv. For each scFv, two different KD values were obtained from 2 binding events against PD-L1 . The Mutant 4 scFv values were not obtainable due to the very slow off-rate.

FIGs. 27A-27G. Mass photometry analysis of the monomer/dimer/multimer state of PD-L1 and WT/Mutant 4 scFv. Each panel image has approximate MW (x-axis) and counts (y-axis). PD-L1 MW (theoretical) = 35-38 kDa (range due to glycosylation) and scFv = 31.4 kDa. Batches 1 and 2 are scFv from 2 separate purifications. The sigma symbol represents the standard deviation of the MW estimate as calculated from the software, the top number above each peak is apparent MW.

FIGs. 28A-28B. SDS-PAGE gels (4-20%) of purified scFv-Hise fusions. M4 is mutant 4. The two gel images represent proteins from 2 separate purifications (2 separate SDS-PAGE gels from 2 separate purifications).

FIGs. 29A-29F. Binding curves between PD-L1 and yeast surface scFv. MFI means fluorescent intensity, M4 means Mutant 4. Each experiment is a separate/independent assay performed on a separate day. FIGs. 30A-30B. Side-by-side comparison of amino acids close to R98C in the wild type (PDB structure 5XXY, FIG. 30A, top) and the IgFold predicted structure (FIG. 30B, bottom) for the top scoring mutant.

FIG. 31 . Overlay of PDB structure 5XXY with the IgFold predicted structure for the top-scoring mutant. Note the large shift in conformation in CDR H3 and the more subtle shift in conformation in the other region that makes contact with the antigen.

FIGs. 32A-32C. Model assigned scores for sequences from Mason et aL, Nat. Biomed. Eng. 5, 600-612, 2021. FIG. 32A. The distribution of scores assigned by the model for binding and nonbinding sequences from the training set. 0 is the cutoff at which a sequence is predicted to be a binder. FIG. 32B. The distribution of scores for sequences selected by Mason et al. as likely binders. Experimentally Mason et al. demonstrated that all of these sequences bind with KD < 20 nM, and all of them are indeed predicted to be binders by the model disclosed Example 2. FIG. 32C. The model’s uncertainty regarding its assigned score for predictions on the test set. Incorrect predictions have higher associated uncertainty than correct predictions, indicating that as anticipated, model assigned uncertainty can assist in determining whether a prediction should be considered reliable. There are 3,158 unique sequences in this test set. Significance was assessed using a two-sided Mann-Whitney II test; the resulting p-value was 1 e-23. The following conventions apply for each boxplot. The upper and lower bounds of the box are the 25th and 75th percentile of the data, and the whiskers are drawn at 1 .5x the interquartile range (the distance from the 25th percentile to the 75th percentile). The center is drawn at the median of the data, and the “notch” represents the 95% confidence interval on the median (as determined by nonparametric bootstrap). The diamonds represent “flier” points which lie outside 1 .5x the interquartile range. Four asterisks indicate the p-value is < 0.0001 .

FIG. 33. Clustering of simulated annealing results. The dendrogram for clustering of sequences harvested from the simulated annealing procedure after sequences with scores associated with wide confidence intervals had been removed. The dendrogram suggests the presence of at least two main regions of sequence space identified by the modified simulated annealing algorithm.

FIG. 34. Example of FACS Gating Strategy. 10⁶ cells (“hi” gate) collected out of 6x10⁶ sorted yeasts. Y-axis is V5 epitope expression and X-axis is PD-L1 binding. V5 tag monoclonal mAb detected with secondary-AF647 polyclonal IgG. Biotin-PD-L1 detected with SA-PE reagent.

Additional advantages of the present disclosure are set forth in part in the description which follows, and in part could be obvious from the description, or can be learned by practice of the disclosure. The advantages of the disclosure could be realized and attained by means of the elements and combinations particularly pointed out in the appended claims. It is to be understood that both the foregoing general description and the following detailed description are exemplary and explanatory only and are not restrictive of the invention, as claimed.

DETAILED DESCRIPTION

Gaussian process (GP) is a Bayesian model which provides several advantages for regression tasks in machine learning such as reliable quantitation of uncertainty and improved interpretability. Their adoption has been precluded by their excessive computational cost and by the difficulty in adapting them for analyzing sequences (e.g., amino acid sequences) and graphs (e.g., small molecules). In certain embodiments, the present disclosure provides a group of random feature-approximated kernels for sequences and graphs that exhibit linear scaling with both the size of the training set and the size of the sequences or graphs. These new kernels were incorporated into the new Python library for GP regression, xGPR, and an efficient and scalable algorithm was developed for fitting GPs equipped with these kernels to large datasets. The performance of xGPR was compared on 17 different benchmarks with both standard and state of the art deep learning models and it was found that GP regression achieves highly competitive accuracy for these tasks while providing with well-calibrated uncertainty quantitation and improved interpretability. Further, the present disclosure provides how xGPR can be used as part of an active learning strategy to engineer a protein with a desired property in an automated way without human intervention.

The present disclosure also provides a novel pipeline called RESP that is composed of direct evolution, sequencing, and machine learning analysis for efficient identification of high affinity antibodies. In certain embodiments, a learned representation trained on over 3 million human B-cell receptor sequences to encode antibody sequences was developed. A variational Bayesian neural network was also developed to perform ordinal regression on a set of the directed evolution sequences binned by off-rate and quantify their likelihood to be tight binders against an antigen. Importantly, this model can assess sequences that are not present in the directed evolution library and thus greatly expands the search space to uncover the best sequences for evaluations.

The potential of RESP in facilitating general antibody development is further illustrated. In certain embodiments, this pipeline was applied to improve the affinity of a well-known antibody to a well-studied antigen by engineering mutants of the heavy chain of Atezolizumab (brand name Tecentriq) with improved affinity for programmed death ligand 1 (PD-L1 ). To do so, a large portion of the Atezolizumab heavy chain sequence was randomized and the mutants were binned using yeast display and FACS, and each bin was sampled for sequencing. The mutant Atezolizumab sequence data was converted to a low-dimensional representation using the autoencoder model, and the encoded mutant Atezolizumab sequences were used to train a Bayesian ordinal regression model that scores each sequence on the probability of its strong binders to PD-L1 . Finally, a modified simulated annealing algorithm was used to select sequences for testing.

In certain embodiments, the present disclosure provides a set of scFv sequences comprising a fusion sequence of the heavy and light chain variable regions of Atezolizumab mAb. The scFvs disclosed herein are discovered using the machine learning-assisted pipeline, discussed herein, and show tight binding to PD-L1 . One of the scFvs disclosed herein shows a 17-fold improvement in the KD of the PD-L1 antibody Atezolizumab and a 10-fold decrease in off-rate. This scFvs disclosed herein could be more potent than any existing antibodies for PD-L1 and can be used in career immunotherapy as a potent antibody against PD-L1 , as well as in CAR-T cell therapy.

Many modifications and other embodiments disclosed herein will come to mind to one skilled in the art to which the disclosed compositions and methods pertain having the benefit of the teachings presented in the foregoing descriptions and the associated drawings. Therefore, it is to be understood that the disclosures are not to be limited to the specific embodiments disclosed and that modifications and other embodiments are intended to be included within the scope of the appended claims. The skilled artisan will recognize many variants and adaptations of the aspects described herein. These variants and adaptations are intended to be included in the teachings of this disclosure and to be encompassed by the claims herein.

Although specific terms are employed herein, they are used in a generic and descriptive sense only and not for purposes of limitation.

As will be apparent to those of skill in the art upon reading this disclosure, each of the individual embodiments described and illustrated herein has discrete components and features which may be readily separated from or combined with the features of any of the other several embodiments without departing from the scope or spirit of the present disclosure.

Any recited method can be carried out in the order of events recited or in any other order that is logically possible. That is, unless otherwise expressly stated, it is in no way intended that any method or aspect set forth herein be construed as requiring that its steps be performed in a specific order. Accordingly, where a method claim does not specifically state in the claims or descriptions that the steps are to be limited to a specific order, it is no way intended that an order be inferred, in any respect. This holds for any possible non-express basis for interpretation, including matters of logic with respect to arrangement of steps or operational flow, plain meaning derived from grammatical organization or punctuation, or the number or type of aspects described in the specification.

All publications and patents cited in this specification are cited to disclose and describe the methods and/or materials in connection with which the publications are cited. All such publications and patents are herein incorporated by references as if each individual publication or patent were specifically and individually indicated to be incorporated by reference. Such incorporation by reference is expressly limited to the methods and/or materials described in the cited publications and patents and does not extend to any lexicographical definitions from the cited publications and patents. Any lexicographical definition in the publications and patents cited that is also not expressly repeated in the instant application should not be treated as such and should not be read as defining any terms appearing in the accompanying claims. The citation of any publication is for its disclosure prior to the filing date and should not be construed as an admission that the present disclosure is not entitled to antedate such publication by virtue of prior disclosure. Further, the dates of publication provided could be different from the actual publication dates that may need to be independently confirmed.

While aspects of the present disclosure can be described and claimed in a particular statutory class, such as the system statutory class, this is for convenience only and one of skill in the art will understand that each aspect of the present disclosure can be described and claimed in any statutory class.

It is also to be understood that the terminology used herein is for the purpose of describing particular aspects only and is not intended to be limiting. Unless defined otherwise, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which the disclosed compositions and methods belong. It will be further understood that terms, such as those defined in commonly used dictionaries, should be interpreted as having a meaning that is consistent with their meaning in the context of the specification and relevant art and should not be interpreted in an idealized or overly formal sense unless expressly defined herein.

Aspects of the present disclosure will employ, unless otherwise indicated, techniques of molecular biology, microbiology, organic chemistry, biochemistry, physiology, cell biology, blood vessel biology, and the like, which are within the skill of the art. Such techniques are explained fully in the literature.

Prior to describing the various aspects of the present disclosure, the following definitions are provided and should be used unless otherwise indicated. Additional terms may be defined elsewhere in the present disclosure.

Definitions

As used herein, “comprising” is to be interpreted as specifying the presence of the stated features, integers, steps, or components as referred to, but does not preclude the presence or addition of one or more features, integers, steps, or components, or groups thereof. Moreover, each of the terms “by”, “comprising,” “comprises”, “comprised of,” “including,” “includes,” “included,” “involving,” “involves,” “involved,” and “such as” are used in their open, non-limiting sense and may be used interchangeably. Further, the term “comprising” is intended to include examples and aspects encompassed by the terms “consisting essentially of” and “consisting of.” Similarly, the term “consisting essentially of” is intended to include examples encompassed by the term “consisting of.

As used herein, the term “and/or” includes any and all combinations of one or more of the associated listed items. Expressions such as “at least one of,” when preceding a list of elements, modify the entire list of elements and do not modify the individual elements of the list.

As used in the specification and the appended claims, the singular forms “a,” “an” and “the” include plural referents unless the context clearly dictates otherwise. Reference to "a/an" chemical compound, therapeutic agent, and pharmaceutical composition each refers to one or more molecules of the chemical compound, therapeutic agent, and pharmaceutical composition rather than being limited to a chemical compound, therapeutic agent, and pharmaceutical composition, the one or more molecules may or may not be identical, so long as they fall under the category of the chemical compound, therapeutic agent, and pharmaceutical composition. Thus, for example, "a" therapeutic agent is interpreted to include one or more molecules of the therapeutic agent, where the therapeutic agent molecules may or may not be identical (e.g., comprising different isotope abundances and/or different degrees of hydration or in equilibrium with different conjugate base or conjugate acid forms).

It should be noted that ratios, concentrations, amounts, and other numerical data can be expressed herein in a range format. It will be further understood that the endpoints of each of the ranges are significant both in relation to the other endpoint, and independently of the other endpoint. It is also understood that there are a number of values disclosed herein, and that each value is also herein disclosed as "about" that particular value in addition to the value itself. For example, if the value "10" is disclosed, then "about 10" is also disclosed. Ranges can be expressed herein as from "about" one particular value, and/or to "about" another particular value. Similarly, when values are expressed as approximations, by use of the antecedent "about," it will be understood that the particular value forms a further aspect. For example, if the value "about 10" is disclosed, then "10" is also disclosed.

Where a range is expressed, a further aspect includes from the one particular value and/or to the other particular value. Where a range of values is provided, it is understood that each intervening value, to the tenth of the unit of the lower limit unless the context clearly dictates otherwise, between the upper and lower limit of that range and any other stated or intervening value in that stated range, is encompassed within the disclosure. The upper and lower limits of these smaller ranges may independently be included in the smaller ranges and are also encompassed within the disclosure, subject to any specifically excluded limit in the stated range. Where the stated range includes one or both of the limits, ranges excluding either or both of those included limits are also included in the disclosure. For example, where the stated range includes one or both of the limits, ranges excluding either or both of those included limits are also included in the disclosure, e.g. the phrase "x to y" includes the range from 'x' to 'y' as well as the range greater than 'x' and less than 'y'. The range can also be expressed as an upper limit, e.g. 'about x, y, z, or less' and should be interpreted to include the specific ranges of 'about x', 'about y', and 'about z' as well as the ranges of 'less than x', less than y', and 'less than z'. Likewise, the phrase 'about x, y, z, or greater' should be interpreted to include the specific ranges of 'about x', 'about y', and 'about z' as well as the ranges of 'greater than x', greater than y', and 'greater than z'. In addition, the phrase "about 'x' to 'y'", where 'x' and 'y' are numerical values, includes "about 'x' to about 'y'".

It is to be understood that such a range format is used for convenience and brevity, and thus, should be interpreted in a flexible manner to include not only the numerical values explicitly recited as the limits of the range, but also to include all the individual numerical values or sub-ranges encompassed within that range as if each numerical value and sub-range is explicitly recited. To illustrate, a numerical range of "about 0.1% to 5%" should be interpreted to include not only the explicitly recited values of about 0.1% to about 5%, but also include individual values (e.g., about 1%, about 2%, about 3%, and about 4%) and the sub-ranges (e.g., about 0.5% to about 1 .1 %; about 5% to about 2.4%; about 0.5% to about 3.2%, and about 0.5% to about 4.4%, and other possible sub-ranges) within the indicated range.

As used herein, "about," "approximately," "substantially," and the like, when used in connection with a numerical variable, can generally refers to the value of the variable and to all values of the variable that are within the experimental error (e.g., within the 95% confidence interval for the mean) or within +/- 10% of the indicated value, whichever is greater. As used herein, the terms "about," "approximate," "at or about," and "substantially" can mean that the amount or value in question can be the exact value or a value that provides equivalent results or effects as recited in the claims or taught herein. That is, it is understood that amounts, sizes, formulations, parameters, and other quantities and characteristics are not and need not be exact but may be approximate and/or larger or smaller, as desired, reflecting tolerances, conversion factors, rounding off, measurement error and the like, and other factors known to those of skill in the art such that equivalent results or effects are obtained. In some circumstances, the value that provides equivalent results or effects cannot be reasonably determined. In general, an amount, size, formulation, parameter or other quantity or characteristic is "about," "approximate," or "at or about" whether or not expressly stated to be such. It is understood that where "about," "approximate," or "at or about" is used before a quantitative value, the parameter also includes the specific quantitative value itself, unless specifically stated otherwise.

As used herein, the terms "optional" or "optionally" means that the subsequently described event or circumstance can or cannot occur, and that the description includes instances where said event or circumstance occurs and instances where it does not.

As used interchangeably herein, "subject," "individual," or "patient" can refer to a vertebrate organism, such as a mammal (e.g., human). "Subject" can also refer to a cell, a population of cells, a tissue, an organ, or an organism, preferably to human and constituents thereof.

As used herein, the terms "treating" and "treatment" can refer generally to obtaining a desired pharmacological and/or physiological effect. The effect can be, but does not necessarily have to be, prophylactic in terms of preventing or partially preventing a disease, symptom or condition thereof. The effect can be therapeutic in terms of a partial or complete cure of a disease, condition, symptom or adverse effect attributed to the disease, disorder, or condition. The term "treatment" as used herein can include any treatment of inflammation associated with any disease in a subject, particularly a human and can include any one or more of the following: (a) preventing the disease from occurring in a subject which may be predisposed to the disease but has not yet been diagnosed as having it; (b) inhibiting the disease, i.e., arresting its development; and (c) relieving the disease, i.e., mitigating or ameliorating the disease and/or its symptoms or conditions. The term "treatment" as used herein can refer to both therapeutic treatment alone, prophylactic treatment alone, or both therapeutic and prophylactic treatment. Those in need of treatment (subjects in need thereof) can include those already with the disorder and/or those in which the disorder is to be prevented. As used herein, the term "treating", can include inhibiting the disease, disorder or condition, e.g., impeding its progress; and relieving the disease, disorder, or condition, e.g., causing regression of the disease, disorder and/or condition. Treating the disease, disorder, or condition can include ameliorating at least one symptom of the particular disease, disorder, or condition, even if the underlying pathophysiology is not affected, e.g., such as treating the pain of a subject by administration of an analgesic agent even though such agent does not treat the cause of the pain.

As used herein, the term "therapeutically effective amount" refers to an amount that is sufficient to achieve the desired therapeutic result or to have an effect on undesired symptoms but is generally insufficient to cause adverse side effects. The specific therapeutically effective dose level for any particular patient will depend upon a variety of factors including the disorder being treated and the severity of the disorder; the specific composition employed; the age, body weight, general health, sex and diet of the patient; the time of administration; the route of administration; the rate of excretion of the specific compound employed; the duration of the treatment; drugs used in combination or coincidental with the specific compound employed and like factors within the knowledge and expertise of the health practitioner and which may be well known in the medical arts. In the case of treating a particular disease or condition, in some instances, the desired response can inhibit the progression of the disease or condition. This may involve only slowing the progression of the disease temporarily. However, in other instances, it may be desirable to halt the progression of the disease permanently. This can be monitored by routine diagnostic methods known to one of ordinary skilled in the art for any disease. The desired response to treatment of the disease or condition also can be delaying the onset or even preventing the onset of the disease or condition.

It is understood that unless otherwise specified, temperatures referred to herein are based on atmospheric pressure (i.e., one atmosphere).

Now having described the aspects of the present disclosure, in general, the following Examples describe some additional and/or more detailed aspects of the present disclosure. While aspects of the present disclosure are described in connection with the following examples and the corresponding text and figures, there is no intent to limit aspects of the present disclosure to this description. On the contrary, the intent is to cover all alternatives, modifications, and equivalents included within the spirit and scope of the present disclosure.

EXAMPLES

The following examples are put forth to provide those of ordinary skill in the art with a complete disclosure and description of how the compounds, compositions, articles, devices and/or methods claimed herein are made and evaluated and are intended to be purely exemplary of the disclosure and are not intended to limit the scope of what the inventors regard as their disclosure. Efforts have been made to ensure accuracy with respect to numbers (e.g., amounts, temperature, etc.), but some errors and deviations should be accounted for. Unless indicated otherwise, parts are parts by weight, temperature is in °C or is at ambient temperature, and pressure is at or near atmospheric.

EXAMPLE 1

Linear-scaling kernels for protein sequences and small molecules outperform deep learning while providing uncertainty quantitation and improved interpretability

In this Example, approximate Gaussian process models were disclosed that can achieve good scalability and performance for sequences and graphs. First, efficient methods were introduced for fitting random features-approximated Gaussian processes using preconditioned conjugate gradients (PCG) with a new preconditioner. The scalability of the resulting approximate GP to larger datasets was demonstrated, showing that it outperforms an alternative method for approximating Gaussian processes on 9 different datasets, and find that it achieves the same test set performance as an exact (non-approximated) GP on a protein engineering task.

Next, a group of efficient, scalable, random feature-approximated kernels were introduced for graphs and sequences. The test set performance of GP regression was compared using these kernels with deep learning methods on 17 different benchmarks. The GP regression models achieve highly competitive performance.

The possibility of combining deep learning with GPs was further explored by using embeddings as input to the approximate GP and it showed that this is frequently beneficial. It was further demonstrated that uncertainty estimates provided by an approximate GP are significantly more well-calibrated than those provided by a variety of uncertainty estimation techniques for deep learning. In a simple study, it was shown how an approximate GP can be used to “engineer” a protein with improved fitness with no human intervention required. Further, how the approximate Gaussian process can be used to cluster the training set and retrieve the datapoints in the training set most similar to a test point was also illustrated.

METHODS

The random features approximation: background

In a typical Gaussian process model, given a training set X with associated y- values, predictions were made for the y-value associated with a new datapoint x,: where is a Gaussian distribution with mean and variance of it,

Here k* is a vector of length N for N training datapoints formed by evaluating the kernel function k(x„x_z) for each x_z in the training set, and K is an N x N matrix where element ij is formed by evaluating the kernel function for kfa, xj points in the training set. A is a kernel hyperparameter which denotes the level of expected noise in the data, k* is the kernel function evaluated on the new datapoint against all training datapoints, and k„ is the kernel function evaluated on the new datapoint against itself. When fitting, no need to “learn” any parameters, but do need to construct and invert an Nx N matrix, which is expensive if N is large.

To avoid this cost, kernel functions were approximated using random Fourier features²⁷. This well-known approximation specifies a random map for a given kernel such that fc(x₁,x₂) ~ z(x₁)^Tz(x₂) (for details, see Supporting Information Section S1 ). The mean and variance are now predicted for new datapoints using:

Where Z is the transformed input matrix, such that row / is z(x_;) and z, is z(x.) where x* is a new datapoint, A is a kernel hyperparameter. For common stationary kernels, the error of the approximation decreases exponentially with an increasing number of random features, so that there are diminishing returns (increasing the number of random features from 1 ,000 to 2,000 achieves a larger relative improvement than going from 2,000 to 4,000 and so forth).

Further, vanilla random features described above were replaced with the structured orthogonal random features (SORF) procedure. Briefly (for details, see Supporting Information Section S2), with the RBF kernel as an example, the following transformation is applied to each input vector:

SHD₁HD₂HD₃ (4)

Where S is a diagonal matrix with elements drawn from a j-distribution with d degrees of freedom, H is the normalized Hadamard matrix, and D_r, D₂, D₃ are diagonal matrices with elements drawn from the Rademacher distribution. The Hadamard matrix multiplication can be replaced with a fast Hadamard transform, reducing the cost of generating random features from 0(NM²) to 0(NM logM). The resulting model is very lightweight, since only four diagonal matrices were needed to store. Generating predictions with a trained model is fast, since only a series of transforms & diagonal matrix multiplications needed to be performed, then a dot product was taken. This alteration is not only beneficial for speed but also for performance. It was demonstrated that this modification improves the accuracy of kernel approximation with random features for the RBF kernel. The fast Hadamard transform for both CPU and GPU was implemented in the library (see Supporting Information section S3 for benchmarking and other details).

Faster fitting with a new preconditioner

To fit the approximate GP model, \Z^TZ + Z²I]~¹Z^Ty - w must be found. If the number of random features is large (e.g. 5000 - 30000), directly forming and inverting the matrix \Z^TZ + A²/] requires expensive matrix multiplications and decompositions. This problem was avoided by instead using the well-known method of conjugate gradients (CG), which achieves subquadratic scaling in the number of random features and iteratively solves Z^Ty = \Z^TZ + A²/] w for w. Briefly (for details, see Supporting Information section S4), rather than computing \Z^TZ + A²/]^-1, this algorithm requires that on each iteration, the matrix-vector product \Z^TZ + A²/]w was computed. No matrix multiplications or decompositions are required at any time.

Note that \Z^TZ + A²/]w = (^t ^zi^Tzi + A²/)w for N row vectors in Z; in other words, the matrix vector product can be formed as a sum of products over individual rows. The dataset was then processed in minibatches. For each minibatch, random features were generated using the Hadamard transform-based procedure and a matrix-vector multiplication was performed with w, which was then added to the result. This implementation also allows easy parallelization. On each iteration of conjugate gradients, each worker can be assigned its own copy of the kernel, its own copy of w, and its own subset of the data to process. At the end of the iteration, the results from all the workers are summed and used to update w for the next iteration.

Note that unlike with stochastic gradient descent-based fitting methods, minibatch size has no effect on model performance or on the number of iterations required for convergence, because CG evaluates the full gradient on each iteration. Larger minibatches can slightly reduce the time required to fit but do increase memory footprint. Thus, the memory footprint is independent of the dataset size, and can be determined by the user by choosing minibatch size. Since CG is iterative, it is important to minimize the number of iterations to make model fitting as fast as possible. xGPR uses the randomized Nystrom approximation of the matrix \Z^TZ + A²/] in equation (1 ) as a preconditioner (for details, see Supporting Information Section S4). While this is a known method for matrix approximation, it has not previously been used as a preconditioner for Gaussian processes.

A modification was also added to the randomized Nystrom approximation that reduces the time needed to construct the preconditioner (see Supporting Information Section S4 for details and for benchmarks). Importantly, in this implementation the preconditioner is constructed by loading only one minibatch of the data into memory at a time, minimizing the memory footprint. The details on how to use random features and more specifically preconditioned conjugate gradients for hyperparameter tuning are provided in the below Supporting Information S5 and S6.

In FIGs. 1A & 1 B, the preconditioned conjugate gradient (PCG) approach was compared with the Scipy library’s implementation of L-BFGS across 7 datasets; L-BFGS (limited memory BFGS) is a popular algorithm for minimizing the loss function with an iteratively constructed approximation to the Hessian. The PCG strategy reduces the number of iterations required for convergence and hence the time to fit by as much as two orders of magnitude compared with L-BFGS or with non-preconditioned CG, enabling to achieve a tight fit in as few as 10-20 iterations. In all cases, a relative tolerance of 1 e- 6 was used as the threshold for convergence. For details on preconditioner rank, see below Supporting Information S4.

In FIGs. 2A-2D, the PCG approach was further compared with two stochastic gradient descent methods, stochastic variance reduced gradient descent (SVRG) and AMSGrad, an improved variant of the Adam algorithm which is popular for training neural networks. It was illustrated how the residual error decreases with the number of epochs (full passes over the dataset). For PCG, each iteration is equivalent to one epoch, while for stochastic gradient descent, there are many iterations per epoch. As illustrated, PCG is the clear winner, converging faster than stochastic gradient descent by orders of magnitude and thereby greatly reducing the number of epochs required to fit. Unlike stochastic gradient descent, CG does not require tuning of a learning rate to achieve good results. The PCG approach coupled with the random feature generation method enables to fit a > 400,000 datapoint dataset in less than 7 minutes including all hyperparameter tuning. If using a minibatch size of 2,000 datapoints, the process requires less than 1.5 GB of RAM at all times.

AMSGrad was used in these studies since in initial experiments, the closely related Adam stochastic gradient descent method popular for fitting neural networks was even worse (failed to reduce error significantly in less than 100 epochs).

Note that fitting half a million datapoints using an exact Gaussian process would require constructing a 400,000 by 400,000 kernel matrix, which would require 1.3 terabytes of memory, and the entire kernel matrix would be required to fit the model. xGPR is thus much more scalable than an exact Gaussian process. Indeed, as described under “Comparison with stochastic variational inference and exact Gaussian processes”, when using an exact Gaussian process as implemented in the GPyTorch library, an out of memory error was immediately encountered when trying to fit a dataset with merely 21 ,000 training datapoints, while xGPR was able to fit half a million datapoints without any issues.

RESULTS

Comparison with stochastic variational inference and exact Gaussian processes

Random features are not the only method for approximating a Gaussian process; one alternative is stochastic variational inference (SVI), which uses a set of “inducing points” much smaller than the number of training datapoints to approximate the GP. As described in the below Supporting Information section S8, both xGPR and a GP approximated were fit using stochastic variational inference (as implemented in GPyTorch) to 9 different datasets ranging in size from 400 to half a million datapoints. The same kernel was used in all cases to ensure a fair comparison (See Supporting Information section S7 for details on these and all other datasets used for benchmarking in this Example).

These studies were performed on a GTX1070 GPU with 8 GB of RAM. Two of the datasets are sufficiently small (< 5,000 datapoints) that can be fit in an exact Gaussian process, and so an exact GP for those two datasets was fit, as well. An exact GP was tried to fit to a third dataset with about 21 ,000 datapoints, but the GPyTorch library immediately threw an out of memory error, demonstrating that an exact GP is already expensive for datasets of a fairly modest size.

It was found that xGPR achieves the same or better performance as SVI in every single case; there is no case where SVI outperforms xGPR. For the two datasets where xGPR can be compared with an exact GP, xGPR achieves very similar performance, while SVI is similar in one case and in the other significantly worse. xGPR was also compared with an exact Gaussian process on an active learning / Bayesian optimization task from the literature. For example, data were generated for 52,910 mutants of the green fluorescent protein (GFP), the fluorescence of each mutant was measured. An exact Gaussian process was fitted with an RBF kernel to this dataset, using 1 ,1 15 sequences with only one mutation as their training set and representing the input datapoints using learned embeddings. In Bayesian optimization, an acquisition function is used to select new datapoints for experimental evaluation. The upper confidence bound (UCB) was used, which is the predicted mean value of the GP, plus a confidence interval. UCB prioritizes both datapoints where the model predicts the mutant sequence will have a high fitness value and datapoints where the model has high uncertainty, thus achieving a balance between exploitation of existing knowledge and exploration of sequence space. It was shown that the UCB values for the test set are reasonably well-correlated with the actual measured fluorescence values. xGPR was fit to the same dataset, using the same training set, kernel and representations. It was found the UCB acqusition function exhibited a Spearman’s-r correlation coefficient of 0.78 with measured fluorescence in the test set. Surprisingly, the result for xGPR is slightly better at 0.8. Next, the 50 mutants from the test set with the largest UCB values (to simulate the effect of experimentally assessing the model predictions) were picked. The average log fluorescence of the top fifty for xGPR is about the same as the result reported in the literature (both are about 3.75). These results suggest that when using the same kernel, xGPR compares well with an exact Gaussian process while achieving greatly improved scalability. A first new set of efficient kernels for sequences that cannot be aligned

Up to this point, kernels for fixed-length inputs were only considered. In protein engineering problems, it is often possible to form a multiple sequence alignment from all the input sequences. Alternatively, if a learned embedding generated by a large language model is used to represent the input, then the learned embedding can be averaged over the sequence to generate a fixed-length input. In many other cases, however, the sequences are of different lengths and cannot easily be aligned - if the input sequences are from different families, for example. Two types of random features approximated kernels that can efficiently handle this problem were introduced: the first type is referred to as FHT-Conv-1d and the second is referred to as Fast-Conv-1d.

The first type of kernel can be constructed as follows. Consider extracting all k- mers of a given length, e.g., 9, from two sequences, and calculating the similarity of each possible pairing of k-mers from the two sequences using any common fixed-vector kernel of interest. There are a range of different kernels that could be used for such a comparison, including the RBF, Matern, neural network and arc-cosine kernels, all of which can be approximated using random features. The “convolutional kitchen sinks” kernel is a special instance of this type of kernel. In this case, the comparison between the two sequences is as follows:

Where k is the fixed-vector kernel that is selected to measure the similarity of any two subsequences. As written, this kernel scales as the square of the sequence length. It is therefore impractical if used in an exact Gaussian process or in stochastic variational inference, because in those cases it will need to be implemented as written, resulting in a very large number of pairwise k-mer comparisons. If comparing two protein sequences of length 250 with subsequence length 9, for example, all 242 k-mers from both sequences would need to be extracted and then over 29,000 pairwise k-mer comparisons would need to be performed.

How to use the random features approximation to implement this kernel in a way that scales linearly with sequence length is shown below. Using the random features approximation, this kernel becomes:

Where the two z-vectors are the random features generated for each pair of k- mers in sequences x and y. The last expression shown above is equivalent to:

Random features can therefore be generated for each subsequence of length k in the two input sequences and be summed, resulting in a single vector representing sequence xand sequence y, then these representations were used as inputs for Bayesian linear regression. The computational cost of this approach exhibits linear scaling with both dataset size and sequence length. This can be considered a convolution kernel since generating the random features for the k-mers in each input sequence can be accomplished through convolution with filters populated with random weights drawn from the distribution appropriate for the kernel we want to approximate. This approximated kernel still has, however, limited practicality since convolution with thousands of filters is slow.

A second innovation was then introduced that greatly accelerates the calculation of random features for these kernels. To understand this implementation, first note that in deep learning libraries convolution is often implemented using matrix multiplication as follows. Consider an input K x P matrix with P sequence elements (e.g., P amino acids) and /("features per element (or amino acid), and let us assume the convolution filter is of width L. The input matrix can be used to populate an LK x (P - L + 1) matrix where each column / is the concatenated features associated with subsequence /. The m convolution filters meanwhile can be stacked to form an mx LK matrix; multiplying this matrix against the LK x (P - L + 1) input matrix is equivalent to performing 1d convolution with m filters of length L.

In this implementation, the m x LK matrix of convolution filters was replaced with the SORF operation (equation 4), and matrix multiplications were replaced with fast diagonal matrix multiplications and fast transforms, analogous to replacing a discrete Fourier transform with FFT. For simplicity, let the number of random features desired, m, be the same as LK. Then the implementation of this kernel with vanilla random features scales as O(N (P — L + l) m²). The proposed modification by contrast scales as O(N (P — L + 1) m log (m)).

Depending on the fixed-vector kernel used to measure distances between k-mers, how encoding the amino acids in the two input sequences to this class of kernels can be chosen so that distance between k-mers is measured in a particular way. If one-hot encodes amino acids, for example, the distance between any two k-mers is twice the Hamming distance. Alternatively, each amino acid can be represented using a learned representation, e.g., from a large language model. These possibilities are described below under “Evaluating convolution kernels for sequence data”, and this first type of kernel is referred to as FHT-Conv-1d.

A second efficient convolution kernel for sequences which cannot be aligned

A second efficient group of convolution kernels were introduced which mimic a three-layer neural network, which is called Fast-Conv-1d. The following simple lemma was described (see Supporting Information Section S9 for a proof):

Lemma 1.1. Let K-. IR^D% IR^D-> IR be a positive definite kernel on IR^D, and let /: A'-> ]R^D be any mapping from X to IR^D, where X is some non-empty set. Then f<(/(u),/(w)) is a positive definite kernel for all u,w E X.

Consider a three-layer neural network with random weights where the first layer is a convolutional layer with ReLU activation followed by global max pooling and the second layer is fully connected with a selected activation function. This is analogous to a fixed vector kernel (e.g., RBF or Linear + RBF) whose input is a random feature convolution with ReLU activation and global max pooling (FIG. 3). The fixed vector kernel is approximated using a fast Hadamard transform-based random features procedure. This combination kernel is positive definite following lemma 1.1. This kernel has the property that the convolutions can all be performed exactly once before training and the results saved on disk. The fast Hadamard transform-based procedure outlined above can be used to perform the random feature convolutions in an efficient way.

To understand how this type of kernel quantifies the similarity of sequences, consider performing a convolution of a one-hot encoded sequence with a width-9 filter containing weights drawn from a normal distribution. For each subsequence of length 9, the convolution is large and positive //the filter happens to contain large positive values at the same positions that are one-hot in the subsequence. Since ReLU activation and max-pooling across the sequence were performed, for each convolution filter, the resulting feature is the largest positive value encountered - in other words, the best “match” for that filter found in the sequence. Thus, this kernel essentially “profiles” each sequence input to quantify the best match in the input for each filter in a set of random filters, then compares the profiles of any two input sequences using a specified fixed vector kernel (e.g., RBF or Linear + RBF; RBF was used in all studies for simplicity). As before, sequences can be encoded using learned representations, one-hot encoding, or some other scheme as desired.

Evaluating xGPR and convolution kernels for sequence data xGPR and the FHT-Conv-1 d and Fast-Conv-1d kernels were evaluated by fitting protein engineering benchmark data (the TAPE benchmarks and FLIP benchmarks, respectively). These datasets involve situations where a pre-constructed training set and test set consisting of protein sequences were provided, and the goal is to predict a real- valued outcome for each sequence. These datasets are briefly below (for more details on these and all other datasets used in this Example, see below Supporting Information Section S7).

The AAV dataset involves mutations including insertions and deletions to a 28- amino acid window of an adeno-associated virus capsid; some of these mutations were randomly generated while others were “designed” (selected intentionally). The data were split up in various ways and these various partitions were referred to as “splits”. For example, they train a model on the randomly mutated proteins and test on the designed ones (“mutant vs designed split”), or train on the designed proteins and test on the mutated ones (“designed vs mutant split”), or train on proteins with seven mutations only and test on the rest (“seven vs rest split”), etc. The goal is to predict the fitness of the mutant, where the higher the better.

The thermostability dataset measures the thermal stability (higher is better) of 48,000 proteins from 13 species using a mass spectrometry approach; the resulting proteins are highly diverse, two splits were used, “mixed” (including all species) or “human” (a smaller subset including human data only). The GB1 dataset measures the fitness of mutants of the binding domain of protein G from a species of bacteria, with mutations at just four sites in the protein. The dataset includes almost complete coverage, i.e. , most of the possible four-site mutants were generated and tested. The “three vs rest” split uses mutants with three mutations only as the training set and all others as test, while “two vs rest” uses mutants with two mutations only as training.

The fluorescence dataset from TAPE involves mutations to the green fluorescent protein; the goal is to predict log fluorescence for test set proteins based on a training set. This dataset is the same as the one used in the literature, discussed under “Comparison with stochastic variational inference and exact Gaussian processes”, although the TAPE benchmark splits the data into training and test differently. Finally, the stability dataset from TAPE measures the stability of a diverse set of proteins in a protease susceptibility assay, where a larger value indicates better stability.

It was fit in triplicate using different random seeds, the resulting error bar measures variability across different random seeds. Since the focus is on Gaussian process regression, only benchmarks which involve regression tasks were used.

As shown in Table 1 , the FHTConv-1d kernel and FastConv-1d kernels trained on one-hot encoded data were compared with convolutional neural networks trained on one- hot encoded data. As shown in Table 2, learned embeddings - representations of the input sequence generated by the FAIR ESM-1 v pretrained language model for protein sequences - to represent the input sequences, were used, and these results were compared with those of fine-tuned language models reported in the literature. The best result from any fine-tuned large language model described in the literature for each benchmark is reported. Table 1. Spearman’s-r correlation coefficient of predictions vs ground truth values for convolutional neural networks and xGPR trained on one-hot encoded data a. In many cases, there are multiple “splits” for the same dataset; for the AAV dataset, for example, the “designed vs mutant” split uses designed sequences as a training set and mutants as test, while the “mutant vs designed” split is the reverse. ^b All convolution kernels use convolution filters of width 9 as a simple default. The best performer for each benchmark is highlighted in bold. For all the models, the predictions were repeated using three different seeds for the random number generator, then the mean and standard error of the mean are reported. For the scores shown here, higher is better. ^c. For the thermostability dataset, some of the sequences are quite long (thousands of amino acids). The Fast-Conv-ld kernel is much faster for these kinds of sequences since it only needs to perform convolutions once during hyperparameter tuning, only this kernel was considered here. ^d. All of the results for the one-hot encoded CNN are either from Dallago et al. (Proc. Neural Inf. Process. Syst. Track Dataset Benchmarks, 2021) or Rao et al. (ArXiv 190608230 Cs Q-Bio Stat 2019). ^e. For the GB1 dataset only, the input sequence is four amino acids in length, so for this dataset only an RBF kernel was used instead (equivalent to a convolution kernel of width 4). Table 2. Spearman’s-r correlation coefficient of predictions vs ground truth values for fine-tuned large language models and for xGPR trained on ESM-1 v embeddings ^a- In many cases, there are multiple “splits” for the same dataset; for the AAV dataset, for example, the “designed vs mutant” split uses designed sequences as a training set and mutants as test, while the “mutant vs designed” split is the reverse. ^b All convolution kernels use convolution filters of width 9 as a simple default. The best performer for each benchmark is highlighted in bold. For all the models, the predictions were repeated using three different seeds for the random number generator, then the mean and standard error of the mean are reported. For the scores shown here, higher is better. ^D. All of the results reported for the pretrained language model and one-hot encoded CNN are either from Dallago et al. or Rao et al. These authors report results for multiple pretrained language models; the best result they obtained in each case was reported. ^e. For the GB1 dataset only, the input sequence is four amino acids in length, so for this dataset only the FastConv-ld kernel was not used and no average over the sequence. For the learned embeddings, there are at least two ways that can be used as input to a Gaussian process, and both were evaluated. First, the representation can be averaged for each amino acid across all amino acids in the sequence, among others, then an RBF kernel can be used. Alternatively, the learned embeddings of size D can be arranged for each sequence of length M as an M x D array and FHTConv-1 d or FastConv- 1d kernels can be used to perform convolution across the sequence. It is easier to train using the FastConv-1d kernel since the convolutions can all be performed once at the start of training, thus eliminating the need to save the learned embeddings to disk (which can be quite large). FastConv-1 d was therefore used here.

Three interesting findings were noted. First, xGPR when trained on one-hot encoded input with a modest number of random features (< 10,000) outperforms a CNN on 8 out of 11 benchmarks and ties on one more. When trained on learned embeddings, the approximate GP outperforms the fine-tuned language models on 6 benchmarks and ties on one more. Note that the pretrained language models required up to 50 GPU-days training time to fit, while the FastConv-1d kernel takes less than 6 minutes per fit to fit every dataset shown here (including hyperparameter tuning), with a maximum memory footprint of 1.7 GB when using a minibatch size of 2000 datapoints. (The FHTConv-1 d kernel is slower but always takes less than 15 minutes including hyperparameter tuning). Using learned embeddings generally improves the performance of the GP, especially on diverse datasets with sequences from many different protein families (e.g., thermostability), but not always; on several benchmarks the GP performs better with one-

hot encoding.

Second, using a larger number of random features in all cases slightly improves the performance of the Gaussian process; it improves the accuracy of the kernel approximation. This feature provides the user with a simple way to improve performance if desired. Third, neither the FHT-Conv-1d kernel nor the Fast-Conv-1d kernel is consistently superior to the other; rather, each outperforms the other on some of the datasets.

Uncertainty calibration

A major motivation for using a GP is to quantify uncertainty. Various methods for uncertainty estimation for deep learning have been reported in the literature. xGPR provides better-calibrated estimates of uncertainty than these alternatives.

For an uncertainty estimate to be useful, uncertainty needs to be increased as moving far away from the training set - this is true for a GP with a stationary kernel by construction. Additionally, an uncertainty estimate that is well-calibrated is desired, in other words, the uncertainty reflects how often a predicted outcome actually happens. If a weather forecast claims there is a 10% chance of rain, for example, there should be rain on ten out of a hundred days when it makes this prediction. Likewise, for a regression model, roughly 10% of the test datapoints is expected to lie within a 10% confidence interval on its predictions (or a “credible interval” in Bayesian inference), roughly 20% of datapoints is expected to lie within a 20% confidence interval, 90% within a 90% confidence interval, and so on.

Calibration can be quantified as an area under the calibration error curve (ALICE). For 100 values in the range from 0 to 1 (0.01 , 0.02, etc.) it was determined what fraction of the test set predictions lie within the corresponding confidence interval (1 %, 2% etc.). The observed fraction was subtracted within the interval from the expected fraction and the absolute value was taken - for example, if the confidence interval is 10% and the observed fraction within that interval is 90%, the difference is |0.9 - 1|. Finally, the area under this curve was calculated; the larger the area, the more poorly calibrated the uncertainty figures provided by the model. A perfectly calibrated model achieves a score of 0. Five different methods were described for uncertainty estimation (dropout uncertainty, last-layer stochastic variational inference, model ensembling, deep evidential regression and mean-variance estimation) and used in combination with deep learning models on five of the benchmarks described in this Example. Both one-hot and ESM embeddings were used as input to the models and AUCE was calculated. As shown in Table 3, the AUCE for their uncertainty estimation methods was compared with the AUCE for xGPR for the same datasets using the same encodings. To give the deep learning models the benefit of the doubt, the best AUCE reported for each uncertainty estimation technique was presented for each dataset, and the xGPR model from Tables 1 and 2 that displayed the best performance for that dataset (as assessed by correlation between predicted and measured fitness) was used.

Table 3. Area under the calibration error curve (AUCE) for the best xGPR and for four deep learning uncertainty estimation methods using the same input representation

Note that in all cases, to give the deep learning models the benefit of the doubt, the best performance for that uncertainty estimation technique was used, regardless of whether using one-hot or ESM embeddings as input.³. In many cases, there are multiple “splits” for the same dataset; for the AAV dataset, for example, the “designed vs mutant” split uses designed sequences as a training set and mutants as test, while the “mutant vs designed” split is the reverse. Where uncertainty calibration is concerned, on all five benchmarks, xGPR performs substantially better, with an AUCE from threefold to tenfold better than the best deep learning uncertainty estimation method. The AUCE for xGPR is < 0.15 in every case, while the deep learning methods often (and for some methods always) achieve AUCE scores of 0.3 or above.

Using uncertainty for protein engineering

Uncertainty can be used to determine whether a prediction generated by a model should be trusted; if the uncertainty is high compared to uncertainty for the training set, it may consider a prediction less likely to be reliable. Alternatively, the uncertainty estimates generated by xGPR can be used as part of a Bayesian optimization strategy. In this framework, we first fit the model to the available data. A new batch of sequences were chosen to experimentally evaluate using an acquisition function that takes into account both the model predictions and the uncertainty, thereby achieving a balance between exploration and exploitation. The upper confidence bound (UCB), for example, adds the upper half of the confidence interval to the model prediction for each test point. After experimentally evaluating the selected points, the model was re-fitted, and this process was repeated until a desired protein function is achieved. For an overview of other common acquisition functions, see Frazier 2018.

In this study, the GB1 dataset discussed under “Evaluating convolution kernels for sequence data” was used in which mutations were introduced to four sites in the binding region of protein G. This dataset is convenient for an experiment of this kind because close to 150,000 of the 160,000 possible mutants have been evaluated. The fitness of all mutants in the dataset was normalized so that it lies between 0 (worst possible) and 1 (best possible). Note that > 99.95% of the sequences have a fitness < 0.6; only one has a fitness of 1 . 384 variants were randomly selected and an xGPR model was trained with an RBF kernel, then the UCB acquisition function score was generated for all the remaining mutants and used to select a batch of 96 to “experimentally” evaluate. These 96 was then added to the training set, the model was refit and the process was repeated for up to five iterations. At each iteration, the best fitness achieved was tracking. Finally, this “experiment7”process” was repeated 50 times with different random seeds.

By the fifth iteration, all repeats achieved a best fitness of 0.6 or above. 40 of the 50 repeats had found one of the four best sequences, and 30 of the repeats had achieved a score of 1 - i.e., had found the best sequence in the dataset - despite experimentally testing < 1000 sequences. Half of the repeats were able to find the best sequence (fitness 1 ) in four iterations or less. These results suggest that if combined with a high-throughput experimental system, this approach using xGPR might be able to “discover” proteins with desirable properties with no human intervention required.

Efficient, numbering-invariant kernels for small molecules and other graphs

Two groups of linear-scaling kernels were introduced for sequences under “A first new set of efficient kernels for sequences that cannot be aligned” above. These kernels can easily be extended to graphs, if each node was encoded in the graph as a feature vector containing some information about that node, its edges and adjacent neighbors. For the Fast-Conv-1d kernel, the random feature convolution and max-pooling operation across all feature vectors associated with the graph can be applied. For the FHT-Conv- 1 d kernel, nodes and their associated features were merely used as inputs rather than k- mers, and a width-1 convolution across the nodes of the graph was performed: this is referred to as a Graph-Conv-1d kernel.

Several variants of the Graph-Conv-1d kernel were implemented in xGPR. The simplest variant, Graph-RBF, was focused on here, which approximates a pairwise comparison across the nodes of two graphs using an RBF kernel. To evaluate the performance of Graph-RBF, the prediction of energetic and thermodynamic properties of small molecules were considered. The QM9 dataset which contains 134,000 small molecules whose geometry has been optimized using density functional theory (DFT) was benchmarked on. The kernels described here seem especially likely to be a good fit for energetic/thermodynamic properties based on prior knowledge, since these properties can in principle be decomposed as a sum of per-atom contributions.

Two different modeling approaches were used. In approach (1 ), each atom / was one-hot encoded, then concatenated this to a one-hot encoding of each other atom j in the molecule divided by the distance of atom j to atom / to the 6th power; the one-hot encodings are sorted by distance. Because of the 6th power weighting, out past four or five angstroms the features associated with more distant atoms become negligibly small and can be disregarded. This representation was then used as input to a Graph-RBF kernel. This simple representation is easy to construct but is not very informative (it discards all information about bond angles) and thus unlikely to achieve good performance. Alternatively, in approach (2), each atom was represented using the SOAP features. In this scheme, the atomic density surrounding each atom is expanded in an orthonormal basis set, and the power spectrum of these coefficients is used as a representation. The distance weighting scheme was used with only small modifications to their selected hyperparameters. For more details on the hyperparameters used and on hyperparameter tuning, see Supporting Information section S10.

Remarkably, approach (1 ), despite using uninformative features, achieves chemical accuracy of about 1 kcal/mol. (2) is much better; results are shown in Table 4. For internal energy, enthalpy and Gibbs free energy, xGPR can achieve an MAE < 0.2 kcal/mol, and for one of the six properties we achieve a new state of the art. Note that 1 kcal/mol is typically considered to be “chemical accuracy”, i.e., the accuracy required to make correct predictions for the behavior of chemical systems, so that both xGPR and several of the other models are well above the level of accuracy that is actually needed. Indeed, the density functional theory (DFT) calculations whose outcome we are trying to predict are often only accurate to 2-3 kcal/mol. Consequently, a difference of 0.01 - 0.05 kcal/mol between models for predicting DFT-calculated energies is likely equivalent for practical purposes.

Table 4. Comparison with state-of-the-art deep learning models for prediction of energetic and thermodynamic properties of small molecules, QM9 dataset

Retrieving similar datapoints from the training set

The feature importance can be measured using the SHAP or LIME techniques for a Gaussian process just as for any other ML model. Unlike for other ML models, for a small dataset the kernel matrix can additionally be constructed and used to determine exactly how much each datapoint in the training dataset contributed to a given prediction. This kind of decomposition is not possible with a deep learning model.

Constructing the kernel matrix explicitly, however, is too expensive if the dataset is larger than 5,000 datapoints. Also, for large datasets assigning a precise contribution to each point in the training set may be more information than is required or helpful. Often, practitioners may merely want to know what are the mosts/m/7a/'datapoints in the training set - the ones which contributed the most to a prediction. This can be determined as follows. Let z, be the random features representation of a test datapoint and let z_t be the random features representation of training datapoint; then the kernel function for the corresponding test and training datapoints x* and x_t is approximately zjz_z (see Supporting Information section S1 )

In practice, it may be computationally burdensome to take the dot product of a test datapoint representation with all the training data every time to retrieve the most similar training datapoints. This problem can be resolved by simply clustering the training data. The same kernel used to fit a Gaussian process regression model can be used to cluster it through random features-approximated kernel k-means clustering. Note that in the same way that random features approximate a Gaussian process as linear regression, clustering the random features-represented training set with k-means approximates applying kernel k-means to the original dataset. Once the dataset has been clustered, to retrieve the most similar training points, there no longer needs to measure the similarity of a test point against every training data point, only the similarity against the datapoints in the nearest cluster needs to be measured.

Tools to perform these functions are included in the xGPR library. This procedure is illustrated for the QM9 dataset, with additional details in Supporting Information section S11 . The QM9 dataset was first fitted with energy at 298 K as the target using the Graph- RBF kernel with atoms encoded using the one-hot encoding approach (1 ), then k-means cluster the training set. The resulting elbow plot used to select the number of clusters is described in Supporting Information and FIGs. 9A-9B. The elbow plot suggests the training dataset can be broken into five clusters. Five molecules from the test set were then randomly selected. For each of them, the procedure described in the above paragraph was used to find the fifty most similar molecules from the training set, where similarity is measured by the kernel function. The structures of the query test molecules and the most similar training molecules (as quantified by the kernel using one-hot encoded input) are described in the Supporting Information section S12.

Supporting Information

Supporting information provided below contains: (S1 ): more detailed background on the random features approximation, (S2): more detailed background on orthogonal random features, (S3): performance metrics for our implementation of the fast Hadamard transform, (S4): the preconditioning algorithm and some benchmarks for modifications to this algorithm, (S5): the algorithm for approximating log determinants for random features- approximated Gaussian processes, (S6): a discussion of hyperparameter tuning strategies, (S7): details on benchmark datasets used in this Example, (S8): benchmarking of random features approximations vs stochastic variational inference, (S9): a proof of Lemma 1.1 , (S10): details on the QM9 benchmark, (S1 1 ): clustering and principal components using random features, and (S12): retrieval of similar molecules from a training set using random features.

S1. Random features

Only positive definite kernel functions and initially only stationary kernels, i.e. kernels for which Kfjx^Xj') = f(x_t - x₇) were considered, and a transformation </>(x) such that the kernel (x^xf) ~ 0(x₍)0(x₇) was sought. From Bochner’s theorem, any of the positive definite kernels considered is the Fourier transform of a corresponding positive measure, i.e.: this expectation can be approximated using Monte Carlo sampling, i.e.:

Since Monte Carlo approximation is an unbiased estimator of the mean, the estimate converges on the true value as the number of random features is increased. The variance of this estimator is proportional to (1 ! M).

Although this procedure was originally derived in the context of stationary kernels, it can be applied to a similar Monte Carlo sampling approach to certain nonstationary kernels, most importantly a class of kernels dubbed “neural network kernels” and also called “pointwise nonlinear Gaussian kernels” in more recent literature. Consider a neural network with one hidden layer and one output layer: where x is an input vector, W is a learned weight matrix of dimensions D x N_H where D is the dimensionality of the input, N_H is the number of units in the hidden layer, b is an offset or bias vector of length N_H , h is some activation function and V[ is a learned weight. Assume that all weights in W are iid drawn from 7V(O, , the weights v_t are iid drawn from and that the bias terms are iid drawn from N(0, <7_b ²). If each x is augmented with an additional element that is 1 for all data points, this can be rewritten more simply as:

Where each row of W is drawn iid from N(0,Z_w), when E_w is a diagonal covariance matrix whose first diagonal element is a_b ² and all remaining diagonal elements are a_w ².

Let S denote the collection of all weights; it follows that

Once again, the expectation term can be approximated on the right-hand side using Monte Carlo sampling and thereby a formulation very similar to the one for stationary kernels is obtained. This kernel is only positive definite for certain choices of activation function h. One such choice is the ReLU activation function which yields the order 1 arccosine kernel, although the error function has also been suggested as a useful activation function. Interestingly, the order 1 arccosine kernel as formulated in the literature omits the bias terms. During the implementation of the library, it was found that the arccosine kernel performs very poorly without the bias term and requires the inclusion of the bias term to achieve acceptable performance.

For both stationary and neural network kernels, then, the random Fourier features approach can be implemented as follows. If X is an N x D matrix for N data points and D input features, a D x M matrix is populated with M samples from is the Fourier transform of the kernel for stationary kernels and the normal distribution for neural network kernels. This sample matrix is multiplied by X, and for neural network kernels is added an additional bias vector to the result. The appropriate activation function (e.g., sine and cosine for stationary kernels) is then applied and multiplied by a constant to yield an N x M matrix Z (or, for stationary kernels, N x 2M). The approximate kernel matrix is then given by ZZ^T.

The distribution across predicted values /, at new datapoint for a mean-zero Gaussian process given training dataset X with corresponding ground truth values y is:

Since the kernel is approximated as ZZ^T, the matrix inversion lemma and the Woodbury matrix identity can be used to rewrite these. The mean for example becomes: where Z is the training set transformed as described above, I is the identity matrix, A is a hyperparameter which ensures the kernel matrix is positive definite and z is the new datapoint transformed as described above. If consider

[Z^TZ + A²I]~¹Z^Ty to be the weights, the predictive mean corresponds to linear regression in the feature space corresponding to the random feature map. It is easy to see the correspondence between a neural network with a single hidden layer plus an output layer and the above kernels approximated using this formulation, although there are several key differences. The approximate GP “learns” by changing the distribution of the hidden layer weights, not individual values, and can be fitted by maximizing the marginal log-likelihood rather than the likelihood; also, the marginal likelihood and the posterior predictive distribution are available in closed form. It is also possible to construct kernels which are analogous to neural networks with several hidden layers. In the main text, the Fast-Conv-1 d kernel was introduced which corresponds to a three-layer neural network with a convolution layer, a fully connected layer, and an output layer, all populated with random weights (See FIG. 3)-

Note that since the error in the approximation decays exponentially with an increasing number of random features, it does not generally need a large number to obtain a good approximation; indeed, increasing random features past 16,000 or so yields sharply diminishing returns, and 8,192 is sufficient for many tasks. For predicting the mean /i*, since this is the most important quantity, practitioners may want to use a large number of random features to “buy” small improvements in accuracy, but for quantifying variance 512 - 1 ,024 is generally sufficient. Consequently, xGPR allows users to use a smaller number of random features to quantify variance than for predicting the mean (to reduce the computational expense of predictions for variance).

S2. Orthogonal random features

Consider the case where the number of random features D is the same as the dimensionality of the input, d, and where the input vector x is multiplied by a matrix M of random features drawn iid from a normal distribution.

If the rows of M are likely to be orthogonal, the QR decomposition of M can be first taken, then the resulting Q is multiplied by a diagonal matrix S whose diagonal elements are drawn from a /-distribution with d degrees of freedom. The input vector x is now transformed by computing the product kSQx. If the number of random features is greater than d, this procedure can be repeated m times where m is the result of ceiling division of D by d, and then any features in excess of D are discarded. If D is less than d, d random features can be generated and any in excess can be discarded.

For certain kernels, in particular the squared exponential kernel, it has been shown that fewer orthogonal features are required to accurately approximate the kernel, so a more accurate approximation for the same number of random features is achieved by using orthogonal random features. In addition to this benefit, the orthogonal feature approach can be further modified to realize a large reduction in computational cost. It was shown that nearly equivalent results can be achieved by replacing the matrix Q with the following product:

HDI HD₂HD₃ where H is the normalized Hadamard matrix, and each D is a separate diagonal matrix whose entries are drawn from the Rademacher distribution. The memory footprint of a fitted model is greatly reduced, because only the diagonal elements of the D-matrices need to be stored. More importantly, similar to the substitution of the fast Fourier transform for multiplication with a matrix in the discrete Fourier transform, the fast Hadamard transform is substituted for the Hadamard matrix multiplication, thereby reducing the cost of this procedure to O (D log D) - a dramatic benefit for large M. The general procedure for transforming the input is shown for the RBF kernel in Algorithm S1 (other kernels similar).

Note that it was suggested replacing the S-matrix used here with multiplication by a constant, Vd; this was failed on 1 d data, where the resulting kernel is unable to provide an acceptable fit for simple toy problems, and it provides no benefit on higher dimensional data. The version from above does not suffer from this limitation.

Algorithm S1 Orthogonal random features transform for RBF kernel Kernel object initialization

Input: Expected dimensionality c/ of input; expected number of sampled frequencies M

I Q _ 2^ce^"^lS^F(^°S'2(^mclx(^ci,2)))

2. If D < M: a. P = ceiling(M I D)

3. Else: a. P = 1

4. Create, store diagonal matrix S with size M x Mwith elements drawn from the /- distribution with degrees of freedom D

5. For i in 1 ...P: a. Create, store diagonal matrices D_ai,D_bi,D_ci with size D with elements drawn from the Rademacher distribution

Kernel applied to transform input data

Input: N x d matrix x_h which is a chunk of the input data; kernel hyperparameters 0, a

Output: kernel-transformed matrix Z_t

1 . For i in 1 ...P: a. Create empty N x D array Y_t b. Copy xt into the first d columns of Y_t; set remaining elements to 0 c. For k in [a, b, c]: i. Yt = H(YiD_ki) , where H is a normalized in-place fast Hadamard transform

2. Concatenate Y₁.. . Y_P to form N x D*P array Y

3. \f D*P > M: a. Discard columns M...D*P of Y

4. Y = aYS

5. Y_r = cos(Y), Y₂ = sin(Y) where sin and cos are applied elementwise

6. Concatenate Y_t and Y₂ to form N x 2M matrix Z

8. Return Z

S3. Performance of the fast Hadamard transform

A Python-wrapped C / Cuda implementation of the fast Hadamard transform for CPU and GPU is implemented. The speed of the implementation with matrix multiplication against a matrix of random features as implemented in the Numpy library (v 1 .22.3) and Cupy library for GPU (v 10.4.0), was compared (see Tables 6 and 7 for the results). Briefly, the implementation is faster than matrix multiplication for modest numbers of random features, and the performance advantage of the implementation increases with larger numbers of random features as expected, enabling to quickly scale to large numbers of random features. To compare with other fast transforms, the discrete cosine transform (DCT) on CPU as implemented in the Scipy library (v 1 .8.0) was also compared (Table 1 for the results).

The results in milliseconds are shown in Table 5 and Table 6. In Table 5, the time needed to generate random features if the input matrix has as many columns as there are features was evaluated - this is a somewhat unusual situation but provides a dramatic demonstration of the speed increase for the fast Hadamard transform-based procedure. In Table 6, the more realistic situation was considered where there are a fixed number of 1024 input features, and an increasing number of random features must be generated.

An input matrix with 2000 rows was used since this chunk size is common when work on real datasets. For the GPU implementation, 32-bit (single precision) floats were used due to the substantially greater performance for single precision matrix multiplication on most GPUs, while for the CPU tests 64-bit (double precision) were used. In general, however, unless the user selects otherwise, xGPR performs random feature generation using single precision for both CPU and GPU.

Table 5: Speed of the SORF operation implemented in this study (highlighted in the table) compared with vanilla random features and with a generic Scipy fast transform The results shown here are rounded to reflect precision, include both CPU time and GPU time (where applicable) and are calculated as the average time per operation for 100 repeats of the specified operation. All studies were performed on an Intel I5-7500 CPU and a GeForce GTX1070 GPU. The methods used in this study are in bold.

Table 6: Speed of the SORF operation implemented in this study (highlighted in the table) compared with vanilla random features and with a generic Scipy fast transform

The results shown here are rounded to reflect precision, include both CPU time and GPU time (where applicable) and are calculated as the average time per operation for 100 repeats of the specified operation. All studies were performed on an Intel i5-7500 CPU and a GeForce GTX1070 GPU. The methods used in this study are in bold.

S4. Preconditioning for conjugate gradients and stochastic gradient descent

Consider the well-known conjugate gradients algorithm, which solves the system Aw = b for w. It can be shown that the residual error on iteration n has as its upper bound: where K is the condition number of A, and | |e₀| | is the error of the initial guess. If K is large, the algorithm will therefore converge very slowly. A preconditioner whose inverse M^-1 is a close approximation to the inverse of A will reduce the condition number of A and accelerate convergence. It is critical to have a preconditioner that can be constructed without explicitly forming the matrix Z^TZ , which requires expensive matrix multiplications and for large numbers of random features is too large to store in memory.

Kernel matrices for Gaussian processes have frequently been approximated using the Nystrom method, which provides the following low-rank approximation: where m is the number of rows in A, and q is a randomly selected subset of the columns of A. The Nystrom approximation is cheap to construct and has worked well for many problems in practice. Although sampling from a uniform distribution has traditionally been the most popular approach, the performance of this method is sensitive to the sampling method used to select columns and rows in A.

The randomized Nystrom approximation was found to provide better performance. In this scheme, for the system:

A low-rank approximation of Z^TZ = A is provided by: where for /2, an iid matrix drawn from a zero mean, unit-variance normal distribution was used; is of size M x L, where M is the number of random features. This approach provides a rank-L approximation to A, and a numerically stable approach for generating it (which was implemented for the library) and appears in Algorithm S2. Importantly, this can be constructed as shown without ever explicitly forming the matrix Z^TZ and using a single loop over the dataset, loading one chunk of the data into memory at a time.

Algorithm S2 Constructing a randomized Nystrom preconditioner

Input: Dataset X stored in chunks x₁...x_n on disk, desired preconditioner size L, GP kernel, number of random features M

Output: Factored preconditioner in the form UAU^T, where U is orthonormal and A is diagonal

6. Generate M x L matrix fl by sampling from N(0,l)

7. fl = qr(fl) QR decomposition

8. Generate M x L matrix Q and initialize to all zeros

9. For i in 1 ...n, load x_t a. Use Algorithm S1 to generate z_t For numerical stability

A novel modification was introduced here. The M x L matrix is replaced with the subsampled randomized Hadamard transform, which is defined as follows:

Where D is a diagonal matrix with entries drawn from the Rademacher distribution, H is the normalized Hadamard matrix, and S is a subset of L rows drawn with equal probability from the identity matrix. In practice, the fast Hadamard transform rather than the Hadamard matrix is used. For step 9b, it substitutes:

In step 12, for the product s SRHT(Q_V) is substituted. Finally, to ensure numerical stability, the shift is replaced -then- procedure above (steps 1 1 and 15) is subtracted, wherein rather than taking the Cholesky decomposition of SRHT(Q_V) in step 12, the self-adjoint square root of this matrix is found using SVD, then used to solve for B in the following step.

Using the SRHT to construct low-rank matrix approximations has been suggested before in the literature, but not to construct a preconditioner for CG I stochastic gradient descent as suggested here. It significantly reduces the cost of preconditioner construction, since a matrix multiplication with O(NML) scaling was replaced with a fast transform operation; this is especially advantageous if working on CPU. Empirically, it is found the SRHT-based preconditioner construction routine to be substantially faster even for datasets of only 40,000 datapoints, and the advantage increases with increasing dataset size and number of random features.

In every case encountered, it is found empirically that a preconditioner constructed using SRHT as suggested here performed either as well as or better than a preconditioner constructed using the unmodified algorithm S2. As shown in FIG. 4, the number of iterations required for conjugate gradients to converge was compared to a pre-specified threshold using a preconditioner constructed either with SRHT or with the unmodified algorithm S2, using a variety of different values for L In every case, the SRHT-built preconditioner achieves the same or greater acceleration, and indeed this has been consistently true throughout the studies. For the time being, both the modified and unmodified preconditioner construction algorithms are provided as options the user can select in xGPR and default to the SRHT-based algorithm in light of its superior speed and empirical performance.

Additionally, the preconditioner can be further modified to improve its performance. In step 9 of Algorithm S2, the product Z^T(Z/2) is formed where Z is the matrix of random features for the training data of shape N x M for M random features and N datapoints. The matrix Q can be retrieved from the QR decomposition of this product and used to form a second product Z^r(ZQ) using a second pass over the dataset. This product can then be used in Algorithm S2 in place of the matrix Q from step 9. (this process can even be repeated a third time, although this offers a negligible improvement in preconditioner performance). The resulting preconditioner reduces the number of iterations required for conjugate gradients to converge in all the studies, in most cases by 20-25%, as illustrated FIG. 5. It is, however, more expensive to construct since it requires two passes over the dataset. Whether this is beneficial depends on the size of the dataset and on available hardware. The user is allowed to select “srht_2” as an option for preconditioner construction in xGPR and guidance is provided in the documentation on when to choose this option.

Irrespective of which option is used to construct it, the inverse of the preconditioning matrix takes the form:

Where 0_L is the smallest eigenvalue of the preconditioner and A is a hyperparameter of the GP kernel. The matrix M^-1 is not explicitly formed, but since all the needed values above have been generated, the matrix-vector product M^-1 v can be supplied on demand, which is all that is required for preconditioned conjugate gradients or preconditioned stochastic gradient descent.

The larger L, the better the approximation; if ft > ft > ft... are the eigenvalues of A, the approximation is also better the more rapidly the eigenvalues decay, and the smaller the eigenvalues for ft>_L. Let UAU^T be the eigendecomposition of A_nystrom of rank L, let the preconditioner P be P = and let ft be the smallest eigenvalue of A_nystrom. It can be shown that if K₂ is the condition number of the preconditioned system, the probability that: is > 1 - 6, where 6 is some failure tolerance.

How to choose L? The value for a given L can be estimated by constructing a preconditioner, and the value of L can then be doubled until suggests the number of iterations required to fit the model is likely to be acceptable. As illustrated in FIG. 1 A, the log of the number of iterations often exhibits a roughly linear relationship with for a given desired tolerance. Consequently, the number of iterations required to fit for a given can estimated and if the resulting number of iterations is considered too high, L can be adjusted, and the preconditioner can be reconstructed. “.

The relationship between the value of L and the number of iterations required to fit using CG is illustrated in FIG. 1 B below.

The preconditioner enables to quickly fit the model but does not address the issue of hyperparameter tuning. Calculating the marginal likelihood of a Gaussian process requires calculating the log determinant of an M x M matrix for M random features, which is very expensive if M is large. Therefore, an approach for estimating marginal likelihood for large numbers of random features was developed.

S5. Log determinant approximation and conjugate gradients

Hyperparameters can be tuned for a GP using performance on a validation set (as for any other model), and this approach is implemented as an option in xGPR, but it is also possible to tune hyperparameters by maximizing the marginal likelihood. This procedure is more robust to overfitting and also means no need for a validation set, so use this approach where possible. Nonetheless, tuning using marginal likelihood does present one particular challenge. Recall that for a Gaussian process the following is found for the negative marginal log likelihood:

Where K_y = K(X,X) + A²1. In the random features scheme, approximation Z^TZ is substituted for K(X,X) and rearranged to yield:

-log p(y | X)

As shown above under S4, it can efficiently be computed (Z^TZ + A²l)~¹Z^Ty without ever forming the matrix (Z^TZ + A²/) explicitly using preconditioned conjugate gradients. The log\Z^TZ + A²11 , however, is more challenging to compute if the number of random features is large. Under section S6 highly efficient ways were developed to tune hyperparameters with a small number of random features, e.g., < 4000, which frequently do in practice. First, however, a method to use preconditioned conjugate gradients to approximate log\Z^TZ + >l²/| without ever forming the Z^TZ + A²l matrix explicitly was introduced. The same tool is used for fitting to tune hyperparameters in an efficient way even if the number of random features is large.

A variety of approaches have been outlined in the literature for approximating the log determinant of a large matrix A, including stochastic Lanczos quadrature, stochastic trace estimation combined with Chebyshev polynomials, and subspace iteration. The subspace iteration method is only appropriate if a small subset of the eigenvalues of A are dominant, i.e., there is a large eigengap with eigenvalue A_t » A_i+1,A_i+2...A_N, which is not always true during tuning.

Stochastic Lanczos quadrature (SLQ) is especially attractive in this context, since for large numbers of random features preconditioned conjugate gradients were used to fit the model, and the tridiagonal matrices required for stochastic Lanczos quadrature during the course of the conjugate gradients optimization at negligible additional cost can be automatically generated. It was already demonstrated that the use of this approach to accurately estimate log determinants for exact (non-approximated) Gaussian processes using an incomplete Cholesky decomposition to build the preconditioner.

In Algorithm S5 below, it shows how to use this approach in combination with the random features approximation and the randomized Nystrom approximation to accurately estimate the log determinant term, log\Z^TZ + A²/ 1, while simultaneously calculating the performance penalty — y^TZ(Z^TZ + A²/)^-1Z^Ty (the remaining terms in the marginal likelihood are either constants or not required to retrieve the data and hence are easily calculated). The resulting algorithm exhibits linear scaling with dataset size, subquadratic scaling with number of random features, and never requires loading more than one chunk of data into memory at a time.

As a preliminary, the log determinant was first rewritten and Hutchison’s trace estimator²¹ was used:

Where p_t are probe vectors drawn (typically) from either a Rademacher or a Gaussian distribution.

Next, when use preconditioned conjugate gradients with preconditioner P, it is in fact solving the preconditioned system log\ P~^1/2 Z^TZ + A²1 p-^² | . The same approach adopted for exact (un-approximated) GPs can be adopted and this is rewritten as: The preconditioner therefore needs to be applied to the probe vectors. This is easily done if the probe vectors are drawn from a Gaussian; they can be merely drawn from a normal distribution whose covariance matrix is the preconditioner. These preconditioned probe vectors are now referred to as Q . The term is evaluated as:

This is done using the quadrature-based procedure but using the preconditioned conjugate gradients algorithm; by storing the alpha and beta coefficients produced on each iteration of conjugate gradients, the tridiagonal matrix constructed via the Lanczos algorithm was generated. In this case, merely need to run conjugate gradients with multiple target vectors - both the Z^Ty vector so that the performance penalty (Z^TZ + A²l}~^rZ^Ty can be calculated; and the preconditioned probe vectors so that n_v tridiagonal matrices T₁,T₂... T_nv can de generated. A simplified version of the full procedure is outlined in Algorithm S5 below.

Algorithm S5 Preconditioned conjugate gradients for simultaneous estimation of (Z^TZ + A²/)-¹Z^ry and log\Z^TZ + A²/|

Input: Dataset X with N datapoints stored in chunks x_r.. . x_n and y_±. .. y_n on disk, initialized kernel with fixed hyperparameters that generates length-/W random feature representation k(Xj) for input datapoint x_t, preconditioner P constructed using Algorithm S2 or appropriate variant, number of probe vectors n_v, convergence threshold to/

1 . Create empty array Z^Ty of shape M

3. For i in 1 ,2...n_v: a. Draw probe vector c_t from N(0, P)

4. Create empty array of shape (M,n_v+1 )

5. Populate the first n_v columns of rj with c_lt c₂, c₃... c_nv. Copy Z^Ty into the remaining column. Let L₁ = [||r₇₁| |J |r₇₂ ||J |r_j3 | |_; .. J |r_7.np+1 ||] where | |r_7fc| | is the Euclidean norm of column k in r₇, such that L_} is a row vector of length n_v + 1 Create initial weight matrix w₇ of shape (M, n_v+1 ); set all entries to zero. Create empty array B of shape (M, n_v+1 ). Let 1_M be a length-M vector of ones. Let err be a length n_r+1 vector of ones. pj = Zj = P~^ rj,' j = 1 While max(err) > tol and j < maxiter: a. Set all entries of B to zero b. For i in range 1 ,2.../?: Let a_ik be the element of a_k corresponding to iteration k and probe vector / (corresponding to column / + 1 of p₇). Let /3_ik be the element of p_k corresponding to iteration R and probe vector / (corresponding to column / + 1 of p₇).\ F = 0 For i in 1 ,2...n_v: a. Create empty matrix T of shape Q,j). Set all entries to zero.

15. Return estimated log determinant r , vector Z^Ty and column 1 of w₇ ; this last corresponds to (Z^TZ + A²/)^-1Z^ry and can be used to calculate the performance term of the negative marginal log likelihood

To evaluate the accuracy of this procedure, 8 hyperparameter combinations were randomly selected, each for 7 different datasets ranging in size from a few hundred datapoints to half a million. For each, the negative marginal log likelihood with 4,096 random features were calculated either using Algorithm S5 or using matrix decompositions, using a preconditioner rank of 512, a convergence threshold of 1 e-5 and the “srht_2” preconditioner construction algorithm in all cases. Either 40 probe vectors or 25 probe vectors were used; the distribution of the absolute percent error is plotted as shown in FIG. 7, showing that 40 probe vectors do not provide a substantial improvement on 25 probe vectors. 25 probe vectors are sufficient in this study to achieve high accuracy.

To further assess the accuracy of this procedure, approximate marginal likelihood in combination with Bayesian optimization was used to tune hyperparameters for two of the tutorials/example experiments in the user documentation. The first such experiment involves the UCI proteins dataset, the second involves the QM9 dataset with one-hot encoded inputs. In both cases, the tuning procedure improves on the best validation set performance achieved by tuning with a smaller number of random features with a matrix decomposition-based approach. The same approach was further used to tune hyperparameters for the tabular data and the small molecule experiments on the QM9 dataset described above to achieve highly competitive performance. These results are consistent with the data shown in FIG. 7 and suggest that given an adequate preconditioner rank, a sufficiently small convergence threshold, and a sufficient number of probe vectors, for the datasets considered here, Algorithm S5 is able to estimate the negative marginal log likelihood with high accuracy.

S6. Improved strategies for hyperparameter tuning

S5 above shows how to use preconditioned conjugate gradients to calculate the marginal likelihood in an efficient way if the number of random features is large. Often, however for noisy data, only need to use 1000 - 3000 random features to tune hyperparameters and get a result that provides sufficient performance. If do indeed decide to use a smaller number of random features for tuning than for fitting, a fast and highly efficient strategy can be introduced. Recall that the marginal likelihood in the random features approximation is given by:

-log p(y | X)

Z^TZ is Hermitian so that as long as it is full rank, an eigendecomposition PAP^T can be generated, such that P is orthonormal, i.e. PP^T - I. The above expression can then be rearranged:

-log p(y | X)

For three and four hyperparameter kernels, this simplification can be used in a variety of ways, for example Algorithm S3 and Algorithm S4.

Algorithm S3 Efficient grid search Input: Dataset X with N datapoints stored in chunks x₁...x_n and y₁...y_n on disk, initialized kernel with two, three or four hyperparameters (A,/?, if applicable a₂) for M random features, selected values for hyperparameters (if applicable) a₁ and a₂

1 . Initialize Mx M array l/IZ and M x 1 array S, set to zero

2. a = 0

3. For j in 1 ...N: a. Load x₇, y₇and use Algorithm S1 to form Qj = Z₇ / ft using 8_t b. W+= Q^TQ ; S+= Q^Ty_t

C. a+= yj^Tyj

4. U,A = eigh(W + bl) add small constant delta for numerical stability

5. A—= 8

6. S = U^TS

7. Define loss_fn(/?, A): a. Return ^S (A + X²I)~¹S + ^y^Ty + (n - m)Zo^(A) +

1 ₉ n

-log\A + A²7| + -log(2n)

Ret

Algorithm S4 Bayesian hyperparameter optimization (Thompson sampling)

Input: Dataset X with N datapoints stored in chunks x_x...x_n and y_x...y_3l on disk, initialized kernel with three or four hyperparameters (A,/?,^, if applicable a₂) for M random features, maxiter maximum iterations, n candidates per iteration, convergence criteria tol

1. Evaluate score, A_best,p_best for 5 - 10 randomly selected kernel-specific hyperparameters <H₁ . . . <D_{1 0} using Algorithm S3

2. Let 12 be the set of all kernel-specific hyperparameter combinations evaluated thus far; let (jo_best be the best evaluated so far

3. For niter in 10...maxiter: a. Fit exact Gaussian process Swith Matern kernel, v = 5/2 b. Draw n kernel-specific hyperparameter combinations a>_x... a>₃₃ from Uniform(a_max > &min) c. Draw m samples from S at each of <z>_x.. . a)_n; select oj_new = ar grain (m) d- If | |ro_Jiew O^best 112 < tol. i. break e. Evaluate the loss function from Algorithm S3 at m_new

4. Return

Note that Thompson sampling was used in S4, but other acquisition functions (e.g., expected improvement) can easily be substituted as well. A simple grid search across the kernel-specific hyperparameter can also be used, which is also implemented in xGPR. The above algorithm can also be modified by calculating the eigendecomposition of Z^TZ by taking the singular value decomposition of Z, which is the most stable approach.

The polynomial kernel has only two hyperparameters (the lambda or noise and beta or amplitude hyperparameters shared with all other kernels). Using the strategy above, these can be tuned using a single pass over the data, so that hyperparameter tuning with the polynomial kernel is exceptionally fast.

S7. Datasets used for experiments described in the main text

All datasets used in this Example are publicly available. Details of each dataset are described in the following table.

Table 7. Details of datasets used for the studies described herein

S8. Comparison of xGPR with stochastic variational inference and exact Gaussian processes

The stochastic variational approximation to Gaussian process models and variants has appeared frequently in the literature. From a theoretical standpoint it has many attractive properties, but how does it perform in practice? As shown in Table 8, the performance of SVGP as implemented in the GPyTorch library was compared with the random feature approximation as implemented in xGPR; for sufficiently small datasets, and also with an exact Gaussian process as implemented in GPyTorch (when an exact Gaussian process was trained for larger datasets, an out of memory error occurred for the UCI proteins dataset with about 8,000 datapoints, so that for even larger datasets it is clearly impractical). Stochastic variational inference was trained in GPyTorch for 40 epochs in all cases, and the RBF kernel was used in all cases, since it was not aware that an SVGP-compatible convolution kernel for sequence data is sufficiently efficient that it can be used on a > 100,000 sequence dataset in a short period of time, nor are any such kernels implemented in GPyTorch. Training for 20 epochs yielded worse performance in every case.

For xGPR, hyperparameters were tuned using 3,000 random features using the “fast” strategy in Algorithm S3 described in S6 above, and then fit using either 16,384 or 32,768 random features for tabular data and 8,192 or 16,384 random features for sequence data (a smaller number of random features was used for sequence data to

correspond to the number of features used for comparison with CNN deep learning models on sequence data). Note that several of the datasets used as benchmarks here are protein sequence datasets also used as benchmarks for comparisons with CNNs. While sequence kernels that do not require the input to be aligned were used, since no sequence kernels are available in GPyTorch or scikit-learn, an RBF kernel was used. Thus, in this comparison, for the protein sequence datasets, each was converted into a multiple sequence alignment and one-hot encoded before it was used as input to an RBF kernel, as contrast to how the sequence datasets are modeled.

The data in Table 8 clearly show that xGPR outperforms SVGP on 8 out of 9 datasets, often by large margins. In the case where SVGP is competitive (Song UCI), they achieve essentially equivalent performance, and SVGP takes eighty times longer to train. xGPR generally performs better using the same kernel and at greatly reduced cost; moreover, performance for xGPR can often be further improved by either using more random features for fitting, or by fine-tuning hyperparameters with a larger number of random features. On the Kin40K dataset, for example, “fine-tuning” the hyperparameters using the marginal likelihood approximation strategy from section S5 above requires an additional 10 minutes and reduces the mean absolute error by an additional 5-8% (depending on the number of random features used to fit). Improving the performance achieved by SVGP, by contrast, would be hard, since if the number of inducing points for SVGP were increased, training time on the larger datasets would become unacceptable.

It is possible that performance for SVGP could be improved by modifying the Adam optimizer settings for stochastic gradient descent, e.g., by dataset-specific tuning of the hyperparameters of the Adam algorithm, but this is in one sense merely a further drawback for SVGP: xGPR does not require extensive experimentation with optimizer settings to achieve good outcomes. Given these results, stochastic variational inference was not implemented for the library, preferring the approach described above instead. Table 8. Performance of SVGP (as implemented in the GPyTorch library), random features (as implemented in xGPR) and exact GPs (as implemented in scikit- learn) for various benchmark datasets

S9. Lemma 1.1

Lemma 1.1. Let K-. ]R^Dx IR^D- > IR be a positive definite kernel on IR^D , and let /: X— > IR^D be any mapping from X to IR^D , where X is some non-empty set. Then ^(/(^U)</(X)) is ^a positive definite kernel on X for u,w E X.

Proof: Let f-.X-> IR^D be any mapping from X to IR^D, where X is a non-empty set. Let K-. [R^DxR^D-> IR be a positive definite kernel on IR^D , i.e.: for any choice of x₁,x₂,x₃...x_n e and for any choice of e K . Since this is true for any arbitrary selection ^ar|d since / maps its input into

K^D, it is also true that

for any choice

In other words - this is merely another way to say that any feature engineering process or function which constructs a representation in IR²³ can be “plugged into” a kernel which is positive-definite on and retain the positive definite property. This can be used to construct “plug-in” kernels, e.g., the FastConvI d kernel described above.

S10. Details for QM9 molecular modeling

The full QM9 dataset consists of approximately 133,000 molecules. First, those molecules that the geometry-optimized structure failed consistency checks were removed. The remaining data were then randomly split into a training set of 1 10,000, a validation set of 10,000, and a test set of 10,831 using random seed 123.

For the SOAP descriptors as implemented in the describe library, there are several hyperparameters that were tuned:

1 ) The standard deviation of the Gaussians used to represent the atoms in the neighborhood of the central atom (sigma),

2) Rcut, or the radius at which an atom is considered “outside the neighborhood” and no longer included in the representation,

3) Weighting, i.e., reduced weighting for atoms further from the central atom (several different schemes provided).

There are also the number of radial basis functions, n_max, and the number of spherical harmonics, /_max.These improve the accuracy of the representation if increased, but also dramatically increase the size of the descriptor vector. 12 for n_max and 9 for l_max were settled on the same values, since a larger descriptor vector size was not used for this study. Even this configuration, however, generates roughly 18,000 features per atom.

The weighting scheme and weights identified were described as optimal by Willatt et al. {Phys. Chem. Chem. Phys., 20 (47), 29661 -29668, 2018), who used power weighting in which an atom within the distance rcut is weighted based on its distance using the expression:

For all of the SOAP feature hyperparameters the values selected by Willatt et al. based on their grid search were used with two exceptions: rcut, or the cutoff distance, and rO from the weighting scheme shown above was optimized using the values from Willatt et al. as a starting point. rO and rcut were tuned by evaluating the marginal likelihood on a randomly selected subset of the training data consisting of 25,000 datapoints. Initially, it was found that almost the same accuracy could be achieved by generating features for heavy atoms (i.e., non-hydrogen) only and omitting the hydrogens. This is advantageous since otherwise the amount of disk space required to store the training set data is quite large. After tuning, the following settings were adopted and described:

N_max: 12

L_max: 9

Sigma: 0.25

R_cut: 3.25

Weighting: Power

R0: 1.5

M: 9

/Wand rO determine how rapidly the weighting falls off from the central atom. SOAP features were generated for each atom individually; since each molecule in this dataset has up to 9 heavy atoms, each molecule has up to 9 descriptor vectors corresponding to these atoms. These can be used as inputs to either the graph convolution kernel or graph polynomial kernel described above. The SOAP descriptors for each atom were divided by the norm of the descriptor vector so that each descriptor vector is unit-norm; this ensures that if the Graph Poly kernel is used, it will have the same interpretation as the original SOAP kernel.

To tune the kernel hyperparameters, the workflow described under S5 and S6 above was used. It began by finding a starting point using matrix decomposition based marginal likelihood with a small number of random features, then optimized with a larger number of random features (16,384) using Bayesian optimization. This procedure took 4 hours on a single A6000 GPU. All hyperparameter tuning was performed using internal energy at 298K as a target. The same hyperparameters that work well for internal energy at 298K also work well for the other objectives as illustrated. The final model was then fitted using either 16,384 random features, 32,768 random features or 66,536 random features. The results are shown in Table 4 above.

Small improvements in performance can be “bought” by using a larger number of random features, but even 16,384 random features are sufficient to achieve a mean absolute error 5x smaller than chemical accuracy. Consequently, in practice using larger numbers of random features has limited usefulness and it would be perfectly fine to use 16,384 random features for this problem. Moreover, the most important target for further improvement is likely the SOAP features themselves. The large number of SOAP features generated per atom greatly increases computational expense, since many features must be generated, saved, and loaded from disk while processing each minibatch of data when iterating across the dataset. Furthermore, it constrains the ability to improve the accuracy of the model by increasing the number of basis functions. Multiple authors have suggested more compact representations that achieve similar or better accuracy. A more compact representation would increase the number of basis functions, thereby improving the accuracy of the model, while decreasing the length of the feature vector, thereby reducing its computational expense.

It is interesting to note the advantages provided by the fast Hadamard transform based Convl d kernel (FHTConvI d) here as well. If vanilla random features were used instead, to generate 66,536 random features with 18,000 input features, it would require a matrix with roughly 1.2 billion floats. The FHTConv-1d kernel, by contrast, needs to store roughly 288,000 8-bit integers and 72,000 floats - a dramatic reduction in memory footprint.

S11. Clustering and kernel PCA using representations generated by xGPR

Here, the QM9 dataset was used. Out of the approximately 130,000 molecules in the dataset, 20,000 were randomly selected as a test set. Each atom and its neighbors were encoded using one-hot encoded, out to the 15th most distant neighbor (an arbitrary cutoff). Each one-hot encoded neighbor is weighted by 1 I distance**6; the 6th power used here corresponds to the weighting used in the London matrix sometimes used as a representation for machine learning on small molecules. The standard workflow was followed by tuning hyperparameters first using a smaller number of random features (2048) with a matrix-decomposition based procedure, then “fine-tuned” using a larger number of random features (8192) with an approximate marginal likelihood procedure. As illustrated, depending on the number of random features employed for fitting, this model can achieve a mean absolute error < 1.1 kcal/mol.

Next, the kernel_xpca tool provided in the visualization_toolkit module of xGPR version 0.0.2.1 was used to retrieve the top 500 principal components. The top 2 principal components were used to generate the kernel PCA plot as shown in FIG. 8, for the test set data only.

As illustrated, the top principal component correlates nicely with the ground-truth label predicted (energy at 298K). This will not always be the case, because a great deal of information is discarded by using a 2d representation in place of the 16,384 dimensional representation generated by xGPR. To further explore the data, k-means clustering was then performed using the scikit-learn library with the top 500 principal components as input, using the sum squared error to assess the performance of a given number of clusters. This operation approximates performing kernel k-means clustering on the original training data using a GraphConv-1 d kernel.

The elbow plot is shown in FIG. 9A; this result suggests 5 clusters may be appropriate, and thus the final k-means model was fitted using 5 clusters. The distribution of energy at 298 K in each resulting cluster is shown in FIG. 9B. In many drug discovery scenarios this ability to cluster the data using the same model used to fit it might in itself be useful since it can be used for efficient retrieval of similar molecules from the training as discussed above.

S12. Molecules from the test set and the most similar molecules from the training set. Five random molecules from the QM9 test set are depicted in FIGs. 10-14, together with the 50 most similar molecules from the training set, where similarity is quantified using the (approximated) GraphRBF kernel function with molecules represented using one-hot encoding as described above. Note that this representation of molecules is relatively uninformative, so it is to some extent surprising that the model performs as well as it does (MAE of about 1 kcal/mol) and can retrieve molecules that in many cases do show many points of similarity.

CONCLUSION

Gaussian processes represent a principled Bayesian approach to machine learning that provides straightforward quantitation of uncertainty. Their most important disadvantages are their unacceptable scaling and the lack of efficient kernels for sequences and graphs. This Example described efficient approaches for fitting GP models with random features-approximated kernels that achieve linear scaling with training set size and linear scaling with molecule size or sequence length. These tools were implemented in a publicly available Python library, xGPR, and it was demonstrated that this library achieves competitive performance on a variety of protein and small molecule property prediction tasks.

As discussed herein, a GP equipped with these kernels generates well-calibrated uncertainty estimates on its predictions, unlike any of the deep learning models used as comparators. This quality is of considerable importance for protein engineering and drug discovery; since evaluating a prediction is often very expensive, high-confidence predictions certainly are preferred to use, if possible.

Additionally, a GP equipped with these kernels is substantially more interpretable than the comparator deep learning models. For protein sequences and small molecules alike, it is able for each kernel to explain in detail how the model determines the similarity of any two datapoints, which in turn determines how it makes predictions. This enables understanding of the model’s limitations and how it might be improved. It is not possible to generate a detailed explanation of this kind for any of the deep learning architectures that are currently popular for sequence or graph data; they are all “black box”. Further, as illustrated for small molecules, the kernel from the trained GP regression model can be used to cluster the input data or perform kPCA, and this clustering provides additional insights about the distribution of the input data. This feature is particularly useful for drug design to catalog small molecules into subgroups.

This combination of competitive accuracy, uncertainty quantitation and improved interpretability suggests the xGPR library is potentially useful for a range of protein and small molecule property prediction tasks.

The kernels discussed in this Example all accept a single protein or small molecule as input. For some tasks (e.g., protein-small molecule interaction), it is desirable to use a “pair” kernel that accepts two inputs (a small molecule and a protein sequence for example) and compare them to other pairs. The kernels discussed here can easily be modified to generate several “pair” kernels that may be appropriate for some problems. For example, the first layer of the FastConv-1 d kernel discussed here can be used to generate a representation of both the small molecule and the protein; these representations are then concatenated and used as input to an RBF kernel. Alternatively, a Graph RBF kernel could be used to generate random features for the small molecule while the FHTConv-1 d or FastConv-1 d kernel is used to generate random features for the protein; the random features are then concatenated into a single vector to generate the prediction. This last arrangement corresponds to using the sum of the two kernels. It is likely possible to design other efficient kernels specific to pair comparisons of this kind.

EXAMPLE 2

The Rapid Engineering System for Proteins (RESP) Al Model Accelerates the Identification of Tight-binding Antibodies

This Example describes a development of an easy to implement machine learning- assisted pipeline for the identification of high-affinity antibodies that addresses these challenges. A new autoencoder model was trained on over 3 million B-cell receptor sequences, and it shows that this learned representation provides better results for a task of interest than state-of-the-art embedding schemes. Next, a Bayesian neural network was also developed and trained to perform ordinal regression to model the relationship between sequence and binding affinity or off-rate using the directed evolution data. Importantly, this model provides an estimate of the uncertainty in its predictions that cannot be achieved by the current deep learning methods. To benchmark this approach, it achieves competitive accuracy on a literature dataset, while providing useful uncertainty information not provided by the deep learning model in the original associated study. Further, in silico mutagenesis was performed using a simulated annealing strategy to explore sequences that are not present in the mutation libraries and their binding affinities were assessed for experimental evaluation. Together, the autoencoder, the Bayesian network for ordinal regression, and the search strategy form the key computational components of this RESP pipeline.

To demonstrate the power of RESP, it is attempted to improve the affinity of a well- known antibody to a well-studied antigen by engineering mutants of the heavy chain of Atezolizumab (brand name Tecentriq) with improved affinity for programmed death ligand 1 (PD-L1 ). A large portion of the Atezolizumab heavy chain sequence was randomized, the mutants were binned using yeast display and FACS, and each bin was sampled for sequencing. The mutant Atezolizumab sequence data was converted to a low- dimensional representation using the autoencoder model, and the encoded mutant Atezolizumab sequences were used to train a Bayesian ordinal regression model that scores each sequence on the probability that it is a strong binder to PD-L1 . Further, a modified simulated annealing algorithm was used to select sequences for testing. It is shown here that this pipeline discovered a panel of Atezolizumab scFv mutants with improved off-rates towards PD-L1 , and one characterized mutant displayed about a 10- fold decrease in off-rate and 17-fold improvement in KD between human PD-L1 and Atezolizumab. This method described here is useful as a general approach for improving antibody-antigen interactions while reducing the experimental effort to do so. The mutant discovered by this method could be a useful reagent for treating PD-L1 positive tumors.

METHODS

Software

Analysis and modeling were conducted using Python 3.9 with the PyTorch library version 1.8.1 , the Numpy library version 1.19.5, the Scipy library version 1.5.4 and the scikit-learn library version 0.24.2. 0 was used as a random seed for model weight initialization, train-test splitting etc. Autoencoder model for antibody encoding

The training set used for the autoencoder is derived from the cAb-Rep database, consisting of sequenced antibody repertoires from 121 human donors (Yicheng Guo, Kevin Chen, Peter D. Kwong, Lawrence Shaprio, and Zizhang Sheng. cAb-Rep: a database of curated antibody repertoires for exploring B cell response and predicting antibody prevalence. Front. Immunol.). For these studies, the high-depth repertoire dataset was used. This dataset was further filtered by numbering all sequences using the ANARCI software with Chothia numbering. During this process, the ANARCI software aligns the input sequences to domain-specific hidden Markov model profiles for human antibodies using the HMMer software. Any sequences with low bitscores resulting from this alignment are rejected, thereby minimizing the chance the dataset will contain proteins that are not actually antibodies. After numbering and filtering using ANARCI, 2,725,492 sequences remained. Since these sequences are from a database of antibody sequences and were filtered to select those that had a high probability of being generated by antibody MSAs, they are antibodies. It should be noted that the sequences used here are VH heavy chain only to match the sequences used in the experimental work.

The autoencoder accepts one-hot encoded sequences as input. To number the positions of each antibody, the Chothia numbering scheme was used. A heavy chain may have as many as 132 amino acids in this scheme, although our mutants have fewer since the size of the complementarity-determining regions or CDRs varies between antibodies. Consequently, each one-hot encoded sequence is a matrix with 132 rows and 21 columns. While there are only 20 amino acids present, the 21 st position indicates a blank, since the heavy chains do not contain the full 132 amino acids present and therefore have blanks at some Chothia-numbered positions. It is also possible to merely leave blank positions as all zeros, although explicitly encode blanks are preferred as such for clarity and adopt this approach here. Each row contains a 1 at one position to indicate which amino acid (or a blank if no amino acid) is present.

Frequently in selecting and designing sequences it is important to determine which positions are most important for affinity or other desired properties. This consideration suggests the learned representation generated by the autoencoder should have the same number of rows (positions) as the input. Moreover, in order to ensure a specific row of the encoding contains information relevant to that position, the decoder should be able to reconstruct the amino acid present at each position using the information at that position (and possibly the neighboring positions).

In keeping with these constraints, the encoder portion of the model was designed to compress the input from a 132x21 matrix to a 132x3 matrix. The encoder module consists of a convolutional network with two convolutional layers with appropriate padding and a third linear layer. The first convolutional layer contains 40 kernels, each of width 21 , while the second contains 20 kernels of width 11 . Each convolutional layer uses gated convolutional activation to improve performance of convolutional neural networks on language modeling tasks. For this activation function, the sigmoid function is applied to the first half of the columns in the output from that layer and these are then multiplied elementwise with the second half. In other words, if the output of a layer for a given sequence is a 132x40 matrix, for gated convolutional activation the sigmoid function is applied to the first 20 columns, and these are then multiplied elementwise with the remaining columns to yield a final output of dimensions 132x20.

The final output of the linear layer in the autoencoder is a 132x3 matrix, which is an encoding of the original sequence. This encoding becomes the input both to a decoder module and to a prediction module. The decoder module consists of a single layer of the following form:

Here a is one row of the output of the encoder. Since the output of the encoder is a 132x3 matrix, a is then a 3-dimensional vector, b is a learned bias vector and IV is a learned 3x21 weight matrix. Softmax is the softmax function:

Where z is the 21 -element vector resulting from a ■ W + b. The decoder layer is applied to each row of the encoder output. This design imposes a strong constraint on the autoencoder: the model is required to reconstruct the input using a decoder function with only a relatively small number of parameters shared across all positions. The output of the encoder is also supplied to a prediction module that differentiates sequences of human B-cell receptors from junk sequences. The prediction module consists of a simple logistic regression model, where the probability that the sequence is an antibody is given by: where a is the full output of the encoder flattened from a 132x3 matrix to a length 396 vector and w is a length 396 learned weight vector plus a learned bias term. Since both the prediction and decoder modules have few parameters (by comparison with typical deep learning models), neither can learn a complicated mapping from encoded sequence to input, thereby placing the burden on the encoder to generate as informative and relevant an encoding as possible.

To generate the junk sequences, a mutant version of each of the 2,725,492 sequences selected from the cAb-Rep database was generated. The goal here is to force the model to incorporate information about the relative abundance of specific amino acids at specific positions into its encoding by requiring it to be able to distinguish true antibody sequences from closely related sequences. It is possible to simply generate random sequences, but this would not force the model to generate an informative encoding; random sequences are so different from antibodies that the logistic regression model which predicts whether a sequence is or is not an antibody would be able to distinguish them without any modifications to the representation generated by the encoder. It is therefore preferable to generate decoy sequences that are just similar enough they will be difficult to distinguish and yet different enough to exhibit modifications rare in true antibody sequences.

The experiments suggested 7 mutations provided a good balance; consequently, each copy was altered to a randomly chosen amino acid at seven randomly selected positions. The end result of this process was thus a library of roughly 6 million sequences, half of which are human B-cell receptors and the other half of which are not. The autoencoder model is thus trained both to encode an input sequence and to embed information about typical features observed in true antibody sequences. The autoencoder was implemented using the PyTorch library in Python 3.6.9 and trained on the full 6 million sequence dataset until convergence. The code for this and all other steps described in this paper is available online at: (10.5281 /zenodo.7508853).

Accuracy was assessed separately for both the prediction task and the reconstruction task using a held-out test set. These metrics are used only as diagnostics because they assess the ability of the autoencoder to reconstruct its input. The true test of the autoencoder is the degree to which the encoding it generates affects predictive accuracy of a model trained with that learned representation as input. To evaluate this more critical metric, the WT and mutant Atezolizumab library (construction described below) were encoded using the autoencoder described above, one-hot encoding, the ProtVec encoding scheme, the UniRep encoding scheme, the FAIR-ESM encoding scheme, the AbLang seq-coding scheme and the AntiBertY embeddings. For AbLang, the seq-codings option was used. For the AntiBertY embeddings, cross all the residue- specific representations in the sequence were averaged. For AntiBertY, initially, all the residue-specific embeddings were used without averaging across them, but this led to poor performance. A Bayesian neural network (construction described below), a fully connected neural network with the same number of layers and weights as the Bayesian network, was trained, and a random forest model was used as a baseline using these available encodings. Classification performance of each model on each encoding was evaluated using 5x cross validation on the training set and a held-out test set, with Matthews correlation coefficient and AUC-ROC for identification of RH03 vs rest as metrics. A good encoding should improve or at least not damage performance relative to one-hot encoding across all three model types and especially the Bayesian network, which is the model of most importance for this study.

Generation of the WT and Mutant Atezolizumab scFv Library

The WT Atezolizumab scFv was first cloned into the pYD1 yeast display vector to test its function on the yeast surface (see FIGs. 24A-24B). For the 1 st Atezolizumab library, the WT plasmid was used as a template for PCR to prepare either the WT light chain with Q5 hotstart DNA polymerase (NEB) or the mutated heavy chain by error-prone PCR with Taq Polymerase (Invitrogen) as previously described. The 2 PCR products were assembled into one product by overlap extension PCR and co-electroporated into EBY100 yeasts along with linearized pYD1 vector as described. The library transformation resulted in 7.8 x 10⁷ transformants, based on colony counts after serial dilution onto selection plates.

Atezolizumab scFv Library Screening by Yeast Surface Display

The optimal competition time for the off-rate screens was determined as previously described. The WT or mutant library yeasts were thawed and inoculated into selective growth media and grown at 30°C for 22 hrs. The library/WT were induced at 20°C in galactose induction media 42 hrs. Afterwards, WT or library was labeled with biotin-PD- L1 (Sino Biological 10084-H08H-B) 3 hrs at RT in TBS-BSA, followed by incubation with excess non-biotin PD-L1 (Sino Biological 10084-H08H) for the determined competition period (at RT, in TBS-BSA). For FACS, yeasts were labeled with anti-V5 (Thermo Fisher R960-25, previously known as 46-0705) at 4°C in TBS-BSA, followed by SA-PE (BD 554061 ) and goat anti-mouse lgG2a AF647 (Thermo Fisher A21241 ) on ice for 30 min in TBS-BSA. Cells were sorted for faster, moderate, and slower off-rates (see FIGs. 16Ai 2 & FIG. 34). Hits were grown up to high density at 30°C and made into frozen stocks at - 80°C.

Preparation of the DNA Libraries of Mutants with Faster, WT-Level, and Slower Off- Rates to PD-L1

To isolate plasmids from each binding group from the initial library screen, plasmids were harvested from yeasts by Zymoprep Yeast Plasmid Miniprep II kit (Zymo Research), eluted with ddH20, concentrated with the DCC-5 (DNA Clean and Concentrator 5) kit (Zymo Research), then subjected to the 1 st round of PCR using primers SeqF/SeqR (Table 1 1 , designed to only amplify the heavy chain regions). The 1 st PCR was performed with Q5 Hotstart DNA Polymerase (NEB) using a moderate number of PCR cycles (17, determined to be optimal by qPCR) to avoid over-amplification of the library and the PCR product concentrated with DCC-5 kit and purified by agarose gel extraction. The second PCR was carried out with various primers (Low/Medium/High Binder NGS F, NGS R) to barcode each binding group (see Table 1 1 ) for MiSeq Nano PE250 using KAPA HiFi Hotstart Ready Mix (KK2601 ) for 5 cycles (determined optimal by qPCR). The PCR products were purified with Ampure XP beads (A63880) and submitted for QC using TapeStation analysis (Agilent HS D1000) before submission for MiSeq sequencing at the Institute of Genomic Medicine, UC San Diego.

Generation and Screening of the Focused 21 -Mutant Library (see also S2.4 below)

21 geneblock fragments (IDT) for the 21 mutant heavy chains were fused with the light chain by overlap extension PGR (Q5 hotstart), followed by co-electroporation into EBY100 yeasts with the linearized pYD1 vector as with the 1 st library (>10⁷ transformants). The resulting library was screened essentially as for the 1 st library, except the final sort involved a very stringent competition time (39 hrs at RT). The hits were harvested by yeast plasmid miniprep, transformed into GC10 competent cells, harvested by bacterial miniprep, and sequenced using standard Sanger sequencing.

Cloning WT Atezolizumab, Mutant 4, Durvalumab, and Avelumab scFv into the pYD1 Vector

The geneblocks for Durvalumab and Avelumab were purchased from IDT with yeast-optimized codons and PGR amplified using the following primers: Durvalumab with Atez LE F & Durv LE R, Avelumab with Atez LE F & Avel LE R (Table 11 ). The WT/Mutant 4 genes were PCR amplified from their plasmids using Atez LE F/R primers. The PCR products were double digested using Xhol/Nhel-HF (NEB) and ligated into linearized/dephosphorylated pYD1 (using the same enzymes to digest and rSAP/T4 DNA ligase (NEB) to dephosphorylate/ligate) and the ligation product transformed into GC10 competent cells (42-658, Genesee Scientific). Plasmids were isolated using a Zippy Plasmid Miniprep Kit (Zymo Research) and sequence verified before transformation into EBY100 yeasts. These constructs were used for the yeast koft and KD determinations.

Determination of the WT/Mutant 4/Durvalumab/Avelumab KD and k₀« Values on the Surface of Yeast

The koft was determined at RT essentially as described and in section S2.3 below in TBS-BSA. The resulting data was fit to a one phase decay model with GraphPad Prism 9.3.0 software using the following equation (Y is fraction of yeast bound to bioti n-PD L1 , X is time, Y0 = 1 (fraction of yeast bound to biotin- PDL1 at the Os competition time point), plateau is a constant based on nonspecifically bound yeast, and K is Koff): (4) Y = (YO - Plateau)*exp(-K*X) + Plateau

The KD values on the surface of yeasts were determined basically as described⁴⁹ and the data fit to the following equation:

(Bmax is maximum MFI value, X the concentration of PD-L1 ). Due to the very slow off- rate of Mutant 4, it was necessary to incubate yeasts with PD-L1 for 6 days at RT.

Cloning and Purification of WT/Mutant 4 Atezolizumab scFv (also see S2.6 below)

The scFv sequences were PCR amplified from the pYD1 vector and cloned into the pET27b(+) (69863-3, MilliporeSigma) vector for bacterial expression. The vectors were transformed into Rosetta (DE3) cells (70954-3, MilliporeSigma) and the scFv-Hise fusions purified essentially as described. The scFv stocks were frozen in PBS + 10% glycerol at -80°C. SDS-PAGE was used to assess purity (FIGs. 28A & 28B) and protein concentration determined by BOA assay (Pierce #23227).

Determination of the KD Values by Bio-Layer Interferometry (BLI)

BLI was performed at the Biophysics and Biochemistry Core at The Scripps Research institute on an Octet Red96 at 25 °C. Binding reactions were performed in 1 X kinetic buffer (Sartorius, 18-1105) consisting of 20 mM phosphate buffer, pH 7.6, 2 mM KCI, 150 mM NaCI, and .02% Tween 20, 0.1 % BSA, .05% sodium azide. Biotinylated PDL1 (Sino Biological 10084-H08H-B) was immobilized on SA Biosensors (Sartorius, 18- 5019) by dipping the sensor into 100 nM b-PDL1 until the signal was saturated. A blank loading channel was used as a negative control. Kinetic experiments were performed with both scFv WT and scFv Mutant 4. Experiments were performed using a gradient of concentrations, with the scFv WT ranging from 475 nM to 0.6 nM, and scFv mutant ranging from 160 nM to 0.2 nM. Binding kinetics were assessed via Octet Data Analysis HT Software Version 12 using a 2:2 binding model.

Mass Photometry

Mass Photometry was performed on a Refeyn TwoMP at the Scripps Research Institute, at room temperature. The PD-L1 used was biotin-PD-L1 (SinoBiological 10084- H08H-B) and purified scFv of the WT or Mutant 4. Each characterization was performed in 1 X PBS (phosphate buffer from Cytiva BR100672). For each experiment, the scFv or PD-L1 (or both) were diluted with PBS to the final concentration. Each experiment resulted in a 60s movie, and mass analysis performed using Refeyn DiscoverMP analysis software v2.3.0. Mass calibration was performed using Urease (Sigma, U7752) and Thyroglobulin (Millipore, 609310).

Sequence Processing

Raw paired end reads from the Atezolizumab dataset were checked for quality (for details of the filtering criteria, see the Sequence Quality section). After the sequences had been processed, they were split into an 80% training set and a 20% test set. All training and cross-validation was performed on the 80% training set only.

In many cases, sequences occurred in more than one category but with a different frequency in each. Clearly, the level of confidence in category assignment is reduced when the frequency of the sequence in the assigned category is not much greater than its frequency in others. To encode the level of confidence, each sequence was weighted with the frequency in the assigned category plus one divided by the total number of occurrences plus three. This is the posterior probability using a multinomial likelihood and a uniform Dirichlet prior (a Dirichlet distribution with a = [1 ,1 ,1 ]).

The ordinal regression model used here employs binary cross-entropy loss: log(l - p₀)) where p for datapoint / is the model-assigned probability that sequence belongs to a more stringent/tighter-binding category than category /, out of M possible categories and for N data points, and y(i,j) is either 0 or 1 and indicates whether the sequence does or does not belong to a tighter-binding category. (For details of how p is calculated, see below). This loss function for sequence / is multiplied by the weight for sequence / so that the model is weakly penalized for misclassifying sequences where our confidence in the category assignment is low. Construction of the Ordinal Regression Model

The Atezolizumab mutant library data was analyzed using Bayesian neural network-based ordinal regression. Each mutant Atezolizumab sequence is first encoded by the autoencoder, yielding a 132x3 matrix which becomes the input to the Atezolizumab model. The model architecture is illustrated in FIGs. 23A-23B. At a high level it is similar to the Bayes by Backprop architecture described by Blundell et al. (ArXiv 150505424 Cs Stat, 2015, except this model has been adapted to perform ordinal regression as described below.

In a traditional neural network, each parameter is a learned parameter that is fixed once the model is trained. In the variational network, by contrast, each parameter has an associated Gaussian distribution described by a mean and a standard deviation whose optimal values were learned during training. To generate predictions, rather than using the fixed learned parameters as in a traditional network, the weight distributions A/ times to generate A/ predicted values were sampled from the distributions over the weights. The variance in these predictions provides a measure of our uncertainty around the final prediction. This approach requires a more complicated training procedure than that associated with traditional neural nets discussed below.

Briefly, parameters 0for a distribution on the weights wis sought, that will minimize the evidence lower bound or ELBO given by:

(7) argmin_e KL[q(w\0) \\ p(w^)] - E_qMe)[logQp(iD\w^] where q(w|0) is a distribution over the weights that approximates the true Bayesian posterior, p(w) is the prior and p(P\w~) is the likelihood of the data. The parameters 0 are the parameters of the normal distributions for all the weights in each hidden layer as illustrated in FIGs. 23AI-2 & 23B. The second term is the negative log likelihood, while the first term is a regularization term that measures the divergence between the weight distribution learned by the model and the prior. A Cauchy prior with unit scale and location zero were chosen, which adds some additional flexibility since some of the weights are now expected to be outliers.

The cost was approximated via a Monte Carlo sampling procedure. Using the current set of parameters 6 (the means and standard deviations of the weight distributions), for each minibatch n sample sets of weights were drawn to yield the following approximate cost function: where M is the number of minibatches and the other terms are as above. In other words, for each minibatch of training data, n sets of sample weights were drawn and then the approximate cost function across these and across the minibatches in the training set were averaged. Backpropagation for this cost function is made tractable by using the reparameterization trick of Kingma and Welling. A standard normal distribution was sampled with mean zero, standard deviation 1 , then the mean was added and multiplied by the standard deviation of the distribution for weight j. In order to ensure that the standard deviation is always positive, instead of using the standard deviation itself as a parameter, each distribution was parameterized with a parameter p that is converted to the standard deviation using the softplus function:

By using this reparameterization trick, the gradient of the approximate cost function with respect to the mean and p of the distribution for each weight and bias term in each hidden layer is easily calculated. The advantages of this model structure are two-fold. First, it imposes strong regularization on the model parameters that can provide improved performance for some tasks. Second, it estimates the uncertainty in the predictions and thereby ascribes greater weight to the most confident predictions when selecting sequences for testing.

All these features of this model are shared in common with most Bayesian neural network architectures. Thus, this model was used to perform ordinal regression. As illustrated in FIGs. 23A-23B, the last hidden layer of the Bayesian network outputs a single latent score value. This score is added to M - 1 learned threshold values for the case where there are M categories. The sigmoid function is then applied to each of the M - 1 outputs to generate an output vector of M - 1 probability values. Each element / of this output vector is the model-assigned probability that the input sequence belongs to a binding category more stringent than /. Since there are three binding categories (RH01 , RH02 and RH03 or weak, moderate and strong), the output vector is 2-dimensional; the first element indicates the probability that the sequence belongs to either RH02 or RH03, while the second element indicates the probability the sequence belongs to RH03. The model is trained by minimizing the binary cross-entropy loss described above.

This arrangement treats the categories as ranked: sequences with a higher latent score are thereby assigned to higher categories. Note that it does not predict the actual off rate or binding affinity. Rather, the score reflects the confidence the sequence is a strong binder relative to others in the training set. This approach has been used in literature to rank sequences for a protein engineering task and select them for experimental evaluation.

Sequence Scoring and Selection

The classic simulated annealing algorithm was modified and equipped with trained models to perform in silica directed evolution as illustrated in FIG. 17. First, the frequency of each amino acid at each position across the entire dataset was computed, added 1 to all values and divided by the total number of sequences in the dataset plus 20 to retrieve a marginal probability for each amino acid at each position. By adding 1 to all frequencies, there is a small but nonzero marginal probability for amino acids not observed in the dataset.

Next, the 500 highest-scoring sequences were selected in the Atezolizumab dataset and the top 10 most frequently mutated positions in these sequences were identified. The selection of 10 sites here is arbitrary, since more or fewer could be used, if desired. On each iteration, equal probability any of these top ten sites was selected. The selected site is randomly reassigned to a new amino acid; the probability for the selection of any given new amino acid is determined by the marginal probabilities calculated as described above. Assume, for example, that 80% of all sequences observed in the dataset carry an arginine at position 100, another 10% carry lysine and so on. If position 100 is selected on a given iteration, this position will be mutated with an 80% chance of being converted to an arginine, a 10% chance of being converted to a lysine and so forth.

The current sequence and the mutated proposed sequence are both encoded using the autoencoder and the trained ordinal regression model is used to assign a score to each. The proposed sequence is accepted with a probability given by:

~t^sbest~^sproposed)

(10) p (acceptance') = e T

Where Tis the temperature and S_best and S_proposed are the best score to date and the score of the proposed sequence respectively.

This algorithm is further illustrated as follows:

Modified simulated annealing

This procedure is a simple approach for exploring the sequence space while ensuring it is not venture too far from the training set. To ensure reproducibility and avoid stochastic fluctuations, for this stage of modeling the Bayesian neural network generated predictions using the means of the distribution for each weight as the weight value.

S2.1 Sequence Quality

Sequence read pairs that contained one or more base pairs with a phred quality score < 10 were discarded since the sequence read in this case may be unreliable. In the event the overlapping region of the paired ends did not match, both reads were discarded, so that no mismatches between the paired ends were allowed. The reads that met these quality criteria were merged and translated to yield mutant atezolizumab sequences of 1 18 amino acids in length. Each sequence that occurred more than once in any given category was assigned a frequency for that category indicating the number of times it was found. Sequences with mutations in the first 30 or last 8 positions were excluded from further consideration, since these positions were not targeted for mutation and mutations at these positions were very rare, so it was difficult to assess their importance.

S2.2 Generation of the WT and Mutant Atezolizumab scFv Library

To test the function of the WT scFv in the yeast display format, the gene for WT Atezolizumab scFv was purchased as a geneblock (using codons optimized for yeast) from Integrated DNA Technologies and cloned into a modified pYD1 (Addgene #73447) yeast display vector by PCR amplification of the geneblock with AtezF/AtezR primers (Table 1 1 ) and then using double digestion/ligation into a digested vector by Golden Gate Assembly using Esp3l (Thermo Scientific) and T4 DNA Ligase (New England Biolabs (NEB)) and sequence verified before transformation into EBY100 yeasts via heat shock. For the Atezolizumab library, the WT scFv in the pYD1 vector was used as template for either high fidelity PCR of the light chain region with Q5 hotstart DNA polymerase (NEB) (primers IF1 F/IF1 R) or error-prone PCR of the heavy chain region with Taq DNA polymerase (Invitrogen 18038-018, IF2F/IF2R) using 8-oxo-dGTP (TriLink N-2034-1 ) and dPTP (TriLink N-2037-1 ), except 30 cycles of error-prone PCR with 200 pM each dPTP and 8-oxo-dGTP (final concentration 20 uM each) were used to increase the mutation rate. The PCR products of the light/heavy chains were concentrated by the DNA Clean and Concentrator-5 kit (DCC-5, Zymo Research) and purified by an agarose gel extraction kit (Zymo Research gel extraction kit). The 2 DNA fragments (WT light chain, mutated heavy chain) were assembled into a single DNA molecule by overlap extension PCR using the OF/OR primers (all primers in Table 11 ) with Q5 hotstart DNA polymerase, the resulting product was again concentrated by DCC-5 kit and purified by gel extraction. TA cloning (Invitrogen TA Cloning Kit, K202020) by adding 3’A to the product with Taq was used to examine the mutation rate in the heavy chain, random sequencing of clones revealed 0-8 mutations per gene. !n preparation for yeast electroporation, the pYD1 vector was double digested with Xhol and EcoRI-HF, concentrated by DCC-5 kit, and gel extracted/eluted into ddH20.

The scFv library was transformed into yeast by electroporation (3 pg digested pYD1 , 9 pg scFv gene per electroporation) and assembled in the EBY100 cells (ATCC MYA-4941 ) by homologous recombination⁶, giving 7.8 x 10⁷ transformants as determined by serial dilution onto plates with selective growth media (16.7 g/L BD bactoagar for solid media). The library was passaged several times in selective growth media to ensure 1 plasmid/cell before being frozen in aliquots at -80°C in 85% ddH20, 10% glycerol, 5% DMSO (3.6 x 10⁸ cells/tube). As a control for sorting, the WT scFv gene/pYD1 was also electroporated into EBY100 so as to be in the exact same position in the vector as the library.

S2.3 Atezolizumab scFv Library Screening by Yeast Surface Display

First, to determine the optimal competition time, the WT scFv koff (off-rate) towards PD-L1 was determined basically as previously described. Yeasts expressing the WT scFv were washed with TBS-BSA and resuspended to 16 x 10⁶ cells in 400 uL volume. Then, 12 uL b-PDL1 (final concentration 178 nM, Sino Biological 10084-H08H-B) was added to the yeasts and labeled at RT for 3hrs. Afterwards, the yeasts were washed with TBS-BSA and resuspended to 800 uL total volume with 174 nM unlabeled PDL1 (Sino Biological 10084-H08H). 2 x 10⁶ cells were added to each eppendorf tube and rotated at RT for x hrs (0-21 hrs). For each timepoint, the cells were washed 3x with TBS-BSA (cold) and kept on ice in 198 uL TBS-BSA until the experiment was completed. Once all data points were collected, 1 :100 V5 antibody (R960-25, previously known as 46-0705) was added and the yeasts labeled on ice for 1 hr. The yeasts were washed with TBS-BSA (cold) and resuspended in 200 uL with 1 :100 SA-PE (BD 554061 ) and goat anti-mouse AF647 (Thermo Fisher A21241 ) and labeled on ice for 30 minutes. Afterwards the yeasts were washed with TBS-BSA (cold) and resuspended in TBS-BSA (cold) for flow cytometry. The optimal competition time was determined to be 29 hours at RT. To screen the library against PD-L1 , either the WT or library was thawed and added to selective growth media (20 x 10⁶ cells into 50 mL growth media + ampicillin, grown with shaking at 30°C for 22 hrs, followed by induction by pelleting the cells, washing them once with sterile ddH20, resuspension into the same type of media but with 20 g/L galactose instead of glucose, pH adjusted to 5.6 and ODeoo = 0.2 with ampicillin). The scFv were induced for 42 hr at 20°C with shaking. After 42 hours, the flasks were placed on ice and 40 x 10⁶ yeasts (WT or library) were pelleted at 17,000xG for 30s, resuspended in 1 mL TBS-BSA (25 mM Tris, 150 mM NaCI, pH 7.5, 5 mg/mL BSA), re-pelleted, resuspended in TBS-BSA buffer at 40 million cells/mL. WT or library was then labeled with biotin-PD-L1 (Sino Biological 10084-H08H-B) at a final concentration of 178 nM for 3 hrs at RT with rotation. Afterwards, the cells were pelleted/washed with TBS-BSA, then resuspended in non- biotin PD-L1 (Sino Biological 10084-H08H, 174 nM) and rotated for various amounts of time (mostly 29 hours, but one screen was for 7 hours) at RT. For sorting by FACS, cells were pelleted/washed 3x with cold TBS-BSA and labeled with anti-V5 antibody (Thermo Fisher R960-25, validated for specificity by the vendor, lots 2001339/2249078/2212258) at 1 :100 dilution on ice for 60 minutes (at 10⁷ cells/mL), pelleted/washed with cold buffer, then labeled with 1 :100 each secondary detection reagent (PE Streptavidin BD 554061 , goat anti-mouse lgG2a AF647 (Thermo Fisher A21241 ) on ice for 30 minutes at 10⁷ cells/mL. Cells were then pelleted/washed with cold buffer and resuspended in cold TBS- BSA buffer for sorting on a BD FACSAria II (Moores Cancer Center, UC San Diego) using BD FACSDiva v. 8.0.1 software, see FIG. 34 for example FACS plot. Sort gates were set to collect mutants with faster, WT-level, or slower off-rates to PD-L1 , see FIG. 34 for example of gating. Hits were collected in selective growth media and grown up to high density at 30°C and aliquoted to frozen stocks at -80°C. The process was repeated twice for binders with faster or WT-level off-rates and four times for the mutants with slower off- rates.

S2.4 Generation and Screening of the Focused 21-Mutant Library

Twenty-one geneblocks coding the 21 heavy chain mutants were purchased from IDT, followed by PCR amplification, and overlap extension PCR (Q5 Hotstart DNA Polymerase) to generate 21 full-length Atezolizumab scFv mutant genes. The sequences of the 21 heavy chain mutants appear in Table 12. The genes were electroporated (along with linearized pYD1 ) into EBY100 yeasts to assemble the focused library and the library was passaged/frozen as described for the initial library. Titration onto selective media plates revealed >10⁷ transformants.

In order to screen the library for the most improved (slowest off-rate) mutants, increasingly long competition times were used in the presence of excess unlabelled PD- L1 (after initially saturating the WT or mutant library with 200 nM biotin-PD-L1 , the cells were pelleted/washed to remove unbound biotin-PD-L1 and resuspended in excess (128 nM unlabelled PD-L1 ) for 7.5h, 17h, then 39h at RT). The most intense clones were collected by FACS (using the WT scFv as a reference) in growth media and propagated to the next round of sorting. After 3 rounds of selection, plasmids from the hits were harvested by yeast miniprep and transformed into GC10 competent cells, which were picked as individual colonies and subject to bacterial miniprep. Plasmids were sequenced using Eton Bioscience sequencing services.

S2.5 Determination of the WT/Mutant 4/Durvalumab/Avelumab KD and koit Values on the Surface of Yeast

The RT Kott (off-rate) of the scFv-PD-L1 complex was determined essentially as described using biotin/non-biotin PD-L1 (same versions used during library screening), with TBS-BSA as the buffer and the same antibody reagents used for library screening. The yeasts were grown/induced in the same manner as with the screening protocol and labeled with 200 nM b-PD-L1 at 4 x 10⁷ cells/mL, pelleted/washed, followed by 64 nM unlabelled PD-L1 (at 10⁷cells/mL) and rotated at RT. The mean fluorescence intensity (MFI) of the V5 positive yeasts was recorded at each time point and the data fit to a one phase decay model with GraphPad Prism 9.3.0 software ( Y= (Y0 - Plateau)*exp(-K*X) + Plateau, Y is fraction of yeast bound to biotin-PD L1 , X is time, Y0 = 1 (fraction of yeast bound to biotin-PDL1 at the 0s competition time point), plateau is a constant based on nonspecifically bound yeast, and K is Koff). The KD on the yeast surface was determined as previously described, where the MFI of the V5+ population was plotted vs. the antigen concentration and the biotin-PD-L1 the same as used for the library screening. Important notes on the KD determination include using 10⁵ cells per data point, using 1 :200 anti-V5 antibody after resuspension of each data point in 100 uL TBS-BSA (cold) for 1 hr, washing cells with cold TBS-BSA after V5 mAb incubation, then labeling cells with 1 :100 SA- PE/lgG-AF647 for 20 minutes on ice, and washing cells in cold TBS-BSA and resuspension in cold TBS-BSA for flow cytometry. The data was fit to the following equation: y = Bmax*X/(KD+X) (Bmax was maximum MFI value, X the antigen concentration, the MFI without PD-L1 was subtracted from all values). Due to the very slow off-rate of Mutant 4, it was necessary to incubate the samples at RT for 6 days to allow the system to get to equilibrium.

S 2.6 Cloning and Purification of WT/Mutant 4 Atezolizumab scFv

The genes encoding the scFv sequences were amplified by PCR out of the pYD1 vector with primers containing Ncol-HF and Xhol cut sites (scFv F/R, Table 11 ) using Q5 hotstart DNA polymerase (NEB). The PCR products were double digested (along with the pET27b(+) vector (69863-3, MilliporeSigma) (the vector was also dephosphorylated with Quick CIP (NEB)). The products of the digestion were ligated using T4 DNA Ligase (NEB) and the ligation transformed into NEB 5-alpha cells. The plasmids were purified using the Zippy Plasmid Miniprep Kit and sequence verified.

The scFv plasmids were transformed into Rosetta (DE3) cells (70954-3, MilliporeSigma) onto LB-0.5% glucose plates with 30 pg/mL kanamycin, 34 pg/mL chloramphenicol and grown in LB-0.5% glucose media supplemented with kanamycin/chloramphenicol (30 pg/mL kanamycin, 34 pg/mL chloramphenicol) (10 mL) for 5 hours at 37°C with shaking until ODsoo = 0.8-1 . Afterwards, 4 mL of the media was added to 400 mL autoinduction media (71300-4, MilliporeSigma, with 100 pg/mL kanamycin and 34 pg/mL chloramphenicol) and shaken at 200 RPM at 25°C for 24 hr with added antibiotics. Afterwards, the scFv-Hise fusions were purified essentially as described, except 5 mM MgSO4 was used instead of MgSOs, supernatant was mixed with 0.3M NaCI (not 0.5M), 5 mL of Ni-NTA agarose (Qiagen 30210) was used instead of 0.5 mL, and protein was exchanged into PBS buffer (140 mM NaCI, 2.7 mM KCI, 10 mM NazHPO^, 1 .8 mM KH2PO4, pH7.4) after elution using Vivaspin-20 (10,000 MWCO, Cytiva 28932360) filters by centrifugation. Also, the Ni-NTA agarose resin was washed with 50 mM HEPES, 300 mM NaCI, 50 mM imidazole, pH 7.5, then the same buffer but with 70 mM imidazole, followed by elution with the same buffer but with 250 mM imidazole. The concentrated scFv was diluted to have 10% glycerol and flash frozen and stored at -80°C. Protein purity was assessed by SDS-PAGE (FIGs. 28A & 28B) and concentration determined by BCA assay.

Statistics and Reproducibility

No statistical method was used to predetermine sample size. When processing raw sequence data, unreliable sequence reads (reads containing one or more bases with a phred quality score < 10 or where the paired end reads did not match in the overlap region) were discarded before any further analysis or processing was conducted. These steps were taken to ensure that only reliable reads were used for analysis. No data was otherwise excluded from any subsequent analysis or model training. The test set for evaluating model performance was constructed by randomly selecting 20% of the assembled sequences and assigning these to test. The random partition was generated using the Mersenne Twister random number generation algorithm as implemented in Python’s numpy library version 1 .19.5 with a seed value of 0. When model performance was assessed using cross-validations, the cross-validation splits were generated by randomly partitioning the dataset into 5 splits of equal size using the KFold function in Python’s scikit-learn library version 0.24.2.

The final evaluation of model performance was conducted “blind” by generating predictions for sequences that are not present in these data and experimentally evaluating these predictions as described above.

DATA AVAILABILITY

The construction of the cAbRep database is described in Guo et al. {Front. Immunol. 10, 2365, 2019). Source data are provided with this paper.

CODE AVAILABILITY

The code used in this study is available online (DOI 10.5281 /zenodo.7508853), together with instructions on how to reproduce all key computational experiments.

RESULTS

Overview of the RESP pipeline

The RESP pipeline comprises four key components. First, a simple new encoding scheme was developed in which an autoencoder is designed to learn representations that incorporate features distinguishing human B-cell receptor (BCR) sequences from closely related sequences. This learned representation enables more efficient and accurate modeling of trends in fluorescence-activated cell sorting (FACS) data than provided by other popular learned representations of protein and antibody sequences (see Table 9). This part of the pipeline is general to any antibody sequence and this representation can be reused for any project.

Table 9 provides comparisons of both different encoding types (one-hot, the autoencoder, Uni Rep, ProtVec etc) and different models (a random forest model, a Bayesian network and a traditional fully connected network) based on classification accuracy.

Table 9 Performance comparisons across different encoding types and model architectures for classification performance on the Atezolizumab dataset

Second, a yeast surface display library of mutants was constructed for a starting antibody sequence (e.g., Atezolizumab here, where residues were randomly mutated in the antibody heavy chain). The mutant library was incubated with the labeled target antigen (e.g., PD-L1 here) and screened for mutants with slower dissociation kinetics (slower off-rate - i.e., likely tighter binders) by incubation in the presence of excess unlabelled target antigen. The flow cytometry experiment collects mutants with lower, moderate, and faster off-rates than the WT antibody, thereby binning the population into a series of ranked groups (FIG. 15 & FIGs. I 6A1-2 & I 6B1-4). The sequences are then determined by sequencing and their group identities are distinguished using an encoding scheme.

Third, the sequence off-rate was predicted by developing a variational Bayesian neural network to model the experimental data through ordinal regression. This affinity model takes as input the encoded sequences from the yeast surface display library and tries to predict the likelihood that a given sequence has a slow off-rate (i.e., is a tight binder). The variational Bayesian architecture provides strong regularization that minimizes the risk of overfitting and estimates the model’s uncertainty on each prediction. By using ordinal regression, each sequence was mapped to a one-dimensional latent score that indicates the model’s level of confidence to which sorting group the sequence belongs. This approach naturally takes into account the ranked ordering of the groups and is a useful approach as demonstrated previously for protein engineering problems⁴¹. The affinity model should be retrained for a new antigen or a different starting antibody sequence, while its architecture can be kept unchanged.

Further, the simulated annealing algorithm was modified to develop an in silica directed evolution algorithm that harnesses the three previous pieces of the pipeline to efficiently explore the sequence space surrounding the training set and sequences likely to exhibit off-rates significantly lower than that of the parent sequence identified. The resulting analysis pipeline is illustrated in FIG. 17. This component provides an efficient approach for exploring the sequence space and can be performed using a Bayesian neural network trained for any antibody-antigen pairing of interest.

In this study, the RESP pipeline was applied to improve the binding affinity of Atezolizumab, an antibody targeting PD-L1. Twenty-one (21 ) candidate sequences with predicted slower off-rates than the WT sequence were selected. One of them has an off- rate about 10-fold slower and KD 17-fold tighter than the WT Atezolizumab scFv.

Furthermore, the core computational components of RESP were evaluated on data acquired by Mason et al. {Nat. Biomed. Eng. 5 (6), 600-612, 2021 ) and it was demonstrated that the models described herein showed the same accuracy as the CNN model in the original study while providing additional useful uncertainty information. Each component of the pipeline is further discussed in more detail below.

Encoding antibody sequences using a new autoencoder model

First, a new autoencoder model trained to represent antibody sequences and distinguish true antibody sequences from other closely related sequences was developed. The requirement to distinguish antibodies from closely related sequences might force the encoder to embed information about typical preferences at specific positions into its learned representation. To this end, a convolutional autoencoder consisting of three main modules or components was built, that is illustrated in FIG. 18, and its performance was compared with one-hot encoding, UniRep, ESM-1 b, ProtVec, AntiBertY and AbLang. This component of the pipeline is general and does not need to be retrained for a new antigen or wild type.

The training set for the autoencoder from the cAb-Rep database was drawn, which contains sequenced B-cell receptor repertoires for 121 human donors representing true antibody sequences. Specifically, the high-depth repertoire dataset was filtered using the ANARCI software package to remove incomplete sequences and the surviving sequences were numbered using the Chothia numbering scheme, resulting in a dataset of 2,725,492 sequences. Then, this dataset was augmented with an equal number of decoys, generated by making a copy of each true antibody sequence and randomly modified at 7 positions. This number was chosen to strike a balance: increasing the number of mutations makes it more unlikely that any mutant will coincide with actual human B-cell receptor sequences occurring in nature. Too many mutations, however, makes it too easy for the classifier to distinguish the decoys from the rest, so that the encoder unit will no longer be forced to learn an information-rich representation of the input. Different numbers of mutations were experimented, and 7 mutations provided a good balance.

To ensure that the decoys were in fact different from typical human antibody sequences, a random sample of approximately 50,000 original sequences was scored and decoyed with three different models of BioPhi, AbLSTM, and ANARCI (see FIGs. 19A-19C). Using any of the three models I tools, the decoys exhibit significantly different human-ness scores from the original sequences (two-sided Mann -Whitney U test). For all three rating methods, using Python 3.9, Scipy 1 .5.4 and Pingouin 0.5.2, both Scipy and Pingouin indicate the p-value is too small to calculate reliably given floating point error on double precision arithmetic and return » 0.0). The decoys are therefore clearly rated by the models as less human, demonstrating that they are suitable to serve as decoys.

The autoencoder was trained on the full one-hot encoded cAb-Rep plus decoy dataset with a test set of 200,000 sequences set aside for 8 epochs at which point the training loss had converged. The reconstruction accuracy on the test set was > 99.99%, while the prediction accuracy for the B-cell receptor vs mutant task was 97.4%. These results suggest the autoencoder can compress the input sequence while retaining all the information needed to reconstruct or classify it. The performance of the autoencoder is assessed as an input for models to predict binding affinity below under the Atezolizumab modeling section.

Generation and Screening of the Atezolizumab scFv Heavy Chain Library for Improved Off-Rates

To use the RESP pipeline to develop antibodies against a target antigen, training data specific to that target antigen must first be generated, so that a model can be trained to recognize the types of sequences that will bind it. Here, the test case is the antigen PD-L1 , and the starting point is the single-chain fragment variable (scFv) heavy chain of the Atezolizumab antibody against PD-L1 , hereafter termed the wild-type (WT).

The WT Atezolizumab scFv was first tested for functionality in the yeast display format by testing binding to human PD-L1 (FIGs/ 24A-24B), and robust binding with 6.1 nM of antigen present was demonstrated. The variable heavy (VH) domain of Atezolizumab was chosen for mutagenesis because the structures of the Fab of Atezolizumab bound to PD-L1 show that the heavy chain is primarily involved with the binding interaction (PDB codes 5XXY, 5X8L). Also, mutating only the heavy chain facilitates deep sequencing because of the shorter region needed to be read by MiSeq (in this case, a region of approximately 316 bp). Error-prone PGR of the heavy chain region was used to randomly generate mutations of the gene, followed by transformation into EBY100 yeast to create a library of up to 78 million Atezolizumab scFv variants. Because the binding interaction and off-rate between the WT scFv and PD-L1 was very strong/slow (KD = 1 .75 nM, kotr = 1 .56 x 10^-4 s^-1), yeast display is better suited under these conditions at decreasing the off-rate rather than directly improving the KD because of the large volumes necessary to maintain a large molar excess of antigen over the scFv on the yeast surface along with the very long incubation times needed to screen for improved Kowhen the off-rate is already very slow. Using the WT scFv as a control during sorting, variants with faster, WT-level, and slower off-rates were isolated from FACS sorts (FIGs. I 6A1-2) and had their plasmids harvested and subjected to multiple rounds of PCR to barcode each group for MiSeq Nano PE250. Notably, the clones selected to have slower off-rates seemed to split into two populations, one more intense than the WT and less intense than the WT (FIGs. I 6B1-4, “HB” panel).

The Atezolizumab dataset consists of sequences derived from mutants in 3 different bins of decreasing off-rates: RH01 , RH02 and RH03. RH01 contains weak binders with a faster off-rate, RH02 contains moderate binders with an off-rate similar to the wild type, and RH03 contains stronger binders with slower off-rates. After filtering for quality, 92,553 unique sequences (550,215 total sequences, since a sequence can appear more than once in a bin or more than once in different bins) were identified, among which 15,004 sequences exhibiting mutations outside the region of the protein targeted for randomization were discarded. There were 15,070, 26,122 and 34,439 unique sequences in RH03, RH02 and RH01 , respectively (Table 10).

Table 10: Read counts by sort category

Each sequence was assigned to the category in which it occurred with the greatest frequency and assigned a weight given by its frequency in the assigned category plus one divided by the total frequency plus three. If the frequency for a given sequence was equal in two categories, it was discarded since it could not be unambiguously assigned. As a result of these filters, 75,631 unique sequences remained. The weight of each sequence corresponds to the posterior probability that it belongs to a given category using a multinomial likelihood and uniform Dirichlet prior. This weighting is important since it enables a model to distinguish between sequences that can and cannot be reliably assigned to a single category. A sequence that appears 5 times in both RH02 and RH03, for example, is clearly less likely to be a strong binder than a sequence that occurs 5 times in RH03 only.

Atezolizumab modeling

The next component of the pipeline is an affinity model which takes as input a representation of a candidate sequence and predicts its sort category (RH01 , RH02 or RH03) to which the sequence should belong. This component requires experimental data specific to the antigen of interest and will need to be retrained with fresh experimental data if a new antigen of interest is selected.

Unlike most classifiers for protein engineering, the affinity model is trained to perform ordinal regression and the last layer of the network outputs a latent score value. A traditional classification model treats the categories as nominal -- they have no particular ordering - and thus the model does not learn any way to rank the categories. This is clearly inappropriate for ordinal sort data. Ordinal regression solves this problem by imposing an ordering on the categories and by using a latent score to determine into which category the sequence should fall. This approach provides a straightforward means to rank sequences and select them for experimental evaluation. In past experiments, it was demonstrated that given protein data with > 3 binding categories, ordinal regression provides improved performance for correctly predicting which sequences will occur in future, more stringent sorts, and in particular outperforms the sequence ranking approach suggested by Liu et al. (Bioinformatics 36 (7), 2126-2133, 2020). Those experiments are difficult to reproduce here for only 3 binding categories but given the past results ordinal regression was used here as well.

The affinity model output score quantifies the extent of the model’s belief that the sequence is likely to be a strong binder - higher scores indicate the model is more certain the sequence should belong to a more stringent sort. The score is next added to M - 1 learned threshold values for M categories, followed by application of the sigmoid to generate an output vector of M - 1 probability values. Each element / of this output vector is the model-assigned probability that the input sequence belongs to a binding category more stringent than i.

Briefly, the affinity model described herein employs an architecture similar to the Bayes by Backprop algorithm, with the difference that this model is adapted to perform ordinal regression. A traditional fully connected neural network learns a specific value for each weight and bias term in each hidden layer during training. A Bayesian neural network, by contrast, treats each weight as a (usually Gaussian) distribution and learns the parameters for each weight’s distribution during training. Instead of fitting using maximum likelihood as is typical for neural networks, the Bayesian architecture is fitted by approximating the posterior probability distribution using a variational method. By sampling over the weight distributions repeatedly, the uncertainty in a prediction can be estimated - or if preferred use the mean of each weight’s Gaussian distribution to generate a point estimate prediction. In addition to generating an estimate of uncertainty, the Bayesian neural network also provides strong regularization on the model parameters by penalizing deviations from the prior.

The uncertainty information provided by the model can clearly be used in one of two ways. Predictions with high associated uncertainty correspond to sequences that lie in relatively unexplored regions of the input space. In practice, only predictions with relatively small associated uncertainty should be selected for experimental evaluation considering the cost and time, and that is the strategy pursued here.

The results generated by a vanilla fully connected neural network trained to perform simple classification were compared with a Bayesian net trained to perform ordinal regression across various encoding types (Table 9). Note that both neural networks compared have the same number of hidden layers and weights. Although the non-Bayesian network offers a modest improvement in performance, the Bayesian network offers additional information about the uncertainty associated with a given prediction that is crucial for analyzing the sort data, since it enables to gauge the relative reliability of the model’s predictions for candidate sequences.

The affinity model can use one-hot encoded sequences as input or another representation. To determine whether the representation generated by the autoencoder is useful as an encoding for the affinity model, the Atezolizumab sequences were encoded using 1 ) the fully trained autoencoder, 2) one-hot encoding, 3) the ProtVec encoding, 4) the UniRep embedding, 5) the ESM-1 b encoding, 6) the AbLang embedding and 7) the AntiBERTy embedding. The encoded datasets were split into a training (80%) and test (20%) portion, and a 5-fold cross-validation was performed on the training portion for each encoding type, using a Bayesian neural network trained to perform ordinal regression (Table 9).

While different encodings may well prove preferable for specific tasks, the autoencoder is the only representation equivalent or superior to one-hot encoding of antibodies for any model type (Table 9). It offers equivalent performance despite a significant reduction in model size and computational expense compared to the one-hot encoding. Note that the autoencoder has roughly 20,000 parameters, while the FAIR- ESM model has 750 million parameters and the AntiBERTy model has 26 million.

The autoencoder therefore provides at least two concrete benefits. First, it reduces model size and complexity significantly compared with simple one-hot encoding, and it is substantially cheaper than language models, since it has a small fraction of the number of parameters. Second, it provides a real-valued representation of each sequence that can be used to cluster sequences or determine similarity between them in later stages of the pipeline. Given these advantages, the autoencoder to one-hot encoding was preferred even though here they achieve similar performance.

In silico directed evolution to select the most promising candidates

The Atezolizumab dataset was encoded using the autoencoder and then the Bayesian neural network (the affinity model) on the full dataset was trained for off-rates for 30 epochs. At this point, the pipeline uses the trained affinity model to search for the sequence space surrounding the training set to generate new candidate antibody sequences that are not present in the training set.

At this stage it is desirable to reduce the size of the search space. If a 1 18 amino acid sequence is mutated in silico, the search space to be covered is impossibly vast. Therefore, it was focused on a smaller subset of positions that contribute significantly to binding affinity - those most frequently mutated positions in the top-scoring sequences, i.e., the most promising sequences considered by the model. For this dataset, it was found that, using the top 500, 1000, 2000 or 4,000 highest scoring sequences to select the top ten most frequently mutated positions, the same set of positions that are used in the subsequent search steps was retrieved. This step in the pipeline is flexible and positions can be selected using other criteria (e.g., only frequently mutated positions present in CDRs).

In the following search steps, on each iteration the model randomly selects a position in the wild-type sequence and randomly mutates it to any amino acid. Note that all mutants generated through this procedure have the same length as the wild-type sequence. The probability of selecting any specific amino acid at a given position is given by the frequency of that amino acid in the training set plus one, divided by the total number of sequences in the training set plus 20. The new candidate sequence is scored by encoding it using the encoder and scoring it using the trained affinity model. The score is assessed using the classic simulated annealing criterion, whereby the candidate is accepted with a probability determined by its score, the score of the last accepted candidate and the temperature. The temperature begins at a high value to allow exploration of a large sequence space and decreases on each iteration, so that the probability of accepting a proposed sequence without an improved score decreases.

Since the average sequence in the training set contains just 7 mutations and the average high-scoring sequence contains even fewer (3 mutations on average), at most positions the most common amino acid is the one present in the wild type. Consequently, the algorithm will tend to heavily sample sequences that are similar to the wild type. Given the stochastic nature of the algorithm, however, it can explore combinations that are not present in the training set and thereby enables finding of new high-scoring sequences. FIGs. 2OA1-2, 20B & 20C illustrate how the best achieved score to date evolves over the course of this optimization.

Ten (10) simulated annealing chains were run; all converged in under 1000 iterations. The accepted candidates from each chain with scores > 90th percentile were harvested and duplicates were removed. The uncertainty was evaluated on the scores by sampling 1000x from the Bayesian neural net for each sequence. The top 50% of sequences with the largest standard deviation on assigned scores were removed. The 50% threshold here is arbitrary and this process yielded 127 sequences.

These sequences were then clustered using median hierarchical clustering (the resulting dendrogram is provided in FIG. 33). The results clearly suggest the selected sequences can be divided into two main subgroups. The marginal distributions (see FIGs. 2OA1-2, 20B and 20C) indicate the two major groups differ primarily at position 77 but also contain subgroups with some other interesting differences. Most selected sequences, for example, exhibit an R98C mutation, but a subset of cluster 1 is unchanged from the wild type at that position.

Cutting the tree at a lower height yields 1 1 subclusters and the two highest scoring sequences (or one if only one sequence was present) in each cluster were selected to yield 21 final candidates. This threshold was selected to yield a manageable number of sequences for experimental evaluation.

Examination of the mutations within the 21 sequences revealed mutations at residues A40, K43, T58, I70, N77, A79, S85, A97, and R98 (in various combinations, see Table 12 and FIG. 21 ).

Table 12. The amino acid sequences of the twenty-one (21 ) mutants were selected for further evaluation

The amino acid sequence of each Full-length IgG scFv sequence of each of the twenty 21 mutants are presented below:

Mutant 1

DIQMTQSPSSLSASVGDRVTITCRASQDVSTAVAWYQQKPGKAPKLLIYSASFLYSGVP

SRFSGSGSGTDFTLTISSLQPEDFATYYCQQYLYHPATFGQGTKVEIKRGGGGSGGGG

SGGGGSGGGGSEVQLVESGGGLVQPGGSLRLSCAASGFTFSDSWIHWVRQAPGQG LEWVAWISPYGGSTYYADSVKGRFTISADTSKNTTYLQMNSLRAEDTAVYYCVCRHWP GGFDYWGQGTLVTVSS (SEQ ID NO: 22)

Mutant 2

DIQMTQSPSSLSASVGDRVTITCRASQDVSTAVAWYQQKPGKAPKLLIYSASFLYSGVP SRFSGSGSGTDFTLTISSLQPEDFATYYCQQYLYHPATFGQGTKVEIKRGGGGSGGGG SGGGGSGGGGSEVQLVESGGGLVQPGGSLRLSCAASGFTFSDSWIHWVRQAPGQG

LEWVAWISPYGGSTYYADSVKGRFTASADTSKNTTYLQMNSLRAEDTAVYYCVCRHW PGGFDYWGQGTLVTVSS (SEQ ID NO: 23)

Mutant 3

DIQMTQSPSSLSASVGDRVTITCRASQDVSTAVAWYQQKPGKAPKLLIYSASFLYSGVP SRFSGSGSGTDFTLTISSLQPEDFATYYCQQYLYHPATFGQGTKVEIKRGGGGSGGGG SGGGGSGGGGSEVQLVESGGGLVQPGGSLRLSCAASGFTFSDSWIHWVRQAPGRG

LEWVAWISPYGGSTYYADSVKGRFTTSADTSKNTTYLQMNSLRAEDTAVYYCVRRHW PGGFDYWGQGTLVTVSS (SEQ ID NO: 24)

Mutant 4

DIQMTQSPSSLSASVGDRVTITCRASQDVSTAVAWYQQKPGKAPKLLIYSASFLYSGVP SRFSGSGSGTDFTLTISSLQPEDFATYYCQQYLYHPATFGQGTKVEIKRGGGGSGGGG SGGGGSGGGGSEVQLVESGGGLVQPGGSLRLSCAASGFTFSDSWIHWVRQAPGKG

LEWVAWISPYGGSTYYADSVKGRFTASADTSKNTTYLQMNSLRAEDTAVYYCVRRHW PGGFDYWGQGTLVTVSS (SEQ ID NO: 25) Mutant 5

DIQMTQSPSSLSASVGDRVTITCRASQDVSTAVAWYQQKPGKAPKLLIYSASFLYSGVP

SRFSGSGSGTDFTLTISSLQPEDFATYYCQQYLYHPATFGQGTKVEIKRGGGGSGGGG

SGGGGSGGGGSEVQLVESGGGLVQPGGSLRLSCAASGFTFSDSWIHWVRQAPGEG

LEWVAWISPYGGSAYYADSVKGRFTTSADTSKNTTYLQMNSLRAEDTAVYYCVCRHW

PGGFDYWGQGTLVTVSS (SEQ ID NO: 26)

Mutant 6

DIQMTQSPSSLSASVGDRVTITCRASQDVSTAVAWYQQKPGKAPKLLIYSASFLYSGVP

SRFSGSGSGTDFTLTISSLQPEDFATYYCQQYLYHPATFGQGTKVEIKRGGGGSGGGG

SGGGGSGGGGSEVQLVESGGGLVQPGGSLRLSCAASGFTFSDSWIHWVRQAPGRG

LEWVAWISPYGGSAYYADSVKGRFTTSADTSKNTTYLQMNSLRAEDTAVYYCVCRHW

PGGFDYWGQGTLVTVSS (SEQ ID NO: 27)

Mutant 7

DIQMTQSPSSLSASVGDRVTITCRASQDVSTAVAWYQQKPGKAPKLLIYSASFLYSGVP

SRFSGSGSGTDFTLTISSLQPEDFATYYCQQYLYHPATFGQGTKVEIKRGGGGSGGGG

SGGGGSGGGGSEVQLVESGGGLVQPGGSLRLSCAASGFTFSDSWIHWVRQAPGKG

LEWVAWISPYGGSTYYADSVKGRFTASADTSKNTTYLQMNSLRAEDTAVYYCVCRHW

PGGFDYWGQGTLVTVSS (SEQ ID NO: 28)

Mutant 8

DIQMTQSPSSLSASVGDRVTITCRASQDVSTAVAWYQQKPGKAPKLLIYSASFLYSGVP

SRFSGSGSGTDFTLTISSLQPEDFATYYCQQYLYHPATFGQGTKVEIKRGGGGSGGGG

SGGGGSGGGGSEVQLVESGGGLVQPGGSLRLSCAASGFTFSDSWIHWVRQAPGKG

LEWVAWISPYGGSAYYADSVKGRFTTSADTSKNTTYLQMNSLRAEDTAVYYCVCRHW

PGGFDYWGQGTLVTVSS (SEQ ID NO: 29)

Mutant 9

DIQMTQSPSSLSASVGDRVTITCRASQDVSTAVAWYQQKPGKAPKLLIYSASFLYSGVP

SRFSGSGSGTDFTLTISSLQPEDFATYYCQQYLYHPATFGQGTKVEIKRGGGGSGGGG

SGGGGSGGGGSEVQLVESGGGLVQPGGSLRLSCAASGFTFSDSWIHWVRQAPGKG

LEWVAWISPYGGSAYYADSVKGRFTISADTSKSTVYLQMNSLRAEDTAVYYCVCRHWP

GGFDYWGQGTLVTVSS (SEQ ID NO: 30)

Mutant 10

DIQMTQSPSSLSASVGDRVTITCRASQDVSTAVAWYQQKPGKAPKLLIYSASFLYSGVP

SRFSGSGSGTDFTLTISSLQPEDFATYYCQQYLYHPATFGQGTKVEIKRGGGGSGGGG

SGGGGSGGGGSEVQLVESGGGLVQPGGSLRLSCAASGFTFSDSWIHWVRQAPGKG

LEWVAWISPYGGSAYYADSVKGRFTISADTSKDTVYLQMNSLRAEDTAVYYCVCRHWP

GGFDYWGQGTLVTVSS (SEQ ID NO: 31 )

Mutant 11 DIQMTQSPSSLSASVGDRVTITCRASQDVSTAVAWYQQKPGKAPKLLIYSASFLYSGVP

SRFSGSGSGTDFTLTISSLQPEDFATYYCQQYLYHPATFGQGTKVEIKRGGGGSGGGG

SGGGGSGGGGSEVQLVESGGGLVQPGGSLRLSCAASGFTFSDSWIHWVRQTPGRG

LEWVAWISPYGGSAYYADSVKGRFTISADTSKDTTYLQMNSLRAEDTAVYYCACRHWP

GGFDYWGQGTLVTVSS (SEQ ID NO: 32)

Mutant 12

DIQMTQSPSSLSASVGDRVTITCRASQDVSTAVAWYQQKPGKAPKLLIYSASFLYSGVP

SRFSGSGSGTDFTLTISSLQPEDFATYYCQQYLYHPATFGQGTKVEIKRGGGGSGGGG

SGGGGSGGGGSEVQLVESGGGLVQPGGSLRLSCAASGFTFSDSWIHWVRQTPGEGL

EWVAWISPYGGSAYYADSVKGRFTISADTSKDTTYLQMNSLRAEDTAVYYCTCRHWP

GGFDYWGQGTLVTVSS (SEQ ID NO: 33)

Mutant 13

DIQMTQSPSSLSASVGDRVTITCRASQDVSTAVAWYQQKPGKAPKLLIYSASFLYSGVP

SRFSGSGSGTDFTLTISSLQPEDFATYYCQQYLYHPATFGQGTKVEIKRGGGGSGGGG

SGGGGSGGGGSEVQLVESGGGLVQPGGSLRLSCAASGFTFSDSWIHWVRQTPGRG

LEWVAWISPYGGSAYYADSVKGRFTASADTSKDTTYLQMNSLRAEDTAVYYCVCRHW

PGGFDYWGQGTLVTVSS (SEQ ID NO: 34)

Mutant 14

DIQMTQSPSSLSASVGDRVTITCRASQDVSTAVAWYQQKPGKAPKLLIYSASFLYSGVP

SRFSGSGSGTDFTLTISSLQPEDFATYYCQQYLYHPATFGQGTKVEIKRGGGGSGGGG

SGGGGSGGGGSEVQLVESGGGLVQPGGSLRLSCAASGFTFSDSWIHWVRQVPGRG

LEWVAWISPYGGSAYYADSVKGRFTASADTSKDTTYLQMNSLRAEDTAVYYCVCRHW

PGGFDYWGQGTLVTVSS (SEQ ID NO: 35)

Mutant 15

DIQMTQSPSSLSASVGDRVTITCRASQDVSTAVAWYQQKPGKAPKLLIYSASFLYSGVP

SRFSGSGSGTDFTLTISSLQPEDFATYYCQQYLYHPATFGQGTKVEIKRGGGGSGGGG

SGGGGSGGGGSEVQLVESGGGLVQPGGSLRLSCAASGFTFSDSWIHWVRQVPGKG

LEWVAWISPYGGSAYYADSVKGRFTASADTSKDTTYLQMNTLRAEDTAVYYCVCRHW

PGGFDYWGQGTLVTVSS (SEQ ID NO: 36)

Mutant 16

DIQMTQSPSSLSASVGDRVTITCRASQDVSTAVAWYQQKPGKAPKLLIYSASFLYSGVP

SRFSGSGSGTDFTLTISSLQPEDFATYYCQQYLYHPATFGQGTKVEIKRGGGGSGGGG

SGGGGSGGGGSEVQLVESGGGLVQPGGSLRLSCAASGFTFSDSWIHWVRQVPGKG

LEWVAWISPYGGSAYYADSVKGRFTASADTSKDTIYLQMNSLRAEDTAVYYCVCRHWP

GGFDYWGQGTLVTVSS (SEQ ID NO: 37)

Mutant 17 DIQMTQSPSSLSASVGDRVTITCRASQDVSTAVAWYQQKPGKAPKLLIYSASFLYSGVP SRFSGSGSGTDFTLTISSLQPEDFATYYCQQYLYHPATFGQGTKVEIKRGGGGSGGGG SGGGGSGGGGSEVQLVESGGGLVQPGGSLRLSCAASGFTFSDSWIHWVRQTPGKGL EWVAWISPYGGSAYYADSVKGRFTTSADTSKDTTYLQMNSLRAEDTAVYYCVCRHWP

GGFDYWGQGTLVTVSS (SEQ ID NO: 38)

Mutant 18

DIQMTQSPSSLSASVGDRVTITCRASQDVSTAVAWYQQKPGKAPKLLIYSASFLYSGVP SRFSGSGSGTDFTLTISSLQPEDFATYYCQQYLYHPATFGQGTKVEIKRGGGGSGGGG SGGGGSGGGGSEVQLVESGGGLVQPGGSLRLSCAASGFTFSDSWIHWVRQVPGKG

LEWVAWISPYGGSAYYADSVKGRFTTSADTSKDTTYLQMNSLRAEDTAVYYCVCRHW PGGFDYWGQGTLVTVSS (SEQ ID NO: 39)

Mutant 19

DIQMTQSPSSLSASVGDRVTITCRASQDVSTAVAWYQQKPGKAPKLLIYSASFLYSGVP SRFSGSGSGTDFTLTISSLQPEDFATYYCQQYLYHPATFGQGTKVEIKRGGGGSGGGG SGGGGSGGGGSEVQLVESGGGLVQPGGSLRLSCAASGFTFSDSWIHWVRQAPGKG

LEWVAWISPYGGSTYYADSVKGRFTVSADTSKDTTYLQMNSLRAEDTAVYYCVRRHW PGGFDYWGQGTLVTVSS (SEQ ID NO: 40)

Mutant 20

DIQMTQSPSSLSASVGDRVTITCRASQDVSTAVAWYQQKPGKAPKLLIYSASFLYSGVP SRFSGSGSGTDFTLTISSLQPEDFATYYCQQYLYHPATFGQGTKVEIKRGGGGSGGGG SGGGGSGGGGSEVQLVESGGGLVQPGGSLRLSCAASGFTFSDSWIHWVRQAPGKG

LEWVAWISPYGGSTYYADSVKGRFTASADTSKDTTYLQMNSLRAEDTAVYYCVRRHW PGGFDYWGQGTLVTVSS (SEQ ID NO: 41 )

Mutant 21

DIQMTQSPSSLSASVGDRVTITCRASQDVSTAVAWYQQKPGKAPKLLIYSASFLYSGVP SRFSGSGSGTDFTLTISSLQPEDFATYYCQQYLYHPATFGQGTKVEIKRGGGGSGGGG SGGGGSGGGGSEVQLVESGGGLVQPGGSLRLSCAASGFTFSDSWIHWVRQAPGKG

LEWVAWISPYGGSAYYADSVKGRFTTSADTSKDTAYLQMNFLRAEDTAVYYCACRHW PGGFDYWGQGTLVTVSS (SEQ ID NO: 42)

Examination of the existing structures of the Atezolizumab Fab and PD-L1 (PDB codes 5XXY and 5X8L) reveals that none of the mutated residues contact PD-L1 in the WT sequence (FIG. 21 ). This suggests these mutations improve the binding affinity through alterations of the heavy chain conformation rather than directly improving binding contacts. Another possibility is that these mutations stabilize the conformation rather than altering it. It was demonstrated, for example, that in some cases stabilizing a specific antibody conformation can be beneficial for affinity.

To explore these possibilities, predicted structures for the top-scoring mutants from the twenty-one (21 ) sequences were generated using IgFold without antigen. The mutations present in high-scoring sequences shift the conformation of both antigen- contacting regions of the protein, especially the CDR H3. FIGs. 30A-30B & FIG. 31 show the conformational change of R98C (a mutation that presents in many of the top-scoring sequences) and its surrounding residues as well as the formation of an apparent hydrophobic pocket in the top mutant compared to the wildtype. Note that the mutations described here are not in direct contact with the antigen. These predicted structures suggest that the mutations selected by this pipeline function through changing the conformation of the contact regions rather than by directly forming new contacts with the antigen by themselves.

The trained affinity model was also used to computationally assess the importance of individual mutations. A) how much a mutation would contribute to binding when only itself is introduced individually to the WT without other mutations and B) how much binding affinity would change if a single mutation is removed from a beneficial mutation combination were scored. The full results are presented in Tables 13 and 14.

Table 13. Change in affinity model-assigned score for each mutation found in the final 21 selected sequences, if that mutation were introduced into the wild type with no other modifications

Table 14. Change in affinity model assigned score (positive = more likely to be a strong binder) if removing the specified mutation from the specified sequence selected for experimental evaluation til SH

The results do suggest that some mutations may be more important than others, and also suggest that 12 of the 19 mutations introduced are predicted to be beneficial even in isolation, i.e., to increase the likelihood that the sequence in question is a tight binder. The remainder is very slightly detrimental in isolation or causes no change, although they are predicted to have beneficial impact in certain contexts. A[79] I or A[79]T, R[98]C, l[70] A, A[97]V and T[58]A are predicted by the model to be mutations with large beneficial impacts even in isolation.

Interestingly, the model is able to predict synergy between individual mutations. For example, K43Q is predicted to have an almost negligible beneficial impact if introduced into the wild type in isolation, but is predicted to have a much larger negative impact on the score if it is removed from the mutant K[43]Q, A[79]T, A[97]V, R[98]C. A[97]T is very slightly detrimental in isolation, but is predicted to be beneficial for A[40]T, K[43]E, T[58]A, N[77]D, A[79]T, A[97]T, R[98]C. While A[79]T, for example, is predicted to be beneficial no matter in what context it is introduced, the size of that impact on the score varies as much as twofold depending on the other mutations present.

Validation of the predicted tight binding antibodies The twenty-one (21 ) mutants (all in the heavy chain) generated by the model were purchased as geneblocks and fused with the WT light chain sequence by PCR, followed by transformation into EBY100 yeasts. The small library of the twenty-one (21 ) mutants was screened in a similar manner as the naive library, except a longer competition time was used in the final screen (39 hrs) in an attempt to separate the mutant with the slowest off-rate from the rest of the 21 mutants (FIGs. 25A-25D). This process did not result in a clearly separated mutant and sequencing of random clones after the most stringent sort showed that out of 17 random sequences, 12 separate mutants (occurring in similar frequencies) were found to be present. This suggests that a significant proportion of the 21 -member pool of mutants had significantly slower off-rates than the WT (which exhibited significantly reduced fluorescence when compared with the pool of mutants in the final, most stringent sort (FIGs. 25C-25D, lower panels). It also appeared that essentially all the 21 mutants selected by the model had slower off-rates than the WT upon examination of the original/unsorted 21 -member library after 8-hours of off-rate competition at RT (FIGs. 25A-25B, top panels).

From the various mutants selected from the 21 -member library after the most stringent sort, the I70A/A79T/A97V mutant (named Mutant 4, see Table 1 1 for sequence) was selected for further characterization.

Table 11. WT Atezolizumab, Mutant 4, Durvalumab, and Avelumab scFv geneblocks, cloning primers, primers used to generate the mutant library and DNA library for next- generation sequencing

The amino acid sequence of the Full-Length IgG Heavy Chain/Constant Region of Mutant 4 is presented below: MGWSCIILFLVATATGVHSEVQLVESGGGLVQPGGSLRLSCAASGFTFSDSWIHWVRQ APGKGLEWVAWISPYGGSTYYADSVKGRFTASADTSKNTTYLQMNSLRAEDTAVYYCV RRHWPGGFDYWGQGTLVTVSSASTKGPSVFPLAPSSKSTSGGTAALGCLVKDYFPEP VTVSWNSGALTSGVHTFPAVLQSSGLYSLSSVVTVPSSSLGTQTYICNVNHKPSNTKV DKKVEPKSCDKTHTCPPCPAPELLGGPSVFLFPPKPKDTLMISRTPEVTCVVVDVSHE DPEVKFNWYVDGVEVHNAKTKPREEQYASTYRVVSVLTVLHQDWLNGKEYKCKVSNK ALPAPIEKTISKAKGQPREPQVYTLPPSREEMTKNQVSLTCLVKGFYPSDIAVEWESNG QPENNYKTTPPVLDSDGSFFLYSKLTVDKSRWQQGNVFSCSVMHEALHNHYTQKSLS LSPGK (SEQ ID NO: 67)

The amino acid sequence of the Full-Length IgG Light Chain/Constant Region of Mutant 4 is presented below:

MGWSCIILFLVATATGVHSDIQMTQSPSSLSASVGDRVTITCRASQDVSTAVAWYQQKP GKAPKLLIYSASFLYSGVPSRFSGSGSGTDFTLTISSLQPEDFATYYCQQYLYHPATFGQ GTKVEIKRTVAAPSVFIFPPSDEQLKSGTASVVCLLNNFYPREAKVQWKVDNALQSGN SQESVTEQDSKDSTYSLSSTLTLSKADYEKHKVYACEVTHQGLSSPVTKSFNRGEC (SEQ ID NO: 68)

This mutant, unlike the other mutants sequenced after the stringent sort, did not have cysteine introduced into its sequence so was more appealing because additional disulfide bonds would not be formed by its mutations. To characterize the improvement in off-rate/binding towards PD-L1 of the isolated mutant, the KD and kotr values of both the WT and Mutant 4 scFv towards PD-L1 were determined using the yeast display format. First, the apparent off-rates at RT on the yeast surface was determined (FIG. 22A), and it was found that Mutant 4 has a kotr approximately 10-fold slower than WT Atezolizumab (6.3 x 10^-5 s^-1, half-life 3.04 hrs for the WT vs. 6.5 x 10’⁶ s^-1, half-life 29.8 hrs for Mutant 4). Then, the kott values of Mutant 4 were compared to the two other FDA approved anti- PD-L1 mAbs (Durvalumab and Avelumab, converted to the scFv format on the yeast surface) (FIG. 22B). Mutant 4 displayed a substantially slower off-rate compared to both WT Atezolizumab and Durvalumab and a noticeably slower value than Avelumab. The order of off-rates of the FDA-approved mAb PD-L1 scFvs (Durvalumab < Atezolizumab < Avelumab) is consistent with a previous study of these scFv values by a surface plasmon resonance (SPR) binding assay.

Next, the binding affinity (KD) values of the WT and Mutant 4 scFv were determined by titrating PD-L1 against each scFv on the yeast surface. The Mutant 4 KD was about 17-fold tighter (on average) than the WT (5.3 pM vs. 92 pM, FIG. 22C). It should be noted that determining the KD on the yeast surface at such low concentrations of antigen is difficult, given the excessive volumes necessary to maintain a significant molar excess of antigen to scFv on the yeast surface at lower antigen concentrations⁴⁹, which also makes performing replicates for each data point difficult. However, three separate KD measurements gave significant (>10-fold) improvements in the binding affinity (FIG. 22C).

In order to assess if the improved binding affinity and off-rate observed on the yeast surface can be repeated after purifying each scFv and reversing the orientation of the binding assay (instead of scFv fused to the yeast surface and binding to soluble PD- L1 , PD-L1 is immobilized and allowed to bind to soluble scFv), a BLI (Bio-Layer Interferometry) binding assay was performed, and it was found that Mutant 4 has a much slower off rate than the WT, consistent with the yeast surface display results. The dissociation of the mutant from PD-L1 in this assay was minimal even after 2 hours at RT while the WT was more significant over about a 60-minute dissociation time. However, BLI gave two KD values for the WT (219pM for the major species and 5.93nM for the minor species) while the Mutant 4 KD values were immeasurable because of its very slow off-rate, leading to an inability to fit/determine the kinetic constant values (see FIGs. 26Ai- 2 and 26B). The minor species may be caused by a minor population of partially unfolded scFv in the protein preparation. This observation leads to examine whether the scFv and PD-L1 were monomeric or formed oligomers, so this was tested by mass photometry (a method that can image oligomeric distributions of a protein in solution). Imaging PD-L1 (residues Met1 -T239, expected MW 35-38 kDa, range in MW due to glycosylation) revealed a predominantly monomeric protein at the concentrations tested (7.4 nM, 18.5 nM, and 74 nM) (FIGs. 27A-27G), consistent with previous studies of human PD-L1 which found it to be monomeric. These data support that the KD measured on the yeast surface is between monomeric PD-L1 and scFv immobilized to the yeast surface. Imaging of the WT or Mutant 4 scFv in solution revealed a mixed population of monomer, dimer, trimer, and possibly higher MW oligomers at 100 nM each scFv. The WT and Mutant 4 scFv were purified twice (using the same exact method) and tested twice by mass photometry. For the WT scFv, in both preparations it was a mixture of monomer and dimer in solution while for an unknown reason the mutant was either predominantly monomeric (1 st batch) or a mixture of monomer/dimer/trimer (batch 2) (FIGs. 27A-27G). Furthermore, it was noticed that the mutant scFv seemed more prone to aggregation than the WT scFv. Taken together, the oligomerization and aggregation of scFv in solution makes binding assays like SPR (surface plasmon resonance) or BLI, where scFv is in solution and PD-L1 is immobilized, inappropriate for KD determination because of the difficulty of knowing the soluble monomeric concentration of scFv.

Evaluation on the Mason et al. (Nat. Biomed. Eng., 5 (6), 600-612, 2021) dataset

To ensure the generality of RESP, it was next tested on the data from the Mason et al. study. Starting from a nonbinding mutant of trastuzumab, they successfully recovered antibodies with similar affinity to the original wild type trastuzumab. Using the same training and test sets, their sequence data were encoded using the autoencoder and a variational Bayesian network with the same structure used for the data was trained to predict binding and nonbinding status.

This model achieves the same accuracy as theirs (Matthews correlation coefficient 0.68, AUC-ROC 0.91 ) for their test set and correctly predicts the sequences they experimentally tested to be binders (see FIGs. 32A-32C). This model likewise assigns high scores to most of these sequences. Trastuzumab, for example, is assigned a higher score than 94% of the sequences in the training set, indicating a high confidence that it is a tight binder. The distribution of scores for the training set and for the experimentally evaluated sequences appear in FIGs. 32A-32C.

Unlike the Mason et al. pipeline, however, this pipeline discussed in this Example provides additional uncertainty information not available from their model. The affinity model’s uncertainty about test set predictions which turn out to be incorrect is significantly higher than its uncertainty about test set predictions which are correct (p < 1 e-23, two- sided Mann-Whitney U test). Consequently, uncertainty can be used to assist in determining whether a prediction is likely to be reliable, which is not possible with the pipeline developed by Mason et al., and no need to sacrifice accuracy to obtain this advantage. DISCUSSION

A pipeline for efficiently developing antibodies is provided with the following innovations. First, a simple learned encoding specific to antibodies is developed. Such encoding contains not only the information in the original sequence in easily recoverable form but also additional encoded information describing key sequence features that differentiate human antibodies from the surrounding sequence space. It was demonstrated that when training a model on antigen-specific experimental data to capture trends in binding affinity, the same model is more accurate if the input sequences are encoded using the autoencoder-generated representation than if the input sequences are instead encoded using popular state-of-the-art encodings like UniRep, ESM-1 b, AntiBertY and AbLang. Remarkably, for this particular task, using all of the models considered, UniRep, ESM-1 b, AntiBertY and AbLang exhibit performance inferior to one- hot encoding. This is consistent with results reported by Makowski et al. (Nat. Commun. 13, 3788, 2022), who found that UniRep or physicochemical properties did not improve performance for antibody affinity prediction compared with simple one-hot encoding.

Next, the training set was fitted using models designed to provide both straightforward and easy to interpret sequence ranking coupled with quantitation of uncertainty. It was shown that the distribution of sequence reads across categories can be incorporated into model fitting as a datapoint weight in a principled and straightforward way. The Bayesian ordinal regression model yields an estimate of the predictive posterior, thereby providing additional information not available from traditional deep learning classifiers, whose predicted probability distribution across categories does not indicate the reliability of a given prediction.

Further, an algorithm was designed to explore the sequence space spanned by the training set. By estimating the reliability of each prediction and by restricting the search to the space spanned by the training set, the time and expense wasted on evaluating poor candidate sequences were minimized.

Moreover, the power of the pipeline disclosed in this Example was also validated. By training the model on a single large library, it was able to select a mutant with an off- rate/binding affinity improvement of 10-fold/17-fold. This is consistent with the previous results for protein engineering of CBX1 , demonstrating a similar strategy yielded improvement in binding affinity equivalent to that achieved by a much lengthier directed evolution process. The Mutant 4 identified and selected here could be useful in cancer therapy as an scFv as was a previously reported high affinity anti-PD-L1 protein (a mutant form of PD-1 , which bound PD-L1 with a KD of 1 10 pM). Like PD-1 , the Mutant 4 scFv is significantly smaller than a monoclonal antibody (30 kDa vs. 150 kDa) and so could possibly be more effective at tumor penetration.

This pipeline disclosed in this Example may also have several important advantages over purely experimental approaches. Phage and yeast display only permit selection of a small population, not of single clones, so that additional experiments (such as ELISA or yeast KD measurements) are needed to assess the clones having the tightest binding affinity. Only strong binders present in the original library can be identified via these techniques, so that often very large libraries and/or multiple libraries are used to maximize coverage of sequence space. It often happens that the best binders identified through this process still do not possess sufficiently strong affinity, so that the desired affinity must often be achieved through so-called affinity maturation. In this process, random mutations are introduced at selected sites and the resulting focused library undergoes further rounds of screening and experimental evaluation. Notably, the antibodies generated through this process are not guaranteed to possess other desirable characteristics like solubility or stability.

Take for example the process by which the Atezolizumab antibody itself was originally discovered (as described in US Patent US8217149B2). Four rounds of panning on a phage display library screened against the PDL1 target were first used to retrieve 96 enriched clones. Two sets of positions were then selected to construct two further libraries for affinity maturation to improve binding. These in turn were used to conduct one plate sort followed by five or six rounds of solution sorting. Finally, enriched clones from the last sort were run through high-throughput ELISA screening to find the best candidate. It is notable that this procedure, while reliable, is expensive in time and cost. It is very common that antibody engineering requires construction of multiple libraries and high- throughput ELISA of thousands of clones; for a couple of other examples involving antibodies now in clinical trials, see US Patent 20180086848A1 and 8313746B2 among others.

In contrast, in the approach disclosed herein, it was able to select a tight-binding antibody after constructing only a single large library followed by FACS-based sorting for off-rate and binning. This approach does not require any high-throughput KD determinations or subsequent focused library screens. Thereby, the time needed to construct multiple libraries, which is considerable, and the time and expense needed for high-throughput ELISA screens/Ko measurements, are all eliminated. Instead, this approach only requires some computational steps which are easy to implement and run quickly on a single computer equipped with a GPU. The Bayesian neural network needs to be retrained any time a new antigen is selected on data specific to that antigen and acquired as described in this Example. The autoencoder, by contrast, can be reused and need not be retrained. Moreover, the approach disclosed in this Example identified tight binders that are not present in the original library, unlike traditional phage and yeast display, in which only sequences present in the library can be screened.

It is also important to note that this approach can easily be modified to introduce in silica screens for stability, solubility, and other desired properties, unlike purely experimental techniques which cannot easily optimize for these other properties simultaneously. It is straightforward to add additional filters to the search step of the pipeline - i.e., to reject candidates suggested by the modified simulated annealing algorithm if predicted solubility or immunogenicity is poor. In a purely experimental approach, by contrast, these properties must be optimized separately at considerable additional expense. The ability to achieve improved affinity while easily introducing additional filters as desired is a key advantage for a machine learning-assisted approach.

There are other strategies for computer-assisted antibody design described in recent literature, however, the approach disclosed herein compares favorably and offers several advantages. For example, Mason et al. achieve a 3x improvement in affinity over the wild-type trastuzumab, and most of the mutants selected by their algorithm as promising in fact show weaker affinity for the target, whereas a 17-fold improvement in affinity was achieved with the approach disclosed herein. Mason et al. performed multiple rounds of mutagenesis and library generation, including a step where positions for further mutagenesis were selected using rational design based on data from a previous selection step. They require additional effort to constrain the search space. In contrast, it was able to constrain in silico using the modified simulated annealing algorithm disclosed herein. Unlike the CNN they propose, the Bayesian neural network disclosed herein provides uncertainty estimates correlated with the accuracy of a prediction so that predictions which are likely to be unreliable can be eliminated.

A rational design approach (as opposed to a machine learning based approach) was described. Their rational design component, however, can only be as accurate as the predictions of the Rosetta modeling software, which results in a low rate of correct predictions, nor is it possible to determine which predictions are most likely to be reliable. They sought for example to improve on the binding of an antibody called G6 to its target, VEGF, using a computational search procedure requiring approximately 250 cpu-days, which is orders of magnitude more expensive than the computational procedure adopted and disclosed herein. Of the mutants selected by this AbLift procedure for experimental evaluation, 60% were worse than the wild type, and most of the remainder were only comparable. Only one of the designs suggested by their approach actually improved the KD, resulting in a fivefold improvement, which is a smaller improvement than the improvement disclosed herein. It is worth noting that most or all the twenty-one (21 ) mutants selected by the model disclosed herein showed slower off-rates (thus likely higher binding affinity) than the WT.

Gaussian process was used to model trends in binding affinity as a function of input sequence. These authors did not experimentally validate their approach, instead using affinity predictions from the Absolut! Software to determine whether a sequence was a strong binder; the Absolut! Software is itself based on docking-generated structures and affinities, so that it is not clear how closely it tracks experimental data. Importantly, their approach suffers from some of the well-known limitations of exact Gaussian processes. Those models exhibit O(N³) scaling in the number of training points (or O(N²) in some more efficient modern implementations) and are thus completely infeasible for datasets larger than 5-10,000 sequences or so. In contrast, the variational Bayesian network-based approach disclosed herein does not suffer from any such limitation. While the approach disclosed herein was validated using Atezolizumab as a starting point and PDL1 as a target, there is nothing in this approach which is specific to the selected target, and thus this approach can easily be adapted to other targets and problems. The pipeline described herein can easily be modified to optimize only a single CDR or a subset of the available positions, and to incorporate other models that rank candidate sequences for other desired properties in addition to affinity. It is anticipated that this pipeline and modified versions of it may therefore prove to be a useful tool for accelerated antibody discovery and development.

Claims

What is claimed is:

1. A method for designing and engineering a protein with a desired property of interest, comprising steps of directed evolution, sequencing, and machine learning analysis.

2. The method of claim 1 , wherein the steps of directed evolution, sequencing, and machine learning analysis comprise using a trained autoencoder computer model to encode directed evolution sequences.

3. The method of claim 2, wherein the steps of directed evolution, sequencing, and machine learning analysis further comprise using a trained variational Bayesian neural network to perform ordinal regression on the encoded directed evolution sequences to quantify the desired property of interest.

4. The method of claim 2, wherein the steps of directed evolution, sequencing, and machine learning analysis further comprise using an approximated Gaussian process (GP) computer model that utilizes random feature-approximated kernels to perform the ordinal regression on the encoded directed evolution sequences.

5. The method of claim 4, wherein the GP computer model reduces a number of computations that need to be performed when utilizing the random feature-approximated kernels by utilizing a method of preconditioned conjugate gradients.

6. The method of claim 3, wherein in silico mutagenesis is performed using a simulated annealing strategy to explore the directed evolution sequences that are not present in mutation libraries and assess their property of interest.

7. The method of claim 1 , wherein the protein is selected from a group consisting of a therapeutics protein, a protein domain, an enzyme, and a protein for enhancing efficiency of CRISPR and enzyme activities.

8. The method of claim 7, wherein the therapeutic protein is an antibody.

9. The method of claim 7, wherein the protein domain is used for an imaging probe.

10. A method for designing and engineering an antibody comprising desired mutations and with an improved binding affinity to or off-rates towards a targeted antigen, comprising: a) establishing an autoencoder model based on encoded mutant sequence data of a known antibody binding to the targeted antigen; b) developing a Bayesian ordinal regression model that scores each mutant sequence on a probability of being a strong binder to the targeted antigen; and c) selecting the antibody comprising the desired mutations using a modified simulated annealing algorithm for testing binding affinity to or off-rates towards the targeted antigen.

1 1 . The method of claim 10, wherein the targeted antigen is PD-L1 .

12. The method of claim 10, wherein the known antibody is Atezolizumab.

13. The method of claim 12, wherein the antibody comprising the desired mutations comprises a scFv amino acid sequence as set forth in SEQ ID Nos:1 -21 .

14. The method of claim 13, wherein the antibody comprising the desired mutations of Mutant 4 which comprises an amino acid sequence of SEQ ID NO:4, wherein a full-length IgG scFv sequence of mutant 4 comprises an amino acid sequence as set forth in SEQ ID NO: 25.

15. A method for improving an antibody-antigen interaction comprising the method of any one of claims 1 -14.

16. An antibody, or a composition thereof, comprising a mutant sequence of a known antibody for a targeted antigen, wherein said antibody is engineered by the method of any one of claims 1 -14.

17. The antibody of claim 16, wherein the antibody comprising the desired mutations comprises a scFv amino acid sequence as set forth in SEQ ID Nos:1 -21 .

18. The antibody of claim 16, wherein the antibody comprising the desired mutations is Mutant 4 which comprises an amino acid sequence of SEQ ID NO:4, wherein a full- length IgG scFv sequence of mutant 4 comprises an amino acid sequence as set forth in SEQ ID NO: 25.

19. The antibody of claim 17 or claim 18, wherein the antibody comprises a mutant sequence of Atezolizumab.

20. The antibody of any one of claims 16-19, wherein the antibody is used for treating a PD-L1 positive tumor.

21. A method for treating a PD-L1 positive tumor comprising administering a subject of interest an effective amount of antibody, or the composition thereof, of any one of claims 17-19.

22. Use of the antibody, or a composition thereof, of any one of claims 16-20 in cancer immunotherapy CAR-T cell therapy or enhancing efficiency of CRISPR and enzyme activities.