Intelligent protein design and engineering

Abstract

The pipeline for protein design and engineering consists of directed evolution, sequencing, and machine learning analysis. This pipeline can explore directed evolution sequences not present in libraries and can be used to facilitate the discovery of therapeutic proteins and imaging probes and / or enhance CRISPR efficiency and enzymatic activity. This pipeline discovers scFv mutants that exhibit tight binding and / or slow off-rates to PD-L1. These scFvs can be used as potent antibodies against PD-L1 in cancer immunotherapy and / or CAR-T cell therapy.

Description

[Technical Field]

[0001] CROSS-REFERENCE TO RELATED APPLICATIONS This application claims the benefit of U.S. Provisional Patent Application No. 63 / 431,556, filed December 9, 2022, and U.S. Provisional Patent Application No. 63 / 432,836, filed December 15, 2023, respectively, the entire contents of each of the above-referenced applications being incorporated herein by reference.

[0002] STATEMENT REGARDING FEDERALLY SPONSORED RESEARCH OR DEVELOPMENT This invention was made with government support under grants R01GM111941 and R21AI158114, respectively, awarded by the National Institutes of Health (NIH). The U.S. Federal Government has certain rights in this invention. [Background technology]

[0003] Identifying proteins and small molecules with desired properties is a very important task for the pharmaceutical and chemical industries. In recent years, machine learning-assisted methods have become more popular. Neural network / deep learning models are widely used due to their flexibility, their ability to learn complex relationships from large datasets, and the development of libraries such as PyTorch, which are easy to implement.

[0004] Despite the success of deep learning, it suffers from several limitations. First, many deep learning architectures do not quantify their uncertainty about their predictions, tending to be "overconfident" when extrapolating and are "black boxes," making it very difficult to determine why a model makes a particular prediction about a particular input. Second, surprisingly, deep learning models have repeatedly been shown to lack robustness against so-called "adversarial attacks," which are small perturbations to their inputs. For example, the addition of subtle noise imperceptible to the human eye can reliably cause a deep learning model 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 costs. Many state-of-the-art models use tens or hundreds of billions of learned parameters. In many cases, it is possible to remove a large fraction of these parameters without compromising performance, suggesting that more efficient approaches may exist. For example, it was reported that 97% of the parameters of the BERT model could be truncated in natural language processing and still achieve comparable performance.

[0005] An alternative model architecture is Gaussian Process (GP) regression. GP is a Bayesian model that prescribes a multivariate normal distribution for possible functions mapping input x variables to output y. In regions with little or no training data, the model exhibits high uncertainty. In regions with substantial training data, the model predicts outcomes with greater confidence. The type of function in the distribution is determined by a kernel function that must be chosen and its hyperparameters (as with any other model) that are "tuned" for the specific problem. The kernel function measures the similarity of any two x inputs.

[0006] GP models have at least four overwhelming advantages. First, like deep learning, when equipped with an appropriate kernel, GPs can approximate arbitrary relationships. Second, GPs can compute marginal likelihoods (probabilities of data averaged over all possible function values) in a closed form. As a result, kernel hyperparameters "learn" by maximizing the marginal likelihood rather than the likelihood, reducing the risk of overfitting and thus making the model more robust. Third, model predictions are generated using the similarity of new data points to data points in the training set, where similarity is quantified by a kernel function. Unlike deep learning models, GPs are not black boxes because the model can measure the similarity between data points and determine how to generate predictions. The same kernel function that GPs use to make predictions can also be used to cluster training data and retrieve the most "similar" data points in the training set, where "similarity" is determined by the kernel function we choose.

[0007] Fourth, GPs robustly quantify uncertainty and assign high uncertainty to data points that are highly different from the training set. Uncertainty is important for protein engineering and drug design, where experimentally evaluating model predictions is expensive. In these scenarios, only model predictions that are likely to be reliable would be used. While GPs detect "distribution shifts" where new data points are highly different from the training set by assigning high uncertainty to their predictions, deep learning models can "silently fail," producing low-accuracy predictions without providing any obvious indication of failure. Uncertainty can also be used for Bayesian optimization / active learning, where the model assists the practitioner in selecting which data points to experimentally evaluate next. Several techniques for estimating uncertainty in deep learning have been introduced in the literature. However, the uncertainty estimates they provide are much more poorly calibrated than those provided by GPs.

[0008] GPs suffer from several limitations. Because designing an appropriate kernel (e.g., for image classification) without sufficient prior knowledge is challenging, the ability to select a kernel function can be both a drawback and an advantage. In some instances, it may be possible to overcome this drawback by combining deep learning with Gaussian processes, and indeed, this possibility has been explored by some practitioners in the literature. For example, a deep learning model can be trained in an unsupervised manner to learn feature representations that can be used as inputs to a GP. In this strategy, the GP essentially serves as the last layer of the deep learning model, thereby hopefully capturing some of the benefits of both approaches.

[0009] The main challenge with GP is its computational cost. A simple implementation runs O(N 3 ) scales with dataset size. In fact, the weak scaling of GPs is often cited as a reason to prefer deep learning. For example, when building models for traffic prediction, it was argued that GPs cannot be used for those problems because "GPs have difficulty scaling to large datasets." A second major challenge is the lack of efficient kernels for array and graph data. Many kernels for graphs have been described, but they typically exhibit quadratic or worse scaling in the size of the graph.

[0010] One approach for approximating a GP is the Random Fourier Feature (RFF) approach of Rahimi and Recht. Briefly, the RFF method converts each input data point into a "random feature" representation and then approximates the GP using Bayesian linear regression in the random feature space, thereby approximating the kernel function by a random map corresponding to the desired kernel such that linear scaling with training set size can be achieved.

[0011] Monoclonal antibodies are one of the most successful biological drugs. Despite their remarkable versatility and specificity, therapeutic antibody development continues to pose various complex challenges. Typically, early hits have insufficient affinity, and their binding must first be improved by in vitro affinity maturation, whereby antibodies with improved affinity are subjected to repeated rounds of mutagenesis and selection. This process is often time-consuming, sometimes taking months to complete, and does not allow for simultaneous optimization of other desirable properties, such as good solubility and low immunogenicity. Computational techniques that can help more quickly identify high-affinity antibodies with desirable properties would likely accelerate this process.

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

[0013] Machine learning-based approaches face at least two major challenges. The first is the lack of estimated uncertainty in predicting binding affinity or other properties. Because training data can cover only a small fraction of sequence space, machine learning models typically perform poorly when asked to extrapolate far beyond the boundaries of their training set. Deep learning models, while flexible and powerful, typically (except for certain architectures) do not provide confidence intervals or estimates of uncertainty in their predictions. While Gaussian process (GP) models, as discussed above, have been proposed as an option to provide well-calibrated confidence intervals, they do not scale well to large datasets without the use of approximations and are often infeasible for datasets of more than 5,000 sequences.

[0014] A second key challenge is the selection of an appropriate representative for the input. Many different encoding schemes for proteins have been described in the literature, including the most recent ones that employ language models; some of these are antibody-specific, and some are general to protein sequences. However, there is little consensus on which of these is most appropriate for a given problem. Classical one-hot encoding schemes are simple to implement, but unnecessarily high-dimensional and therefore useless, because all amino acids are treated as completely different from each other. Summary of the Invention [Problem to be solved by the invention]

[0015] Thus, high-affinity antibodies are often identified by directed evolution, which can require many iterations of mutagenesis and selection to find the best candidates. While deep learning techniques hold the potential to facilitate this process, current methods cannot provide the confidence intervals or uncertainties needed to assess the reliability of predictions.

[0016] Programmed cell death ligand 1 (PD-L1) is a protein that, in humans, binds to another protein called PD-1 as part of a pathway that suppresses the immune response against cells recognized as foreign. Many cancers overexpress PD-L1 to downregulate the immune response and prevent their own destruction by the immune system. Interference with PD-L1-mediated inhibition of the immune response has been shown to be a highly effective approach for treating some cancers, and several monoclonal antibodies (mAbs) that bind to PD-L1, including atezolizumab, durvalumab, and avelumab, have been approved by the FDA for various indications. However, there is a need to develop more potent and higher-affinity antibodies against PD-L1. [Means for solving the problem]

[0017] The present disclosure provides a novel pipeline that overcomes the limitations of the state of the art. 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 identifying high-affinity antibodies, a new autoencoder model trained on over 3 million human B cell receptor sequences is used to efficiently encode antibody sequences. A variational Bayesian neural network is then developed to perform ordinal regression on the directed evolution sequences and quantify their likelihood as strong binders to the antigen. Importantly, the model can also evaluate sequences not present in the directed evolution library, thus greatly expanding the search space for identifying the best sequences for further evaluation.

[0018] In certain embodiments, the power of this pipeline is demonstrated by achieving a 17-fold improvement in the KD of the PD-L1 antibody atezolizumab, a success that demonstrates the pipeline's ability to facilitate general antibody development and other protein design / engineering tasks. In certain embodiments, the present disclosure provides a set of 21 scFv sequences (fusions of heavy and light chain variable regions) that bind tightly to PD-L1. All of these scFvs are engineered variants of the parent atezolizumab mAb and were discovered using the machine-learning-assisted pipeline disclosed herein to discover antibodies that are tight binders to their targets of interest. The measured binding affinity of the top candidate from this set of 21 is approximately 17-fold stronger than atezolizumab, meaning it has one of the tightest affinities of any current anti-PD-L1 mAb. These scFvs are also believed to have slower off-rates than durvalumab and avelumab. Given their exceptionally tight binding, these scFvs may have potential as therapeutic agents or as tools for in vitro tumor research.

[0019] Because the scFvs provided herein function by forming tight, reversible, non-covalent interactions with the PD-L1 protein expressed on the surface of human cells, this may make it possible to treat cancer by inhibiting the PD-L1 / PD-1 interaction that some cancers use to suppress immune responses or, alternatively, serve as tools for in vitro assays involving cancer cells that overexpress PD-L1.

[0020] In addition to their potential use as antibody therapeutics or in vitro tools, the scFvs provided herein may be useful for CAR-T-based tumor therapy. CAR-T cells are immune cells engineered 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 seen significant growth recently, with the first FDA approval in 2017. CAR-T cells can potentially be engineered to express the exceptionally tight-binding anti-PD-L1 antibody scFvs disclosed herein and thereby recognize cancer cells that overexpress PD-L1. Indeed, CAR-T cells expressing anti-PD-L1 nanobodies have been reported to exhibit anti-tumor activity in vitro and in mouse models, suggesting that this may be a viable strategy. Given the exceptionally tight binding exhibited by the scFvs disclosed herein, these scFvs may be exceptionally well-suited for this task.

[0021] Thus, the present disclosure provides for the design and engineering of therapeutic proteins (such as antibodies), imaging probes (such as protein domains), and enhancing the efficiency of CRISPR and enzymatic activity.

[0022] Other systems, methods, features, and advantages of the present disclosure will be or become apparent to one of ordinary skill in the art upon examination of the following drawings and detailed description. All such additional systems, methods, features, and advantages are intended to be included herein, within the scope of the present disclosure, and protected by the accompanying claims. Furthermore, all necessary and preferred features and modifications of the described embodiments can be used in all aspects of the present disclosure taught herein. Furthermore, the individual features of the independent claims and all necessary and preferred features and modifications of the described embodiments can be combined with and interchanged with each other.

[0023] Many aspects of the present disclosure can be better understood by referring to the following drawings (also referred to as "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. [Brief explanation of the drawings]

[0024] [Figure 1A] Figures 1A-1B show the number of epochs (passes through the entire data set) required for fitting. In all cases, a relative error tolerance of 1e-6 is used as the threshold for convergence. [Figure 1B] Same as above. [Figure 2] Figures 2A-2D show the relative error versus epoch number: convergence rate of stochastic gradient descent and CG, with and without preprocessing. In early experiments, the closely related Adam stochastic gradient descent method, popular for fitting neural networks, performed even worse (it failed to significantly reduce error in less than 100 epochs), so we used AMSGrad. [Figure 3] Schematic of the Fast-Conv-1d kernel compared to a standard 3-layer convolutional neural network. [Figure 4] Comparison of the acceleration achieved for conjugate gradients using a preconditioner constructed using SRHT or using the unmodified algorithm S2. [Figure 5] Number of iterations for convergence across different benchmark datasets using different settings for preconditioner ranking and three different variants of the preconditioner construction algorithm. [Figure 6] Figures 6A-6B show the number of epochs / iterations required for convergence using CG for different datasets as a function of beta / lambda^2. [Figure 7] Distribution of absolute error percentages for marginal likelihood estimates across 112 evaluations using 25 or 40 probe vectors. [Figure 8] Kernel PCA on the QM9 test set using the GraphConv-1d kernel. [Figure 9A] 9A-9B show kernel k-means clustering of the QM9 training set using the GraphRBF kernel. [Figure 9B] Same as above. [Figure 10] Query molecule 1 is tested together with the training molecule to which it is most similar. [Figure 11] Query molecule 2 is tested along with the most similar training molecule. [Figure 12] Query molecule 3 is tested along with the most similar training molecule. [Figure 13] Query molecule 4 is tested along with the most similar training molecule. [Figure 14] The query molecule 5 is tested together with the most similar training molecule. [Figure 15] Schematic representation of the RESP computational and experimental pipeline. U-PDL1 and B-PDL1 are unlabeled and biotin-labeled PDL1, respectively. [Figure 16A]Figures 16A1-2 and 16B1-4 show sorting of the naive library by FACS and quality assessment of enriched populations by flow cytometry, respectively. Figure 16A1-2 shows sorting of the naive atezolizumab scFv library on the yeast surface compared to 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 (intermediate) gate was set for WT-level off-rates, and the left (lowest) gate was 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 the percentage of each group of the total number of sorted cells. Figures 16B1-4 show testing of each sorted population for binding intensity relative to WT after 8 hours of off-rate competition. LB is a low binder (faster off-rate), while MB is a medium binder (off-rate of WT), and HB is a high binder (slower off-rate). 105 cells per group were analyzed, and mean fluorescence intensity (MFI) is given on the x-axis per plot, along with V5 expression from AF647 intensity (y-axis) and binding in PE intensity (x-axis). [Figure 16B] Same as above. [Figure 17] In silico directed evolution using approximate Bayesian inference and learned embeddings for efficient candidate sequence evaluation. The 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. [Figure 18] The structure of an autoencoder model adapted to the task. The encoder module generates a training representation for each input sequence; the decoder module attempts to reconstruct the input sequence, while the classifier generates a binary prediction of the decoys for human B cell receptors in the training set. [Figure 19]Figures 19A-19C show the humanity scoring of decoy 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 humanity using the literature-derived AbLSTM model (Figure 19A), the ANARCI tool (Figure 19B), and the BioPhi model (Figure 19C). In all three cases, the model scores for the decoys were significantly different from those for the non-decoys, with the decoys being less human than the original sequences. In all three cases, using a two-tailed Mann-Whitney U test implemented in Python's Scipy library, version 1.5.4, the calculated p-value is 0.0 (meaning it is approximately 0 given floating-point error). The following conventions apply to each boxplot: The upper and lower limits of the box are the first and third quartiles of the data, and the whiskers are drawn at 1.5 times the interquartile range (the distance from the first to the third quartile). The center is drawn at the median of the data, and the "notches" represent the 95% confidence interval around the median (determined by non-parametric bootstrap). Diamonds represent "flyer" points outside 1.5 times the interquartile range. Four asterisks indicate a p-value less than 0.0001. Source data is provided as a source data file. [Figure 20]Figures 20A1-2 to 20C show an analysis of the pipeline results. Figures 20A1-2 show the marginal distributions per position in the two major clusters. Each major cluster has several subclusters, as evidenced by the marginal distributions. Only positions mutated in either cluster are shown. Each position is first numbered using the numbering from the wild-type sequence, and then numbered in parentheses using Chothia's numbering. Figure 20B illustrates the distribution of scores assigned to sequences in each binding category by a variational Bayesian model designed to perform ordinal regression (classification on ranked categories). The experimental data do not directly measure off-rates, but rather classify sequences into three bins: RH01 (high off-rate), RH02 (intermediate off-rate), and RH03 (slow off-rate). A higher score assigned by the ordinal regression model indicates a more reliable model for sequences with low off-rates, while a lower score indicates a more reliable model for sequences with high off-rates. Figure 20C shows the tolerance score versus iterations for a typical simulated annealing strand. The algorithm initially explores sequence space unbiasedly, and as the temperature is decreased, it is gradually forced to focus on the most promising regions it has found so far. Source data is provided as source data files for all panels. [Figure 21] Location of mutations among the 21 mutants in the structure of the atezolizumab heavy chain. Mutated residues are labeled in the heavy chain. Structure from the RCSB Protein Data Bank (5XXY) in Mol*Viewer46 71. [Figure 22]Figures 22A-22C show experimental validation of Koff and KD on the yeast surface. Figure 22A shows WT vs. mutant 4 scFv dissociation after 92 hours at RT on the yeast surface, where T1 / 2 is half-life. Figure 22B shows a comparison of WT atezolizumab, mutant 4, durvalumab, and avelumab scFv dissociation over 28 hours at RT on the yeast surface. Figure 22C shows binding affinity (KD) measurements determined on the yeast surface between the scFvs and PD-L1 (three independent measurements, see also Figures 29A-29F). [Figure 23] Figures 23A1-2 and 23B show variational Bayesian neural network architectures. Figure 23A1-2 shows a comparison of a variational Bayesian neural network with a fully connected architecture. Both networks map vectors of input values ​​x1, x2, ... xn to hidden layer vectors H1_1, H1_2, ... H1_n through matrix multiplication followed by nonlinear activation. However, in a fully connected network, each element of the hidden layer weight matrix is ​​a learned value that is fixed once the network is trained. In contrast, in a Bayesian network, 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, we can sample them from the weight distribution, which provides an estimate of the uncertainty about the prediction. Figure 23B shows the structure of the ordinal regression scoring model used in the pipeline. [Figure 24] Figures 24A-24B show that WT atezolizumab scFv binds to human PD-L1 on the yeast surface. The y-axis of each panel represents the expression of the scFv (V5 epitope / AF647), while the x-axis represents the binding intensity to human PD-L1 / PE (left / 24A: 0 nM PD-L1 vs. right / 24B: 6.1 nM PD-L1). The number of double-positive cells out of the total cells analyzed is listed below each panel. [Figure 25]Figures 25A-25D show the initial screening of the 21-member library after 8 hours of off-rate competition (vs. WT, top plots in Figures 25A and 25B) and the final screening of the library after 39 hours of competition (bottom plots in Figures 25C and 25D). The right triangle represents the sorting gate containing the relevant % of cells collected per sort. The y-axis of each panel represents scFv (V5 epitope / AF647) expression, while the x-axis represents binding intensity to human PD-L1 / PE intensity. [Figure 26] Figures 26A1-2 to 26B show BLI binding analysis of soluble WT and Mutant 4 scFvs to immobilized biotin-PD-L1. Figure 26A1-2 shows the binding curve (association and dissociation) profiles for different scFv concentrations. WT was titrated over a range of 0.65 to 475 nM, while Mutant 4 was titrated over a range of 0.219 to 160 nM. Figure 26B shows the on-rate / off-rate for WT and Mutant 4 scFvs. For each scFv, two different KD values ​​were obtained from two binding events to PD-L1. A value for Mutant 4 scFv was not obtained due to its very slow off-rate. [Figure 27] Figures 27A-27G show mass spectrometric analysis of PD-L1 and WT / mutant 4 scFv in monomeric, dimeric, and multimeric states. The images in each panel have 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 two separate purifications. The sigma symbol represents the standard deviation of the MW estimates calculated by the software, and the numbers above each peak are the apparent MW. [Figure 28A] Figures 28A-28B show SDS-PAGE gels (4-20%) of purified scFv-His6 fusions. M4 is mutant 4. The two gel images represent proteins from two separate purifications (two separate SDS-PAGE gels from two separate purifications). [Figure 28B] Same as above. [Figure 29-1] Figures 29A-29F show the binding curves between PD-L1 and scFv on the yeast surface. MFI means mean fluorescence intensity, and M4 means mutant 4. Each experiment is a separate / independent assay performed on a different day. [Figure 29-2] Same as above. [Figure 29-3] Same as above. [Figure 30] Figures 30A-30B show a side-by-side comparison of amino acids near R98C in the wild type (PDB structure 5XXY, Figure 30A, top) and the IgFold predicted structure for the top-scoring mutant (Figure 30B, bottom). [Figure 31] Superposition of the PDB structure 5XXY with the IgFold predicted structures for the top-scoring mutants. Note the large shift in conformation of CDR H3 and the subtle shifts in conformation of other regions that contact the antigen. [Figure 32]Figures 32A-32C show scores assigned by the model for sequences from Mason et al., Nat. Biomed. Eng. 5, 600-612, 2021. Figure 32A shows the distribution of scores assigned by the model for binding and non-binding sequences from the training set. 0 is the cutoff for predicting a sequence as a binder. Figure 32B shows the distribution of scores for sequences selected by Mason et al. as likely binders. Mason et al. experimentally demonstrated that all of these sequences bind with a KD of less than 20 nM and that all of them are indeed predicted to be binders by the model disclosed in Example 2. Figure 32C shows the model's uncertainty for its assigned scores for predictions on the test set. Inaccurate predictions have higher associated uncertainty than accurate predictions, indicating, as expected, that the uncertainty assigned by the model can be useful in determining whether a prediction should be considered reliable. There are 3,158 unique sequences in this test set. A two-tailed Mann-Whitney U test was used to assess significance; the resulting p-value was 1e-23. The following conventions are applied to each box plot: the upper and lower limits of the box are the first and third quartiles of the data, and the whiskers are depicted at 1.5 times the interquartile range (the distance from the first to the third quartile). The center is depicted at the median of the data, and the "notches" represent the 95% confidence interval above the median (determined by non-parametric bootstrap). Diamonds represent "flyer" points outside 1.5 times the interquartile range. Four asterisks indicate a p-value less than 0.0001. [Figure 33] Clustering of simulated annealing results. Dendrogram of the clustering of sequences harvested from the simulated annealing procedure after removing sequences with scores associated with wide confidence intervals. The tree suggests the presence of at least two major regions of sequences identified by the modified simulated annealing algorithm. [Figure 34]Example of FACS gating strategy. 106 cells ("high" gate) collected from 6 x 106 sorted yeast. Y-axis is V5 epitope expression, X-axis is PD-L1 binding. V5 tag monoclonal mAb detected with secondary AF647 polyclonal IgG. Biotin-PD-L1 detected with SA-PE reagent. DETAILED DESCRIPTION OF THE INVENTION

[0025] Additional advantages of the present disclosure will be set forth in part in the description which follows, and in part will be obvious from the description, or may be learned by practice of the present disclosure. The advantages of the present disclosure may 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.

[0026] Gaussian processes (GPs) are Bayesian models that offer several advantages for regression tasks in machine learning, such as reliable quantification of uncertainty and improved interpretability. Their adoption has been hindered by their excessive computational cost and the difficulty of adapting them to analyze sequences (e.g., amino acid sequences) and graphs (e.g., small molecules). In certain embodiments, the present disclosure provides a family of kernels approximated by random features for sequences and graphs that exhibit linear scaling with both the size of the training set and the size of the sequence or graph. These new kernels were incorporated into a 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 with both standard and state-of-the-art deep learning models on 17 different benchmarks, and GP regression was found to achieve highly competitive accuracy for these tasks while providing well-calibrated quantification of uncertainty and improved interpretability. Furthermore, the present disclosure provides how xGPR can be used as part of an active learning strategy to engineer proteins with desired properties in an automated manner without human intervention.

[0027] The present disclosure also provides a novel pipeline called RESP, which consists of directed evolution, sequencing, and machine learning analysis for the efficient identification of high-affinity antibodies. In one specific embodiment, we developed a learning representation trained on 3 million human B cell receptor sequences encoding antibody sequences. We also developed a variational Bayesian neural network to perform ordinal regression on a set of directed evolution sequences binned by off-rate to quantify their likelihood as strong binders to the antigen. Importantly, this model can evaluate sequences not present in the directed evolution library, thus greatly expanding the search space for identifying the best sequences for evaluation.

[0028] The ability of RESP to facilitate general antibody development is further illustrated. In one specific embodiment, this pipeline was applied to improve the affinity of a well-known antibody to a well-studied antigen by engineering mutants of the atezolizumab (trade name Tecentriq) heavy chain with improved affinity for programmed cell death ligand 1 (PD-L1). To do so, the atezolizumab heavy chain sequence was largely randomized, mutants were binned using yeast display and FACS, and each bin was sampled for sequencing. The sequence data for mutant atezolizumab was converted to a low-dimensional representation using an autoencoder model, and the encoded mutant atezolizumab sequences were used to train a Bayesian ordinal regression model that scored each sequence for its probability of being a strong binder to PD-L1. Finally, a modified simulated annealing algorithm was used to select sequences for testing.

[0029] In certain embodiments, the present disclosure provides a set of scFv sequences comprising fusion sequences of the heavy and light chain variable regions of atezolizumab mAb. The scFvs disclosed herein were discovered using the machine learning-assisted pipeline discussed herein and exhibit strong binding to PD-L1. One of the scFvs disclosed herein is the K of atezolizumab, a PD-L1 antibody. D The scFv disclosed herein is more potent than current antibodies against PD-L1 and can be used as a potent antibody against PD-L1 in cancer immunotherapy as well as in CAR-T cell therapy.

[0030] Those skilled in the art to which the disclosed compositions and methods pertain will envision many modifications and other embodiments of the present disclosure with the benefit of the teachings presented in the foregoing description and the associated drawings. It is therefore to be understood that the present disclosure is not limited to the particular embodiments disclosed, and that modifications and other embodiments are intended to be included within the scope of the appended claims. Those skilled in the art will recognize many variations and adaptations of the aspects described herein. These variations and adaptations are intended to be included within the teachings of the present disclosure and encompassed by the claims of this specification.

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

[0032] As will be apparent to those skilled in the art upon reading this disclosure, each individual embodiment described and illustrated herein has distinct components and features that can 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 disclosure.

[0033] 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 expressly stated otherwise, it is in no way intended that any method or embodiment described herein be construed as requiring that its steps be performed in a specific order. Thus, if a method claim does not specifically recite in the claim or specification that the steps are to be limited to a particular order, no order is intended to be inferred in any respect. This preserves any possible non-expressive basis for interpretation, including matters of logic regarding the arrangement of steps or operational flow, apparent meaning derived from grammatical construction or punctuation, or the number or type of embodiments described in the specification.

[0034] All publications and patents cited herein are cited to disclose and describe the methods and / or materials in connection with which the publication is cited. All such publications and patents are incorporated by reference herein as if each individual publication or patent was 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 lexicographic definitions derived from the cited publications and patents. Lexicographic definitions in cited publications and patents that are not expressly repeated in this application should not be treated as such and should not be read to define terms that appear in the appended claims. The citation of all publications is for their 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 publications provided may be different from the actual publication dates, which may need to be independently confirmed.

[0035] Although aspects of the present disclosure may be described and claimed in particular statutory classifications, such as systems statutory classifications, this is for convenience only, and those skilled in the art will understand that each aspect of the present disclosure may be described and claimed in any statutory classification.

[0036] It should also be understood that the terms used herein are for the purpose of describing specific embodiments only and are not intended to be limiting.Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by those skilled in the art to which the disclosed compositions and methods belong.Terms such as those defined in commonly used dictionaries should be interpreted as meanings consistent with their meanings in the context of this specification and related art, and should not be interpreted in an idealized or overly formal sense unless expressly defined herein.

[0037] Aspects of the present disclosure employ, unless otherwise indicated, techniques of molecular biology, microbiology, organic chemistry, biochemistry, physiology, cell biology, vascular biology, and the like, which are within the skill of the art, and such techniques are fully explained in the literature.

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

[0039] definition As used herein, "comprising" should be interpreted as specifying the presence of the stated features, integers, steps, or components being described, but does not exclude the presence or addition of one or more other features, integers, steps, or components, or groups thereof. Furthermore, the terms "by," "comprising," "comprises," "comprised of," "including," "includes," "included," "involving," "involves," "involved," and "such as" are each used in their open, non-limiting sense and may be used interchangeably. Furthermore, the term "comprising" is intended to include examples and embodiments 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."

[0040] As used herein, the term "and / or" includes any and all combinations of one or more of the associated listed items. When preceding a list of elements, phrases such as "at least one" modify the entire list of elements and not each individual element of the list.

[0041] As used herein and in the appended claims, the singular forms "a," "an," and "the" include plural referents unless the context clearly dictates otherwise. References to "a / an" compound, therapeutic agent, and pharmaceutical composition are not limited to one compound, therapeutic agent, and pharmaceutical composition, respectively, but rather refer to one or more molecules of the compound, therapeutic agent, and pharmaceutical composition, which may or may not be identical, so long as they fall within the category of 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., including different isotope abundances and / or different degrees of hydration or equilibria with different conjugate base or conjugate acid forms).

[0042] 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 both significant in relation to the other endpoint, and independently of the other endpoint. It is also understood that there are several values ​​disclosed herein, and that each value is also disclosed herein 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, the use of the antecedent "about" will understand that the particular value forms a further aspect. For example, if the value "about 10" is disclosed, then "10" is also disclosed.

[0043] When a range is expressed, a further embodiment includes from one particular value and / or to the other particular value. When a range of values ​​is provided, it is understood that each intervening value between the upper and lower limits of that range, to the tenth of the unit of the lower limit, and any other stated or intervening value in that stated range, is encompassed within the disclosure unless the context clearly dictates otherwise. The upper and lower limits of these smaller ranges may be independently included in the smaller ranges, and are also encompassed within the disclosure, subject to any specifically excluded limit in the stated range. When a stated range includes one or both limits, ranges excluding either or both of those included limits are also included within the disclosure. For example, when a stated range includes one or both limits, ranges excluding either or both of those included limits are also included within the disclosure; for example, the phrase "x to y" includes ranges from "x" to "y" as well as ranges greater than "x" and less than "y." Ranges may also be expressed as upper limits, e.g., "about x, y, z, or less," and should be interpreted as including the specific ranges of "about x," "about y," and "about z," as well as the ranges "less than x," "less than y," and "less than z." Similarly, the phrase "about x, y, z, or greater" should be interpreted as including the specific ranges of "about x," "about y," and "about z," as well as the ranges "greater than x," "greater than y," and "greater than z." Furthermore, the phrase "about 'x' to 'y'" (where 'x' and 'y' are numbers) includes "about 'x' to about 'y'."

[0044] Such range formats are used for convenience and brevity and should thus be interpreted in a flexible manner to include not only the numerical values ​​explicitly recited as range boundaries, but also all individual numerical values ​​or subranges encompassed within that range, as if each numerical value and subrange were expressly recited. For example, 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 individual values ​​(e.g., about 1%, about 2%, about 3%, and about 4%) and subranges (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 subranges) within the stated range.

[0045] As used herein, "about," "approximately," "substantially," and the like, when used in connection with a numerical variable, may generally refer to that value of the variable and all values ​​of the variable that are within experimental error (e.g., within a 95% confidence interval about the mean) or within ±10% of the stated value, whichever is greater. As used herein, the terms "about," "approximately," "at or about," and "substantially" may mean that the amount or value in question may be its exact value or a value that will produce an equivalent result or effect 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 necessary, to reflect tolerances, conversion factors, rounding, measurement error, and the like, and other factors known to those skilled in the art that will produce equivalent results or effects. In some circumstances, values ​​that will produce equivalent results or effects cannot be reasonably determined. In general, an amount, size, formulation, parameter, or other quantity or characteristic is "about," "approximately," or "at or about," regardless of whether it is expressly stated as such. When "about," "approximately," or "at or about" is used before a quantitative value, it is understood that the parameter also includes the particular quantitative value itself, unless specifically stated otherwise.

[0046] As used herein, the terms "optional" or "optionally" mean that the subsequently described event or circumstance may or may not occur, and that the description includes instances in which the event or circumstance occurs and instances in which it does not occur.

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

[0048] As used herein, the terms "treating" and "treatment" may generally refer to obtaining a desired pharmacological and / or physiological effect. The effect may be, but need not necessarily be, prophylactic, in terms of preventing or partially preventing a disease, its symptoms, or condition. The effect may also be therapeutic, in terms of partially or completely curing a disease, condition, symptom, or adverse effects resulting from a disease, disorder, or condition. The term "treatment" as used herein may include any treatment of inflammation associated with a disease in a subject, particularly a human, and may include any one or more of the following: (a) preventing the disease from occurring in a subject who may be predisposed to the disease but has not yet been diagnosed as having it; (b) inhibiting the disease, i.e., halting its development; and (c) alleviating the disease, i.e., reducing or ameliorating the disease and / or its symptoms or condition. The term "treatment" as used herein may refer to therapeutic treatment only, prophylactic treatment only, or both therapeutic and prophylactic treatment. Those in need of treatment (subjects in need thereof) may include those already with the disorder and / or those for whom a disorder is to be prevented. As used herein, the term "treating" can include inhibiting a disease, disorder, or condition, e.g., preventing its progression; as well as alleviating a disease, disorder, or condition, e.g., causing regression of the disease, disorder, and / or disorder. Treating a disease, disorder, or condition can also include improving at least one symptom of a particular disease, disorder, or condition, even if the underlying pathophysiology is unaffected, such as treating pain in a subject by administering a painkiller, even if such an agent does not treat the cause of the pain.

[0049] As used herein, the term "therapeutically effective amount" refers to an amount that is sufficient to achieve the desired therapeutic result or 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 depends on various factors, including the disorder being treated and the severity of the disorder; the specific composition used; the patient's age, weight, general health, sex and diet; the time of administration; the route of administration; the excretion rate of the specific compound used; the duration of treatment; the drugs used in combination with or simultaneously with the specific compound used, as well as similar factors that are within the knowledge and expertise of the physician and may be well known in the medical field.When treating a specific disease or condition, in some cases, the desired response can inhibit the progression of the disease or condition.This may only involve temporarily slowing down the progression of the disease.However, in other cases, it may be desirable to permanently stop the progression of the disease.This can be monitored by routine diagnostic methods known to those skilled in the art for any disease.The desired response to the treatment of a disease or condition may be to delay the onset of the disease or condition, or even prevent its onset.

[0050] Unless otherwise specified, temperatures referred to herein are understood to be based on atmospheric pressure (ie, 1 atmosphere).

[0051] Now that the embodiments of the present disclosure have been generally described, the following examples will describe some further and / or more detailed aspects of the present disclosure. Although the embodiments of the present disclosure will be described in relation to the following examples and corresponding text and drawings, it is not intended to limit the embodiments of the present disclosure to this description. On the contrary, it is intended to cover all alternatives, modifications, and equivalents that fall within the spirit and scope of the present disclosure.

[0052] The following examples are presented so as to provide those of ordinary skill in the art with a complete disclosure and description of how to make and evaluate the compounds, compositions, articles, devices and / or methods claimed herein, 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 otherwise indicated, parts are parts by weight, temperature is in ° C. or is at ambient temperature, and pressure is at or near atmospheric. [Example]

[0053] [Example 1]

[0054] Linear scaling kernels for protein sequences and small molecules outperform deep learning while providing improved uncertainty quantification and interpretability In this example, we present an approximate Gaussian process model that can achieve good scalability and performance on sequences and graphs. First, we introduce an efficient method for fitting approximate Gaussian processes to random features using preconditioned conjugate gradients (PCG) with a novel preconditioner. We prove the scalability of the resulting approximate GP to larger datasets, show that it outperforms alternative methods for approximating Gaussian processes on nine different datasets, and find that it achieves the same test set performance as an exact (unapproximated) GP for a protein manipulation task.

[0055] We then introduce a family of efficient, scalable, random feature approximation kernels for graphs and sequences. We compare the test set performance of GP regression with deep learning methods using these kernels on 17 different benchmarks. GP regression models achieve highly competitive performance.

[0056] The possibility of combining deep learning and GP was further explored by using embeddings as inputs to approximate GPs, which was shown to be often beneficial. It was further demonstrated that the uncertainty estimates provided by approximate GPs are significantly better calibrated than those provided by various uncertainty estimation techniques for deep learning. A brief study showed how approximate GPs can be used to "engineer" proteins with improved fitness without the need for human intervention. Furthermore, we demonstrated how approximate Gaussian processes can be used to cluster training sets and recover data points in the training set that are most similar to test points.

[0057] method Random Feature Approximation: Background In a typical Gaussian process model, given a training set X along with associated y values, a new data point x * and the associated y value f * Predictions were made regarding: p(f * |X,y,x * )=N(u * ,σ * ) (1) (where N(u * ,σ * )teeth, u * =k * T (λ 2 I+K) -1 y and σ * =k ** -k * k * T (λ 2 I+K) -1 k * (2) (The distribution is a Gaussian distribution showing the mean and variance of

[0058] where k * is the sum of each x in the training set i The kernel function k(x * ,x i) on the N training data points, and K is the length N vector formed by evaluating the k(x i ,x j ) points. λ is a kernel hyperparameter that indicates the expected level of noise in the data. * is the kernel function evaluated at new data points for all training data points, and k ** is the kernel function evaluated at new data points against itself. When fitting, we do not need to "learn" any parameters, but we do need to construct and invert an NxN matrix, which is expensive when N is large.

[0059] To avoid this cost, we can use the kernel function as a random Fourier feature 27 This well-known approximation gives k(x1,x2)≒z(x1) T Identify the random map for a given kernel such that z(x2) is (see Supporting Information section S1 for details). The mean and variance are

[0060]

number

[0061] Z is the row i of z(x i ) and z * is z(x * ) (wherein x * where λ is the transformed input matrix such that λ is the new data point. For a general invariant kernel, the error of the approximation decreases exponentially with increasing number of random features, such 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, etc.).

[0062] Additionally, we replaced the vanilla random features described above with a Structured Orthogonal Random Features (SORF) procedure. Briefly (see Supporting Information section S2 for details), for an example RBF kernel, the following transformations: SHD1HD2HD3(4) is applied to each input vector.

[0063] where S is a diagonal matrix with elements drawn from the χ distribution with d degrees of freedom, H is a normalized Hadamard matrix, and D1, D2, and D3 are diagonal matrices with elements drawn from the Rademacher distribution. The Hadamard matrix multiplication can be replaced with a fast Hadamard transform to reduce the cost of generating random features to O(NM 2 ) to O(NM logM). The resulting model is very lightweight, since only four diagonal matrices need to be stored. Generating predictions using the trained model is fast, since only a series of transformations and diagonal matrix multiplications need to be performed, followed by dot products. This change is not only beneficial in terms of speed, but also in terms of performance. This modification has been shown to improve the accuracy of kernel approximations, along with random features for RBF kernels. A fast Hadamard transform for both CPU and GPU has been implemented in the library (see Section S3 in the Supporting Information for benchmarking and other details).

[0064] Faster fitting with new preconditioner To fit an approximate GP model, we use [Z T Z+λ 2 I] -1 Z T We need to find y=w. If the number of random features is large (e.g., 5000-30000), we can use the matrix [Z T Z+λ 2 Directly forming and inverting [I] requires expensive matrix multiplications and factorizations. This problem achieves a nearly square scaling in the number of random features, Z with respect to w. T y=[ZT Z+λ 2 This was circumvented by instead using the well-known method of conjugate gradient (CG), which iteratively solves for [I]w. Briefly (see Section S4 of the Supporting Information for details), T Z+λ 2 I] -1 Rather than computing the matrix-vector product [Z T Z+λ 2 I]w. At no point is matrix multiplication or decomposition required.

[0065] For a vector of N rows in Z,

[0066]

number

[0067] Note that because CG evaluates the full gradient for each iteration, unlike fitting methods based on stochastic gradient descent, the mini-batch size has no effect on model performance or on the number of iterations required for convergence. Larger mini-batches can slightly reduce the time required for fitting, but cannot increase the memory footprint. Thus, the memory footprint is independent of the size of the dataset, and the user can determine this by choosing the mini-batch size.

[0068] Since CG is iterative, minimizing the number of iterations is important to make the model fitting as fast as possible. xGPR uses the matrix [Z T Z+λ 2 We use a randomized Nystrom approximation of [I] (see Section S4 of the Supporting Information for details), which is a known method for matrix approximation but has not previously been used as a preconditioner for Gaussian processes.

[0069] We also made modifications to the randomized Nystrom approximation that reduce the time required to build the preconditioner (see Supporting Information, Section S4 for details and benchmarks). Importantly, in this implementation, we build the preconditioner by loading only one mini-batch of data into memory at a time, minimizing memory footprint. Further details on how to use random features and more specifically preconditioned conjugate gradients for hyperparameter tuning are provided in Supporting Information, S5 and S6 below.

[0070] In Figures 1A and 1B, we compare the preconditioned conjugate gradient (PCG) method with the Scipy library's implementation of L-BFGS across seven datasets. L-BFGS (limited-memory BFGS) is a popular algorithm for minimizing a loss function using an iteratively constructed approximation to the Hessian matrix. The PCG strategy reduces the number of iterations required for convergence, and therefore the time to fit, by two orders of magnitude compared to L-BFGS or unpreconditioned CG, achieving a robust fit in as few as 10–20 iterations. In all cases, a relative error tolerance of 1e-6 was used as the threshold for convergence. For more details regarding the preconditioner ranking, see Supporting Information S4 below.

[0071] In Figures 2A–2D, we further compare the PCG method with two stochastic gradient descent methods: stochastic variance-reduced gradient descent (SVRG) and AMSGrad, an improved variant of the Adam algorithm, popular for training neural networks. We demonstrate how the residual error decreases with the number of epochs (full passes through the dataset). For PCG, each iteration is equivalent to one epoch, whereas for stochastic gradient descent, there are many iterations per epoch. As shown, PCG is the clear winner, converging several orders of magnitude faster than stochastic gradient descent, thereby significantly reducing the number of epochs required for fitting. Unlike stochastic gradient descent, CG does not require learning rate tuning to achieve good results. The PCG method combined with random feature generation can fit a dataset of over 400,000 data points in less than 7 minutes, including all hyperparameter tuning. When using a mini-batch size of 2,000 data points, the process consistently requires less than 1.5 GB of RAM.

[0072] In early experiments, the closely related Adam stochastic gradient descent method, which is popular for fitting neural networks, was even worse (it failed to significantly reduce error in less than 100 epochs), so AMSGrad was used in these studies.

[0073] Fitting 500,000 data points using an exact Gaussian process requires constructing 400,000 by 400,000 kernel matrices, which requires 1.3 terabytes of memory, and the entire kernel matrix is ​​required to fit the model. xGPR is thus much more scalable than exact Gaussian processes. In fact, as described under "Comparison with Stochastic Variational Inference and Exact Gaussian Processes," using the exact Gaussian process implemented in the GPyTorch library quickly encountered an out-of-memory error when attempting to fit a dataset with only 21,000 training data points, whereas xGPR was able to fit 500,000 data points without any issues.

[0074] result Stochastic variational inference and comparison with exact Gaussian processes Random features are not the only way to approximate Gaussian processes; one alternative is stochastic variational inference (SVI), which uses a set of "inducing points" much smaller than the number of training data points to approximate a GP. As described in Supporting Information section S8 below, both xGPR and approximate GP were fitted using stochastic variational inference (implemented in GPyTorch) to nine different datasets ranging in size from 400 to 500,000 data points. The same kernel was used in all cases to ensure a fair comparison. (See Supporting Information section S7 for details regarding these and all other datasets used for benchmarking in this example.)

[0075] These studies were performed on a GTX1070 GPU with 8GB of RAM. Two datasets were small enough (fewer than 5,000 data points) that they could be fitted with exact Gaussian processes, and therefore exact GPs were fitted for these two datasets as well. When we attempted to fit exact GPs to a third dataset with approximately 21,000 data points, the GPyTorch library immediately threw an out-of-memory error, indicating that exact GPs are already expensive for datasets of even modest size.

[0076] xGPR was found to perform as well as or better than SVI in all single cases; there were no cases in which SVI outperformed xGPR. For the two datasets in which xGPR could be compared to exact GP, xGPR achieved very similar performance, while SVI was similar in one case and significantly worse in the other.

[0077] We also compared xGPR with exact Gaussian processes on active learning / Bayesian optimization tasks from the literature. For example, data were generated for 52,910 mutants of green fluorescent protein (GFP), and the fluorescence of each mutant was measured. An exact Gaussian process was fitted to this dataset using 1,115 sequences with a single mutation as its training set and a RBF kernel, which uses learned embeddings to represent the input data points. In Bayesian optimization, an acquisition function was used to select new data points for experimental evaluation. The upper confidence bound (UCB), which is the predicted mean value of the GP plus a confidence interval, was used. The UCB prioritizes both data points where the model predicts that the mutant sequence has a high fitness value and data points where the model has high uncertainty, thus achieving a balance between utilizing current knowledge and exploring sequence space. UCB values ​​on the test set were shown to correlate reasonably well with actual fluorescence measurements.

[0078] Using the same training set, kernel, and representation, we fitted xGPR to the same dataset. The UCB acquisition function was found to exhibit a Spearman r correlation coefficient of 0.78 with the measured fluorescence in the test set. Surprisingly, the xGPR results were slightly better, at 0.8. We then selected 50 mutants from the test set with the largest UCB values ​​(to simulate the effect of experimentally evaluating the model predictions). The average log fluorescence of the top 50 for xGPR is nearly identical to results reported in the literature (both approximately 3.75). These results suggest that when using the same kernel, xGPR compares favorably with an exact Gaussian process while achieving significantly improved scalability.

[0079] The first new set of efficient kernels for unalignable sequences So far, only kernels for fixed-length inputs have been considered. In protein engineering problems, it is often possible to form multiple sequence alignments from all input sequences. Alternatively, when representing inputs using learned embeddings generated by large language models, the learned embeddings are averaged across sequences to generate fixed-length inputs. However, in many other cases, for example, when input sequences come from different families, the sequences are of different lengths and cannot be easily aligned. We introduce two types of random feature approximation kernels that can efficiently handle this problem: the first type is called FHT-Conv-1d, and the second type is called Fast-Conv-1d.

[0080] A first type of kernel can be constructed as follows: Consider extracting all k-mers of a given length, say 9, from two sequences and calculating the similarity of each possible pair of k-mers from the two sequences using a kernel of any general fixed vector of interest. There are a variety of different kernels that can be used for such comparisons, including RBF, Matern, neural network, and arccosine kernels, all of which can be approximated using random features. The "convolutional kitchen sink" kernel is a special case of this type of kernel. In this case, the comparison between two sequences is as follows:

[0081]

number

[0082] where k is a fixed vector kernel chosen to measure the similarity of any two subsequences. As written, this kernel scales as the square of the sequence length. Therefore, in these cases, it must be implemented as written, resulting in a very large number of pairwise k-mer comparisons that would be infeasible when used in exact Gaussian process or stochastic variational inference. For example, comparing two protein sequences of length 250 with a subsequence of length 9 would require extracting all 242 k-mers from both sequences and then performing over 29,000 pairwise k-mer comparisons.

[0083] A method for using random feature approximation to implement this kernel in a manner that scales linearly with the length of the sequence is shown below. Using random feature approximation, this kernel is

[0084]

number

[0085] where the two z vectors are random features generated for each pair of k-mers in sequences x and y. The final equation shown above is

[0086]

number

[0087] Thus, random features can be generated for each subsequence of length k in two input sequences and summed to obtain a single vector representing sequences x and y, after which these representations are used as inputs for Bayesian linear regression. The computational cost of this approach scales linearly with both dataset size and sequence length. The generation of random features for k-mers in each input sequence can be achieved by convolution with filters that are randomly weighted and drawn from a distribution appropriate for the kernel we wish to approximate, which can be thought of as a convolution kernel. However, this approximated kernel remains of limited practical use because convolution with thousands of filters is slow.

[0088] We then introduced a second innovation that greatly accelerates the computation of random features for these kernels. To understand this implementation, first note that in deep learning, library convolutions are often performed using matrix multiplication as follows: Consider an input KxP matrix with P sequence elements (e.g., P amino acids) and K features per element (or amino acid), and assume the convolution filter is of width L. Using the input matrix, we can append LKx(P-L+1), where each column i is the concatenated feature associated with subsequence i. Meanwhile, we can stack m convolution filters to form an mxLK matrix; multiplying this matrix by an input matrix of LKx(P-L+1) is equivalent to performing a 1D convolution with m filters of length L.

[0089] In this implementation, we replace the mxLK matrix of convolution filters with SORF operations (Equation 4), replacing matrix multiplications with fast diagonal matrix multiplications and fast transforms, similar to replacing the Discrete Fourier Transform with an FFT. For simplicity, let m be the number of desired random features, which is the same as LK. Then, the implementation of this kernel using vanilla random features takes O(N(P-L+1)m 2 In contrast, the proposed modification scales as O(N(P-L+1)mlog(m)).

[0090] Depending on the kernel of the fixed vector used to measure the distance between k-mers, how this class of kernels encodes the amino acids in two input sequences can be chosen so that the distance between k-mers is measured in a particular way. For example, if one-hot encodes amino acids, the distance between any two k-mers is twice the Hamming distance. Alternatively, each amino acid can be represented using a learned representation derived, for example, from a large language model. These possibilities are described below under "Evaluating Convolutional Kernels on Sequence Data," and this first type of kernel is referred to as FHT-Conv-1d.

[0091] A second efficient convolution kernel for unalignable sequences We introduce a second efficient group convolution kernel, called Fast-Conv-1d, which mimics a three-layer neural network. We state the following simple lemma (see Section S9 in the Supporting Information for the proof):

[0092] Lemma 1.1.

[0093]

number

[0094]

number

[0095]

number

[0096]

number

[0097] 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 similar 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 (Figure 3). The fixed vector kernel is approximated using a random feature procedure based on the fast Hadamard transform. This combined kernel is positive definite according to Lemma 1.1. This kernel has the property that all convolutions can be performed exactly once before training, and the results are saved to disk. The fast Hadamard transform-based procedure outlined above can be used to perform random feature convolutions in an efficient manner.

[0098] To understand how this type of kernel quantifies sequence similarity, consider performing a convolution of a one-hot encoded sequence with a filter of width 9 containing weights drawn from a normal distribution. For each subsequence of length 9, the convolution is large and positive if the filter happens to contain a large positive value at the same position in the subsequence that is one-hot. ReLU activation and max-pooling across the sequence are performed for each convolution filter, so 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” the input of each sequence to quantify the best match in the input for each filter in a set of random filters, and then compares the profiles of any two input sequences using a specific fixed vector kernel (e.g., RBF or linear+RBF; for simplicity, RBF was used in all studies). As before, sequences can be encoded using learned representations, one-hot encoding, or some other scheme, as desired.

[0099] Evaluation of xGPR and convolution kernels on sequence data We evaluated the xGPR and FHT-Conv-1d and Fast-Conv-1d kernels by fitting protein engineering benchmark data (the TAPE and FLIP benchmarks, respectively). These datasets involve situations where participants are provided with a pre-constructed training set and a test set of protein sequences, with the goal of predicting a real-valued outcome for each sequence. Briefly, these datasets are as follows (see Section S7 of the Supporting Information below for further details regarding these and all other datasets used in this example):

[0100] The AAV dataset contains mutations, including insertions and deletions, across a 28-amino acid window of the adeno-associated virus capsid; some of these mutations were randomly generated, while others were "designed" (intentionally selected). The data were divided in various ways, and these various partitions were called "splits." For example, they trained the model on randomly mutated proteins and tested them against designed ones ("mutant vs. designed split"), or trained on designed proteins and tested them against mutated ones ("design vs. mutant split"), or trained on proteins with only seven mutations and tested them against the rest ("7 vs. rest split"). The goal was to predict the fitness of the mutants; the higher the better.

[0101] The thermostability dataset uses mass spectrometry to measure the thermostability (higher is better) of 48,000 proteins from 13 species; the resulting proteins are highly diverse, and two splits were used: "mixed" (including all species) or "human" (a smaller subset including only human data). The GB1 dataset measures the fitness of mutants in the binding domain of a G protein from a bacterial species, with mutations at exactly four sites in the protein. The dataset contains near-complete coverage, i.e., many of the possible four-site mutants were generated and tested. The "3 vs. rest" split uses mutants with only three mutations as the training set and all others as test, while the "2 vs. rest" split uses mutants with two mutations as training.

[0102] The fluorescence dataset derived from TAPE contains mutations to green fluorescent protein; the goal is to predict the log fluorescence for a test set of proteins based on the training set. This dataset is the same as that used in the literature, discussed under "Comparison with Stochastic Variational Inference and Exact Gaussian Processes," but the TAPE benchmark splits the data into separate training and testing sets. Finally, the stability dataset derived from TAPE measures the stability of a diverse set of proteins in a protease susceptibility assay, with larger values ​​indicating better stability.

[0103] It was fitted with three iterations using different random seeds, and the resulting error bars measure the variability across different random seeds. Since our focus is on Gaussian process regression, we only used benchmarks that involve regression tasks.

[0104] As shown in Table 1, we compared the FHTConv-1d and FastConv-1d kernels trained on one-hot encoded data with convolutional neural networks trained on one-hot encoded data. As shown in Table 2, we used learned embeddings to represent input sequences (representations of input sequences generated by the FAIR ESM-1v pre-trained language model for protein sequences) and compared these results with those of fine-tuned language models reported in the literature. For each benchmark, we report the best result from a fine-tuned large language model described in the literature.

[0105] [Table 1]

[0106] [Table 2]

[0107] There are at least two ways that the learned embeddings can be used as input to a Gaussian process, and we evaluated both. First, we can use an RBF kernel, specifically after averaging the representation for each amino acid across all amino acids in the sequence. Alternatively, we can arrange the learned embeddings of size D as an MxD array for each sequence of length M, and perform convolutions across the sequence using the FHTConv-1d or FastConv-1d kernel. Training with the FastConv-1d kernel is easier because all convolutions can be performed once at the beginning of training, thus eliminating the need to save the learned embeddings on disk (which can be very large). Therefore, we used FastConv-1d here.

[0108] We found three interesting findings. First, when trained on one-hot encoded inputs with a small number of random features (<10,000), xGPR outperformed CNNs on 8 of the 11 benchmarks and tied on another. When trained on learned embeddings, approximation GP outperformed fine-tuned language models on 6 benchmarks and tied on another. Note that while the pre-trained language models required up to 50 GPU-days of training time for fitting, the FastConv-1d kernel took less than 6 minutes per fit (including hyperparameter tuning) to fit all datasets shown here, with a maximum memory footprint of 1.7GB when using a mini-batch size of 2000 data points. (The FHTConv-1d kernel is slower, but always takes less than 15 minutes, including hyperparameter tuning.) The use of learned embeddings generally improves the performance of GPs, especially for diverse datasets containing sequences from many different protein families (e.g., thermostability), but not always, and for some benchmarks, GPs perform better than one-hot encoding.

[0109] Second, the use of a large 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 users with a simple way to improve performance if needed. Third, neither the FHT-Conv-1d nor the Fast-Conv-1d kernel consistently outperforms the other; rather, each outperforms the other on some datasets.

[0110] Uncertainty calibration The main motivation for using 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.

[0111] For uncertainty estimates to be useful, they need to increase as you move further away from the training set—this is true for GPs with invariant-by-construction kernels. Furthermore, well-calibrated uncertainty estimates are desirable; in other words, the uncertainty reflects how often the predicted outcome actually occurs. For example, if a weather forecast claims there is a 10% chance of rain, then 10 out of 100 days should rain at the time this prediction is made. Similarly, for a regression model, roughly 10% of the test data points are expected to fall within a 10% confidence interval (or "credible interval" in Bayesian inference) in its prediction, roughly 20% of the data points are expected to fall within a 20% confidence interval, 90% within a 90% confidence interval, etc.

[0112] Calibration can be quantified as the area under the calibration error curve (AUCE). For 100 values ​​ranging from 0 to 1 (0.01, 0.02, etc.), we determined what fraction of the test set predictions fell within the corresponding confidence interval (1%, 2%, etc.). We subtracted the observed fraction within the interval from the predicted fraction and took the absolute value. For example, if the confidence interval is 10% and the observed fraction within that interval is 90%, the difference is |0.9-1|. Finally, we calculated the area under this curve; the larger the area, the weaker the calibration of the uncertainty provided by the model. A perfectly calibrated model would achieve a score of 0.

[0113] Five different methods for uncertainty estimation were described (dropout uncertainty, final-layer stochastic variational inference, model ensemble, deep evidential regression, and mean-variance estimation) and used in combination with deep learning models on the five benchmarks described in this example. Both one-hot and ESM embeddings were used as inputs to the models to calculate AUCE. As shown in Table 3, the AUCE for these uncertainty estimation methods was compared to the AUCE for xGPR on the same dataset using the same encoding. To give the deep learning model 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 showed the best performance for that dataset (as assessed by the correlation between predicted and measured fit) was used.

[0114] [Table 3]

[0115] When uncertainty calibration is performed for all five benchmarks, xGPR performs substantially better, with AUCEs 3-10 times better than the best deep learning uncertainty estimation methods. AUCE for xGPR is below 0.15 in all cases, while deep learning methods often (and for some methods, consistently) achieve AUCE scores of 0.3 or above.

[0116] Using uncertainty for protein engineering Uncertainty can be used to determine whether predictions generated by a model should be trusted; if the uncertainty is high compared to the uncertainty for the training set, the prediction may be considered 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 a model to the available data. A new batch of sequences was selected and experimentally evaluated using an acquisition function that takes into account both the model predictions and the uncertainty, thereby achieving a balance between exploration and exploitation. For example, an upper confidence limit (UCB) adds the upper half of the confidence interval to the model prediction for each test point. After experimental evaluation of the selected points, the model was fit, and this process was repeated until the desired protein function was achieved. For a review of other common acquisition functions, see Frazier 2018.

[0117] In this study, we used the GB1 dataset discussed under "Evaluation of Convolution Kernels on Sequence Data," in which mutations were introduced at four sites in the binding region of a G protein. This dataset is favorable for this type of experiment because approximately 150,000 of the 160,000 possible mutants were evaluated. The fitness of all mutants in the dataset was normalized so that it was between 0 (worst possible) and 1 (best possible). Note that over 99.95% of the sequences had a fitness below 0.6, with only one having a fitness of 1. After randomly selecting 384 variants and training the xGPR model with the RBF kernel, UCB acquisition function scores were generated for all remaining mutants and used to select 96 batches for "experimental" evaluation. These 96 were then added to the training set, and the model was fit; the process was repeated for up to five iterations. The best fitness achieved was tracked in each iteration. Finally, this "experiment" / "process" was repeated 50 times with different random seeds.

[0118] By the fifth iteration, all iterations achieved a best fit of 0.6 or higher. Forty of the 50 iterations found one of the four best sequences, and 30 of the iterations achieved a score of 1, i.e., they found the best sequence in the dataset despite experimentally testing fewer than 1,000 sequences. Half of the iterations were able to find the best sequence (fitness 1) in four iterations or less. These results suggest that when combined with a high-throughput experimental system, this approach using xGPR may be able to "discover" proteins with desirable properties without the need for human intervention.

[0119] 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 unalignable sequences" above. These kernels can be easily extended to graphs if each node is encoded in the graph as a feature vector containing some information about the node, its boundary, and its immediate neighbors. For the Fast-Conv-1d kernel, we can apply random feature convolution and max-pooling operations across all feature vectors associated with the graph. For the FHT-Conv-1d kernel, we simply used the node and its associated features as inputs rather than k-mers and performed a width-1 convolution across the nodes of the graph: this is referred to as the Graph-Conv-1d kernel.

[0120] Several variants of the Graph-Conv-1d kernel were implemented in xGPR. Here, we focused on the simplest variant, Graph-RBF, which uses an RBF kernel to approximate pairwise comparisons across the nodes of two graphs. To evaluate the performance of Graph-RBF, we considered predictions of energetic and thermodynamic properties of small molecules. We benchmarked the QM9 dataset, which contains 134,000 small molecules whose geometry was optimized using density functional theory (DFT). Because these properties can, in principle, be decomposed as a sum of per-atom contributions, the kernel described here is particularly likely to be a good fit for energetic / thermodynamic properties based on prior knowledge.

[0121] Two different modeling approaches were used. In approach (1), each atom I was one-hot encoded, then this was concatenated with the one-hot encoding of each other atom j in the molecule and divided by the sixth power of atom j's distance to atom i; the one-hot encodings were filtered by distance. Due to the sixth power weighting, after 4 or 5 angstroms, features associated with more distant atoms become negligibly small and can be ignored. This representation was then used as input to the Graph-RBF kernel. This simple representation, while easy to construct, is not very informative (it discards any information about bond angles) and thus is unlikely to achieve good performance.

[0122] Alternatively, in method (2), each atom was represented using SOAP features. In this scheme, the atomic density around each atom is expanded with an orthonormal basis set, and the power spectrum of these coefficients is used as the representation. A distance weighting scheme was used with only minor modifications to its selected hyperparameters. For further details on the hyperparameters used and their tuning, see Section S10 in the Supporting Information.

[0123] Notably, method (1) achieves a chemical accuracy of approximately 1 kcal / mol despite using uninformative features. (2) is much better; the results are shown in Table 4. For internal energy, enthalpy, and Gibbs free energy, xGPR is able to achieve an MAE of less than 0.2 kcal / mol, achieving a new state-of-the-art for one of six properties. Note that 1 kcal / mol is typically considered "chemical accuracy," i.e., the accuracy required to make accurate predictions about the behavior of chemical systems, so both xGPR and some other models easily exceed the level of accuracy required in practice. In fact, density functional theory (DFT) calculations of the outcomes we are trying to predict are often only accurate to 2–3 kcal / mol. As a result, a difference of 0.01–0.05 kcal / mol between models for predicting DFT-calculated energies is likely equivalent for practical purposes.

[0124] [Table 4]

[0125] Recovering similar data points from the training set Feature importance can be measured using SHAP or LIME techniques for Gaussian processes just as with any other ML model. Unlike other ML models, for small datasets, a kernel matrix can also be constructed and used to determine how accurately each data point in the training dataset contributed to a given prediction. This kind of decomposition is not possible with deep learning models.

[0126] However, when the dataset is larger than 5,000 data points, constructing the kernel matrix is ​​clearly too expensive. Also, for large datasets, assigning an exact contribution to each point in the training set may be more informative or helpful than necessary. In many cases, a practitioner may simply want to know which are the most similar data points in the training set, which contributed most to the prediction. This can be determined as follows: * Let be the random feature representation of the test data points, and z i Let be a random feature representation of the training data points, then the corresponding test and training data points x * and x i The kernel function for is approximately Z * T z i (See Section S1 of the Supporting Information).

[0127] In practice, taking the dot product of a test data point representation with all training data each time to retrieve the most similar training data point can be computationally intensive. This problem can be solved by simply clustering the training data. The same kernel used to fit a Gaussian process regression model can be used to cluster it with kernel k-means clustering approximated by random features. Note that in the same way that random features approximate a Gaussian process as a linear regression, clustering a random feature-represented training set with k-means approximates applying kernel k-means to the original dataset. Once the dataset is clustered, it is no longer necessary to measure the similarity of a test point to all training data points to retrieve the most similar training points; it is only necessary to measure its similarity to the data points in the nearest cluster.

[0128] The tools for performing these functions are included in the xGPR library. This procedure is illustrated with further details for the QM9 dataset in Section S11 of the Supporting Information. After initially fitting the QM9 dataset using the Graph-RBF kernel with atoms encoded using the one-hot encoding method (1) and an energy of 298 K as the target, k-means clustering of the training set is performed. The resulting elbow plots used to select the number of clusters are shown in the Supporting Information and Figures 9A-9B. The elbow plots suggest that the training dataset can be collapsed into five clusters. Five molecules from the test set were then randomly selected. For each of them, the procedure described in the previous paragraph was used to find the 50 most similar molecules from the training set, whose similarity is measured by the kernel function. The structures of the query test molecule and the most similar training molecule (quantified by the kernel using the one-hot encoded inputs) are shown in Section S12 of the Supporting Information.

[0129] Support Information The supporting information provided below includes: (S1): a more detailed background on random feature approximation; (S2): a more detailed background on orthogonal random features; (S3): performance metrics for our implementation of the fast Hadamard transform; (S4): preprocessing of the algorithm and some benchmarks for modifications to this algorithm; (S5): an algorithm for approximating the logarithmic determinant for random feature approximation Gaussian processes; (S6): a discussion of hyperparameter tuning strategies; (S7): details on the benchmark dataset used in this example; (S8): benchmarking random feature approximation versus stochastic variational inference; (S9): a proof of Lemma 1.1; (S10): details on the QM9 benchmark; (S11): clustering and principal components using random features; and (S12): recovering similar molecules from a training set using random features.

[0130] S1. Random Features Only positive definite kernel functions and only initially invariant kernels, i.e., K(x i ,x j )=f(x i -x j ) and consider the kernel (x i ,x j )≒φ(x i )φ(x j ) From Bochner's theorem, any of the considered positive definite kernels is the Fourier transform of the corresponding positive measurement, i.e.:

[0131]

number

[0132] This expectation can be approximated using Monte Carlo sampling, i.e.:

[0133]

number

[0134] The Monte Carlo approximation is an unbiased estimator of the mean, so as the number of random features increases, the estimate converges to the true value. The variance of this estimator is proportional to (1 / M).

[0135] Although this procedure originally came from the context of invariant kernels, it can be applied to similar Monte Carlo sampling techniques for certain non-invariant kernels, most importantly a class of kernels referred to as "neural network kernels", and in more recent literature also as "pointwise nonlinear Gaussian kernels". A neural network with one hidden layer and one output layer:

[0136]

number

[0137] where x is the input vector and W is the learned weight matrix of dimension DxN H where D is the dimensionality of the input and N H is the number of units in the hidden layer, and b is the length N H is the offset or bias vector of , h is some activation function, and ν i are the learned weights. All weights in W are N(0, σ w 2 ) and the weight ν i is N(0,σ ν 2 ) and the bias term is N(0,σ b 2 ) If each x is augmented with an additional element that is 1 for all data points, this can be written more simply as

[0138]

number

[0139] where each row of W is N(0,Σ w ) and Σ w If is a diagonal covariance matrix, its first diagonal element is σ b 2 and all remaining diagonal elements are σ w 2 is.

[0140] Let S denote the set of all weights, which is

[0141]

number

[0142] Again, the expectation term can be approximated on the right-hand side using Monte Carlo sampling, resulting in a formulation very similar to that for the invariant kernel. This kernel is uniquely positive definite for a particular choice of activation function h. One such choice is the ReLU activation function, which results in an arccosine kernel of degree 1, although the error function has also been suggested as a useful activation function. Interestingly, the arccosine kernel of degree 1 formulated in the literature omits the bias term. During library implementation, it was found that the arccosine kernel does not perform very well without the bias term and requires its inclusion to achieve acceptable performance.

[0143] For both invariant and neural network kernels, the random Fourier feature approach can then be implemented as follows: If X is an NxD matrix for N data points and D input features, then a DxM matrix is ​​added with M samples from p(ω), where p(ω) is the Fourier transform for the invariant kernel and the normal distribution for the neural network kernel. This sample matrix is ​​multiplied by X, and an additional bias vector for the neural network kernel is added to the result. Then, an appropriate activation function (e.g., sine and cosine for the invariant kernel) is applied and multiplied by a constant to obtain an NxM matrix Z (or Nx2M for the invariant kernel). The approximate kernel matrix is ​​then denoted by ZZ. T Give by.

[0144] A new data point x on a training dataset X given by a zero-mean Gaussian process with corresponding ground truth value y. * Predicted value f at * The distribution over p(f * |X,y)~N(μ * ,Σ * ) μ * =K(x * ,X)[K(X,X)+λ 2 I]-1 y, σ * =K(x * ,x * )-K(x * ,X)[K(X,X)+λ 2 I] -1 K(X,x * ) is.

[0145] The kernel is ZZ T These can be rewritten using the matrix inversion lemma and Woodbury's identity, since the mean is approximated as

[0146]

number

[0147] where Z is the training set transformed as above, I is the identity matrix, λ is a hyperparameter that ensures the kernel matrix is ​​positive definite, and z are the new data points transformed as above. [Z T Z+λ 2 I] -1 Z T y If we consider the weights to be , the predicted mean corresponds to a linear regression in feature space that matches the random feature map. It is easy to see the correspondence between a neural network with a single hidden layer plus output layer and the above kernel approximated using this formulation, but there are some key differences. The approximate GP "learns" by varying the distribution of hidden layer weights rather than individual values, and can be fitted by maximizing the marginal log-likelihood rather than the likelihood; and the marginal likelihood and posterior predictive distributions are available in closed form. It is also possible to construct a kernel similar to a neural network with several hidden layers. In this paper, we introduce the Fast-Conv-1d kernel, which corresponds to a three-layer neural network with convolutional, fully connected, and output layers, all with random weights added (see Figure 3).

[0148] Note that the error in the approximation decays exponentially with increasing number of random features, so a large number is generally not required to obtain a good approximation; in fact, increasing the number of random features past 16,000 or so has sharply diminishing returns, and 8,192 is sufficient for many tasks. * For predicting σ, a practitioner may want to use a larger number of random features to "buy" a slight improvement in accuracy, since this is the most important quantity, but for quantifying variance, 512 to 1,024 is generally sufficient. As a result, xGPR allows users to use fewer random features to quantify variance rather than to predict the mean (reducing the computational expense of predicting variance).

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

[0150] If the rows of M are likely to be orthogonal, first take a QR decomposition of M, then multiply the resulting Q by a diagonal matrix S and subtract its diagonal elements from a χ distribution with d degrees of freedom. Here, the input vector x is transformed by calculating the product kSQx. If the number of random features is greater than d, this procedure may be repeated m times, where m is the result of a ceiling division of D by d, and then any features in D in excess are discarded. If D is less than d, d random features can be generated and any in excess discarded.

[0151] It has been shown that for certain kernels, in particular the square exponential kernel, fewer orthogonal features are required to accurately approximate the kernel, and therefore 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 achieve a significant reduction in computational cost. The matrix Q can be further modified by multiplying it by: HD1HD2HD3 It has been shown that nearly equivalent results can be achieved by replacing

[0152] where H is the normalized Hadamard matrix, and each D is a separate diagonal matrix whose entities are drawn from the Rademacher distribution. The memory footprint of the fitted model is greatly reduced because only the diagonal elements of the D matrix need to be stored. More importantly, by replacing the Fast Hadamard Transform with a Hadamard matrix multiplication, similar to the replacement of the Fast Fourier Transform with a matrix in the Discrete Fourier Transform for multiplication, we reduce the cost of this procedure to O(DlogD), a dramatic gain for large M. The general procedure for transforming inputs is shown for the RBF kernel in Algorithm S1 (other kernels are similar).

[0153] The S matrix used here is a constant,

[0154]

number

[0155] [Table 5]

[0156] S3. Performance of the Fast Hadamard Transform A Python-wrapped C / Cuda implementation of the fast Hadamard transform for CPU and GPU has been implemented. We compared the speed of our implementation for matrix multiplication on matrices of random features implemented in the Numpy library (v1.22.3) and the Cupy library for GPU (v10.4.0) (see Tables 6 and 7 for results). Briefly, our implementation is faster than matrix multiplication for a small number of random features, and the performance advantage of our implementation increases with a larger number of random features, as expected, and can be quickly scaled to large numbers of random features. For comparison with other fast transforms, we also compared the discrete cosine transform on CPU implemented in the Scipy library (v1.8.0) (see Table 1 for results).

[0157] The results in milliseconds are shown in Tables 5 and 6. In Table 5, we evaluated the time required to generate random features when the input matrix has the same number of columns as there are features. This is a somewhat unorthodox situation, but provides a dramatic demonstration of the speedup for procedures based on the fast Hadamard transform. In Table 6, we considered a more realistic situation where there are a fixed number of 1024 input features and a large number of random features need to be generated.

[0158] An input matrix with 2000 rows was used because this chunk size is common when working on real-world datasets. For the GPU implementation, 32-bit (single precision) floating-point numbers were used due to the substantially higher performance for single-precision matrix multiplication on many GPUs, while for the CPU tests, 64-bit (double precision) was used. However, in general, xGPR performs random feature generation using single precision for both CPU and GPU unless the user selects otherwise.

[0159] [Table 6]

[0160] [Table 7]

[0161] S4. Preprocessing for conjugate gradients and stochastic gradient descent Consider the well-known conjugate gradient algorithm for solving the system Aw=b for w. The residual error at iteration n is bounded by

[0162]

number

[0163] where κ is the condition number of A and ||e0|| is the error of the initial guess. If κ is large, the algorithm will therefore converge very slowly. The inverse M -1 A preconditioner where A is an approximation to the inverse of A reduces the condition number of A and accelerates convergence. The matrix Z requires expensive matrix multiplications and is too large to store in memory due to the large number of random features. T It is important to have a preconditioner that can be constructed without specifically forming Z.

[0164] The kernel matrix for a Gaussian process is a low-rank approximation of: A≒A m,q A q -1 A m,q T has often been approximated using the Nystrom method, which gives

[0165] where m is the number of rows in A and q is a randomly chosen subset of the columns of A. The Nystrom approximation is inexpensive to construct and has worked well in practice for many problems. Sampling from a uniform distribution has traditionally been the most popular approach, but the performance of this method is sensitive to the sampling method used to select the columns and rows in A.

[0166] The randomized Nystrom approximation was found to provide better performance. In this scheme, for the system: (Z T Z+λ 2 I) -1 Z T y=w is.

[0167] Z T A low-rank approximation of Z=A is A nystrom ≒(AΩ)(Ω T AΩ) † (AΩ) T Provided by.

[0168] where Ω is a zero-mean, iid matrix, unit variance normal distribution; Ω is of size MxL, where M is the number of random features. This approach provides a rank-L approximation to A, and a numerically stable method for generating it (implemented for the library), found in Algorithm S2. Importantly, this is because the matrix Z T It can be constructed as shown without ever explicitly forming Z or using a single loop over the data set, loading the data into memory one chunk at a time.

[0169] [Table 8]

[0170] Here, we introduce a new modification: we define the MxL matrix Ω as follows:

[0171]

number

[0172] where D is a diagonal matrix with entities drawn from the Rademacher distribution, H is a 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 is used rather than the Hadamard matrix. For step 9b, it is Q+SRHT(Z T )Z Replace with.

[0173] In step 12, the product Ω T Q ν Regarding SRHT(Q ν ) is replaced. Finally, to ensure numerical stability, the above procedure (steps 11 and 15) is subtracted after replacing the shift, where SRHT(Q ν ), the self-adjoint square root of this matrix is ​​found using SVD and then used to solve for B in the following steps.

[0174] Rather than constructing a preconditioner for CG / stochastic gradient descent proposed here, constructing a low-rank matrix approximation using SRHT has been previously proposed in the literature. It significantly reduces the cost of preconditioner construction by replacing matrix multiplications with fast transform operations using O(NML) scaling; this is particularly advantageous when running on a CPU. Empirically, the SRHT-based preconditioner construction routine has been found to be substantially faster even for datasets of only 40,000 data points, and the advantage increases with increasing dataset size and number of random features.

[0175] Empirically, in all cases encountered, preconditioners constructed using the proposed SRHT performed as well as or better than preconditioners constructed using the unmodified algorithm S2. As shown in Figure 4, the number of iterations required for the conjugate gradient to converge was compared to a prespecified threshold using preconditioners constructed using either the SRHT or the unmodified algorithm S2, using a variety of different values ​​for L. In all cases, the preconditioners constructed using the SRHT achieved the same or greater acceleration, and indeed, this was consistently true throughout the study. For the time being, both the modified and unmodified preconditioner construction algorithms are offered as options for users to select in xGPR, with the latter becoming the default over the SRHT-based algorithm in light of its superior speed and experimental performance.

[0176] Additionally, the preconditioner can be further modified to improve its performance. Step 9 of Algorithm S2 calculates the product Z T (ZΩ) is formed. The matrix Q is recovered from the QR decomposition of this product and used to compute the second product Z T(ZQ) can be formed. This product can then be used in algorithm S2 instead of matrix Q from step 9 (this process can be repeated three more times, but this results in negligible improvement in preconditioner performance). The resulting preconditioner reduces the number of iterations required for the conjugate gradient to converge in all studies by 20-25%, as shown in Figure 5. However, it is more expensive to build because it requires two passes through the dataset. Whether this is beneficial depends on the size of the dataset and the available hardware. Users can select "srht_2" as an option for preconditioner construction in xGPR, and guidance on when to select this option is provided in the documentation.

[0177] Regardless of which alternative you use to construct it, the inverse of the preconditioner matrix is M -1 =(β L +λ 2 )U(Λ+λ 2 I) -1 U T +(I-UU T ) It takes the form of:

[0178] In the formula, β L is the smallest eigenvalue of the preconditioner, and λ is a hyperparameter of the GP kernel. -1 is not clearly formed, but since all the necessary values ​​above have been generated, the matrix-vector product M -1 We can supply v on demand, which is all that is needed for preconditioned conjugate gradients or preconditioned stochastic gradient descent.

[0179] The larger L, the better the approximation; if β1 ≥ β2 ≥ β3... are eigenvalues ​​of A, the eigenvalues ​​decay more rapidly and β j≧L The smaller the eigenvalues ​​for UΛU, the better the approximation. T Rank L Anystrom Let the preconditioner P be

[0180]

number

[0181]

number

[0182] How to choose L? For a given L

[0183]

number

[0184]

number

[0185]

number

[0186]

number

[0187] The relationship between the value of L and the number of iterations required for fitting using CG is shown in Figure 1B below.

[0188] While preconditioners allow for fast model fitting, they do not address the problem of hyperparameter tuning. Computing the marginal likelihood of a Gaussian process requires computing the logarithm determinant of an MxM matrix over M random features, which is very expensive when M is large. Therefore, we developed a method for estimating the marginal likelihood over a large number of random features.

[0189] S5. Logarithmic Determinant Approximation and Conjugate Gradients While performance on the validation set can be used to tune hyperparameters for the GP (as for any other model), and this approach is implemented as an option in xGPR, hyperparameters can also be tuned by maximizing the marginal likelihood. This procedure is more robust to overfitting and also means that a validation set is not required, so we use this approach when possible. Nevertheless, tuning using marginal likelihoods presents one particular challenge. For a Gaussian process, for the negative marginal log-likelihood,

[0190]

number

[0191] In the formula, K y =K(X,X)+λ 2 I. In the random feature scheme, the approximate Z T By replacing Z with K(X,X) and rearranging, we get

[0192]

number

[0193] As shown in S4 above, we use the preconditioned conjugate gradient to find the matrix (Z T Z+λ 2 Without explicitly forming I, (Z T Z+λ 2 I) -1 Z T y can be calculated efficiently. However, when the number of random features is large, log|Z T Z+λ 2 |I| is more difficult to compute. Under Section S6, a very efficient method is developed for tuning hyperparameters using a small number of random features, e.g., less than 4000, as is often the case in practice. However, we first T Z+λ 2 Without explicitly forming the I matrix, we can approximate log|Z T Z+λ 2 We introduce a method that uses preconditioned conjugate gradients for |I|. The same approach is used for fitting to tune hyperparameters in an efficient way whenever the number of random features is large.

[0194] Various techniques have been outlined in the literature to approximate the logarithmic determinant of a large matrix A, including stochastic Lanczos quadrature, stochastic trace estimation combined with Chebyshev polynomials, and subspace iteration. A small subset of A's eigenvalues ​​dominates, i.e., the eigenvalues ​​λ are not always true during fitting. i >>λ i+1 ,λ i+2 ...λ N If there is a large eigengap for , then only subspace iteration is appropriate.

[0195] Stochastic Lanczos quadrature (SLQ) is particularly attractive in this context because, for a large number of random features, it is possible to fit models using preconditioned conjugate gradients and automatically generate the tridiagonal matrices required for stochastic Lanczos quadrature during the process of conjugate gradient optimization at negligible additional cost. The use of this technique to accurately estimate the logarithmic determinant for exact (non-approximate) Gaussian processes using an incomplete Cholesky decomposition to construct preconditioners has already been demonstrated.

[0196] In the following algorithm S5, we use the random feature approximation and the logarithmic determinant term,

[0197]

number

[0198]

number

[0199] Previously, we first rewrote the logarithmic determinant and used Hutchison's trace estimator 21 Used:

[0200]

number

[0201] In the formula, p iis the probe vector drawn from (typically) either a Rademacher or Gaussian distribution.

[0202] Then, if we use a preconditioned conjugate gradient with a preconditioner P, it is in fact the preconditioned system log|P -1 / 2 Z T Z+λ 2 |P -1 / 2 The goal is to solve |. We can adopt the same approach as for exact (non-approximate) GP, which is log|Z T Z+λ 2 I|=log|P -1 / 2 (Z T Z+λ 2 I)P -1 / 2 |+log|P| log|P -1 / 2 (Z T Z+λ 2 I)P -1 / 2 |=trace(P -1 / 2 log(Z T Z+λ 2 I)P -1 / 2 )≒

[0203]

number

[0204] Therefore, a preconditioner needs to be applied to the probe vectors. This is easily done if the probe vectors are drawn from a Gaussian; we can simply draw them from a normal distribution whose covariance matrix is ​​the preconditioner. These preconditioner probe vectors are then denoted by c i It is called.

[0205] The term is

[0206]

number

[0207] This is done using a quadrature-based procedure, but using a preconditioned conjugate gradient algorithm; by saving the alpha and beta coefficients produced at each iteration of the conjugate gradient, we generated a tridiagonal matrix constructed by the Lanczos algorithm. In this case, we simply multiply the target vector by the performance penalty (Z T Z+λ 2 I) -1 Z T Z such that y can be calculated T y vector and n ν Tridiagonal matrices T1,T2...T nν A conjugate gradient needs to be performed using both the preprocessed probe vector and the vector σ such that σ can be generated. A simplified version of the entire procedure is outlined in Algorithm S5 below.

[0208] [Table 9-1] [Table 9-2]

[0209] To evaluate the accuracy of this procedure, we randomly selected eight hyperparameter combinations for seven different datasets, each ranging in size from a few hundred data points to 500,000 data points. For each, we calculated the negative marginal log-likelihood with 4,096 random features using algorithm S5 or matrix decomposition, using a preconditioner rank of 512, a convergence threshold of 1e-5, and the "srht_2" preconditioner construction algorithm in all cases. We used either 40 or 25 probe vectors; the distribution of absolute percent error is plotted as shown in Figure 7, showing that 40 probe vectors do not offer a substantial improvement over 25 probe vectors. 25 probe vectors are sufficient in this study to achieve high accuracy.

[0210] To further evaluate the accuracy of this procedure, we used approximate marginal likelihood combined with Bayesian optimization to tune hyperparameters for two tutorial / experimental examples in user documentation. The first such experiment involves the UCI protein dataset, and the second involves the QM9 dataset with one-hot encoded inputs. In both cases, the tuning procedure improves on the performance achieved by tuning fewer random features with a matrix decomposition-based approach in the best validation set. We further used the same approach to tune hyperparameters for tabular data and small molecule experiments on the QM9 dataset described above, achieving highly competitive performance. These results are consistent with the data shown in Figure 7 and suggest that for the dataset considered here, given a sufficient preconditioner rank, a sufficiently small convergence threshold, and a sufficient number of probe vectors, algorithm S5 can estimate negative marginal log-likelihood with high accuracy.

[0211] S6. Improved strategies for hyperparameter tuning S5 above shows how to use preprocessed conjugate gradients to compute marginal likelihoods in an efficient way when the number of random features is large. However, for noisy data, it is often necessary to use only 1000-3000 random features to tune the hyperparameters and obtain results that yield sufficient performance. If one actually decides to use fewer random features for tuning than for fitting, a fast and highly efficient strategy can be introduced. The marginal likelihoods for random feature approximations are

[0212]

number

[0213] Z T Z is orthonormal as long as it has full rank, i.e., P-P T =I, the eigenvalue decomposition PΛPT is a Hermitian matrix such that it can generate: Then the above equation can be rearranged:

[0214]

number

[0215] For 3 and 4 hyperparameter kernels, this simplification can be used in various ways, for example, in Algorithm S3 and Algorithm S4.

[0216] [Table 10]

[0217] [Table 11]

[0218] Note that in S4 we used Thompson sampling, but other acquisition functions (e.g., expected improvement) could just as easily be substituted. A simple grid search over kernel-specific hyperparameters may also be used, and is also implemented in xGPR. The above algorithm can be further refined by taking the singular value decomposition of Z, which is the most stable approach. T It can also be modified by computing the eigenvalue decomposition of Z.

[0219] The polynomial kernel has only two hyperparameters (the lambda or noise and beta or amplitude hyperparameters shared with all other kernels), which can be tuned using the above strategy using a single pass through the data, making hyperparameter tuning with polynomial kernels exceptionally fast.

[0220] S7. Datasets used for the experiments described in the text All datasets used in this example are publicly available. Details of each dataset are listed in the table below.

[0221] [Table 12]

[0222] Comparing S8.xGPR with Stochastic Variational Inference and Exact Gaussian Processes Stochastic variational approximation to Gaussian process models and variants has been widely found in the literature. From a theoretical perspective, it has many attractive properties, but how does it perform in practice? As shown in Table 8, we compared the performance of SVGP implemented in the GPyTorch library with the random feature approximation implemented in xGPR; for sufficiently small datasets, we also compared it with an exact Gaussian process implemented in GPyTorch. (Training an exact Gaussian process on larger datasets resulted in an out-of-memory error for the UCI protein dataset, which has approximately 8,000 data points, making it clearly infeasible for even larger datasets.) Because we did not find a convolutional kernel compatible with SVGP on sequence data that could be used on a dataset of over 100,000 sequences in a short period of time, nor did such a kernel implemented in GPyTorch, we trained stochastic variational inference in all cases for 40 epochs in GPyTorch, using the RBF kernel in all cases. Training for 20 epochs resulted in lower performance in all cases.

[0223] For xGPR, we used the "fast" strategy in algorithm S3 described in S6 above to tune hyperparameters using 3,000 random features, followed by fitting using 16,384 or 32,768 random features for the tabular data and 8,192 or 16,384 random features for the sequence data (fewer random features were used for the sequence data to match the number of features used for comparison with the CNN deep learning model on the sequence data). Note that several datasets used as benchmarks here are protein sequence datasets that are also used as benchmarks for comparison with CNN. We used a sequence kernel, which does not require the inputs to be aligned, but we used the RBF kernel because sequence kernels are not available in GPyTorch or scikit-learn. Thus, for the protein sequence datasets in this comparison, we converted each into a multiple sequence alignment and one-hot encoded it before using it as input to the RBF kernel, in contrast to how the sequence dataset is modeled.

[0224] The data in Table 8 clearly show that xGPR outperforms SVGP on eight of the nine datasets, often by a large margin. Where SVGP is competitive (Song UCI), they achieve essentially equivalent performance, with SVGP taking 80 times longer to train. xGPR generally performs better using the same kernel and at a greatly reduced cost; furthermore, performance relative to xGPR can often be further improved by using more random features for fitting or by fine-tuning hyperparameters with a larger number of random features. On the Kin40K dataset, "fine-tuning" hyperparameters using, for example, 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 for fitting). In contrast, performance improvements achieved by SVGP are difficult to achieve when the number of induction points for SVGP is increased, as training times for larger datasets become unacceptable.

[0225] While performance with respect to SVGP can be improved by modifying the Adam optimizer settings for stochastic gradient descent, e.g., by dataset-specific tuning of the Adam algorithm's hyperparameters, this is in some sense just a further drawback with respect to SVGP: xGPR does not require expensive experimentation with optimizer settings to achieve good outcomes. Given these results, we did not implement stochastic variational inference for our library, opting instead for the techniques described above.

[0226] [Table 13]

[0227] S9. Lemma 1.1 Lemma 1.1.

[0228]

number

[0229]

number

[0230]

number

[0231]

number

[0232] evidence:

[0233]

number

[0234]

number

[0235]

number

[0236]

number

[0237]

number

[0238]

number

[0239]

number

[0240]

number

[0241]

number

[0242]

number

[0243] In other words, this is

[0244]

number

[0245]

number

[0246] More information on S10.QM9 molecular modeling The complete QM9 dataset consists of approximately 133,000 molecules. First, we removed those molecules whose shape-optimized structures failed consistency checks. We then randomly split the remaining data into a training set of 110,000, a validation set of 10,000, and a test set of 10,831 using a random seed.

[0247] For the SOAP descriptors implemented in the described library, there are several hyperparameters that were tuned: 1) the standard deviation (sigma) of the Gaussian used to represent the neighboring atoms of the central atom, 2) Rcut, or radius, below which atoms are considered "outside the neighborhood" and are no longer included in the representation; 3) Weighting, i.e., decreasing weighting for atoms farther from the central atom (several different schemes are provided).

[0248] There is also the number of radial basis functions, n_max, and the number of spherical harmonics, I_max. These not only improve the accuracy of the representation when increased, but also dramatically increase the size of the descriptor vector. Since larger descriptor vector sizes were not used in this study, we decided on the same values ​​of 12 for n_max and 9 for I_max. However, even with this configuration, we generate roughly 18,000 features per atom.

[0249] The weighting scheme and identified weights are determined by the following formula:

[0250]

number

[0251] For all SOAP feature hyperparameters, the values ​​selected by Willatt et al. based on their grid search were used, with two exceptions: rcut, or cutoff distance, and r0, derived from the weighting scheme presented above, were optimized using the values ​​derived from Willatt et al. as a starting point. r0 and rcut were tuned by evaluating the marginal likelihood on a randomly selected subset of the training data consisting of 25,000 data points. It was found that by first generating features for only heavy atoms (i.e., non-hydrogens) and excluding hydrogens, nearly the same accuracy could be achieved. This is advantageous because the disk space required to store the training set data would otherwise be very 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

[0252] M and r0 determine how quickly the weight drops off from the central atom. SOAP features were generated for each atom separately; each molecule in this dataset has a maximum of nine heavy atoms, and therefore each molecule has a maximum of nine descriptor vectors corresponding to these atoms. These can be used as input to the graph convolution kernel or graph polynomial kernel described above. The SOAP descriptor for each atom was divided by the norm of the descriptor vector so that each descriptor vector has unit norm; this ensures that when using the GraphPoly kernel, it has the same interpretation as the original SOAP kernel.

[0253] To tune the kernel hyperparameters, we used the workflow described under S5 and S6 above. It began by finding a starting point using marginal likelihood based on matrix decomposition with a small number of random features, and 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 298 K as the target. The same hyperparameters that performed well for internal energy at 298 K also performed well for other targets as shown. The final model was then fitted using 16,384 random features, 32,768 random features, or 66,536 random features. The results are shown in Table 4 above.

[0254] Although small improvements in performance can be achieved by using a larger number of random features, 16,384 random features are sufficient to achieve a mean absolute error five times smaller than the chemical accuracy. As a result, in practice, using a larger number of random features has limited usefulness, and using 16,384 random features for this problem is perfectly fine. Furthermore, the most important target for further improvement is likely the SOAP features themselves. When iterating over a dataset, generating a large number of SOAP features per atom significantly increases the computational cost, since many features must be generated, saved, and loaded from disk while processing each mini-batch of data. Furthermore, it inhibits the ability to improve the model's accuracy by increasing the number of basis functions. Several authors have suggested more compact representations that achieve similar or better accuracy. A more compact representation reduces its computational cost by decreasing the length of the feature vector, while increasing the number of basis functions and thereby improving the model's accuracy.

[0255] It is interesting to note here the advantages offered by the Conv1d kernel based on the Fast Hadamard Transform (FHTConv1d). If vanilla random features were used instead to generate 66,536 random features from 18,000 input features, it would require a matrix with roughly 1.2 billion floating-point numbers. In contrast, the FHTConv-1d kernel requires storing roughly 288,000 8-bit integers and 72,000 floating-point numbers, a dramatic reduction in memory footprint.

[0256] Clustering and Kernel PCA using representations produced by S11.xGPR Here, we used the QM9 dataset. 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 encoding up to the 15th furthest neighbor (an arbitrary cutoff). Each one-hot encoded neighbor was then scaled by 1 / distance. ** The weighting is done by 6; the 6th power used here corresponds to the weighting used in the London matrix, which is sometimes used as a representation for machine learning for small molecules. Following the standard workflow, hyperparameters were first adjusted using a smaller number of random features (2048) with a matrix decomposition-based procedure, and 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 used for fitting, this model can achieve a mean absolute error of less than 1.1 kcal / mol.

[0257] Next, we retrieved the top 500 principal components using the kernel_xpca tool provided in the visualization_toolkit module of xGPR version 0.0.2.1. We used the top two principal components to generate a kernel PCA plot, shown in Figure 8, for the test set data only.

[0258] As illustrated, the top principal components correlate well with the predicted ground truth labels (energy at 298K). This is not always the case, as a large amount of information is discarded by using a 2D representation instead of the 16,384-dimensional representation generated by xGPR. To further explore the data, we then performed k-means clustering using the scikit-learn library with the top 500 principal components as input, using the sum of squared error to evaluate the performance for a given number of clusters. This operation approximates performing kernel k-means clustering on the original training data using the GraphConv-1d kernel.

[0259] The elbow plot is shown in Figure 9A; the results suggest that five clusters may be appropriate, and thus the final k-means model was fitted using the five clusters. The distribution of energies at 298 K in each resulting cluster is shown in Figure 9B. In many drug discovery scenarios, this ability to cluster data using the same model used to fit it can be useful in itself, as it can be used for the efficient recovery of similar molecules from the training discussed above.

[0260] S12. A molecule from the test set and the most similar molecule from the training set. Five random molecules from the QM9 test set, along with the 50 most similar molecules from the training set, are depicted in Figures 10-14, where similarity is quantified using the (approximated) GraphRBF kernel function for molecules represented using one-hot encoding as described above. Note that this representation of the molecules is relatively useless, and therefore it is somewhat surprising that the model works, as well as that it does so with an MAE of about 1 kcal / mol and is able to recover molecules that exhibit many similarities in many cases.

[0261] conclusion Gaussian processes are a sensible Bayesian approach to machine learning that provides easy quantification of uncertainty. Their most significant drawbacks are their unacceptable scaling and lack of efficient kernels for sequences and graphs. This example describes efficient methods for fitting GP models using kernels approximated by random features that achieve linear scaling with training set size and linear scaling with molecule size or sequence length. These methods were implemented in a publicly available Python library, xGPR, which was proven to achieve competitive performance for a variety of protein and small molecule property prediction tasks.

[0262] As discussed herein, GPs equipped with these kernels produce well-calibrated uncertainty estimates for their predictions, unlike either of the deep learning models used as a benchmark. This quality is of considerable importance for protein engineering and drug discovery; evaluation of predictions is often very expensive, so it is certainly preferable to use high-confidence predictions when possible.

[0263] Furthermore, GPs equipped with these kernels are substantially more interpretable than the baseline deep learning models. For protein sequences and small molecules alike, each kernel can provide a detailed description of how the model determines the similarity between any two data points and then makes a prediction. This allows for understanding the limitations of the model and how to improve it. It is impossible to generate this kind of detailed description for any of the currently popular deep learning architectures for sequence or graph data; they are all "black boxes." Furthermore, as shown for small molecules, kernels derived from trained GP regression models can be used to cluster input data or perform kPCA, which provides further insight into the distribution of the input data. This feature is particularly useful for drug design to classify small molecules into subgroups.

[0264] This combination of competitive accuracy, uncertainty quantification and improved interpretability suggests that xGPR libraries are potentially useful for a range of protein and small molecule property prediction tasks.

[0265] All kernels discussed in this example accept a single protein or small molecule as input. For some tasks (e.g., protein-small molecule interactions), it is desirable to use a "pairwise" kernel that accepts two inputs (e.g., a small molecule and a protein sequence) and compares them with other pairs. The kernels discussed herein can be easily modified to generate several "pairwise" kernels that may be appropriate for some problems. For example, the first layer of the FastConv-1d kernel discussed herein can be used to generate representations of both the small molecule and the protein; these representations can then be concatenated and used as input to the RBF kernel. Alternatively, the GraphRBF kernel can be used to generate random features for the small molecule, while the FHTConv-1d or FastConv-1d kernel generates random features for the protein; the random features are then concatenated in a single vector to generate a prediction. This last arrangement is consistent with using the sum of two kernels. It is likely that other efficient kernels can be designed that are specific to this type of pairwise comparison. [Example 2]

[0266] The Rapid Processing System for Proteins (RESP) AI model accelerates the identification of tight-binding antibodies This example describes the development of an easy-to-implement machine-learning-assisted pipeline for identifying high-affinity antibodies that addresses these challenges. A novel autoencoder model was trained on over 3 million B cell receptor sequences, demonstrating that this learned representation yields better results for the desired task than state-of-the-art embedding schemes. Next, a Bayesian neural network was developed and trained to perform ordinal regression to model the relationship between sequence and binding affinity or off-rate using directed evolution data. Importantly, this model provides an estimate of the uncertainty in its predictions, which cannot be achieved by current deep learning methods. To benchmark this approach, it achieves competitive accuracy on a literature dataset while providing useful uncertainty information not provided by deep learning models in original related studies. Furthermore, in silico mutagenesis was performed using a simulated annealing strategy to explore sequences not present in the mutation library and assess their binding affinity 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.

[0267] To demonstrate the power of RESP, we attempted to improve the affinity of a well-known antibody for a well-studied antigen by engineering mutants of the heavy chain of atezolizumab (trade name Tecentriq) with improved affinity for programmed cell death ligand 1 (PD-L1). We randomized a large portion of the atezolizumab heavy chain sequence, binned the mutants using yeast display and FACS, and sampled each bin for sequencing. We converted the sequence data of mutant atezolizumab to a low-dimensional representation using an autoencoder model, and used the encoded mutant atezolizumab sequences to train a Bayesian ordinal regression model that scored each sequence for the probability that it was a strong binder to PD-L1. Furthermore, we used a modified simulated annealing algorithm to select sequences for testing. This pipeline discovered a panel of atezolizumab scFv mutants with improved off-rates for PD-L1, with one characterized mutant demonstrating an approximately 10-fold reduction in off-rate and a K value for atezolizumab interaction with human PD-L1. D It is shown here that the antibody-antigen interaction showed a 17-fold improvement over the previous study. The method described herein is useful as a general approach to improving antibody-antigen interactions while reducing experimental efforts to do so. Mutants discovered by this method may be useful reagents for treating PD-L1-positive tumors.

[0268] method software Analysis and modeling were performed using Python 3.9 with the PyTorch library version 1.8.1, Numpy library version 1.19.5, Scipy library version 1.5.4, and scikit-learn library version 0.24.2. 0 was used as a random seed for initializing model weights, train-test splits, etc.

[0269] Autoencoder model for antibody encoding The training set used for the autoencoder was derived from the cAb-Rep database, which consists 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, a deep repertoire dataset was used. This dataset was further filtered by numbering all sequences using ANARCI software using Chothia's numbering system. During this process, ANARCI software aligns the input sequences to domain-specific hidden Markov model profiles for human antibodies using HMMer software. Sequences with low bit scores resulting from this alignment were rejected, thereby minimizing the possibility that the dataset contains proteins that are not actually antibodies. After numbering and filtering using ANARCI, 2,725,492 sequences remained. These sequences were derived from a database of antibody sequences and filtered to select those with a high probability of being generated by the antibody MSA, and therefore are antibodies. Note that the sequences used here are only the VH heavy chains that match the sequences used in the experimental work.

[0270] The autoencoder accepts one-hot encoded sequences as input. We used the Chothia numbering scheme to number the positions of each antibody. The heavy chain may have 132 amino acids in this scheme, but our mutants have fewer amino acids due to the size of the complementarity-determining regions (CDRs) varying between antibodies. As a result, each one-hot encoded sequence is a matrix with 132 rows and 21 columns. Because the heavy chain does not contain the full 132 amino acids and therefore has blanks at some Chothia numbering positions, only 20 amino acids are present, but the 21st position indicates a blank. While it is also possible to simply fill blank positions with all zeros, it is preferable to leave explicitly encoded blanks as they are for clarity, and this approach is adopted here. Each row contains a 1 at one position, indicating the presence of an amino acid (or a blank, if no amino acid is present).

[0271] In sequence selection and design, it is often important to determine which positions are most important for affinity or other desired properties. This consideration suggests that the training representations generated by an autoencoder should have the same number of rows (positions) as the input. Furthermore, to ensure that a particular row of encoding contains information associated with that position, the decoder should be able to reconstruct the amino acids present at each position using information from that position (and possibly neighboring positions).

[0272] To accommodate 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 and a third linear layer with appropriate padding. The first convolutional layer contains 40 kernels, each 21 wide, and the second contains 20 kernels, each 11 wide. Each convolutional layer uses gated convolutional activation to improve the performance of convolutional neural networks on language modeling tasks. For this activation function, a sigmoid function is applied to the first half of the columns in the output from that layer, and then these are multiplied element-wise by the second half. In other words, if the output of a layer for a given sequence is a 132x40 matrix, for gated convolutional activation, a sigmoid function is applied to the first 20 columns, and then these are multiplied element-wise by the remaining columns to obtain a final output of dimension 132x20.

[0273] The final output of the linear layer in an autoencoder is a 132x3 matrix, which is an encoding of the original sequence. This encoding is the input to both the decoder and prediction modules. The decoder module consists of a single layer of the following form: (1) softmax(a·w+b)

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

[0275]

number

[0276] where z is a 21-element vector obtained 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.

[0277] The output of the encoder is also fed into a prediction module that distinguishes human B cell receptor sequences from junk sequences. The prediction module consists of a simple logistic regression model, where the probability that a sequence is an antibody is given by:

[0278]

number

[0279] where a is the full output of the encoder flattened from a 132x3 matrix to a vector of length 396, and w is the learned weight vector of length 396 plus the learned bias term. Because both the predictor and decoder modules have few parameters (compared to typical deep learning models), neither can learn a complex mapping from the encoded sequence to the input, thereby burdening the encoder to produce as useful and meaningful an encoding as possible.

[0280] To generate junk sequences, we generated mutant versions of each of the 2,725,492 sequences selected from the cAb-Rep database. The goal here is to force the model to incorporate information about the relative abundance of specific amino acids at specific positions during its encoding by requiring it to be able to distinguish true antibody sequences from closely related sequences. While random sequences can be easily generated, this does not force the model to generate useful encodings; random sequences are quite different from antibodies, which a logistic regression model that predicts whether a sequence is an antibody or not can distinguish without any modifications to the representation generated by the encoder. Therefore, it is preferable to generate decoy sequences that are just similar enough to make them difficult to distinguish, but different enough to exhibit modifications that are rare in true antibody sequences.

[0281] Experiments suggested that seven mutations would provide a good balance; consequently, each copy was changed to a randomly selected 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 were human B cell receptors and the other half were not. An autoencoder model was thus trained to encode the input sequences and 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).

[0282] Accuracy was assessed separately for both the prediction and reconstruction tasks using a holdout test set. These metrics assess the autoencoder's ability to reconstruct its inputs, so they are used only as diagnostics. The true test of an autoencoder is the extent to which the encodings it produces affect the predictive accuracy of a model trained using its learned representations as input. To evaluate this more important metric, WT and mutant atezolizumab libraries (constructs 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 encoding scheme, and the AntiBertY embedding. For AbLang, the seq encoding option was used. For the AntiBertY embedding, residue-specific representations in the sequence were averaged across all sequences. For AntiBertY, all residue-specific embeddings were initially used without averaging across them, which yielded good 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. The classification performance of each model on each encoding was evaluated using 5x cross-validation on the training set and holdout test set, along with the Matthews correlation coefficient and AUC-ROC for RH03 vs. remaining identification as metrics. A good encoding should improve, or at least not impair, performance compared to one-hot encoding across all three model types, especially the Bayesian network, which is the most important model for this study.

[0283] Generation of WT and mutant atezolizumab scFv libraries The WT atezolizumab scFv was first cloned into the pYD1 yeast display vector to test its function on the yeast surface (see Figures 24A-24B). For the first atezolizumab library, the WT plasmid was used as a template for PCR to prepare either the WT light chain using Q5 hot-start DNA polymerase (NEB) or the mutant heavy chain by error-prone PCR using Taq polymerase (Invitrogen) as previously described. The two PCR products were assembled into one product by overlap-extension PCR and co-electroporated into EBY100 yeast along with the linearized pYD1 vector as described. Based on the number of colonies after serial dilutions on selective plates, library transformation yielded 7.8 x 10 7 This resulted in transformants.

[0284] Atezolizumab scFv library screening by yeast surface display The optimal competition time for off-rate screening was determined as previously described. WT or mutant library yeast were thawed, inoculated into selective growth medium, and grown at 30°C for 22 hours. The library / WT was induced in galactose induction medium at 20°C for 42 hours. Subsequently, the WT or library was labeled with biotin-PD-L1 (Sino Biological 10084-H08H-B) in TBS-BSA at RT for 3 hours, and then incubated with excess non-biotinylated PD-L1 (Sino Biological 10084-H08H) for the designated competition period (at RT in TBS-BSA). For FACS, yeast were labeled with anti-V5 (Thermo Fisher R960-25, formerly known as 46-0705) in TBS-BSA at 4°C, followed by SA-PE (BD 554061) and goat anti-mouse IgG2a AF647 (Thermo Fisher A21241) in TBS-BSA for 30 min on ice. Cells were sorted for faster, intermediate, and slower off-rates (Figure 16A). 1-2and Figure 34). Hits were grown to high density at 30°C and made into frozen stocks at -80°C.

[0285] Preparation of DNA libraries of mutants with faster, WT-level, and slower off-rates for PD-L1 To isolate plasmids from each binding group from the initial library screening, the plasmids were harvested from yeast using the Zymoprep Yeast Plasmid Miniprep II Kit (Zymo Research), eluted with ddH2O, and concentrated using the DCC-5 (DNA Clean and Concentrator 5) Kit (Zymo Research), followed by a first round of PCR using primers SeqF / SeqR (Table 11, designed to amplify only the heavy chain region). The first round of PCR was performed with Q5 Hot Start DNA Polymerase (NEB) using a moderate number of PCR cycles (17, determined to be optimal by qPCR) to avoid overamplification of the library, and the PCR products were concentrated using the DCC-5 Kit and purified by agarose gel extraction. A second PCR was performed using various primers (Low / Medium / High Binder NGS F, NGS R) to barcode each binding group (see Table 11) for MiSeq Nano PE250 using KAPA HiFi Hotstart Ready Mix (KK2601) for 5 cycles (determined optimal by qPCR). PCR products were purified using Ampure XP beads (A63880) and submitted for QC using TapeStation analysis (Agilent HS D1000) before being submitted for MiSeq sequencing at the Institute of Genomic Medicine, UC San Diego.

[0286] Generation and screening of a focused library of 21 mutants (see also S2.4 below) The 21 geneblock fragments (IDT) for the 21 mutant heavy chains were fused with the light chains by overlap extension PCR (Q5 hot start) and then the first library (10 7 The first library, as well as the first transformants (>100 transformants), was co-electroporated with the linearized pYD1 vector into EBY100 yeast. The resulting library was screened essentially the same as the first library, except that the final sorting included a highly stringent competition period (39 hours at RT). Hits were harvested by yeast plasmid miniprep, transformed into GC10 competent cells, harvested by bacterial miniprep, and sequenced using standard Sanger sequencing.

[0287] Cloning of WT atezolizumab, mutant 4, durvalumab, and avelumab scFv into pYD1 vector Geneblocks for durvalumab and avelumab were purchased from IDT with yeast-optimized codons and PCR-amplified using the following primers: Atez LE F & Durv LE R for durvalumab, and Atez LE F & Avel LE R for avelumab (Table 11). WT / mutant 4 genes were PCR-amplified from their plasmids using Atez LE F / R primers. The PCR products were double-digested with XhoI / NheI-HF (NEB) and ligated into linearized / dephosphorylated pYD1 (digested with the same enzymes and dephosphorylated / ligated with rSAP / T4 DNA ligase (NEB)). The ligation products were transformed into GC10 competent cells (42-658, Genesee Scientific). Plasmids were isolated using the Zippy Plasmid Miniprep Kit (Zymo Research) and sequence-verified before transformation into EBY100 yeast. These constructs were then transformed into yeast kappa cells (Kappa Cells, 2019). off and K. D was used for the determination of

[0288] K of WT / Mutant 4 / Durvalumab / Avelumab on the surface of yeast D and k off Determining Values K off was determined at RT in TBS-BSA essentially as described in section S2.3 below. The resulting data were analyzed using the following equation: Y = 1 (fraction of yeast bound to biotin-PDL1 at the 0 second competition time point), where Y is the fraction of yeast bound to biotin-PDL1, X is time, Y = 1 (fraction of yeast bound to biotin-PDL1 at the 0 second competition time point), plateau is a constant based on nonspecifically bound yeast, and K is K. off (This is): (4) Y = (Y0 - plateau) * exp(-K * X) + plateau was fitted with a one-phase decay model using GraphPad Prism 9.3.0 software.

[0289] K on the yeast surface D The values ​​are basically as described 49 Determine the data using the following formula:

[0290]

number

[0291] Cloning and purification of WT / mutant 4 atezolizumab scFv (see also S2.6 below) The scFv sequence was PCR amplified from the pYD1 vector and cloned into the pET27b(+) (69863-3, MilliporeSigma) vector for bacterial expression. The vector was transformed into Rosetta (DE3) cells (70954-3, MilliporeSigma), and the scFv-His6 fusion was purified essentially as described. The scFv was frozen in PBS + 10% glycerol at -80°C. Purity was assessed using SDS-PAGE (Figures 28A and 28B), and protein concentration was determined by BCA assay (Pierce #23227).

[0292] Biolayer Interferometry (BLI) D Determining Values BLI was performed on Octet Red 96 at 25°C at The Scripps Research Institute's Biophysics and Biochemistry Core. Binding reactions were performed in 1X kinetic buffer (Sartorius, 18-1105) consisting of 20 mM phosphate buffer, pH 7.6, 2 mM KCl, 150 mM NaCl, and 0.02% Tween 20, 0.1% BSA, and 0.05% sodium azide. Biotinylated PDL1 (Sino Biological 10084-H08H-B) was immobilized on the sensor by soaking SA Biosensors (Sartorius, 18-5019) in 100 nM b-PDL1 until the signal was saturated. A blank loading channel was used as a negative control. Kinetic experiments were performed using both scFv WT and scFv mutant 4. Experiments were performed using a concentration gradient of scFv WT ranging from 475 nM to 0.6 nM and scFv mutants ranging from 160 nM to 0.2 nM. Binding kinetics were evaluated by Octet Data Analysis HT Software Version 12 using a 2:2 binding model.

[0293] mass photometry Mass spectrophotometric analysis was performed at room temperature on a Refeyn TwoMP at Scripps Research Institute. The PD-L1 used was biotin-PD-L1 (SinoBiological 10084-H08H-B) and purified scFv of wild-type or mutant 4. Each characterization was performed in 1X PBS (phosphate buffer from Cytiva BR100672). For each experiment, the scFv or PD-L1 (or both) was diluted to the final concentration in PBS. Each experiment generated a 60-second movie, and mass analysis was performed using Refeyn DiscoverMP analysis software v2.3.0. Mass calibration was performed using urease (Sigma, U7752) and thyroglobulin (Millipore, 609310).

[0294] Sequence Processing The raw paired-end reads from the atezolizumab dataset were checked for quality (see the Sequence Quality section for details on filtering criteria). After processing the sequences, they were divided into an 80% training set and a 20% test set. All training and cross-validation were performed on only the 80% training set.

[0295] In many cases, a sequence was present in more than one category, but at different frequencies. Clearly, if the frequency of a sequence in its assigned category is not significantly higher than its frequency in other categories, the level of confidence in the category assignment is low. To encode the level of confidence, each sequence was weighted by its frequency in its assigned category + 1 divided by the total number of occurrences + 3. This is a posterior probability using a multinomial likelihood and a uniform Dirichlet prior (Dirichlet distribution with α = [1,1,1]).

[0296] The ordinal regression model used here is the binary cross-entropy loss:

[0297]

number

[0298] where p for data point i is the probability that the model assigned the sequence to belong to a more stringent / tight-binding category than category j, out of M possible categories, and for N data points, and y(i,j) is either 0 or 1, indicating whether the sequence belongs to the tighter-binding category (see below for details on how p is calculated). This loss function for sequence i is multiplied by the weight for sequence i, so that the model is weakly penalized for misclassifying sequences when our confidence in the category assignment is low.

[0299] Building an Ordinal Regression Model The atezolizumab mutant library data was analyzed using ordinal regression based on a Bayesian neural network. Each mutant atezolizumab sequence was first encoded by an autoencoder to obtain a 132x3 matrix, which served as input to the atezolizumab model. The model architecture is shown in Figures 23A-23B. At a high level, it is similar to the Bayesian backpropagation architecture described by Blundell et al. (ArXiv150505424 Cs Stat, 2015), except that this model was adapted to perform ordinal regression as described below.

[0300] In traditional neural networks, each parameter is a learned parameter that is fixed once the model is trained. In variational networks, by contrast, each parameter has an associated Gaussian distribution whose optimal value is described by a mean and standard deviation learned during training. Rather than using fixed learned parameters as in traditional networks, to generate predictions, we sampled from the distribution across the weights, multiplied by N to generate N predicted values. The variance in these predictions provides a measure of our uncertainty around the final prediction. This approach requires a more complex training procedure than that associated with traditional neural networks, discussed below.

[0301] In short, (7)argmin θ KL[q(w|θ||p(w)]-E q(w|θ) [log(p(D|w))] Find the parameters θ for the distribution over the weights w that minimize the evidence lower bound or ELBO given by

[0302] where q(w|θ) is the distribution over the weights that approximates the true Bayesian posterior distribution, p(w) is the prior distribution, and p(D|w) is the likelihood of the data. The parameter θ is calculated using the Bayesian Bayesian model shown in Figure 23A. 1-2 and are the parameters of the normal distribution for all weights in each hidden layer shown in 23B. The first term is a regularization term that measures the discrepancy between the weight distribution learned by the model and the prior distribution, while the second term is the negative log-likelihood. Because some weights are expected to be outliers here, we chose a Cauchy prior with unit scale and zero location, which adds some additional flexibility.

[0303] The cost was approximated by a Monte Carlo sampling procedure. Using the current set of parameters θ (mean and standard deviation of the weight distribution), for each mini-batch, we subtracted n sample sets of weights to obtain the following approximate loss function:

[0304]

number

[0305] where M is the number of mini-batches, and the other terms are as above. In other words, for each mini-batch of training data, we subtracted n sets of sample weights, then averaged an approximate loss function over these in the training set and over the mini-batch. We make backpropagation for this loss function tractable by using Kingma and Welling's reparameterization trick. We sampled a standard normal distribution with mean zero and standard deviation 1, then added the mean and multiplied by the standard deviation of the distribution for weight j. To ensure that the standard deviation is always positive, instead of using the standard deviation itself as a parameter, we parameterize each distribution using a softplus function: (9)σ = log(1 + e ρ ) It is parameterized by the parameter ρ, which is converted to standard deviation using

[0306] Using this reparameterization trick, the gradient of the approximate loss function with respect to the mean and ρ of the distribution for each weight and bias term in each hidden layer is easily calculated. The advantage of this model structure is two-fold. First, it imposes a strong regularization on the model parameters that can result in improved performance for some tasks. Second, it estimates the uncertainty in the predictions, thereby attributing more weight to the most confident predictions when selecting sequences for testing.

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

[0308] This arrangement processes ranked categories, whereby sequences with higher potential scores are assigned to higher categories. Note that it does not predict the actual off-rate or binding affinity. Rather, the score reflects the confidence that the sequence is a strong binder relative to others in the training set. This approach has been used in the literature to rank sequences for protein engineering tasks and select them for experimental evaluation.

[0309] Sequence scoring and selection To perform in silico directed evolution as illustrated in Figure 17, a classical simulated annealing algorithm was modified and equipped with a trained model. First, the frequency of each amino acid at each position across the entire dataset was calculated, and the marginal probability of each amino acid at each position was retrieved by adding 1 to all values ​​and dividing by the total number of sequences in the dataset plus 20. By adding 1 to all frequencies, there is a small, but non-zero, marginal probability of an amino acid not being observed in the dataset.

[0310] 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, as more or fewer sites can be used as needed. In each iteration, there is an equal probability of selecting any of these top 10 sites. The selected sites are randomly reassigned to new amino acids; the probability of selecting any given new amino acid is determined by the marginal probability calculated above. For example, assuming that 80% of all sequences observed in the dataset carry arginine at position 100, another 10% carry lysine, etc. If position 100 is selected for a given iteration, this position is mutated with an 80% chance of being converted to arginine, a 10% chance of being converted to lysine, etc.

[0311] Both the current sequence and the mutated proposed sequence are encoded using an autoencoder, and a score is assigned to each using a trained ordinal regression model.

[0312]

number

[0313] where T is the temperature and S 最良 and S 提唱are the best score to date and the score of the proposed sequence, respectively.

[0314] This algorithm is further illustrated as follows.

[0315] Modified and simulated annealing

[0316] [Table 14]

[0317] This procedure is a simple approach to explore sequence space while ensuring that it does not risk straying too far from the training set. For this stage of modeling, to ensure reproducibility and avoid stochastic fluctuations, the Bayesian neural network generated predictions using the mean of the distribution for each weight as the weight value.

[0318] S2.1 Sequence quality Sequence reads containing one or more base pairs with a Phred quality score of less than 10 were discarded because the sequence read was unreliable. In the event that the overlapping regions of the paired ends did not match, both reads were discarded, so mismatches between the paired ends were not allowed. Reads that met these quality criteria were combined and translated to obtain a mutant atezolizumab sequence of 118 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. The first 30 or last 8 positions were not targeted for mutation, and because mutations at these positions are very rare and therefore difficult to assess their significance, sequences with mutations at these positions were excluded from further consideration.

[0319] S2.2 Generation of WT and mutant atezolizumab scFv libraries To test the function of the WT scFv in a yeast display format, the gene for the 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 11), followed by transformation into EBY100 yeast by heat shock using double digestion / ligation in the digested vector by Golden Gate Assembly using Esp3I (Thermo Scientific) and T4 DNA ligase (New England Biolabs (NEB)), which was sequence verified. For the atezolizumab library, WT scFv in the pYD1 vector was used as a template for high-fidelity PCR of the light chain region with Q5 Hot Start DNA Polymerase (NEB) (primers IF1F / IF1R) or error-prone PCR of the heavy chain region with Taq DNA Polymerase (Invitrogen 18038-018, IF2F / IF2R) using 8-oxodGTP (TriLink N-2034-1) and dPTP (TriLink N-2037-1), except that 30 cycles of error-prone PCR with 200 μM dPTP and 8-oxodGTP (final concentration 20 μM each) were used to increase the mutation rate. The light and heavy chain PCR products were concentrated using a DNA Clean and Concentrator-5 kit (DCC-5, Zymo Research) and purified using an agarose gel extraction kit (Zymo Research Gel Extraction Kit). The two DNA fragments (WT light chain, mutated heavy chain) were assembled into a single DNA molecule by overlap extension PCR using OF / OR primers (all primers in Table 11) with Q5 hot-start DNA polymerase, and the resulting product was concentrated again with the DCC-5 kit and purified by gel extraction.Mutation rates in the heavy chain were examined using TA cloning (Invitrogen TA Cloning Kit, K202020) by adding a 3'A to the product with Taq; random sequencing of clones revealed 0-8 mutations per gene. In preparation for yeast electroporation, the pYD1 vector was double-digested with XhoI and EcoRI-HF, concentrated using the DCC-5 kit, and gel-extracted / eluted in ddH2O.

[0320] The scFv library was transformed into yeast by electroporation (3 μg digested pYD1, 9 μg scFv gene per electroporation) and subjected to homologous recombination. 6 and 7.8 x 10 as determined by serial dilution on plates containing selective growth medium (BD Bacto Agar at 16.7 g / L for solid medium). 7 Transformants were obtained. The library was passaged several times in selective growth medium to ensure 1 plasmid / cell and then frozen in aliquots (3.6 x 10) at -80°C in 85% ddH2O, 10% glycerol, 5% DMSO. 8 As a control for selection, the WT scFv gene / pYD1 was also electroporated into EBY100 to have the exact same location in the vector as the library.

[0321] S2.3 Atezolizumab scFv library screening by yeast surface display First, the k of WT scFv against PD-L1 was measured to determine the optimal competition time. off The off-rate was determined essentially as previously described. Yeast expressing WT scFv were washed with TBS-BSA and diluted with 16x10 6The yeast were resuspended at 2 x 10 cells. Then, 12 μL of b-PDL1 (final concentration 178 nM, Sino Biological 10084-H08H-B) was added to the yeast and labeled for 3 hours at room temperature. The yeast were then washed with TBS-BSA and resuspended in 174 nM unlabeled PDL1 (Sino Biological 10084-H08H) in a total volume of 800 μL. 2 x 10 cells were resuspended at 2 x 10 cells. 6 Cells were added to each Eppendorf tube and rotated at RT for x hours (0–21 hours). For each time point, cells were washed three times with cold TBS-BSA and kept on ice in 198 μL of TBS-BSA until the experiment was completed. Once all data points were collected, a 1:100 V5 antibody (R960-25, formerly known as 46-0705) was added and the yeast were labeled on ice for 1 hour. The yeast were washed with cold TBS-BSA and resuspended in 200 μL of 1:100 SA-PE (BD 554061) and goat anti-mouse AF647 (Thermo Fisher A21241) for 30 minutes on ice. The yeast were then washed with cold TBS-BSA and resuspended in cold TBS-BSA for flow cytometry. The optimal competition time was determined to be 29 hours at RT. To screen the library against PD-L1, WT or library was thawed and added to selective growth medium (20x10 6 Cells were added to 50 mL of growth medium plus ampicillin and grown with shaking at 30°C for 22 hours, after which they were induced by settling the cells, washing them once with sterile ddH2O, adjusting the pH to 5.6, and OD with ampicillin. 600 The scFv were resuspended in the same type of medium but with 20 g / L galactose instead of glucose, adjusted to a 0.2 % RI. The scFv were induced for 42 hours at 20°C with shaking. After 42 hours, the flasks were placed on ice and 40 x 10 6Yeast cells (WT or library) were sedimented at 17,000xG for 30 seconds, resuspended in 1 mL of TBS-BSA (25 mM Tris, 150 mM NaCl, pH 7.5, 5 mg / mL BSA), re-sedimented, and resuspended at 40 million cells / mL in TBS-BSA buffer. The WT or library was then labeled with biotin-PD-L1 (Sino Biological 10084-H08H-B) at a final concentration of 178 nM for 3 hours at room temperature with rotation. The cells were then sedimented / washed with TBS-BSA, then resuspended in non-biotinylated PD-L1 (Sino Biological 10084-H08H, 174 nM) and rotated at room temperature for various amounts of time (mostly 29 hours, but one screen was 7 hours). For FACS sorting, cells were precipitated / washed three times with cold TBS-BSA and labeled with anti-V5 antibody (R960-25 from Thermo Fisher, lot 2001339 / 2249078 / 2212258, verified for specificity by the supplier) at a dilution of 1:100 for 60 min on ice (10 7 cells / mL), precipitated / washed with cold buffer, and then 7Cells were labeled with 1:100 of each secondary detection reagent (PE streptavidin BD 554061, goat anti-mouse IgG2a AF647 (Thermo Fisher A21241) at 1:100 cells / mL for 30 minutes on ice. Cells were then precipitated / washed with cold buffer and analyzed on a BD FACSAria II (Moores Cancer Center, UC San Diego, CA) using BD FACSDiva v.8.0.1 software. The hits were collected in selective growth medium, grown to high density at 30°C, and aliquoted into frozen stocks at -80°C. The process was repeated twice for binders with faster or WT-level off-rates and four times for mutants with slower off-rates.

[0322] S2.4 Generation and screening of a focused library of 21 mutants Twenty-one geneblocks encoding the 21 heavy chain mutants were purchased from IDT and then PCR amplified and subjected to overlap extension PCR (Q5 hot start DNA polymerase) to generate 21 full-length atezolizumab scFv mutant genes. The sequences of the 21 heavy chain mutants are found in Table 12. The genes were electroporated (with linearized pYD1) into EBY100 yeast to assemble a focused library, which was passaged / frozen as described for the initial library. Titration on selective media plates yielded 10 7 More than 10 transformants were shown.

[0323] To screen the libraries for the most improved (slowest off-rate) mutants, increasingly longer competition times were used in the presence of excess unlabeled PD-L1 (WT or mutant libraries were first saturated with 200 nM biotin-PD-L1, followed by sedimentation / washing of cells to remove unbound biotin-PD-L1 and resuspending them in excess (128 nM unlabeled PD-L1) for 7.5 hours (17 hours at RT, then 39 hours). The most potent clones were collected by FACS in growth medium (using WT scFv as a reference) and carried through to the next round of selection. After three rounds of selection, plasmids derived from hits were collected by yeast miniprep and transformed into GC10 competent cells, which were picked as individual colonies and subjected to bacterial miniprep. Plasmids were sequenced using Eton Bioscience sequencing service.

[0324] S2.5 K of WT / Mutant 4 / Durvalumab / Avelumab on the yeast surface D and k off Determining Values RT K of scFv-PD-L1 complex off The off-rates were determined essentially as described using biotin / non-biotin PD-L1 (same version used during library screening) with TBS-BSA as the buffer and the same antibody reagents used for library screening. Yeast were grown / induced in the same manner as the screening protocol, and 4x10 7 Cells / mL were labeled with 200 nM b-PD-L1, precipitated / washed, and then treated with 64 nM unlabeled PD-L1 (10 7 Cells / mL) and rotated at RT. The mean fluorescence intensity (MFI) of V5-positive yeast was recorded at each time point, and the data were fit to a one-phase decay model (Y = (Y - plateau)) using GraphPad Prism 9.3.0 software. * exp(-K *where Y is the fraction of yeast bound to biotin-PDL1, X is time, Y = 1 (the fraction of yeast bound to biotin-PDL1 at the 0 second competition time point), plateau is a constant based on non-specifically bound yeast, and K is Koff). K on the yeast surface D was determined as previously described, and the MFI of the V5+ population was plotted against the antigen concentration and biotin-PD-L1 as used for library screening. D An important note regarding the determination is that 10 5 The data included using cells, anti-V5 antibody at 1:200 after resuspension of each data point in 100 μL of cold TBS-BSA for 1 hour, washing of cells with cold TBS-BSA after V5 mAb incubation, then labeling of cells with SA-PE / IgG-AF647 at 1:100 for 20 minutes on ice, washing of cells in cold TBS-BSA, and resuspension in cold TBS-BSA for flow cytometry. Data were calculated using the following equation: y=B max * X / (K D +X) (Bmax is the maximum MFI value, X is the antigen concentration, and 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 equilibrate.

[0325] S2.6 Cloning and purification of WT / mutant 4 atezolizumab scFv The gene encoding the scFv sequence was amplified by PCR from the pYD1 vector using primers (scFv F / R, Table 11) containing NcoI-HF and XhoI cleavage sites using Q5 Hot Start DNA Polymerase (NEB). The PCR product was double-digested (together with the pET27b(+) vector (69863-3, MilliporeSigma) (the vector was also dephosphorylated with Quick CIP (NEB)). The digested products were ligated using T4 DNA ligase (NEB), and the ligation was transformed into NEB 5-alpha cells. The plasmid was purified using a Zippy Plasmid Miniprep Kit, and the sequence was verified.

[0326] The scFv plasmid was transformed into Rosetta (DE3) cells (70954-3, MilliporeSigma) on LB-0.5% glucose plates containing 30 μg / mL kanamycin, 34 μg / mL chloramphenicol, and grown at OD in LB-0.5% glucose medium supplemented with kanamycin / chloramphenicol (30 μg / mL kanamycin, 34 μg / mL chloramphenicol) (10 mL). 600The cells were grown at 37°C for 5 hours with shaking until the RI was 0.8 to 1. Then, 4 mL of the medium was added to 400 mL of autoinduction medium (73100-4, MilliporeSigma, containing 100 μg / mL kanamycin and 34 μg / mL chloramphenicol), and the antibiotics were added and the cells were incubated at 25°C for 24 hours with shaking at 200 RPM. The scFv-His6 fusion was then purified essentially as described, except that 5 mM MgSO was used instead of MgSO, the supernatant was mixed with 0.3 M NaCl (instead of 0.5 M), 5 mL of Ni-NTA agarose (Qiagen 30210) was used instead of 0.5 mL, and the protein was exchanged into PBS buffer (140 mM NaCl, 2.7 mM KCl, 10 mM NaHPO, 1.8 mM KHPO, pH 7.4) after elution using a Vivaspin-20 (10,000 MWCO, Cytiva 28932360) by centrifugation. Ni-NTA agarose resin was also washed with 50 mM HEPES, 300 mM NaCl, 50 mM imidazole, pH 7.5, followed by the same buffer but with 70 mM imidazole, and then eluted with the same buffer but with 250 mM imidazole. The concentrated scFv was diluted with 10% glycerol, snap-frozen, and stored at -80°C. Protein purity was assessed by SDS-PAGE (Figures 28A and 28B), and concentration was determined by BCA assay.

[0327] Statistics and Reproducibility No statistical methods were used to predetermine sample size. When processing raw sequence data, low-confidence sequence reads (reads containing one or more bases with a Phred quality score below 10 or paired-end reads that did not match in overlapping regions) were discarded before further analysis or processing. These steps were performed to ensure that only highly reliable reads were used in the analysis. Otherwise, data were not excluded from any subsequent analysis or model training. A test set for evaluating model performance was constructed by randomly selecting 20% ​​of the assembled sequences and assigning them to the test set. Random partitions were generated using the Mersenne Twister random number generation algorithm implemented in Python's numpy library, version 1.19.5, with a seed value of 0. When cross-validation was used to evaluate model performance, cross-validation splits were generated by randomly dividing the dataset into five equal-sized splits using the KFold function in Python's scikit-learn library, version 0.24.2.

[0328] A final assessment of model performance was performed "blind" by generating predictions for sequences not present in these data and experimentally evaluating these predictions as described above.

[0329] Data availability Construction of the cAbRep database is described in Guo et al. (Front. Immunol. 10, 2365, 2019). The source data are provided with this paper.

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

[0331] result Overview of the RESP Pipeline The RESP pipeline contains four key components. First, we developed a simple new encoding scheme in which an autoencoder is designed to learn a representation that incorporates features that distinguish 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 those 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.

[0332] Table 9 provides a comparison of both different encoding types (one-hot, autoencoder, UniRep, ProtVec, etc.) and different models (random forest model, Bayesian network and traditional fully connected network) based on classification accuracy.

[0333] [Table 15]

[0334] Second, a yeast surface display library of mutants was constructed for a starting antibody sequence (e.g., atezolizumab, where residues were randomly mutated in the antibody heavy chain). The mutant library was incubated with a labeled target antigen (e.g., PD-L1) by incubation in the presence of excess unlabeled target antigen and screened for mutants with slower off-rates (slower off-rates, i.e., more likely to be tight binders). Flow cytometry experiments collected mutants with slower, intermediate, and faster off-rates than the WT antibody, thereby binning the population into a series of ranked groups (Figure 15 and Figure 16A). 1-2 and 16B 1-4The sequence is then determined by sequencing and the identity of the group is identified using an encoding scheme.

[0335] Third, a variational Bayesian neural network was developed to predict the off-rate of sequences by modeling experimental data via ordinal regression. This affinity model takes encoded sequences from a yeast surface display library as input and attempts to predict the likelihood that a given sequence has a low off-rate (i.e., is a tight binder). The variational Bayesian architecture minimizes the risk of overfitting and provides powerful regularization to estimate the model's uncertainty in each prediction. By using ordinal regression, each sequence was mapped to a one-dimensional latent score that indicates the model's level of confidence in the sorting group to which the sequence belongs. This method naturally takes into account the ranked order of the groups, a useful technique as previously shown for protein engineering problems. 41 The affinity model should hold for new antigens or different starting antibody sequences, but the architecture may remain unchanged.

[0336] Additionally, we modified the simulated annealing algorithm to develop an in silico directed evolution algorithm that utilizes the three previous pieces of the pipeline to efficiently explore the sequence space surrounding the training set and sequences likely to exhibit significantly lower off-rates than those of the identified parent sequences. The resulting analysis pipeline is illustrated in Figure 17. This component provides an efficient approach to exploring sequence space and can be implemented using a Bayesian neural network trained on any antibody-antigen pairing of interest.

[0337] In this study, we applied the RESP pipeline to improve the binding affinity of atezolizumab, an antibody targeting PD-L1. We selected 21 candidate sequences predicted to have slower off-rates than the WT sequence. One of them showed an approximately 10-fold slower off-rate and 17-fold tighter K than the WT atezolizumab scFv. D It has.

[0338] Furthermore, the core computational components of RESP were evaluated against data acquired by Mason et al. (Nat. Biomed. Eng. 5 (6), 600-612, 2021), demonstrating that the model described herein exhibits 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.

[0339] Encoding antibody sequences using a novel autoencoder model First, we developed a novel autoencoder model trained to represent antibody sequences and distinguish true antibody sequences from other closely related sequences. The requirement to distinguish antibodies from closely related sequences requires the encoder to embed information about typical selectivity at specific positions into its trained representation. To this end, we constructed a convolutional autoencoder consisting of three main modules or components, illustrated in Figure 18, and compared its performance with one-hot encoding, UniRep, ESM-1b, ProtVec, AntiBertY, and AbLang. The components of the pipeline are general and do not need to be retained for new antigens or wild-type sequences.

[0340] The training set for the autoencoder was drawn from the cAb-Rep database, which contains sequenced B cell receptor repertoires for 121 human donors representing true antibody sequences. Specifically, the deep 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. This dataset was then augmented with an equal number of decoys, generated by creating copies of each true antibody sequence and randomly altering them at seven positions. This number was chosen to achieve a balance: increasing the number of mutations makes it less likely that any mutant will match an actual human B cell receptor sequence present in nature. However, too many mutations makes it too easy for the classifier to distinguish the decoys from the rest, and the encoder unit is no longer forced to learn an informative representation of the input. Various numbers of mutations were experimented, and seven mutations provided a good balance.

[0341] To ensure that the decoys were indeed distinct from typical human antibody sequences, we scored a random sample of approximately 50,000 original sequences and decoyed them using three different models: BioPhi, AbLSTM, and ANARCI (see Figures 19A-19C). Using any of the three models / tools, the decoys exhibit significantly different humanness scores from the original sequences (two-tailed Mann-Whitney U test). For all three scoring methods using Python 3.9, Scipy 1.5.4, and Pingouin 0.5.2, both Scipy and Pingouin show p-values ​​that are too small to reliably calculate given floating-point error in double-precision arithmetic, returning near 0.0. Thus, the decoys are clearly scored by the models as less than human, proving their suitability to serve as decoys.

[0342] The autoencoder was trained on a decoy dataset containing a test set of 200,000 sequences, eight epochs removed from the point at which the full one-hot encoded cAb-Rep+ training loss converged. Reconstruction accuracy on the test set was over 99.99%, while prediction accuracy on the B cell receptor task for mutants was 97.4%. These results suggest that the autoencoder can compress the input sequence while retaining all the information needed to reconstruct or classify it. The performance of the autoencoder is evaluated as input to a model for predicting binding affinity in the atezolizumab modeling section below.

[0343] Generation and screening of an 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 the model can be trained to recognize the type of sequence it binds to. Here, the test case is the antigen PD-L1, and the starting point is the single-chain variable fragment (scFv) heavy chain of the atezolizumab antibody against PD-L1, hereafter referred to as wild-type (WT).

[0344] The WT atezolizumab scFv was first functionally tested in a yeast display format by testing binding to human PD-L1 (Figures 24A-24B), demonstrating robust binding to 6.1 nM of antigen present. Because the structure of atezolizumab Fab bound to PD-L1 indicates that the heavy chain is primarily involved in the binding interaction (PDB codes 5XXY, 5X8L), the variable heavy chain (VH) domain of atezolizumab was selected for mutagenesis. Furthermore, due to the shorter region required for MiSeq sequencing (approximately 316 bp in this case), mutating only the heavy chain facilitates deep sequencing. Using error-prone PCR of the heavy chain region, random mutations were generated, followed by transformation into EBY100 yeast to generate a library of up to 78 million atezolizumab scFv variants. The binding interaction and off-rate between WT scFv and PD-L1 were very strong / slow (K D = 1.75 nM, k off =1.56x10 -4 s -1 ), when the off-rate is already very low, D The large volumes required to maintain a large molar excess of antigen over scFv on the yeast surface, along with the very long incubation times required to screen for improvements in K D Yeast display is better suited under these conditions to reduce the off-rate rather than directly improving the activity. Using the WT scFv as a control during the selection, variants with faster, WT-level, and slower off-rates were isolated from the FACS sort (Figure 16A). 1-2 ), the plasmids were collected and subjected to multiple rounds of PCR, and each group was barcoded for MiSeq Nano PE250. Notably, clones selected for slower off-rates appeared to be divided into two populations, one more potent than WT and one less potent than WT (Figure 16B 1-4 , "HB" panel).

[0345] The atezolizumab dataset consists of sequences derived from mutants in three different bins with low off-rates: RH01, RH02, and RH03. RH01 contains weak binders with faster off-rates, RH02 contains moderate binders with off-rates similar to wild-type, and RH03 contains stronger binders with slower off-rates. After filtering for quality, 92,553 unique sequences were identified (550,215 sequences in total, since sequences can appear more than once in a bin or more than once in different bins), of 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).

[0346] [Table 16]

[0347] Each sequence was assigned to the category in which it was most frequently present and assigned a weight given by its frequency in the assigned category + 1 divided by the total frequency + 3. If a given sequence had equal frequency in two categories, it could not be unambiguously assigned and was discarded. As a result of these filters, 75,631 unique sequences remained. Each sequence's weight corresponds to the posterior probability that it belongs to a given category using multinomial likelihood and a uniform Dirichlet prior. This weighting is important because it allows the model to distinguish between sequences that can be reliably assigned to a single category and those that cannot. For example, a sequence that appears five times in both RH02 and RH03 is clearly less likely to be a strong binder than a sequence that appears five times only in RH03.

[0348] Atezolizumab modeling The next component in the pipeline is an affinity model that takes as input a representation of a candidate sequence and predicts which sorting category (RH01, RH02, or RH03) the sequence should belong to. This component requires experimental data specific to the antigen of interest and needs to be kept with fresh experimental data when new antigens of interest are selected.

[0349] Unlike many classifiers for protein engineering, affinity models are trained to perform ordinal regression, with the final layer of the network outputting a latent score. Traditional classification models treat categories as nominal—they have no particular order—and thus the model does not learn to rank the categories in any way. This is clearly inappropriate for typical classification data. Ordinal regression solves this problem by imposing an ordering on the categories and using latent scores to determine which category a sequence should fall into. This approach provides a simple means for ranking sequences and selecting them for experimental evaluation. Previous experiments, considering protein data in fewer than three binding categories, have shown that ordinal regression provides improved performance for accurately predicting future occurrences of sequences, more stringent classification, and is particularly superior to the sequence ranking method proposed by Liu et al. (Bioinformatics 36 (7), 2126-2133, 2020). These experiments are difficult to replicate here with only three binding categories, but given the previous results, we used ordinal regression here as well.

[0350] The affinity model outputs a score that quantifies the model's confidence that a sequence is likely to be a strong binder—a higher score indicates the model is more confident that the sequence should belong to a more stringent classification. The score is then added to M−1 learned thresholds for the M categories, after which a sigmoid is applied to generate an output vector of M−1 probability values, where each element i of this output vector is the probability assigned to the model that the input sequence belongs to a binding category more stringent than i.

[0351] Briefly, the affinity model described herein uses an architecture similar to Bayes' backpropagation algorithm, differing in that this model is adapted to perform ordinal regression. Traditional fully connected neural networks learn specific values ​​for each weight and a bias term in each hidden layer during training. In contrast, Bayesian neural networks treat each weight as a (usually Gaussian) distribution and learn parameters for each weight distribution during training. Instead of fitting using maximum likelihood, as is typical for neural networks, Bayesian architectures are fitted by approximating the posterior probability distribution using variational methods. Uncertainty in predictions can be estimated by repeatedly sampling the weight distribution. Alternatively, if preferred, the mean of the Gaussian distribution for each weight is used to generate a point estimate. In addition to generating uncertainty estimates, Bayesian neural networks also provide strong regularization for model parameters by penalizing deviations from the prior distribution.

[0352] The uncertainty information provided by the model can be clearly used in one of two ways: predictions with high uncertainty will match sequences in relatively unexplored regions of the input space; indeed, only predictions with relatively low uncertainty should be selected for cost- and time-considered experimental evaluation, i.e., the strategy pursued here.

[0353] The results produced by a vanilla fully connected neural network trained to perform simple classification were compared to a Bayesian network trained to perform ordinal regression across various encoding types (Table 9). Note that both compared neural networks have the same number of hidden layers and weights. While the non-Bayesian network provides a modest improvement in performance, the Bayesian network provides additional information about the uncertainty associated with a given prediction, which is important for analyzing classification data, as it allows us to judge the relative reliability of the model's predictions for candidate sequences.

[0354] Affinity models can use one-hot encoded sequences as input or as a separate representation. To determine whether the representations generated by autoencoders are useful as encodings for affinity models, we encoded the atezolizumab sequence using 1) a fully trained autoencoder, 2) one-hot encoding, 3) ProtVec encoding, 4) UniRep embedding, 5) ESM-1b encoding, 6) AbLang embedding, and 7) AntiBERTy embedding. The encoded dataset was split into training (80%) and testing (20%) parts, and 5-fold cross-validation was performed on the training part for each encoding type using a Bayesian neural network trained to perform ordinal regression (Table 9).

[0355] Although different encodings may prove preferable for specific tasks, the autoencoder is the only representation that is comparable to or better than one-hot encoding of antibodies for any model type (Table 9). It provides comparable performance despite a significant reduction in model size and computational expense compared to 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 parameters.

[0356] Thus, autoencoders offer at least two distinct benefits. First, they have only a few parameters, significantly reducing model size and complexity compared to simple one-hot encoding and are substantially cheaper than language models. Second, they provide a real-valued representation of each sequence that can be used to cluster sequences or determine similarities between them at later stages in the pipeline. Given these advantages, autoencoders were preferred over one-hot encoding here, even though they achieve similar performance.

[0357] In silico directed evolution to select the most promising candidates After encoding the atezolizumab dataset using an autoencoder, a Bayesian neural network (affinity model) on the entire dataset was trained for 30 epochs for off-rates. At this point, the pipeline uses the trained affinity model to search the sequence space surrounding the training set and generate new candidate antibody sequences not present in the training set.

[0358] At this stage, it is desirable to reduce the size of the search space. When 118 amino acid sequences are mutated in silico, the search space covered is extremely large. Therefore, we focused on positions that contribute significantly to binding affinity—top-scoring sequences, i.e., a smaller subset of the most frequently mutated positions in the most promising sequences considered by the model. For this dataset, we found that using the top 500, 1000, 2000, or 4,000 highest-scoring sequences to select the top 10 most frequently mutated positions recovered the same set of positions used in subsequent search steps. This step in the pipeline is flexible, and other criteria (e.g., only frequently mutated positions present in CDRs) could be used to select positions.

[0359] In the following search step, for 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 by this procedure have the same length as the wild-type sequence. The probability of selecting any particular amino acid at a given position is given by the frequency of that amino acid in the training set + 1 divided by the total number of sequences in the training set + 20. New candidate sequences are scored by encoding them using an encoder and scoring them using the trained affinity model. The score is evaluated using a classical simulated annealing criterion, whereby candidates are accepted with a probability determined by their score, the score of the last accepted candidate, and the temperature. The temperature starts at a high value to allow for exploration of a large sequence space and is decreased with each iteration, so that the probability of accepting a proposed sequence without improving the score decreases.

[0360] Because the average sequence in the training set contains just seven mutations and the average high-scoring sequence contains even fewer mutations (three mutations on average), at many positions, many common amino acids are those present in the wild type. As a result, the algorithm tends to oversample sequences similar to the wild type. However, given the stochastic nature of the algorithm, it can explore combinations not present in the training set, thereby enabling the discovery of new high-scoring sequences. 1-2 20B and 20C illustrate how the best achieved scores to date evolve over the course of this optimization.

[0361] Ten simulated annealing strands were run; all converged within 1000 iterations. Accepted candidates from each strand with scores above the 90th percentile were collected, and identical candidates were removed. Uncertainty was assessed on the scores by sampling each sequence from a Bayesian neural network 1000 times. The top 50% of sequences showing the largest standard deviation in terms of assigned scores were removed. Here, the 50% threshold was arbitrary, and this process yielded 127 sequences.

[0362] These sequences were then clustered using median hierarchical clustering (the resulting dendrogram is provided in Figure 33). The results clearly suggest that the selected sequences can be divided into two major subgroups. The marginal distribution (Figure 20A) 1-2 , see Figures 20B and 20C) show that the two major groups differ primarily at position 77, but also contain subgroups that exhibit several other interesting differences. Many selected sequences exhibit, for example, the R98C mutation, while a subset of cluster 1 is unchanged from wild-type at that position.

[0363] Cutting the tree at a lower height resulted in 11 subclusters, and the two highest-scoring sequences in each cluster (or one if only one sequence was present) were selected to yield 21 final candidates. This threshold was chosen to obtain a manageable number of sequences for experimental evaluation.

[0364] Examination of 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 Figure 21).

[0365] [Table 17-1] [Table 17-2]

[0366] The amino acid sequences of the full-length IgG scFv sequences for each of the 21 mutants are provided below: Mutant 1 DIQMTQSPSSLSASVGDRVTITCRASQDVSTAVAWYQQKPGKAPKLLIYSASFLYSGVPSRFSGSGSGTDFTLTISSLQPEDFATYYCQQYLYHPATFGQGTKVEIKRGGGGSGGGGSGGGGSGGGGSEVQLVESGGGLVQPGGSLRLSCAASGFTFSDSWIHWVRQAPGQGLEWVAWISPYGGSTYYADSVKGRFTISADTSKNTTYLQMNSLRAEDTAVYYCVCRHWPGGFDYWGQGTLVTVSS (SEQ ID NO: 22) Mutant 2 DIQMTQSPSSLSASVGDRVTITCRASQDVSTAVAWYQQKPGKAPKLLIYSASFLYSGVPSRFSGSGSGTDFTLTISSLQPEDFATYYCQQYLYHPATFGQGTKVEIKRGGGGSGGGGSGGGGSGGGGSEVQLVESGGGLVQPGGSLRLSCAASGFTFSDSWIHWVRQAPGQGLEWVAWISPYGGSTYYADSVKGRFTASADTSKNTTYLQMNSLRAEDTAVYYCVCRHWPGGFDYWGQGTLVTVSS (SEQ ID NO: 23) Mutant 3 DIQMTQSPSSLSASVGDRVTITCRASQDVSTAVAWYQQKPGKAPKLLIYSASFLYSGVPSRFSGSGSGTDFTLTISSLQPEDFATYYCQQYLYHPATFGQGTKVEIKRGGGGSGGGGSGGGGSGGGGSEVQLVESGGGLVQPGGSLRLSCAASGFTFSDSWIHWVRQAPGRGLEWVAWISPYGGSTYYADSVKGRFTTSADTSKNTTYLQMNSLRAEDTAVYYCVRRHWPGGFDYWGQGTLVTVSS (SEQ ID NO: 24) Mutant 4 DIQMTQSPSSLSASVGDRVTITCRASQDVSTAVAWYQQKPGKAPKLLIYSASFLYSGVPSRFSGSGSGTDFTLTISSLQPEDFATYYCQQYLYHPATFGQGTKVEIKRGGGGSGGGGSGGGGSGGGGSEVQLVESGGGLVQPGGSLRLSCAASGFTFSDSWIHWVRQAPGKGLEWVAWISPYGGSTYYADSVKGRFTASADTSKNTTYLQMNSLRAEDTAVYYCVRRHWPGGFDYWGQGTLVTVSS (SEQ ID NO: 25) Mutant 5 DIQMTQSPSSLSASVGDRVTITCRASQDVSTAVAWYQQKPGKAPKLLIYSASFLYSGVPSRFSGSGSGTDFTLTISSLQPEDFATYYCQQYLYHPATFGQGTKVEIKRGGGGSGGGGSGGGGSGGGGSEVQLVESGGGLVQPGGSLRLSCAASGFTFSDSWIHWVRQAPGEGLEWVAWISPYGGSAYYADSVKGRFTTSADTSKNTTYLQMNSLRAEDTAVYYCVCRHWPGGFDYWGQGTLVTVSS (SEQ ID NO: 26) Mutant 6 DIQMTQSPSSLSASVGDRVTITCRASQDVSTAVAWYQQKPGKAPKLLIYSASFLYSGVPSRFSGSGSGTDFTLTISSLQPEDFATYYCQQYLYHPATFGQGTKVEIKRGGGGSGGGGSGGGGSGGGGSEVQLVESGGGLVQPGGSLRLSCAASGFTFSDSWIHWVRQAPGRGLEWVAWISPYGGSAYYADSVKGRFTTSADTSKNTTYLQMNSLRAEDTAVYYCVCRHWPGGFDYWGQGTLVTVSS (SEQ ID NO: 27) Mutant 7 DIQMTQSPSSLSASVGDRVTITCRASQDVSTAVAWYQQKPGKAPKLLIYSASFLYSGVPSRFSGSGSGTDFTLTISSLQPEDFATYYCQQYLYHPATFGQGTKVEIKRGGGGSGGGGSGGGGSGGGGSEVQLVESGGGLVQPGGSLRLSCAASGFTFSDSWIHWVRQAPGKGLEWVAWISPYGGSTYYADSVKGRFTASADTSKNTTYLQMNSLRAEDTAVYYCVCRHWPGGFDYWGQGTLVTVSS (SEQ ID NO: 28) Mutant 8 DIQMTQSPSSLSASVGDRVTITCRASQDVSTAVAWYQQKPGKAPKLLIYSASFLYSGVPSRFSGSGSGTDFTLTISSLQPEDFATYYCQQYLYHPATFGQGTKVEIKRGGGGSGGGGSGGGGSGGGGSEVQLVESGGGLVQPGGSLRLSCAASGFTFSDSWIHWVRQAPGKGLEWVAWISPYGGSAYYADSVKGRFTTSADTSKNTTYLQMNSLRAEDTAVYYCVCRHWPGGFDYWGQGTLVTVSS (SEQ ID NO: 29) Mutant 9 DIQMTQSPSSLSASVGDRVTITCRASQDVSTAVAWYQQKPGKAPKLLIYSASFLYSGVPSRFSGSGSGTDFTLTISSLQPEDFATYYCQQYLYHPATFGQGTKVEIKRGGGGSGGGGSGGGGSGGGGSEVQLVESGGGLVQPGGSLRLSCAASGFTFSDSWIHWVRQAPGKGLEWVAWISPYGGSAYYADSVKGRFTISADTSKSTVYLQMNSLRAEDTAVYYCVCRHWPGGFDYWGQGTLVTVSS (SEQ ID NO: 30) Mutant 10 DIQMTQSPSSLSASVGDRVTITCRASQDVSTAVAWYQQKPGKAPKLLIYSASFLYSGVPSRFSGSGSGTDFTLTISSLQPEDFATYYCQQYLYHPATFGQGTKVEIKRGGGGSGGGGSGGGGSGGGGSEVQLVESGGGLVQPGGSLRLSCAASGFTFSDSWIHWVRQAPGKGLEWVAWISPYGGSAYYADSVKGRFTISADTSKDTVYLQMNSLRAEDTAVYYCVCRHWPGGFDYWGQGTLVTVSS (SEQ ID NO: 31) Mutant 11 DIQMTQSPSSLSASVGDRVTITCRASQDVSTAVAWYQQKPGKAPKLLIYSASFLYSGVPSRFSGSGSGTDFTLTISSLQPEDFATYYCQQYLYHPATFGQGTKVEIKRGGGGSGGGGSGGGGSGGGGSEVQLVESGGGLVQPGGSLRLSCAASGFTFSDSWIHWVRQTPGRGLEWVAWISPYGGSAYYADSVKGRFTISADTSKDTTYLQMNSLRAEDTAVYYCACRHWPGGFDYWGQGTLVTVSS (SEQ ID NO: 32) Mutant 12 DIQMTQSPSSLSASVGDRVTITCRASQDVSTAVAWYQQKPGKAPKLLIYSASFLYSGVPSRFSGSGSGTDFTLTISSLQPEDFATYYCQQYLYHPATFGQGTKVEIKRGGGGSGGGGSGGGGSGGGGSEVQLVESGGGLVQPGGSLRLSCAASGFTFSDSWIHWVRQTPGEGLEWVAWISPYGGSAYYADSVKGRFTISADTSKDTTYLQMNSLRAEDTAVYYCTCRHWPGGFDYWGQGTLVTVSS (SEQ ID NO: 33) Mutant 13 DIQMTQSPSSLSASVGDRVTITCRASQDVSTAVAWYQQKPGKAPKLLIYSASFLYSGVPSRFSGSGSGTDFTLTISSLQPEDFATYYCQQYLYHPATFGQGTKVEIKRGGGGSGGGGSGGGGSGGGGSEVQLVESGGGLVQPGGSLRLSCAASGFTFSDSWIHWVRQTPGRGLEWVAWISPYGGSAYYADSVKGRFTASADTSKDTTYLQMNSLRAEDTAVYYCVCRHWPGGFDYWGQGTLVTVSS (SEQ ID NO: 34) Mutant 14 DIQMTQSPSSLSASVGDRVTITCRASQDVSTAVAWYQQKPGKAPKLLIYSASFLYSGVPSRFSGSGSGTDFTLTISSLQPEDFATYYCQQYLYHPATFGQGTKVEIKRGGGGSGGGGSGGGGSGGGGSEVQLVESGGGLVQPGGSLRLSCAASGFTFSDSWIHWVRQVPGRGLEWVAWISPYGGSAYYADSVKGRFTASADTSKDTTYLQMNSLRAEDTAVYYCVCRHWPGGFDYWGQGTLVTVSS (SEQ ID NO: 35) Mutant 15 DIQMTQSPSSLSASVGDRVTITCRASQDVSTAVAWYQQKPGKAPKLLIYSASFLYSGVPSRFSGSGSGTDFTLTISSLQPEDFATYYCQQYLYHPATFGQGTKVEIKRGGGGSGGGGSGGGGSGGGGSEVQLVESGGGLVQPGGSLRLSCAASGFTFSDSWIHWVRQVPGKGLEWVAWISPYGGSAYYADSVKGRFTASADTSKDTTYLQMNTLRAEDTAVYYCVCRHWPGGFDYWGQGTLVTVSS (SEQ ID NO: 36) Mutant 16 DIQMTQSPSSLSASVGDRVTITCRASQDVSTAVAWYQQKPGKAPKLLIYSASFLYSGVPSRFSGSGSGTDFTLTISSLQPEDFATYYCQQYLYHPATFGQGTKVEIKRGGGGSGGGGSGGGGSGGGGSEVQLVESGGGLVQPGGSLRLSCAASGFTFSDSWIHWVRQVPGKGLEWVAWISPYGGSAYYADSVKGRFTASADTSKDTIYLQMNSLRAEDTAVYYCVCRHWPGGFDYWGQGTLVTVSS (SEQ ID NO: 37) Mutant 17 DIQMTQSPSSLSASVGDRVTITCRASQDVSTAVAWYQQKPGKAPKLLIYSASFLYSGVPSRFSGSGSGTDFTLTISSLQPEDFATYYCQQYLYHPATFGQGTKVEIKRGGGGSGGGGSGGGGSGGGGSEVQLVESGGGLVQPGGSLRLSCAASGFTFSDSWIHWVRQTPGKGLEWVAWISPYGGSAYYADSVKGRFTTSADTSKDTTYLQMNSLRAEDTAVYYCVCRHWPGGFDYWGQGTLVTVSS (SEQ ID NO: 38) Mutant 18 DIQMTQSPSSLSASVGDRVTITCRASQDVSTAVAWYQQKPGKAPKLLIYSASFLYSGVPSRFSGSGSGTDFTLTISSLQPEDFATYYCQQYLYHPATFGQGTKVEIKRGGGGSGGGGSGGGGSGGGGSEVQLVESGGGLVQPGGSLRLSCAASGFTFSDSWIHWVRQVPGKGLEWVAWISPYGGSAYYADSVKGRFTTSADTSKDTTYLQMNSLRAEDTAVYYCVCRHWPGGFDYWGQGTLVTVSS (SEQ ID NO: 39) Mutant 19 DIQMTQSPSSLSASVGDRVTITCRASQDVSTAVAWYQQKPGKAPKLLIYSASFLYSGVPSRFSGSGSGTDFTLTISSLQPEDFATYYCQQYLYHPATFGQGTKVEIKRGGGGSGGGGSGGGGSGGGGSEVQLVESGGGLVQPGGSLRLSCAASGFTFSDSWIHWVRQAPGKGLEWVAWISPYGGSTYYADSVKGRFTVSADTSKDTTYLQMNSLRAEDTAVYYCVRRHWPGGFDYWGQGTLVTVSS (SEQ ID NO: 40) Mutant 20 DIQMTQSPSSLSASVGDRVTITCRASQDVSTAVAWYQQKPGKAPKLLIYSASFLYSGVPSRFSGSGSGTDFTLTISSLQPEDFATYYCQQYLYHPATFGQGTKVEIKRGGGGSGGGGSGGGGSGGGGSEVQLVESGGGLVQPGGSLRLSCAASGFTFSDSWIHWVRQAPGKGLEWVAWISPYGGSTYYADSVKGRFTASADTSKDTTYLQMNSLRAEDTAVYYCVRRHWPGGFDYWGQGTLVTVSS (SEQ ID NO: 41) Mutant 21 DIQMTQSPSSLSASVGDRVTITCRASQDVSTAVAWYQQKPGKAPKLLIYSASFLYSGVPSRFSGSGSGTDFTLTISSLQPEDFATYYCQQYLYHPATFGQGTKVEIKRGGGGSGGGGSGGGGSGGGGSEVQLVESGGGLVQPGGSLRLSCAASGFTFSDSWIHWVRQAPGKGLEWVAWISPYGGSAYYADSVKGRFTTSADTSKDTAYLQMNFLRAEDTAVYYCACRHWPGGFDYWGQGTLVTVSS (SEQ ID NO: 42)

[0367] Inspection of the existing structures of atezolizumab Fab and PD-L1 (PDB codes 5XXY and 5X8L) shows that none of the mutated residues contact PD-L1 in the WT sequence (Figure 21). This suggests that these mutations improve binding affinity by altering heavy chain conformation rather than directly improving binding contacts. Another possibility is that these mutations stabilize a conformation rather than alter it. For example, it has been demonstrated that stabilizing certain antibody conformations can be beneficial for affinity in some cases.

[0368] To explore these possibilities, predicted structures for the top-scoring mutants from 21 sequences were generated using IgFold without antigen. The mutations present in the high-scoring sequences shift the conformation of both antigen-contacting regions of the protein, particularly CDR H3. Figures 30A-30B and 31 show the conformational changes of R98C (a mutation present in many of the top-scoring sequences) and its surrounding residues, as well as the formation of an apparent hydrophobic pocket in the top mutants compared to the wild type. Note that the mutations described herein do not directly contact the antigen. These predicted structures suggest that the mutations selected by this pipeline function by altering the conformation of contact regions rather than by directly forming new contacts with the antigen.

[0369] The trained affinity model was also used to computationally evaluate the importance of individual mutations. A) How much each mutation contributes to binding when introduced individually into the WT without any other mutations, and B) How much binding affinity changes when a single mutation is removed from a combination of beneficial mutations. The full results are presented in Tables 13 and 14.

[0370] [Table 18]

[0371] [Table 19-1] [Table 19-2] [Table 19-3] [Table 19-4] [Table 19-5]

[0372] The results suggest that some mutations may be more important than others, and that 12 of the 19 introduced mutations are predicted to be beneficial alone, i.e., increase the likelihood that the sequence in question is a tight binder. The remaining mutations are predicted to be only slightly deleterious or no change alone, but to have beneficial effects in certain circumstances. A

[79] I or A

[79] T, R

[98] C, I

[70] A, A

[97] V, and T

[58] A are predicted by the model to be mutations that have a very beneficial effect alone.

[0373] Interestingly, the model can predict synergistic effects between individual mutations. For example, K43Q is predicted to have a nearly negligible beneficial effect when introduced alone into the wild type, but a much larger negative impact on score when it is removed from mutants K

[43] Q, A

[79] T, A

[97] V, and R

[98] C. A

[97] T is predicted to be only slightly deleterious alone but beneficial with respect to A

[40] T, K

[43] E, T

[58] A, N

[77] D, A

[79] T, A

[97] T, and R

[98] C. A

[79] T, for example, is predicted to be beneficial whatever the context in which it is introduced, but the magnitude of its effect on score varies by as much as two-fold depending on the other mutations present.

[0374] Validation of predicted tight binding antibodies The 21 model-generated mutants (all heavy chains) were purchased as geneblocks, fused to the WT light chain sequence by PCR, and then transformed into EBY100 yeast. A small library of 21 mutants was screened in a similar manner to the naive library, except that a longer competition time (39 hours) was used for the final screen in an attempt to isolate the mutant with the slowest off-rate from the remaining 21 mutants (Figures 25A-25D). This process did not result in clearly separated mutants and sequencing of random clones after the most stringent sorting indicated that 12 distinct mutants (occurring at similar frequencies) were found among the 17 random sequences. This suggests that a significant proportion of the 21-member pool of mutants had significantly slower off-rates than the WT, which exhibited significantly lower fluorescence when compared to the pool of mutants in the final, most stringent sorting (Figures 25C-25D, lower panels). Essentially all of the 21 mutants selected by the model also appeared to have slower off-rates than WT when examined with the original / unsorted 21-member library after 8 hours of off-rate competition at RT (Figures 25A-25B, top panels).

[0375] From the various mutants selected from the 21-member library after the most stringent sorting, the I70A / A79T / A97V mutant (designated mutant 4, see Table 11 for sequence) was selected for further characterization.

[0376] [Table 20-1] [Table 20-2]

[0377] The amino acid sequence of the full-length IgG heavy chain / constant region of mutant 4 is presented below: MGWSCIILFLVATATGVHSEVQLVESGGGLVQPGGSLRLSCAASGFTFSDSWIHWVRQAPGKGLEWVAWISPYGGSTYYADSVKGRFTASADTSKNTTYLQMNSLRAEDTAVYYCVRR HWPGGFDYWGQGTLVTVSSASTKGPSVFPLAPSSKSTSGGTAALGCLVKDYFPEPVTVSWNSGALTSGVHTFPAVLQSSGLYSLSSVVTVPSSSLGTQTYICNVNHKPSNTKVDKKVEP KSCDKTHTCPPCPAPELLGGPSVFLFPPKPKDTLMISRTPEVTCVVVDVSHEDPEVKFNWYVDGVEVHNAKTKPREEQYASTYRVVSVLTVLHQDWLNGKEYKCKVSNKALPAPIEKTISKAKGQPREPQVYTLPPSREEMTKNQVSLTCLVKGFYPSDIAVEWESNGQPENNYKTTPPVLDSDGSFFLYSKLTVDKSRWQQGNVFSCSVMHEALHNHYTQKSLSLSPGK (SEQ ID NO: 67)

[0378] The amino acid sequence of the full-length IgG light chain / constant region of mutant 4 is presented below: MGWSCIILFLVATATGVHSDIQMTQSPSSLSASVGDRVTITCRASQDVSTAVAWYQQKPGKAPKLLIYSASFLYSGVPSRFSGSGSGTDFTLTISSLQPEDFATYYCQQYLYHPATFGQGTKVEIKRTVAAPSVFIFPPSDEQLKSGTASVVCLLNNFYPREAKVQWKVDNALQSGNSQESVTEQDSKDSTYSLSSTLTLSKADYEKHKVYACEVTHQGLSSPVTKSFNRGEC (SEQ ID NO: 68)

[0379] This mutant was more preferred because, unlike the other mutants sequenced after stringent sorting, it did not have a cysteine ​​introduced into its sequence, and therefore no additional disulfide bonds were formed by the mutation. To characterize the improved off-rate / binding of the isolated mutants to PD-L1, the K D and K. off The K values ​​were determined using a yeast display format. First, the apparent off-rate was determined on the yeast surface at RT (Figure 22A), revealing that mutant 4 had an approximately 10-fold slower K than WT atezolizumab. off (6.3x10 for WT) -5 s -1 , half-life 3.04 hr vs. 6.5x10 for mutant 4 -6 s -1 (Half-life 29.8 hours). Then, the K off The values ​​were compared to two other FDA-approved anti-PD-L1 mAbs (durvalumab and avelumab, converted to scFv format on the yeast surface) (Figure 22B). Mutant 4 exhibited a substantially slower off-rate compared to both WT atezolizumab and durvalumab, and a significantly lower value than avelumab. The order of off-rates of the FDA-approved mAb PD-L1 scFvs (durvalumab < atezolizumab < avelumab) is consistent with previous studies of these scFv values ​​by surface plasmon resonance (SPR) binding assays.

[0380] Next, the binding affinities (K D ) values ​​were determined by titrating PD-L1 against each scFv on the yeast surface. K for mutant 4 D was approximately 17-fold more robust (on average) than WT (5.3 pM vs. 92 pM, Figure 22C). Excess volume was required to maintain a significant molar excess of antigen to scFv on the yeast surface at lower antigen concentrations. 49Considering that performing an iteration for each data point would also be difficult, the K D It should be noted that it is difficult to determine the K. However, there are three separate K D Measurements showed a significant (more than 10-fold) improvement in binding affinity (Figure 22C).

[0381] To assess whether the improved binding affinity and off-rate observed on the yeast surface could be replicated after purifying each scFv and reversing the orientation of the binding assay (PD-L1 immobilized and bound to soluble scFvn instead of binding to scFv fused to the yeast surface and soluble PD-L1), we performed a biolayer interferometry (BLI) binding assay. Consistent with the yeast surface display results, mutant 4 was found to have a much slower off-rate than WT. Dissociation of the mutant from PD-L1 in this assay was minimal, even after 2 hours at RT, whereas WT was more significant over a dissociation time of approximately 60 minutes. However, BLI revealed a K of 2 for WT. D values ​​(219 pM for the major species and 5.93 nM for the minor species), whereas the K D The value was not measurable due to its very slow off-rate, and no rate constant value could be fitted / determined (Figure 26A). 1-2and 26B). The minor species could be caused by a minor population of partially unfolded scFv in the protein preparation. This observation led us to examine whether the scFv and PD-L1 were monomeric or formed oligomers, and therefore this was tested by mass spectrophotometry, a method that can image the oligomeric distribution of proteins in solution. Imaging of PD-L1 (residues Met1-T239, predicted MW 35-38 kDa, MW range due to glycosylation) showed a predominantly monomeric protein at the concentrations tested (7.4 nM, 18.5 nM, and 74 nM) (Figures 27A-27G), which is consistent with previous studies of human PD-L1, which was found to be monomeric. These data were consistent with the K measured on the yeast surface. D This supports the correlation between monomeric PD-L1 and the scFv immobilized on the yeast surface. Imaging of WT or mutant 4 scFv in solution revealed a mixed population of monomers, dimers, trimers, and possibly higher MW oligomers at 100 nM scFv, respectively. WT and mutant 4 scFv were purified twice (using the same extraction method) and tested twice by mass spectrometry. For the WT scFv, in both preparations, it was a mixture of monomers and dimers in solution, whereas, for unknown reasons, the mutant was either primarily monomeric (first batch) or a mixture of monomers, dimers, and trimers (batch 2) (Figures 27A-27G). Furthermore, the mutant scFv appeared to be more prone to aggregation than the WT scFv. Taken together, because it is difficult to know the soluble monomeric concentration of scFv, oligomerization and aggregation of scFv in solution can be investigated by using binding assays such as SPR (surface plasmon resonance) or BLI, in which scFv and PD-L1 in solution are immobilized. D making it unsuitable for measurement.

[0382] Evaluation against the dataset of Mason et al. (Nat. Biomed. Eng., 5 (6), 600-612, 2021) To ensure the generality of RESP, we next tested it on data from the study by Mason et al. Starting with a non-binding mutant of trastuzumab, we successfully recovered an antibody with similar affinity to the original wild-type trastuzumab. Using the same training and test sets, we encoded the sequence data using an autoencoder and trained a variational Bayesian network with the same structure used for the data to predict the bound and non-bound states.

[0383] This model achieved the same accuracy (Matthews correlation coefficient 0.68, AUC-ROC 0.91) as theirs for the test set, accurately predicting the sequences they experimentally tested to be binders (see Figures 32A-32C). The model also assigns high scores to many of these sequences. Trastuzumab, for example, is assigned a higher score than 94% of the sequences in the training set, indicating high confidence that it is a strong binder. The distribution of scores for the training set and for the experimentally evaluated sequences can be seen in Figures 32A-32C.

[0384] However, unlike the Mason et al. pipeline, the pipeline considered in this example provides additional uncertainty information not available from their model. The affinity model's uncertainty for test set predictions that turn out to be inaccurate is significantly higher than its uncertainty for test sets that are accurate (p<1e-23, two-tailed Mann-Whitney U test). As a result, uncertainty can be used to help determine whether a prediction is likely to be reliable, which is not possible with the pipeline developed by Mason et al., and there is no need to sacrifice accuracy to gain this advantage.

[0385] Consideration A pipeline for efficient antibody development is provided with the following innovations. First, a simple, learned encoding specific to antibodies is developed. Such encoding not only contains the information in the original sequence in an easily recoverable form, but also contains additional encoded information that describes important sequence features that distinguish human antibodies from the surrounding sequence space. When training a model on antigen-specific experimental data to capture trends in binding affinity, the same model proved more accurate when the input sequence was encoded using representations generated by an autoencoder than when the input sequence was encoded using popular state-of-the-art encodings such as UniRep, ESM-1b, AntiBertY, and AbLang. Surprisingly, for this particular task, the use of all possible models, UniRep, ESM-1b, AntiBertY, and AbLang, performed worse than one-hot encoding. This is consistent with the 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 to simple one-hot encoding.

[0386] The training set was then fitted using a model designed to provide simple and easy interpretation of sequence ranks along with quantification of uncertainty. We demonstrated that the distribution of sequence reads across categories can be incorporated into model fitting as data point weights in a principled and straightforward manner. Bayesian ordinal regression models provide estimates of posterior predictions, thereby providing additional information not available from traditional deep learning classifiers, where the predicted probability distribution across categories does not indicate the reliability of a given prediction.

[0387] Furthermore, the algorithm was designed to explore the sequence space spanned by the training set. By estimating the reliability of each prediction and by limiting the search to the space spanned by the training set, the time and expense wasted in evaluating poor candidate sequences was minimized.

[0388] We further validated the power of the pipeline disclosed in this example. By training a model against a single large library, we were able to select mutants that exhibited 10-fold / 17-fold improvements in off-rate / binding affinity. This is consistent with previous results on protein engineering of CBX1, demonstrating that a similar strategy resulted in binding affinity improvements comparable to those achieved by much longer directed evolution processes. Mutant 4 identified and selected here is a previously reported high-affinity anti-PD-L1 protein (K of 110 pM). D As an scFv, it may be useful in cancer therapy, such as a mutant form of PD-1 (which binds to PD-L1 via its cytoplasmic domain). Like PD-1, the scFv of mutant 4 is significantly smaller than the monoclonal antibody (30 kDa vs. 150 kDa) and therefore likely more effective at tumor penetration.

[0389] This pipeline disclosed in this example may also have several important advantages over purely experimental approaches: phage and yeast display only allow the selection of small populations rather than single clones, and therefore additional experiments (ELISA or yeast K) may be required to assess clones with the strongest binding affinities. DMeasurements, etc., are required. Because only strong binders present in the original library can be identified by these techniques, very large libraries and / or multiple libraries are often used to maximize coverage of sequence space. The best binders identified by this process still often do not have sufficiently strong affinities, so the desired affinity often must be achieved by 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, antibodies generated by this process are not guaranteed to have other desirable characteristics such as solubility or stability.

[0390] Take for example the process by which the atezolizumab antibody itself was first discovered (described in US Patent US8217149B2). Four rounds of panning on a phage display library screened against the PDL1 target were first used to recover 96 enriched clones. Two sets of positions were then selected to construct two additional libraries for affinity maturation to improve binding. These clones were then used for one-plate selection, followed by five or six rounds of solution selection. Finally, enriched clones from the final selection were subjected to high-throughput ELISA screening to identify the best candidates. It is worth noting that this procedure, while reliable, is expensive in terms of time and money. It is very common for antibody engineering to require the construction of multiple libraries and high-throughput ELISA of thousands of clones; see, inter alia, US Patents 20180086848A1 and 8313746B2 for several other examples, including antibodies currently in clinical trials.

[0391] In contrast, the approach disclosed herein allowed for the construction of a single large library followed by selection of tightly binding antibodies after FACS-based sorting for off-rate and binding. This approach is suitable for high-throughput KD No determination of the number of copies of a library or subsequent focused library screening is required, thereby eliminating the time required to construct multiple libraries, as well as high-throughput ELISA screening / K screening, which are considerable. D The time and expense required for measurement are entirely eliminated. Instead, this method is easy to implement on a single computer with a GPU and requires only a few computational steps that are quickly executed. The Bayesian neural network needs to be maintained at any time, and new antigens are selected on data specific to that antigen and obtained as described in this example. In contrast, autoencoders can be reused and do not need to be maintained. Furthermore, unlike traditional phage and yeast display, which can only screen sequences present in the library, the method disclosed in this example identified strong binders that were not present in the original library.

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

[0393] While other strategies for computer-aided antibody design have been described in recent literature, the approach disclosed herein compares favorably and offers several advantages. For example, Mason et al. achieved a 3-fold improvement in affinity over wild-type trastuzumab; indeed, many of the mutants selected by their algorithm as promising exhibited weaker affinity for the target, but a 17-fold improvement in affinity was achieved using the approach disclosed herein. Mason et al. performed multiple rounds of mutagenesis and library generation, including a step in which positions for further mutagenesis were selected using rational design based on data from previous selection steps. They required additional efforts to constrain the search space. In contrast, this could be achieved in silico using the modified simulated annealing algorithm disclosed herein. Unlike the CNN they propose, the Bayesian neural network disclosed herein provides an uncertainty estimate that correlates with the accuracy of the prediction, allowing for the elimination of predictions that are likely to be unreliable.

[0394] Rational design approaches (as opposed to machine learning-based approaches) have been described. However, their rational design components result in a low percentage of accurate predictions and are sometimes only as accurate as the predictions of Rosetta modeling software, which cannot determine which predictions are most likely to be reliable. For example, they sought to improve the binding of an antibody called G6 to its target, VEGF, using a computational search procedure requiring approximately 250 cpu-days, several orders of magnitude more expensive than the computational procedure employed and disclosed herein. Of the mutants selected by this AbLift procedure for experimental evaluation, 60% were inferior to the wild type, and many of the remaining were only comparable. Only one of the designs recommended by their approach achieved a K Dindeed improved the binding affinity of WT, resulting in a 5-fold improvement, a smaller improvement than that disclosed herein. It is noteworthy that many or all of the 21 mutants selected by the model disclosed herein exhibited slower off-rates (and thus presumably higher binding affinities) than WT.

[0395] Gaussian processes were 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 Absolut! Software to determine whether a sequence is a strong binder; since Absolut! Software is itself based on docking-generated structures and affinities, it is unclear how closely it tracks experimental data. Importantly, their approach suffers from some well-known limitations of exact Gaussian processes. These models require O(N) of training points. 3 ) scaling (or in some more efficient modern implementations, O(N 2 ) and is thus completely infeasible for datasets larger than about 5-10,000 sequences. In contrast, the variational Bayesian network-based approach disclosed herein does not suffer from such limitations.

[0396] Although the method disclosed herein has been validated using atezolizumab as a starting point and PDL1 as a target, there is nothing specific about this method to the selected target, and thus this method can be easily adapted to other targets and problems.The pipeline described herein can be easily modified to optimize only a single CDR or a subset of available positions, and to incorporate other models that rank candidate sequences for other desirable properties in addition to affinity.It is therefore expected that this pipeline and its modified versions will prove to be a useful tool for facilitating antibody discovery and development.

Claims

1. A method for designing and engineering proteins with desired properties of interest, including the steps of directed evolution, sequencing, and machine learning analysis.

2. 10. 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 the directed evolution sequences.

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

4. 3. The method of claim 2, wherein the steps of directed evolution, sequencing, and machine learning analysis further comprise performing ordinal regression on the encoded directed evolution sequences using an approximated Gaussian Process (GP) computer model that uses a random feature approximation kernel.

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

6. 4. The method of claim 3, wherein the in silico mutagenesis is performed using a simulated annealing strategy to search for directed evolution sequences not present in the mutation library and evaluate their properties of interest.

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

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. 1. A method for designing and engineering antibodies containing desired mutations and having improved binding affinity or off-rate for a target antigen, comprising: a) establishing an autoencoder model based on encoded mutant sequence data of known antibodies binding to a target antigen; b) developing a Bayesian ordinal regression model that scores each mutant sequence for its probability of being a strong binder to the target antigen; and c) selecting antibodies containing desired mutations using a modified simulated annealing algorithm to test binding affinity or off-rate to the target antigen; A method comprising:

11. 11. The method of claim 10, wherein the target 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 containing the desired mutation comprises an scFv amino acid sequence set forth in SEQ ID NOs: 1-21.

14. The method of claim 13, wherein the antibody containing the desired mutation is mutant 4, which comprises the amino acid sequence of SEQ ID NO: 4, and the full-length IgG scFv sequence of mutant 4 comprises the amino acid sequence set forth in SEQ ID NO:

25.

15. A method for improving antibody-antigen interactions, comprising the method of any one of claims 1 to 14.

16. 15. An antibody, or composition thereof, comprising a mutant sequence of a known antibody against a target antigen, said antibody being engineered by the method of any one of claims 1 to 14.

17. The antibody of claim 16, wherein the antibody containing the desired mutation comprises an scFv amino acid sequence set forth in SEQ ID NOs: 1-21.

18. The antibody of claim 16, wherein the antibody containing the desired mutation is mutant 4, which comprises the amino acid sequence of SEQ ID NO: 4, and the full-length IgG scFv sequence of mutant 4 comprises the amino acid sequence set forth in SEQ ID NO:

25.

19. 19. The antibody of claim 17 or claim 18, comprising a mutant sequence of atezolizumab.

20. 20. The antibody of any one of claims 16 to 19, for use in treating a PD-L1 positive tumor.

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

22. 21. Use of an antibody, or composition thereof, according to any one of claims 16 to 20 in cancer immunotherapy CAR-T cell therapy or in enhancing the efficiency and enzymatic activity of CRISPR.

Images