Intelligent design and engineering modification of proteins

Abstract

Protein design and engineering processes include directed evolution, sequencing, and machine learning analysis. The process can explore directed evolution sequences that are not present in libraries, and can be used to facilitate the discovery of therapeutic proteins and imaging probes, and / or to enhance the efficiency of CRISPR and enzymatic activity. Through the flow, the scFv mutant which is tightly combined with the PD-L1 and / or has a slow dissociation rate on the PD-L1 is found. These scFvs can be used as strong anti-PD-L1 antibodies in cancer immunotherapy and / or CAR-T cell therapy.

Description

[0001] Cross-references to related applications

[0002] This application claims the benefits of U.S. Provisional Application No. 63 / 431,556, filed December 9, 2022, and U.S. Provisional Application No. 63 / 432,836, filed December 15, 2023, the entire contents of which are incorporated herein by reference.

[0003] Statement on Federally Funded Development Research

[0004] This invention was developed with government support from the National Institutes of Health (NIH) under license numbers R01GM111941 and R21AI158114. The government holds certain rights to this invention. Background Technology

[0005] Identifying proteins and small molecules with desired properties is crucial for the pharmaceutical and chemical industries. In recent years, machine learning-assisted methods have become increasingly popular. Neural network / deep learning models are widely used due to their flexibility, ability to learn complex relationships from massive datasets, and ease of implementation (thanks to the development of libraries such as PyTorch).

[0006] Despite its success, deep learning still has some limitations. First, most deep learning architectures do not quantify the uncertainty of predictions, tend to be "overconfident" during inference, and are "black boxes," making it difficult to determine why the model makes a specific prediction for a particular input. Second, and this is noteworthy, deep learning models have repeatedly demonstrated a lack of robustness to so-called "adversarial attacks" (i.e., tiny perturbations to the input data). For example, adding subtle noise imperceptible to the human eye can cause a deep learning model to misidentify a picture of a cat as guacamole or a lionfish as eggnog. In biology, protein structure prediction models have also been shown to be vulnerable to adversarial attacks. Finally, deep learning typically requires enormous computational resources. Many state-of-the-art models use hundreds of billions or even trillions of learning parameters. Often, most of these parameters can be removed without compromising performance, suggesting that more efficient methods are possible. For example, in natural language processing, it has been reported that 97% of the parameters of the BERT model were pruned, yet equivalent performance was still achieved.

[0007] Another model architecture is Gaussian process (GP) regression. A Gaussian process is a Bayesian model that defines a multivariate normal distribution on the possible functions that map input variables x to output variables y. In regions where training data is scarce or lacking, this model exhibits high uncertainty. In regions with abundant training data, the model is able to make predictions with higher confidence. The type of function in this distribution is determined by the kernel function and its hyperparameters (as with any other model) that need to be “fine-tuned” for the specific problem. The kernel function measures the similarity between any two input variables x.

[0008] Gaussian Processes (GPs) have at least four core advantages. First, similar to deep learning, a Gaussian process can approximate any relationship if equipped with a suitable kernel. Second, Gaussian processes can compute marginal likelihoods in closed-form, which is the probability that the data averages the value of the function over all possible values. Therefore, kernel hyperparameters can be "learned" by maximizing marginal and non-likelihoods, thus reducing the risk of overfitting and enhancing the model's robustness. Third, the model's predictions are generated by using the similarity between new data points and data points in the training set, where similarity is quantified by the kernel function. Unlike deep learning models, GPs are not black boxes because one can determine how the model measures the similarity between data points and generates predictions. The same kernel function used by the GP model for predictions can also be used to cluster the training data and retrieve the most "similar" data points in the training set, where "similarity" is determined by the kernel function we choose.

[0009] Fourth, GP can accurately quantify uncertainty and assign high uncertainty to data points that differ significantly from their training set. Uncertainty is crucial for protein engineering and drug design because evaluating model predictions experimentally is extremely costly in these situations. In such cases, it's preferable to rely solely on more reliable model predictions. GP detects "distribution shifts"—new data points that differ significantly from the training set—by assigning high uncertainty to its predictions; in contrast, deep learning models may "fail silently," generating low-precision predictions without explicit failure signals. Deterministic estimation can also be used for Bayesian optimization / active learning, where the model guides the operator in selecting data points for next experimental evaluation. Several techniques for estimating uncertainty in deep learning have been introduced in the literature. However, their uncertainty estimates are far inferior to the calibration provided by GP.

[0010] Gaussian Processes (GPs) have several limitations. The ability to choose a kernel function is both a disadvantage and an advantage, as designing a suitable kernel (e.g., for image classification) presents numerous challenges when sufficient prior knowledge is lacking. In some cases, this limitation might be overcome by combining deep learning with Gaussian processes; in fact, some practitioners have already explored this possibility in the literature. For example, deep learning models can be trained unsupervised to learn feature representations that can be used as inputs to GPs. In this strategy, the GP essentially acts as the final layer of the deep learning model, thus potentially combining some of the advantages of both approaches.

[0011] The main challenge facing GP is computational overhead. A simple and straightforward implementation is to scale to O(N) as the dataset size increases. 3 In fact, the poor scaling performance of general GPs is generally considered a reason why deep learning is preferred. For example, when building traffic prediction models, some people believe it's impossible to solve their problems using general GPs because "general GPs are difficult to scale to large datasets." The second major challenge is the lack of efficient kernels for sequential and graphical data. Many kernels for graphs have been described, but they typically exhibit quadratic scaling or worse scaling with respect to graph size.

[0012] One approach to approximate GP is the Random Fourier Feature (RFF) method by Rahimi and Recht. In short, the RFF method approximates the kernel function by using a random mapping corresponding to the target kernel function, thus converting each input data point into a "random feature" representation; then, Bayesian linear regression can be used in the random feature space to approximate GP, thereby achieving linear scaling of the training set size.

[0013] Monoclonal antibodies are among the most successful biotherapies. Despite their impressive versatility and specificity, the development of therapeutic antibodies still faces a variety of complex challenges. Typically, the initial target antibody has insufficient affinity, necessitating in vitro affinity maturation to enhance its binding force—a process involving multiple rounds of mutagenesis and screening to obtain antibodies with higher affinity. This process is often time-consuming and labor-intensive, potentially taking months to complete, and cannot simultaneously optimize other desired properties, such as good solubility and low immunogenicity. Computational techniques that can help identify high-affinity antibodies with desired properties more quickly may accelerate this process.

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

[0015] Machine learning-based approaches face at least two major challenges. First, there is a lack of estimation for the uncertainty in predictions incorporating affinity or other properties. Since training data only covers a small fraction of the sequence space, machine learning models typically perform poorly when asked to make inferences far beyond their training set. While deep learning models are flexible and powerful, they often (outside of their specific architecture) do not provide estimates of confidence intervals or uncertainty in predictions. As mentioned above, Gaussian process (GP) models have been suggested as an alternative, which do provide well-calibrated confidence intervals, but they are difficult to scale to large datasets without approximations and are often infeasible for datasets exceeding 5000 sequences.

[0016] The second key challenge is choosing a suitable input representation. Many different protein-coding schemes are described in the literature, including those recently employing language models; some are antibody-specific, while others are general to protein sequences. However, there is no consensus on which scheme is best suited for a specific problem. Classical one-hot encoding schemes are easy to implement, but because they treat each amino acid as completely different, they are too high-dimensional and lack sufficient information.

[0017] Therefore, high-affinity antibodies are typically identified through directed evolution, which may require multiple iterations of mutagenesis and screening to find the optimal candidate antibody. Deep learning techniques promise to accelerate this process, but existing methods cannot provide the confidence intervals or uncertainties needed to assess the reliability of predictions.

[0018] Programmed death-ligand 1 (PD-L1) is a protein that binds to another protein called PD-1, and this binding pathway suppresses the immune response against cells recognized as foreign. Many cancers overexpress PD-L1 to downregulate the immune response and prevent the immune system from self-destructing. Interfering with PD-L1-mediated immunosuppression has proven to be a highly effective approach for treating certain cancers, and several PD-L1-binding monoclonal antibodies (mAbs) have been approved by the FDA for various indications, including atezolizumab, durvalumab, and avelumab. However, there is a need to develop more effective anti-PD-L1 antibodies with higher affinity. Summary of the Invention

[0019] This disclosure provides a novel process that overcomes the limitations of the prior art described above. In some embodiments, this disclosure provides a novel process that includes novel machine learning / artificial intelligence methods for designing and engineering proteins with desired properties. More specifically, in the application of identifying high-affinity antibodies, a novel autoencoder model is used, trained on over 3 million human B-cell receptor sequences to efficiently encode antibody sequences. A variational Bayesian neural network is then developed for ordinal regression of directed evolution sequences and quantification of their likelihood of tight binding to antigens. Importantly, the model can also evaluate sequences not present in the directed evolution library, thereby greatly expanding the search space to discover optimal sequences for further evaluation.

[0020] In some embodiments, a 17-fold increase in KD of atezolizumab demonstrates the robustness of the procedure, and this success illustrates its potential to facilitate antibody development and other protein design / engineering tasks in general. In some embodiments, this disclosure provides a set of twenty-one (21) tightly bound PD-L1 scFv sequences (fusions of heavy and light chain variable regions). These scFvs are modified variants of the parent atezolizumab (mAb), discovered via a machine learning-assisted procedure disclosed herein for discovering antibodies that bind tightly to targets of interest. The best candidate measured from these twenty-one (21) candidate antibodies exhibited approximately 17-fold higher binding affinity than atezolizumab, implying one of the highest affinities among all currently available anti-PD-L1 mAbs. The dissociation rates (dissociation ratios) of these scFvs also appear to be lower than those of durvalumab and averulimumab. Given their exceptionally tight binding, these scFvs may have the potential to serve as tools for therapeutic or in vitro oncology studies.

[0021] Since the scFvs provided in this article function by forming a tight, reversible, non-covalent interaction with the PD-L1 protein expressed on the surface of human cells, this may enable them to treat cancer by interfering with the PD-L1 / PD-1 interaction that some cancers use to suppress the immune response, or alternatively, as an in vitro assay tool for cancer cells that overexpress PD-L1.

[0022] In addition to its potential use as an antibody therapeutic agent or in vitro tool, the scFvs described herein can also be used for CAR-T cell-based tumor therapy. CAR-T cells are modified immune cells that utilize proteins overexpressed on the surface of specific cancers to recognize and attack cancer cells. The field of CAR-T cell-based cancer therapy has seen significant development in recent years, and was first approved by the FDA in 2017. It is possible to engineer CAR-T cells to express the highly binding anti-PD-L1 antibody scFvs described herein, thereby recognizing cancer cells overexpressing PD-L1. In fact, CAR-T cells expressing anti-PD-L1 nanobodies have been reported to exhibit anti-tumor activity in vitro and in mouse models, suggesting this may be a viable strategy. Given the highly binding affinity exhibited by the scFvs of this disclosure, these scFvs may be well-suited for this task.

[0023] Therefore, this disclosure provides for the design and engineering of proteins for therapeutic agents (e.g., antibodies), imaging probes (e.g., protein domains), enhancing enzyme activity, and improving the efficiency of CRISPR.

[0024] Other systems, methods, features, and advantages of this disclosure will become apparent to those skilled in the art upon reading the following figures and detailed description. All such other / additional systems, methods, features, and advantages are to be included in this specification and fall within the scope of this disclosure, and are protected by the appended claims. Furthermore, all optional and preferred features and modifications of the described embodiments are applicable to all aspects of the disclosure herein. Moreover, the various features of the dependent claims and all optional and preferred features and modifications of the described embodiments are interchangeable.

[0025] Brief description of the attached figures

[0026] A better understanding of many aspects of this disclosure can be achieved by referring to the following accompanying drawings (also referred to as the “Figures”). The components in the Figures are not necessarily drawn to scale; the focus is on clearly illustrating the principles of this disclosure. Furthermore, in the Figures, the same reference numerals denote corresponding components in each view.

[0027] Figure 1A-1B : The number of training epochs required for fitting (the number of times the entire dataset is traversed). In all cases, a relative tolerance of 1e-6 is used as the convergence threshold.

[0028] Figures 2A-2DRelationship between relative error and number of training epochs: convergence speed of conjugate gradient (CG) with and without preconditions and stochastic gradient descent. AMSGrad was used here because in preliminary experiments, the closely related and commonly used Adam stochastic gradient descent method for training neural networks performed worse (failing to significantly reduce error within fewer than 100 training epochs).

[0029] Figure 3 : A schematic diagram of the Fast-Conv-1d kernel compared to a standard 3-layer convolutional neural network.

[0030] Figure 4 Comparison of the acceleration effect of conjugate gradients implemented using preconditioners constructed with SRHT and those constructed with preconditioners constructed with the unmodified algorithm S2.

[0031] Figure 5 The number of iterations to converge on various benchmark datasets using various settings of the preconditioner rank and three different variants of the preconditioner construction algorithm.

[0032] Figures 6A-6B : The number of training epochs / iterations required to achieve convergence using CG on different datasets, as a function of β / λ2.

[0033] Figure 7 The absolute percentage error distribution of marginal likelihood estimation in 112 evaluations using 25 or 40 probe vectors.

[0034] Figure 8 Kernel PCA was performed on the QM9 training set using the GraphConv-1d kernel.

[0035] Figures 9A-9B Kernel k-means clustering of the QM9 training set using GraphRBF kernel.

[0036] Figure 10 : Test query molecule 1, and the most similar training molecule.

[0037] Figure 11 : Test query molecule 2, and the most similar training molecule.

[0038] Figure 12 : Test query molecule 3, and the most similar training molecule.

[0039] Figure 13 : Test query molecule 4, and the most similar training molecule.

[0040] Figure 14 : Test query molecule 5, and the most similar training molecule.

[0041] Figure 15 : Schematic diagram of RESP calculation and experimental process, where U-PDL1 and B-PDL1 are unlabeled and biotin-labeled PDL1, respectively.

[0042] Figure 16A 1-2 and 16B 1-4 Initial library was screened using FACS and the quality of each rich cluster was assessed by flow cytometry. Figure 16A 1-2 Comparison of primary atezolizumab library and wild-type (WT) scFv on yeast cell surface screening (after 29 hours of dissociation rate competition). The right (highest) gate was used to collect mutants with slower dissociation rates, the middle gate was used to collect mutants with dissociation rates comparable to wild-type, and the left (lowest) gate was used to collect mutants with faster dissociation rates. The y-axis represents scFv expression (measured by V5 epitope / AF647 level), and the x-axis is PD-L1 binding strength / PE strength. % values ​​represent the percentage of total cells selected in each group. 16B 1-4 The binding strength of each sorted population relative to the wild type was tested after 8 hours of dissociation rate competition. LB represents low binding strength (faster dissociation rate), MB represents medium binding strength (dissociation rate comparable to wild type), and HB represents high binding strength (slower dissociation rate). 10 groups were analyzed for each population. 5 For each cell, the mean fluorescence intensity (MFI) on the x-axis is given in each figure, the y-axis represents V5 expression in terms of AF647 intensity, and the x-axis represents binding in terms of PE intensity.

[0043] Figure 17 Approximate Bayesian inference and learned embeddings are used for computer-directed evolution to efficiently evaluate candidate sequences. 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.

[0044] Figure 18 The structure of a task-adaptive autoencoder is as follows: The encoder module generates learned representations for each input sequence; the decoder module attempts to reconstruct the input sequence; and the classifier generates binary predictions for the decoys and human B-cell receptors in its training set.

[0045] Figures 19A-19CThe "human-ness" of the decoy and original sequences was scored using the autoencoder training set. 47,772 sequences were randomly selected from the autoencoder training set (half decoy, half human). The human-ness of these sequences was then scored using the AbLSTM model (Fig. 19A), the ANACRCI tool (Fig. 19B), and the BioPhi model (Fig. 19C) from the literature. In all three cases, the model scores for the decoy sequences were significantly different from those for the non-decoy sequences, and the "human-ness" of the decoy sequences was lower than that of the original sequences. In all three cases, the p-value was calculated to be 0.0 (meaning it is approximately 0 given a floating-point error) using the two-sided Man-Whitney U test implemented in Python Scipy library version 1.5.4. The following conventions apply to each box plot. The upper and lower boundaries of the boxes are the 25th and 75th percentiles of the data, respectively, and the beard is drawn at 1.5 times the interquartile range (the distance from the 25th percentile to the 75th percentile). Centered on the median of the data, a "gap" represents the 95% confidence interval of the median (determined by a nonparametric guided method). Diamonds represent "outliers" located beyond 1.5 interquartile ranges. Four asterisks indicate p-values ​​< 0.0001. Source data is provided as a source data file.

[0046] Figure 20A 1-2 -20C: Process Result Analysis. Figure 20A 1-2 Marginal distribution of positions within the two main clusters. The marginal distribution shows that each main cluster contains multiple sub-clusters. Only the locations of mutations within each cluster are shown in the figure. Positions are first numbered using the wild-type sequence number, then followed by the Chothia number in parentheses. Figure 20B shows the score distribution assigned to sequences in each binding category by a variational Bayesian model designed for ordinal regression (classification by ordinal category). Instead of directly measuring dissociation rates, the experimental data categorizes sequences into three classes: RH01 (high dissociation rate), RH02 (medium dissociation rate), and RH03 (low dissociation rate). Higher scores assigned by the ordinal regression model indicate higher confidence in the sequence having a low dissociation rate, while lower scores indicate higher confidence in a high dissociation rate. Figure 20C: Acceptable score versus iteration for a typical simulated annealing chain. The algorithm initially explores the sequence space unbiasedly, gradually becoming forced to focus on the most promising regions discovered so far as the temperature decreases. Source data is provided as source data files for all subgraphs.

[0047] Figure 21 The locations of 21 mutant mutations in the atezolizumab heavy chain structure. The mutated residues have been tagged in the heavy chain. Structure from Mol*Viewer. 46 Chinese RCSB Protein Database (5XXY) 71 .

[0048] Figures 22A-22C : Yeast surface K off and K D Experimental verification. Figure 22A: Dissociation of wild-type (WT) and mutant 4scFv on yeast surface at room temperature for 92 hours, T 1 / 2 The half-life is shown in Figure 22B. Figure 22C shows the dissociation of WT atezolizumab, mutant 4, durvalumab, and averumumab scFv on yeast surfaces after 28 hours at room temperature. Figure 22C shows the binding affinity (Kb) between scFv and PD-L1 on yeast surfaces. D ) Measurements (3 independent measurements, see also Figures 29A-29F).

[0049] Figure 23A 1-2 Figure 23A and 23B: Variational Bayesian Neural Network Architecture. 1-2 A comparison of variational Bayesian neural networks and fully connected architectures. Both networks map a vector of input values ​​x1, x2, ... xn to hidden layer vectors H1_1, H1_2, ... H1_n through matrix multiplication and subsequent nonlinear activation. However, in a fully connected network, each element of the hidden layer weight matrix is ​​a learned value that becomes a fixed value 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 the learned mean and learned standard deviation. To generate predictions, samples can be taken from the weight distribution, thus providing an estimate of the prediction uncertainty. Figure 23B: Structure of the ordinal regression scoring model used in the process.

[0050] Figures 24A-24B WT atezolizumab scFv binds to human PD-L1 on the surface of yeast. The y-axis of each subplot represents the expression of scFv (V5 epitope / AF647), and the x-axis represents the binding strength 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 in the total analysis is shown below each subplot.

[0051] Figures 25A-25D The figures show the initial screening of 21 member libraries after 8 hours of dissociation rate competition (relative to WT, top figures 25A and 25B) and the final screening after 39 hours of competition (bottom figures 25C and 25D). Right-angled triangles represent sorting gating and the percentage of cells collected at each sorting stage. The y-axis of each subplot represents scFv (V5 epitope / AF647) expression, and the x-axis represents the binding strength with human PD-L1 / PE.

[0052] Figure 26A 1-2-26B: BLI binding analysis of soluble WT and mutant 4scFv with immobilized biotin PD-L1. Figure 26A 1-2 Binding (binding and dissociation) curves for different scFv concentrations. The titration range for WT was 0.65 to 475 nM, while that for mutant 4 was 0.219 to 160 nM. Figure 26B: Binding / dissociation rates and KB of WT and mutant 4 scFv. D For each scFv, two different K were obtained in the two binding events targeting PD-L1. D Value. Because the dissociation rate of the mutant 4scFv is very slow, its K value cannot be obtained. D value.

[0053] Figure 27A-27G Mass spectrometry analysis of PD-L1 and WT / mutant 4scFv monomeric / dimeric / multimeric states. Each sub-plot includes approximate molecular weight (MW) (x-axis) and count (y-axis). PD-L1 molecular weight (theoretical value) = 35-38 kDa (range varies depending on glycosylation), scFv = 31.4 kDa. Batch 1 and batch 2 were derived from two separate purified scFv samples. σ (sigma) indicates the standard deviation of the molecular weight estimate calculated by the software, and the numbers above each peak represent the apparent MV.

[0054] Figures 28A-28B SDS-PAGE gel (4-20%) of purified scFv-His6 fusion protein. M4 is mutant 4. The two gel images represent proteins from two different purifications (two different SDS-PAGE gels from two different purifications).

[0055] Figures 29A-29F: Binding curves of PD-L1 to yeast surface scFv. MFI represents fluorescence intensity, and M4 represents mutant 4. Each experiment was conducted individually / independently on different dates.

[0056] Figures 30A-30B Side-by-side comparison of amino acids close to R98C in the wild type (PDB structure 5XXY, Fig. 30A, top) and the IgFold predicted structure of the highest-scoring mutant (Fig. 30B, bottom).

[0057] Figure 31 : Overlay of the PDB structure 5XXY with the IgFold predicted structure of the highest-scoring mutant. Note that the conformational shift in CDRH3 is relatively large, while the conformational shifts in other regions in contact with the antigen are more subtle.

[0058] Figures 32A-32C: Scores for sequence assignments in the model as described in Mason et al., Nat. Biomed. Eng. 5, 600-612, 2021. Figure 32A: Score distribution for binding and non-binding sequence assignments in the training set. 0 is the cutoff value for predicting a sequence as a binding. Figure 32B: Score distribution for sequences selected by Mason et al. as possible binding sequences. Experiments show that the binding KD of all these sequences by Mason et al. is <20 nM, and they are indeed predicted as binding sequences according to the model disclosed in Example 2. Figure 32C: Uncertainty of the model's score for its predicted assignments on the test set. The associated uncertainty of incorrect predictions is higher than that of correct predictions, indicating that, as expected, the uncertainty of model assignments can help determine whether the predictions are reliable. There are 3,158 unique sequences in the test set. Significance was assessed using a two-sided Mann-Whitney U test; the p-value obtained was 1e-23. The following conventions apply to each box plot. The upper and lower boundaries of the box represent the 25th and 75th percentiles of the data, respectively. The stratigraphic line is drawn at 1.5 times the interquartile range (the distance from the 25th to the 75th percentile). The center is drawn at the median of the data, and a "gap" represents the 95% confidence interval of the median (determined by a nonparametric guide method). Diamonds represent outliers located beyond 1.5 times the interquartile range. Four asterisks indicate a p-value < 0.0001.

[0059] Figure 33 Clustering of simulated annealing results. A dendrogram of the clustering results obtained from the sequences obtained after the simulated annealing process, after removing sequences with scores correlated with wide confidence intervals. This dendrogram shows that the modified simulated annealing algorithm identified at least two main spatial regions of the sequences.

[0060] Figure 34 Example of FACS gating policy. From 6x10 6 106 cells were collected from the sorted yeast cells (“high (hi)” phylum). The Y-axis represents V5 epitope expression, and the X-axis represents PD-L1 binding. V5-tagged monoclonal antibody was detected using the secondary antibody AF647 polyclonal IgG. Biotin-PD-L1 was detected using SA-PE reagent.

[0061] Other advantages of this disclosure are set forth in part in the following description, and in part will be apparent from the description or may be learned by practice of this disclosure. The advantages of this disclosure can be realized and obtained by means of the elements and combinations thereof particularly pointed out in the appended claims. It should be understood that the foregoing general description and the following detailed description are exemplary and explanatory only, and do not constitute a limitation on the invention claimed in the claims. Detailed Implementation

[0062] Gaussian processes (GPs) are Bayesian models that offer numerous advantages for regression tasks in machine learning, such as reliable quantification of uncertainty and higher interpretability. Their adoption has been hampered by their high computational cost and difficulty in applying them to analyze sequences (e.g., amino acid sequences) and graphs (e.g., small molecules). In some implementations, this disclosure provides a set of random feature approximation kernels that scale linearly with both the size of the training set and the size of the sequence or graph. These new kernels are integrated into the new Python library xGPR for GP regression, and an efficient and scalable algorithm is developed for fitting GPs equipped with these kernels to large datasets. The performance of xGPR was compared with standard and state-of-the-art deep learning models across 17 different benchmarks, demonstrating that GP regression achieves highly competitive accuracy on these tasks while providing well-calibrated uncertainty quantification and higher interpretability. Furthermore, this disclosure provides how xGPR can be used as part of an active learning strategy to automate the design of proteins with desired properties without human intervention.

[0063] This disclosure also provides a novel workflow called RESP, which comprises directed evolution, sequencing, and machine learning analysis for the efficient identification of high-affinity antibodies. In some embodiments, a learned characterization based on training with over 3 million human B-cell receptor sequences is developed to encode antibody sequences. Furthermore, a variational Bayesian neural network is developed for ordinal regression on a set of directed evolution sequences binned according to dissociation rate, quantifying their likelihood of tight binding to the antigen. Importantly, this model can evaluate sequences not present in the directed evolution library, thereby significantly expanding the search space to discover optimal evaluation sequences.

[0064] The potential of RESP in facilitating the development of universal antibodies is further elaborated. In some implementations, the affinity of atezolizumab (trade name Tecentriq) heavy chain mutants for programmed death-ligand 1 (PD-L1) is enhanced through engineering, a process applied to improve the affinity of known antibodies for well-studied antigens. To this end, we randomized most of the atezolizumab heavy chain sequence and binned the mutants using yeast display and FACS, then sampled and sequenced each bin. The mutant atezolizumab sequence data were converted to a low-dimensional representation using an autoencoder model, and a Bayesian ordinal regression model was trained using the encoded mutant atezolizumab sequences, scoring each sequence based on the probability of strong binding to PD-L1. Finally, a modified simulated annealing algorithm was used to select sequences for testing.

[0065] In some embodiments, this 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 workflow discussed herein and have shown tight binding to PD-L1. One of the scFvs disclosed herein shows the K-type of the PD-L1 antibody atezolizumab. D The dissociation rate was increased by 17-fold, while the dissociation rate was reduced by 10-fold. This scFv disclosed in this paper may be more effective than any existing PD-L1 antibody and could serve as an effective anti-PD-L1 antibody for cancer immunotherapy and CAR-T cell therapy.

[0066] To those skilled in the art, the teachings presented in the foregoing description and drawings will enable them to conceive of numerous modifications and other embodiments of the compositions and methods disclosed herein. Therefore, it should be understood that this disclosure is not limited to the specific embodiments disclosed herein, and various modifications and other embodiments are intended to be included within the scope defined by the appended claims. Those skilled in the art will recognize many variations and modifications to the aspects described herein. These variations and modifications are intended to be incorporated into the teachings of this disclosure and are included in the claims herein.

[0067] Although specific terms are used herein, they are used in a general and descriptive sense only and are not intended to limit the invention.

[0068] As will be apparent to those skilled in the art upon reading this disclosure, each individual embodiment described and illustrated herein has discrete components and features that can be readily separated from or combined with features of any of the remaining plurality of embodiments without departing from the scope and spirit of this disclosure.

[0069] Any method described may be performed in the order of the events described or in any other logically feasible order. That is, unless expressly stated otherwise, no method or aspect described herein should be construed as requiring its steps to be performed in a particular order. This applies to any possible non-express interpretation, including logical questions relating to the arrangement of steps or operational procedures, plain meanings derived from grammatical organization or punctuation, or the number or type of aspects described in the specification.

[0070] All publications and patents referenced in this specification are intended to disclose and describe methods and / or materials in relation to the referenced publications. All such publications and patents are incorporated herein by reference as if each publication or patent were expressly and individually indicated to be incorporated by reference. Such inclusion by reference is expressly limited to the methods and / or materials described in the referenced publications and patents and does not extend to any dictionary definitions in the referenced publications and patents. Any dictionary definitions in the referenced publications and patents, unless expressly repeated in this application, should not be considered definitions, nor should they be construed as defining any terms appearing in the appended claims. References to any publication refer to its disclosure prior to the filing date and should not be construed as an admission that this disclosure is not entitled to a prior disclosure prior to that publication. Furthermore, the publication date provided may differ from the actual publication date, which may require separate confirmation.

[0071] Although aspects of this disclosure may be described and claimed in a particular statutory category (e.g., the system statutory category), this is merely for convenience, and those skilled in the art will understand that aspects of this disclosure may be described and claimed in any statutory category.

[0072] It should also be understood that the terminology used herein is for descriptive purposes only and is not intended to be limiting. Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which the disclosed compositions and methods pertain. It should also be understood that terms (e.g., those defined in common dictionaries) should be interpreted as having the same meaning as they have in the context of the specification and the relevant field, and should not be interpreted in an idealized or overly formal sense unless explicitly defined herein.

[0073] Unless otherwise stated, all aspects of this disclosure will employ techniques within the scope of the art, including molecular biology, microbiology, organic chemistry, biochemistry, physiology, cell biology, and vascular biology. These techniques have been described in detail in the literature.

[0074] Before describing various aspects of this disclosure, the following definitions are provided and should be used unless otherwise stated. Other terms may be defined in other parts of this disclosure.

[0075] definition

[0076] As used herein, the term “comprising” should be interpreted as specifying the presence of the stated feature, integer, step, or component, but does not preclude the presence or addition of one or more features, integers, steps, or components, or combinations thereof. Furthermore, terms such as “comprising,” “including,” “consisting of,” “containing,” “involving,” and “for example” are used in their open, non-limiting sense and may be used interchangeably. Additionally, the term “comprising” is intended to cover examples and aspects covered by “substantially constitutes” and “consisting of.” Similarly, the term “substantially constitutes” is intended to cover examples covered by “consisting of.”

[0077] As used herein, the term "and / or" includes any one or more of the associated listed items and all combinations thereof. Expressions such as "at least one / a" modify the entire list of elements, rather than individual elements within the list, when used before a list of elements.

[0078] Unless the context clearly specifies otherwise, the singular forms “a,” “an,” and “the” used in the specification and appended claims include the plural forms. When referring to “a” compound, therapeutic agent, or pharmaceutical composition, it means one or more molecules of that compound, therapeutic agent, or pharmaceutical composition, and not limited to a single compound, therapeutic agent, or pharmaceutical composition. These one or more molecules may be the same or different, provided they belong to the class of that compound, therapeutic agent, or pharmaceutical composition. Thus, for example, “a” therapeutic agent should be interpreted as including one or more molecules of that therapeutic agent, wherein the therapeutic agent molecules may be the same or different (e.g., containing different isotopic abundances and / or different degrees of hydration, or in equilibrium with different conjugate base or conjugate acid forms).

[0079] It is important to note that proportions / ratios, concentrations, quantities, and other numerical data may be expressed in ranges. It should also be understood that the endpoints of each range are significant, both related to and independent of the other endpoint. Furthermore, it should be understood that multiple values ​​are disclosed herein, and each value is disclosed not only as the value itself but also as “about” that specific value. For example, if the numerical value “10” is disclosed, then “about 10” is also disclosed. In this document, a range may be expressed as beginning approximately to one specific value and / or ending approximately to another specific value. Similarly, when a numerical value is expressed as an approximation using the antecedent “about”, it should be understood that the specific value constitutes another aspect. For example, if the numerical value “about 10” is disclosed, then “10” is also disclosed.

[0080] When describing a range, this also includes a range from one specific value and / or to another specific value. It should be understood that when a numerical range is given, unless otherwise expressly stated herein, all intermediate values ​​between the upper and lower limits of the range, up to one-tenth of the lower limit, and any other notations or intermediate values ​​of the range are included within the scope of this disclosure. Smaller ranges may independently contain the upper and lower limits of these smaller ranges, which are also within the scope of this disclosure and are subject to limits explicitly excluded by the range. When a range includes one or two limits, this disclosure also includes ranges excluding one or both of those limits. For example, when a range includes one or two limits, ranges excluding one or both of those limits are also included within this disclosure; for example, the phrase "x to y" includes a range from "x" to "y" and a range greater than "x" and less than "y". The range may also be expressed as an upper limit, such as "about x, y, z or less", and should be understood to include specific ranges of "about x", "about y", and "about z" as well as ranges of "less than x", "less than y", and "less than z". Similarly, the phrase “about x, y, z or greater” should be understood to include the specific ranges of “about x,” “about y,” and “about z,” as well as the ranges of “greater than x,” “greater than y,” and “greater than z.” Furthermore, the phrase “about x to y” (where “x” and “y” are numerical values) includes “about x to y.”

[0081] It should be understood that this range format is for the sake of brevity and should therefore be flexibly interpreted to include not only the values ​​explicitly listed as range limits, but also all individual values ​​or subranges within that range, as if each value and subrange had been explicitly listed. For example, the range "about 0.1% to 5%" should be understood to include not only the explicitly listed values ​​of about 0.1% to about 5%, but also the individual values ​​within that range (e.g., about 1%, about 2%, about 3%, and about 4%) and their subranges (e.g., about 0.5% to about 1.1%; about 5% to about 2.4%; about 0.5% to about 3.2%; about 0.5% to about 4.4%, and other possible subranges).

[0082] As used herein, the terms "about," "approximately," and "substantially," when used with a numerical variable, generally refer to the numerical value of that variable, and all values ​​of that variable within the experimental error range (e.g., within a 95% confidence interval of the mean) or within + / - 10% of a specified value (whichever is greater). The terms "about," "approximately," "about," and "substantially" as used herein can indicate that the quantity or value discussed can be an exact value or a value that provides an equivalent result or effect as described in the claims or taught herein. That is, it should be understood that quantities, sizes, formulations, parameters, and other quantities and characteristics are not and need not be exact, but can be approximate and / or larger or smaller as needed, reflecting tolerances, conversion factors, rounding, measurement errors, and other factors known to those skilled in the art to obtain an equivalent result or effect. In some cases, it is not reasonable to determine a numerical value that provides an equivalent result or effect. Generally, whether explicitly stated or not, quantities, sizes, formulations, parameters, or other quantities or characteristics should be described using "about," "approximately," or "about." Unless otherwise explicitly stated, when “about,” “approximately,” or “approximately” is used before a quantitative value, the parameter also includes the specific quantitative value itself.

[0083] As used herein, the terms “optional” or “optionally” indicate that the event or situation described below may or may not occur, and the description includes both the occurrence and non-occurrence of the event or situation.

[0084] The terms “object,” “individual,” or “patient,” used interchangeably in this document, can refer to a vertebrate, such as a mammal (e.g., a human). “Object” can also refer to a cell, cell group, tissue, organ, or organism, preferably a human being and its constituent parts.

[0085] As used herein, the terms “treatment” and “management” generally refer to achieving a desired pharmacological and / or physiological effect. This effect may be (but is not necessarily) preventative, i.e., prevention or partial prevention of a disease, symptom, or condition. This effect may be therapeutic, i.e., partial or complete cure of a disease, condition, symptom, or adverse reaction caused by that disease, disorder, or condition. As used herein, the term “treatment” can include any treatment of any disease-related inflammation in a subject (particularly a person) and may include one or more of the following: (a) prevention of disease occurrence in subjects who may be susceptible to the disease but have not yet been diagnosed; (b) suppression of disease, i.e., halting its development; and (c) alleviation of disease, i.e., reducing or improving the disease and / or its symptoms or condition. As used herein, the term “treatment” can refer to therapeutic treatment alone, preventative treatment alone, or both therapeutic and preventative treatment. An individual requiring treatment (the subject requiring treatment) can include an individual who already has the disease and / or an individual who needs to prevent the disease. As used herein, the term “treatment” can include suppression of a disease, disorder, or condition, e.g., halting its progression; and alleviation of a disease, disorder, or condition, e.g., causing the disease, disorder, and / or condition to subside. Treatment of a disease, disorder, or condition may include improving at least one symptom of the particular disease, disorder, or condition, even if its underlying pathophysiology is not affected. For example, treating an individual’s pain by administering an analgesic, even if the analgesic does not treat the cause of the pain.

[0086] As used herein, the term "therapeutic effective dose" refers to a dose sufficient to achieve the desired therapeutic effect or to act on adverse symptoms, but generally insufficient to cause adverse side effects. For any given patient, the specific therapeutic effective dose level will depend on a variety of factors, including the disease being treated and its severity; the specific composition used; the patient's age, weight, general health condition, sex, and diet; the timing of administration; the route of administration; the excretion rate of the specific compound used; the duration of treatment; medications used in combination with or concurrently with the specific compound used; and similar factors that may be well known in the medical field and fall within the knowledge and expertise of a medical practitioner. In the treatment of a specific disease or condition, in some cases, the expected response may suppress the progression of the disease or condition. This may involve only a temporary slowing of disease progression. However, in other cases, it may be necessary to permanently halt disease progression. This can be monitored using routine diagnostic methods known to those skilled in the art for any disease. The expected response to treatment of a disease or condition may also be to delay the onset of the disease or condition or even to prevent its onset.

[0087] It should be understood that, unless otherwise stated, all temperatures referred to in this article are based on atmospheric pressure (i.e., one atmosphere).

[0088] The various aspects of this disclosure have been described above. In general, the following embodiments will further illustrate some other and / or more detailed aspects of this disclosure. While various aspects of this disclosure have been described in conjunction with the following embodiments and corresponding text and drawings, this does not mean that the various aspects of this disclosure will be limited to these descriptions. Rather, all alternatives, modifications, and equivalents are intended to cover the spirit and scope of this disclosure.

[0089] Example

[0090] The following examples are intended to provide those skilled in the art with a complete disclosure and description of the preparation and evaluation of the compounds, compositions, articles, devices, and / or methods claimed herein. These embodiments are merely illustrative of the invention and are not intended to limit the scope of the disclosure as perceived by the inventors. We have made every effort to ensure the accuracy of numerical values ​​(e.g., content, temperature, etc.), but some errors and deviations should still be considered. Unless otherwise stated, parts are by weight, temperatures are in °C or ambient temperature, and pressures are atm or near atm.

[0091] Example 1

[0092] Linear scaling kernels for protein sequences and small molecules outperform deep learning, while also providing quantification of uncertainty and Improved interpretability

[0093] This embodiment discloses an approximate Gaussian process model that achieves good scalability and performance for both sequences and graphs. First, an efficient method is introduced for fitting a stochastic approximate Gaussian process to a preprocessed, conjugate gradient (PCG) method with novel preconditioners. The scalability of the resulting approximate Gaussian process to larger datasets is demonstrated, showing that it outperforms another approximate Gaussian process method on nine different datasets, and achieving the same test set performance as the exact (non-approximate) Gaussian process on protein engineering tasks.

[0094] Then, a set of efficient, scalable, and stochastic feature approximations of kernels were introduced for graphs and sequences. We used these kernels to compare the performance of the Gaussian Process Regression (GP) model with deep learning methods on 17 different benchmark sets. The results show that the GP regression model achieves highly competitive performance.

[0095] By using embeddings as input to an approximate Gaussian process, the possibility of combining deep learning with Gaussian processes is further explored, and it is shown that this is generally beneficial. It is further demonstrated that the uncertainty estimate provided by the approximate Gaussian process has higher calibration than that provided by various deep learning uncertainty estimation techniques. A simple study demonstrates how to "engineer" a more fit protein using an approximate Gaussian process without human intervention. Furthermore, it illustrates how to use an approximate Gaussian process to cluster the training set and retrieve the data points in the training set most similar to the test points.

[0096] method

[0097] Random Feature Approximation: Background

[0098] In a typical Gaussian process model, given a training set X with relevant y-values, for a new data point x... * The relevant y-value f * Make a prediction:

[0099] p(f * |X, y, x * )=N(u * , σ * (1)

[0100] Where N(u) * , σ * The distribution is Gaussian, with mean and variance as follows:

[0101] u * =k * T (λ 2 I+K) -1 y and σ * =k ** -k * T (λ 2 I+K) -1 k * (2)

[0102] Here, k * It is a vector of length N given N training data points, which is obtained by analyzing each x in the training set. i Calculate the kernel function k(x) * x i ) is derived; K is an N×N matrix, where the elements ij are obtained by applying k(x) to the training set. i y j The kernel function is calculated at points. λ is a kernel hyperparameter representing the expected noise level in the data. *The kernel function, k, is evaluated over all training data points on the new data point. ** It is a kernel function that evaluates itself on new data points. During fitting, it does not need to "learn" any parameters, but it requires constructing and reversing an N×N matrix, which can be computationally expensive if N is large.

[0103] To avoid this cost, use random Fourier features. 27 This approximates the kernel function. This well-known approximation assigns a random mapping to a given kernel such that k(x1, x2) ≈ z(x1). T z(x2) (See Supporting Information Section S1 for details). Now use the following formula to predict the mean and variance of the new data points:

[0104]

[0105] Where Z is the transformed input matrix, and thus the i-th row is z(x i ), z * It is z(x) * ), where x * Here are new data points, and λ is the kernel hyperparameter. For common stationary kernels, the approximation error decreases exponentially with the increase of the number of random features, thus resulting in diminishing returns (increasing the number of random features from 1,000 to 2,000 provides a greater relative improvement than increasing it from 2,000 to 4,000, etc.).

[0106] Furthermore, the original random features described above are replaced by the Structured Orthogonal Random Features (SORF) procedure. In short (see Supporting Information Section S2 for details), taking the RBF kernel as an example, the following transformation is applied to each input vector:

[0107] SHD1HD2HD3 (4)

[0108] Where S is a diagonal matrix whose elements are derived from a χ² distribution with d degrees of freedom, H is a normalized Hadamard matrix, and D1, D2, and D3 are diagonal matrices whose elements are derived from a Rademacher distribution. Hadamard matrix multiplication can be replaced by the Fast Hadamard Transform, thereby reducing the cost of generating random features from O(NM) to O(NM). 2 The scaling factor is reduced to O(NM logM). The resulting model is very lightweight, requiring only four diagonal matrices to be stored. Generating predictions using the trained model is fast, as it only requires a series of transformations and diagonal matrix multiplications followed by a dot product. This change benefits not only speed but also performance. This modification has proven to improve the accuracy of the stochastic feature kernel approximation of the RBF kernel. Fast Hadamard transforms for CPU and GPU are implemented in the library (see Support Information section S3 for benchmarks and other details).

[0109] Faster fitting using new preconditioners

[0110] To fit an approximate GP model, [Z] must be found. T Z+λ 2 I] -1 Z T y = w. If the number of random features is large (e.g., 5000-30000), directly generate and invert the matrix [Z]. T Z+λ 2 [I] This requires expensive matrix multiplication and decomposition. This problem can be avoided by using the well-known conjugate gradient method (CG), which achieves sub-quadratic scaling of the random number of features and iteratively solves for Z of w. T y = [Z T Z+λ 2 In short (see Supporting Information section S4 for details), this algorithm does not compute [Z]. T Z+λ 2 I] -1 Instead, it requires calculating the matrix-vector product [Z] in each iteration. T Z+λ 2 Matrix multiplication or decomposition is never required.

[0111] Note that for the Nth row vector in Z... In other words, the matrix-vector product can be formed as the sum of the products of individual rows. The dataset is then processed in mini-batches. For each mini-batch, random features are generated using a Hadamard transform-based procedure, and matrix-vector multiplication is performed using w, which is then added to the result. This implementation also facilitates parallelization. In each iteration of the conjugate gradient, each worker node can be assigned its own copy of the kernel, its own copy of w, and its own subset of data to process. At the end of the iteration, the results of all worker nodes are summed and used to update w for the next iteration.

[0112] Please note that unlike fitting methods based on stochastic gradient descent, mini-batch size does not affect model performance or the number of iterations required for convergence because stochastic gradient descent evaluates the full gradient in each iteration. Larger mini-batches can slightly reduce the fitting time but increase memory usage. Therefore, memory usage is independent of dataset size, and users can decide this by choosing the appropriate mini-batch size.

[0113] Since CG is iterative, it is important to minimize the number of iterations to fit the model as quickly as possible. xGPR uses the matrix [Z] in Equation (1). T Z+λ 2The stochastic Nystrom approximation of [I] is used as a preconditioner (see Supporting Information Section S4 for details). Although this is a known matrix approximation method, it has never been used before as a preconditioner for Gaussian processes.

[0114] We also modified the stochastic Nystrom approximation to reduce the time required to construct the preconditioner (see Support Information S4 for details and benchmarks). Importantly, in this implementation, the preconditioner is constructed by loading only a small batch of data into memory at a time, thus minimizing memory footprint. See Support Information S5 and S6 below for details on how to use stochastic features, and more specifically, preconditioned conjugate gradients, for hyperparameter fine-tuning.

[0115] exist Figure 1A and 1B In this section, we compare the preconditioned conjugate gradient (PCG) method with L-BFGS implemented on seven datasets in the SciPy library. L-BFGS (Limited-Memory BFGS) is a commonly used algorithm that minimizes the loss function by iteratively constructing an approximation of the Hessian matrix. Compared to L-BFGS or unconditional CG, the PCG strategy reduces the number of iterations required for convergence, thereby shortening the fitting time by two orders of magnitude, achieving a tight fit in just 10-20 iterations. In all cases, a relative tolerance of 1e-6 is used as the convergence threshold. For details on preconditioned subranks, please refer to Support Information S4 below.

[0116] exist Figures 2A-2D In this paper, the PCG method is further compared with two stochastic gradient descent methods: stochastic variance-reducing gradient descent (SVRG) and AMSGrad, the latter being an improved variant of the Adam algorithm commonly used to train neural networks. The figure shows how the residuals decrease as the number of training epochs (a complete traversal of the dataset) increases. For PCG, each iteration corresponds to one training epoch, while for stochastic gradient descent, each training epoch contains multiple iterations. As shown, PCG clearly outperforms stochastic gradient descent, converging several orders of magnitude faster, thus significantly reducing the number of training epochs required for fitting. Unlike stochastic gradient descent, PCG does not require fine-tuning the learning rate to achieve good results. Combining the PCG method with a stochastic feature generation method can fit datasets of over 400,000 data points in less than 7 minutes, including all hyperparameter fine-tuning. With mini-batches of 2,000 data points, the process consistently requires less than 1.5GB of RAM.

[0117] These studies used AMSGrad because, in preliminary experiments, the closely related and commonly used Adam stochastic gradient descent method for training neural networks performed worse (failing to significantly reduce error within fewer than 100 training epochs).

[0118] It's important to note that fitting 500,000 data points using the exact Gaussian process requires constructing a 400,000 × 400,000 kernel matrix, which would require 1.3TB of memory and the entire kernel matrix to fit the model. Therefore, xGPR is more scalable than the exact Gaussian process. In fact, as described in "Comparison with Stochastic Variational Inference and Exact Gaussian Processes" below, when using the exact Gaussian process implemented in the PyTorch library, attempting to fit a dataset containing only 21,000 training data points immediately results in an out-of-memory error, while xGPR can fit 500,000 data points without any problems.

[0119] result

[0120] Comparison with stochastic variational inference and exact Gaussian processes

[0121] Random features are not the only way to approximate a Gaussian process; another approach is stochastic variational inference (SVI), which uses a set of “induced points” much smaller than the number of training data points to approximate GP. As described in Supporting Information S8 below, both xGPR and the approximate GP use stochastic variational inference (implemented in GPyTorch) to fit nine different datasets ranging in size from 400 to 500,000 data points. To ensure a fair comparison, the same kernel function is used in all cases (see Supporting Information S7 for details on these datasets and all other datasets used for benchmarking in this example).

[0122] These studies were conducted on a GTX 1070 GPU equipped with 8GB of RAM. Two of the datasets were small enough (<5,000 data points) to be fitted with an exact Gaussian process, so the exact GP was also applicable to these two datasets. An attempt was made to fit a third dataset containing approximately 21,000 data points with the exact GP, but the GPyTorch library immediately threw an out-of-memory error, indicating that the exact GP is quite costly for moderately sized datasets.

[0123] The study found that xGPR achieves the same or better performance as SVI in every case; SVI outperforms xGPR in all cases. For two datasets where xGPR can be compared with exact GP, xGPR's performance is very similar, while SVI performs similarly in one case but significantly worse in another.

[0124] xGPR was also compared with exact Gaussian processes in the literature for active learning / Bayesian optimization tasks. For example, data on 52,910 green fluorescent protein (GFP) mutants were generated, and the fluorescence of each mutant was measured. An exact Gaussian process was fitted to this dataset using the RBF kernel, with 1,115 sequences containing only one mutation as the training set, and the learned embeddings representing the input data points. In Bayesian optimization, a collection function was used to select new data points for experimental evaluation. A confidence upper bound (UCB), which is the predicted mean of GP plus a confidence interval, was used. The UCB prioritizes data points where the model predicts high fitness values ​​for mutant sequences and data points where the model has high uncertainty, thus striking a balance between utilizing existing knowledge and exploring the sequence space. The results showed that the UCB values ​​on the test set correlated quite well with the actual measured fluorescence values.

[0125] Using the same training set, kernel function, and representation, xGPR was fitted to the same dataset. The results showed that the Spearman's-r correlation coefficient between the UCB acquisition function and the fluorescence values ​​measured in the test set was 0.78. Surprisingly, xGPR performed slightly better, with a coefficient of 0.8. Next, the 50 mutants with the highest UCB values ​​were selected from the test set (to simulate experiments and evaluate the model's predictive performance). The mean log fluorescence value of the top 50 xGPR mutants was roughly the same as reported in the literature (around 3.75). These results indicate that when using the same kernel function, xGPR performs comparably to the exact Gaussian process while achieving significantly improved scalability.

[0126] The first group of novel, highly efficient kernels for sequences that cannot be compared.

[0127] So far, we have only considered kernels for fixed-length inputs. In protein engineering problems, multiple sequence alignments can typically be formed from all input sequences. Alternatively, if the input is represented using learned embeddings generated by a large language model, the learned embeddings can be averaged over the entire sequence to generate a fixed-length input. However, in many other cases, sequence lengths vary and alignments are not easy, for example, if the input sequences come from different families. Two stochastic feature approximation kernels are introduced that can efficiently handle such problems: the first is called FHT-Conv-1d, and the second is called Fast-Conv-1d.

[0128] The first type of kernel can be constructed as follows. Consider extracting all k-mers of a given length (e.g., 9) from two sequences and calculating the similarity of each possible k-mer pairing in the two sequences using any common fixed vector kernel of interest. Several different kernels are available for this type of comparison, including the RBF kernel, Matern kernel, neural network kernel, and inverse cosine kernel, all of which can be approximated using random features. The "convolutional kitchen sinks (ConvKS)" kernel is a special case of this type. In this case, the comparison between the two sequences is as follows:

[0129]

[0130] Here, k is a fixed vector kernel used to measure the similarity between any two subsequences. As described, the scaling factor of this kernel is proportional to the square of the sequence length. Therefore, it is impractical to use it for exact Gaussian processes or stochastic variational inference, because in these cases it needs to be implemented in the manner described, resulting in a large number of pairwise k-mer comparisons. For example, if comparing two protein sequences with a length of 250 and a subsequence length of 9, all 242 k-mers in both sequences need to be extracted, requiring more than 29,000 pairwise k-mer comparisons.

[0131] The following section demonstrates how to implement this kernel using a random feature approximation that scales linearly with the sequence length. Using the random feature approximation, the kernel becomes:

[0132]

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

[0134]

[0135] Therefore, 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. These representations are then used as inputs for Bayesian linear regression. The computational cost of this method scales linearly with the dataset size and sequence length. This can be considered a convolutional kernel, as generating k-mer random features for each input sequence can be achieved by convolving with filters padded with random weights extracted from a distribution that fits the kernel we want to approximate. However, this approximate kernel still has limited practicality because convolving with thousands of filters is slow.

[0136] Subsequently, we introduce a second innovation that significantly accelerates the computation of random features for these kernels. To understand this implementation, it's important to first note that in deep learning libraries, convolution is typically implemented using matrix multiplication as shown below. Consider an input K×P matrix containing P sequence elements (e.g., P amino acids), each element (or amino acid) having K features, and assume the width of the convolutional filter is L. This input matrix can be used to fill an LK×(P-L+1) matrix, where each column i is a connection feature associated with subsequence i. Simultaneously, m convolutional filters can be stacked to form an m×LK matrix; multiplying this matrix by the LK×(P-L+1) input matrix is ​​equivalent to performing a one-dimensional (1d) convolution using m filters of length L.

[0137] In this implementation, the m×LK matrix of the convolution filter is replaced by SORF operations (Equation 4), and matrix multiplication is replaced by fast diagonal matrix multiplication and fast transform, similar to replacing the discrete Fourier transform with FFT. For simplicity, assume the required number of random features m is the same as LK. Thus, the size of the original random feature implementation of this kernel is O(N(P-L+1)m). 2 In contrast, the proposed modification has a scale of O(N(P-L+1)m log(m)).

[0138] Based on a fixed vector kernel used to measure the distance between k-mers, one can choose how to encode amino acids in two input sequences to measure the distance between k-mers in a specific way. For example, if amino acids are encoded one-hot, the distance between any two k-mers is twice the Hamming distance. Alternatively, learned representations (e.g., representations learned from a large language model) can be used to represent each amino acid. These possibilities will be described in the "Evaluating Convolutional Kernels for Sequence Data" section below, and the first type of kernel is called FHT-Conv-1d.

[0139] For the second valid convolutional kernel of the misaligned sequence

[0140] A second set of efficient convolutional kernels is introduced, which simulates a three-layer neural network, called Fast-Conv-1d. Its lemma is as follows (proof is given in Supporting Information S9):

[0141] Lemma 1.1: set up for On the positive definite kernel, and set From X to Let X be any mapping of X, where X is some non-empty set. Then, K(f(u), f(w)) is the positive definite kernel for all u, w∈X.

[0142] 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; the second layer is fully connected to a selected activation function. This is analogous to a fixed vector kernel (e.g., RBF or linear+RBF) whose input is a convolution of random features with ReLU activation and global max pooling. Figure 3 The fixed vector kernel is approximated using a stochastic feature procedure based on the Fast Hadamard Transform. According to Lemma 1.1, this combined kernel is a positive definite kernel. This kernel has the property that all convolutions can be executed exactly once before training, and the results are stored on disk. The aforementioned Fast Hadamard Transform-based procedure can be used to efficiently perform stochastic feature convolutions.

[0143] To understand how such kernels quantify sequence similarity, consider convolving a one-hot encoded sequence with a filter of width 9, whose weights are derived from a normal distribution. For each subsequence of length 9, the convolution kernel is large and positive if the filter contains large positive values ​​at the same position in the subsequence, and these positive values ​​are one-hot encoded. Since ReLU activation and max pooling are applied to the entire sequence, the final feature for each convolutional filter is the largest positive value encountered—in other words, the best "match" found in the sequence with that filter. Therefore, the kernel essentially "analyzes" each sequence input to quantify the best match of each filter in a set of random filters within the input, and then compares the analysis results of any two input sequences using a specified fixed-vector kernel (e.g., RBF or linear + RBF; for simplicity, RBF is used in all studies). As before, the learned representations, one-hot encoding, or other schemes can be used to encode the sequence as needed.

[0144] Evaluation of xGPR and convolution kernels for sequence data

[0145] xGPR, along with FHT-Conv-1d and Fast-Conv-1d kernels, were evaluated by fitting protein engineering benchmark datasets (TAPE and FLIP, respectively). These datasets involve providing pre-constructed training and test sets consisting of protein sequences, with the goal of predicting real-valued results for each sequence. A brief description of these datasets is provided below (for more details on these datasets and all other datasets used in this example, please see Supporting Information section S7 below).

[0146] The AAV dataset involves mutations across a 28-amino acid window of the adeno-associated virus (AAV) capsid, including insertions and deletions; some of these mutations are randomly generated, while others are “designed” (intentionally selected). The data is split in various ways, these different partitions are called “segments.” For example, they train the model using proteins with random mutations and test it on designed proteins (“mutant-to-designed segment”), or they train using designed proteins and test it on mutant proteins (“designed-to-mutant segment”), or they train using only proteins with seven mutations and test it on the remaining proteins (“seven-to-remains segment”), and so on. The goal is to predict the fitness of mutants, with higher fitness being better.

[0147] The thermal stability dataset uses mass spectrometry to measure the thermal stability (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" (including only a smaller subset of human data). The GB1 dataset measures the fitness of a bacterial protein G-binding domain mutant located at only four sites on the protein. This dataset is almost fully coverive, generating and testing most possible four-site mutants. The "three relative remainder" split uses mutants with only three mutations as the training set and all other mutations as the test set, while the "two relative remainder" split uses mutants with only two mutations as the training set.

[0148] The TAPE fluorescence dataset involves mutations in green fluorescent protein; the goal is to predict the log fluorescence of proteins in the test set based on the training set. This dataset is the same as the dataset used in the literature, discussed in "Comparison with Stochastic Variational Inference and Exact Gaussian Processes," although the TAPE benchmark divides the data into training and test sets differently. Finally, the TAPE stability dataset measures the stability of various proteins in protease sensitivity assays, with higher values ​​indicating better stability.

[0149] Three fits were performed using different random seeds, and the resulting error bars measured the variability among the different random seeds. Since the focus is on Gaussian process regression, only the benchmarks involved in the regression task were used.

[0150] As shown in Table 1, we compare 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 represent the input sequences using learned embeddings (input sequence representations generated by a protein sequence language model pre-trained by FAIR ESM-1v) and compare the results with those of fine-tuned language models reported in the literature. The best results of any fine-tuned large language models described in the literature for each benchmark are reported.

[0151] Table 1. Spearman's-r correlation coefficients of the predicted values ​​of convolutional neural networks trained on one-hot encoded data and xGPR relative to the true values.

[0152]

[0153]

[0154] a. In many cases, there are multiple “segments” for the same dataset; for example, for the AAV dataset, the “designed relative mutation” segment uses designed sequences as the training set and mutants as the test set, while the “mutant relative designed” segment does the opposite. b.. All convolutional kernels use a convolutional filter with a width of 9 as a simple default. The best performer for each benchmark is highlighted in bold. For all models, repeated predictions are made using three different seeds from a random number generator, and the mean and standard error of the mean are then reported. For the scores shown here, higher scores are better. c. For thermally stable datasets, some sequences are very long (thousands of amino acids). The Fast-Conv-1d kernel is much faster for these types of sequences because it only needs to perform one convolution during hyperparameter fine-tuning, and only this kernel is considered here. d. Results for all one-hot encoded CNNs are from Dallago et al. (Proc. Neural Inf. Process. Syst. Track DatasetBenchmarks, 2021) or Rao et al. (ArXiv190608230 Cs Q-Bio Stat 2019). e For the GB1 dataset only, the input sequence length is 4 amino acids, so only the RBF kernel (equivalent to a convolution kernel with a width of 4) was used for this dataset.

[0155] Table 2. Spearman rank correlation coefficients of xGPR predictions relative to true values ​​trained on fine-tuned large language models and ESM-1v embeddings.

[0156]

[0157] a. In many cases, there are multiple "segments" for the same dataset; for example, for the AAV dataset, the "designed relative mutation" segment uses designed sequences as the training set and mutants as the test set, while the "mutant relative designed" segment does the opposite. b. All convolutional kernels use a convolutional filter with a width of 9 as a simple default. The best performer for each benchmark is highlighted in bold. For all models, repeated predictions are made using three different seeds from a random number generator, and then the mean and the standard error of the mean are reported. For the scores shown here, higher scores are better. D.. All results for the pre-trained language models and one-hot encoded CNNs reported are from Dallago et al. or Rao et al. These authors report results for multiple pre-trained language models; they also report their best results obtained in each case. e For the GB1 dataset only, the input sequence length is 4 amino acids. Therefore, for this dataset, the FastConv-1d kernel was not used and the sequence was not averaged.

[0158] For the learned embeddings, at least two methods can be used as input to the Gaussian process, and both methods were evaluated. First, the characterization of each amino acid in the sequence can be averaged, and then an RBF kernel can be used. Alternatively, the learned embeddings of size D can be arranged into an M×D array for each sequence of length M, and the entire sequence can be convolved using either the FHTConv-1d or FastConv-1d kernel. Using the FastConv-1d kernel makes training easier because all convolution operations can be performed once at the start of training, eliminating the need to save the learned embeddings to disk (which can be very large). Therefore, FastConv-1d was used here.

[0159] Three interesting findings are noted. First, when trained with one-hot encoded inputs containing a moderate amount of random features (<10,000), xGPR outperforms CNNs on 8 out of 11 benchmarks and matches them on another. When trained with learned embeddings, the approximate GP outperforms fine-tuned language models on 6 benchmarks and matches them on another. Notably, the pre-trained language model requires up to 50 GPU days of training time to fit, while the FastConv-1d kernel fits the datasets shown here (including hyperparameter fine-tuning) in less than 6 minutes per fit, with a maximum memory footprint of 1.7GB when using mini-batch sizes of 2000 data points (the FastConv-1d kernel is slower, but always takes less than 15 minutes including hyperparameter fine-tuning). Using learned embeddings generally improves GP performance, especially on diverse datasets with sequences from many different protein families (e.g., thermal stability), but this is not always the case; in several benchmarks, GP performs better with one-hot encodings.

[0160] Secondly, using more random features in all cases slightly improves the performance of the Gaussian process; it increases the accuracy of the kernel approximation. This feature provides users with a simple way to improve performance as needed. Thirdly, the FHT-Conv-1d kernel and Fast-Conv-1d are not always superior to each other; instead, they perform slightly differently on certain datasets.

[0161] Uncertainty calibration

[0162] The primary motivation for using Gaussian processes (GP) is to quantify uncertainty. Various methods for uncertainty estimation in deep learning have been reported in the literature. xGPR provides more accurate uncertainty estimates than these alternatives.

[0163] For uncertainty estimates to be useful, uncertainty needs to increase with distance from the training set—this is true for generalizations (GPs) that build stationary kernels. Furthermore, a well-calibrated uncertainty estimate is needed; in other words, the uncertainty reflects how often the predicted outcome actually occurs. For example, if a weather forecast says there is a 10% probability of rain, then it should rain on 10 out of the 100 days that the forecast was made. Similarly, for regression models, approximately 10% of the test data points are expected to fall within their predicted 10% confidence interval (or the “confidence interval” in Bayesian inference), approximately 20% within the 20% confidence interval, 90% within the 90% confidence interval, and so on.

[0164] Calibration can be quantified as the area under the calibration error curve (AUCE). For 100 values ​​(0.01, 0.02, etc.) in the range of 0 to 1, determine what percentage of the predicted values ​​in the test set fall within the corresponding confidence interval (1%, 2%, etc.). Subtract the observed percentage within the interval from the expected percentage and take the absolute value; for example, if the confidence interval is 10% and the observed percentage within that interval is 90%, the difference is |0.9-1|. Finally, calculate the area under the curve; the larger the area, the worse the numerical calibration of the uncertainty provided by the model. A perfectly calibrated model scores 0.

[0165] Five different uncertainty estimation methods (dropout uncertainty, last-layer stochastic variational inference, model ensemble, deep evidence regression, and mean-variance estimation) are described and used in conjunction with deep learning models for the five benchmarks described in this example. One-hot embeddings and ESM embeddings are both used as inputs to the models, and AUCEs are computed. As shown in Table 3, the AUCEs of their uncertainty estimation methods are compared with the AUCEs of xGPR using the same encoding on the same dataset. To allow the deep learning models to benefit from skepticism, the best AUCEs reported for each uncertainty estimation technique are provided for each dataset, and the xGPR model that demonstrates the best performance on that dataset (evaluated by the correlation between the fitness of predictions and measurements) from Tables 1 and 2 is used.

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

[0167]

[0168] Note: In all cases, to benefit the deep learning model, the best performance of uncertainty estimation techniques was used, regardless of whether one-hot embeddings or ESM embeddings were used as inputs. a.. In many cases, there are multiple “segments” for the same dataset; for example, for the AAV dataset, the “designed relative mutation” segment uses designed sequences as the training set and mutants as the test set, while the “mutant relative designed” segment does the opposite.

[0169] Regarding uncertainty calibration, xGPR significantly outperformed the best deep learning uncertainty estimation methods across all five benchmarks, achieving an AUCE three to ten times better. xGPR's AUCE was less than 0.15 across all scenarios, while deep learning methods typically (and some always) achieved AUCE scores of 0.3 or higher.

[0170] Protein engineering using uncertainty

[0171] Uncertainty can be used to determine the reliability of predictions generated by the model; if the uncertainty is high compared to the uncertainty in the training set, the predictions may be considered less reliable. Alternatively, the uncertainty estimates generated by xGPR can be used as part of a Bayesian optimization strategy. In this framework, we first fit the model to the available data. A new set of sequences is selected, and experimental evaluation is performed using a collection function that considers both model predictions and uncertainty, thus striking a balance between exploration and exploitation. For example, the upper confidence limit (UCB) adds the upper half of the confidence interval to the model predictions at each test point. After experimental evaluation at the selected points, the model is refitted, and this process is repeated until the desired protein function is achieved. For an overview of other common collection functions, see Frazier 2018.

[0172] In this study, the GB1 dataset, discussed under “Evaluating Convolutional Kernels for Sequence Data,” was used, which introduces mutations at four sites in the protein G-binding region. This dataset is convenient for such experiments as nearly 150,000 of the 160,000 possible mutants have been evaluated. The fitness of all mutants in the dataset was normalized to be between 0 (potentially worst) and 1 (potentially best). Note that >99.95% of the sequences had a fitness less than 0.6; only one sequence had a fitness of 1. 384 variants were randomly selected and an xGPR model was trained using an RBF kernel. UCB acquisition function scores were then generated for all remaining mutants and used to select a batch of 96 for “experimental” evaluation. These 96 were then added to the training set, the model was refitted, and this process was repeated for up to five iterations. In each iteration, the best fitness achieved was tracked. Finally, this “experiment” / “process” was repeated 50 times with different random seeds.

[0173] By the fifth iteration, all replicates achieved a fitness of 0.6 or higher. Forty out of fifty replicates found one of the four best sequences, with 30 replicates achieving a score of 1, meaning they found the best sequence in the dataset—despite testing fewer than 1000 sequences. Half of the replicates found the best sequence (fitness 1) within four or fewer iterations. These results suggest that, when combined with high-throughput experimental systems, this xGPR-based method could potentially "discover" proteins with desirable properties without human intervention. It is efficient and number-invariant for small molecules and other graphical applications.

[0174] In the previous section, “A New Set of Efficient Kernels for Unaligned Sequences,” we introduced two sets of linear scaling kernels for sequences. These kernels can be easily scaled to graphs if each node in the graph is encoded as a feature vector containing some information about that node, its edges, and neighboring nodes. For the Fast-Conv-1d kernel, random feature convolutions and max pooling operations can be applied to all feature vectors associated with the graph. For the FHT-Conv-1d kernel, nodes and their associated features are used only as inputs, instead of k-mer, and convolutions of width 1 are performed between nodes in the graph: this is called the Graph-Conv-1d kernel.

[0175] Several variants of the Graph-Conv-1d kernel are implemented in xGPR. This paper focuses on the simplest variant, Graph-RBF, which uses an RBF kernel to perform approximate pairwise comparisons of nodes in two graphs. To evaluate the performance of Graph-RBF, we consider predictions of small molecule energy and thermodynamic properties. We benchmarked using the QM9 dataset, which contains 134,000 small molecules whose geometries have been optimized using density functional theory (DFT). The kernel described in this paper appears particularly well-suited for energy / thermodynamic properties based on prior knowledge, as these properties can, in principle, be decomposed into the sum of contributions from each atom.

[0176] Two different modeling methods were used. In method (1), each atom I is one-hot encoded and then concatenated with the one-hot encodings of all other atoms j in the molecule, which are divided by the sixth power of the distance from atom j to atom i; the one-hot encodings are sorted by distance. Due to the effect of the sixth power weighting, features related to more distant atoms become negligible beyond distances of four or five angstroms. This representation is then used as input to the Graph-RBF kernel. This simple representation is easy to construct but has little information (it discards all information about bond angles), and therefore is unlikely to achieve good performance.

[0177] Alternatively, in method (2), each atom is represented using SOAP features. In this scheme, the atomic density around each atom is expanded in an orthogonal basis set, and the power spectrum using these coefficients is used as the representation. A distance-weighted scheme is used, with only minor modifications to its selected hyperparameters. For more details on the hyperparameters used and their fine-tuning, please refer to the supporting information section S10.

[0178] It is noteworthy that method (1), despite using non-informative features, achieves a chemical precision of approximately 1 kcal / mol. Method (2) is significantly better; the results are shown in Table 4. For internal energy, enthalpy, and Gibbs free energy, xGPR achieves a MAE < 0.2 kcal / mol, while achieving a new level of technology for one of the six properties. Note that 1 kcal / mol is generally considered "chemical precision," i.e., the precision required to correctly predict the behavior of a chemical system; therefore, xGPR and several other models are far above the level of precision actually required. In fact, density functional theory (DFT) calculations for which we attempted to predict results are typically only accurate to 2–3 kcal / mol. Therefore, the 0.01–0.05 kcal / mol difference between models used to predict energies calculated by DFT may be the same in practical applications.

[0179] Table 4. Comparison with state-of-the-art deep learning models for predicting the energy and thermodynamic properties of small molecules, QM9 dataset

[0180]

[0181] Retrieve similar data points from the training set

[0182] For Gaussian processes, SHAP or LIME techniques can be used to measure feature importance, just like any other machine learning model. Unlike other machine learning models, for small datasets, a kernel matrix can be additionally constructed and used to precisely determine the contribution of each data point in the training dataset to a given prediction. This decomposition is not possible in deep learning models.

[0183] However, if the dataset exceeds 5,000 data points, directly constructing the kernel matrix is ​​too costly. Furthermore, for large datasets, assigning precise contributions to each point in the training set may provide more information than is actually needed or useful. Often, practitioners may simply want to know which data points in the training set are the most similar, i.e., those data points that contribute the most to the predictions. This can be determined as follows. Let z... * Let z be a random feature representation of the test data points. i Let x be the random feature representation of the training data points; then the kernel function x corresponding to the test data points and training data points is... * Approximate x i For z * T z i (See Supporting Information Section S1)

[0184] In practice, retrieving the most similar training point by calculating the dot product of the test data point representation and all training data each time can be computationally expensive. This problem can be solved simply by clustering the training data. The same kernel used to fit the Gaussian process regression model can be used to cluster it using kernel K-means clustering, which approximates the random features. It is worth noting that, just as random features approximate a Gaussian process in a linear regression manner, clustering the training set represented by random features using K-means is equivalent to applying kernel K-means to the original dataset. Once the dataset is clustered, to retrieve the most similar training point, it is no longer necessary to measure the similarity between the test point and each training data point; it is only necessary to calculate its similarity to the data points within the nearest cluster.

[0185] The tools for performing these functions are included in the xGPR library. This procedure is described for the QM9 dataset; see Support Information section S11 for more details. First, the QM9 dataset is fitted using a Graph-RBF kernel with a target energy of 298K, and atoms are encoded using the one-hot encoding method (1). Then, k-means clustering is performed on the training set. The final elbow plot used to select the number of clusters is shown in the Support Information and Figures 9A-9B The description is as follows. The elbow diagram shows that the training dataset can be divided into five clusters. Then, five molecules are randomly selected from the test set. For each molecule, the 50 most similar molecules from the training set are found using the procedure described above, where similarity is measured by a kernel function. The structure of the query test molecule and the most similar training molecule (quantified by a kernel function using one-hot encoded input) is described in the supporting information section S12.

[0186] Support information

[0187] Supporting information includes: (S1): More detailed background on random feature approximation; (S2): More detailed background on orthogonal random features; (S3): Performance metrics of the Fast Hadamard Transform implementation; (S4): Preconditioning algorithm and some benchmarks for modifications to this algorithm; (S5): Algorithm for approximating random features – logarithmic determinant of approximation Gaussian process; (S6): Discussion of hyperparameter fine-tuning strategies; (S7): Details of the benchmark dataset used in this embodiment; (S8): Benchmarks for random feature approximation and random variational inference; (S9): Proof of Lemma 1.1; (S10): Details of the QM9 benchmark; (S11): Clustering and principal components using random features; and (S12): Retrieving similar molecules from the training set using random features.

[0188] S1 random features

[0189] In our initial study, we considered only positive definite kernels, and limited them to stationary kernels, i.e., those satisfying K(x) i x j)=f(x i -x j The nucleus of ) is sought, and a change in φ(x) is sought to make the nucleus (x) i x j )≈φ(x i )φ(x j This satisfies the condition. According to Bochner's theorem, any considered positive definite kernel is a Fourier transform of the corresponding positive measure, i.e.:

[0190]

[0191] This expectation can be approximated by Monte Carlo sampling, i.e.:

[0192]

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

[0194] Although this process was initially derived in the context of stationary kernels, it can be applied to Monte Carlo sampling methods with certain non-stationary kernels, most importantly a class of kernels called "neural network kernels," and in recent literature also referred to as "pointwise nonlinear Gaussian kernels." Consider a neural network with one hidden layer and one output layer:

[0195]

[0196] Where x is the input vector, and W is a vector of dimension D×N. H The learning weight matrix, where D is the dimension of the input, and N is the weight matrix. H 'b' is the number of hidden layer units, and 'b' is a unit of length N. H The offset or bias vector, where h is a certain activation function, and v i This refers to the learning weights. Assume all weights in W are independently and identically distributed (i.e., i.i.d.) following N(0, σ). w 2 ), weight v i They are independently and identically distributed, following N(0, σ). v 2 The deviation terms are also independently and identically distributed, following N(0, σ). b 2 If we add an element to each x, that element is 1 for all data points, then it can be rewritten more simply as follows:

[0197]

[0198] Where each row of W is independently and identically distributed, following N(0, Σ). W), where when Σ W When it is a diagonal covariance matrix, its first diagonal element is σ. b 2 And all remaining diagonal elements are σ. w 2 .

[0199] Let S denote the set of all weights; therefore

[0200]

[0201] Furthermore, Monte Carlo sampling can be used to approximate the expectation term on the right side, resulting in a formulation very similar to the stationary kernel. This kernel is positive definite only with respect to certain choices of activation function h. One such choice is the ReLU activation function, which produces a first-order arccosine kernel, although the error function is also considered a useful activation function. Interestingly, the first-order arccosine kernel formulated in the literature omits the bias term. In our library implementation, we found that the arccosine kernel performs very poorly without the bias term and requires its inclusion to achieve acceptable performance.

[0202] For stationary kernels and neural network kernels, the stochastic Fourier feature method can be implemented as follows. If X is an N×D matrix with N data points and D input features, then a D×M matrix is ​​filled with M samples from p(ω), where p(ω) is the Fourier transform of the kernel for stationary kernels and the normal distribution for neural network kernels. This sample matrix is ​​multiplied by X, and for neural network kernels, an additional bias vector is added to the result. An appropriate activation function (e.g., sine and cosine for stationary kernels) is then applied and multiplied by a constant to obtain an N×M matrix Z (or, N×2M for stationary kernels). The approximate kernel matrix is ​​then obtained by ZZ. T Provided.

[0203] Given a training dataset X and its corresponding true value y, a zero-mean Gaussian process, its value at a new data point x * Predicted value f * The distribution is as follows:

[0204] p(f * |X,y)~N(μ * , Σ * )

[0205] μ * =K(x) * ,X)[K(X,X)+λ 2 I] -1 y, σ *

[0206] =K(x) * x * )-K(x* ,X)[K(X,X)+λ 2 I] -1 K(X,x * ).

[0207] Since the kernel approximation is ZZ T Therefore, these formulas can be rewritten using the matrix inversion lemma and the Woodbury matrix identity. For example, the mean becomes:

[0208]

[0209] Where Z is the training set after the above transformation, I is the identity matrix, λ is the hyperparameter ensuring the kernel matrix is ​​a positive definite matrix, and z is the new data point after the above transformation. If we consider...

[0210] [Z T Z+λ 2 I] -1 Z T y

[0211] As weights, the predicted mean corresponds to a linear regression in the feature space corresponding to the random feature map. It's easy to see the correspondence between a neural network with a single hidden layer and an output layer and the kernel described above approximated using this formula, although some key differences exist. The approximate Gaussian process "learns" by changing the distribution of the hidden layer weights (rather than a single value) and can be fitted by maximizing the marginal log-likelihood rather than the likelihood; furthermore, both the marginal likelihood and the posterior prediction distribution are provided in closed-form. Kernels similar to those of neural networks with multiple hidden layers can also be constructed. In the main text, the Fast-Conv-1d kernel is introduced, which corresponds to a three-layer neural network with convolutional layers, fully connected layers, and an output layer, all filled with random weights (see...). Figure 3 ).

[0212] Note that since the error in the approximation decays exponentially with the number of random features, a large number is usually not required to obtain a good approximation; in fact, the gains decrease sharply after increasing the number of random features to around 16,000, while 8,192 is sufficient for many tasks. For the predicted mean μ... * Since this is the most important quantity, practitioners might want to use a large number of random features in exchange for a small improvement in accuracy, but for quantizing variance, 512-1,024 is usually sufficient. Therefore, xGPR allows users to use fewer random features than the predicted mean to quantize variance (in order to reduce the computational cost of predicting variance).

[0213] S2. Orthogonal random characteristics

[0214] Consider a case where the number of random features D is the same as the dimension d of the input, and the input vector x is multiplied by an independent and identically distributed matrix of random features M drawn from a normal distribution.

[0215] If the rows of M may be orthogonal, we can first perform a QR decomposition on M, and then multiply the resulting Q by a diagonal matrix S, whose diagonal elements are taken from a χ-distribution with d degrees of freedom. The input vector is then transformed by calculating the product kSQx. If the number of random features is greater than d, this process can be repeated m times, where m is the result of dividing D by d and rounding up, then any features exceeding D are discarded. If D is less than d, d random features can be generated, and any redundant features are discarded.

[0216] For certain kernels, particularly the quadratic exponential kernel, it has been shown that only a small number of orthogonal features are needed to accurately approximate the kernel. Therefore, using orthogonal random features can achieve a more accurate approximation with the same number of random features. In addition to this advantage, the orthogonal feature method can be further modified to significantly reduce computational costs. Results show that almost the same result can be obtained by replacing matrix Q with the following product:

[0217] HD1HD2HD3

[0218] Where H is the normalized Hadamard matrix, and each D is an independent diagonal matrix whose elements are derived from the Lamach distribution. The memory footprint of the fitted model is significantly reduced because only the diagonal elements of the D matrices need to be stored. More importantly, similar to replacing matrix multiplication with Fast Fourier Transform in Discrete Fourier Transform, replacing Hadamard matrix multiplication with Fast Hadamard Transform reduces the cost of the process to O(Dlog D), which is a significant advantage for larger M values. Algorithm S1 illustrates the general steps of the input transformation for the RBF kernel (similar for other kernels).

[0219] Please note that some have suggested replacing the S matrix used here with a constant multiplied by it. However, this fails on one-dimensional data because the generated kernel cannot provide an acceptable fit for simple basic validation problems, and offers no advantage on high-dimensional data. The version mentioned above does not have this limitation.

[0220] Algorithm S1: Orthogonal random eigentransformation of RBF kernel

[0221] Kernel object initialization

[0222] Input: Expected dimension d of the input; expected number of sampling frequencies M

[0223] 1

[0224] 2. If D <M:

[0225] aP = ceiling(M / D)

[0226] 3. Otherwise:

[0227] aP = 1

[0228] 4. Create and store a diagonal matrix S of size M×M, whose elements are taken from a χ-distribution with D degrees of freedom.

[0229] 5. For i from 1…P:

[0230] a. Create and store a diagonal matrix D of size D. ai D bi D ci Its elements are taken from the Lamach distribution.

[0231] kernel used for transforming input data

[0232] Input: N×d matrix x i It is a block of input data; kernel hyperparameters β, σ

[0233] Output: Kernel transformation matrix Z i

[0234] 1. For i from 1…P:

[0235] a. Create an empty N×D array Y i

[0236] b. Copy x i To Y i The first d columns; set the remaining elements to 0.

[0237] c. For k in [a,b,c]:

[0238] iY i =H(Y) i D ki ), where H is the normalized in-situ fast Hadamard transform.

[0239] 2. Y1...Y P Concatenate into an N×D*P array Y

[0240] 3. If D*P>M:

[0241] a. Discard the M…D*Pth columns of Y.

[0242] 4. Y = σYS

[0243] 5. Y1 = cos(Y), Y2 = sin(Y), where sin and cos are element-wise operations.

[0244] 6. Concatenate Y1 and Y2 to form an N×2M matrix Z.

[0245] 7.

[0246] 8. Return to ZS3: Performance of Fast Hadamard Transform

[0247] This paper presents a C / CUDA implementation of the Fast Hadamard Transform in Python for both CPUs and GPUs. We compared the speed of matrix multiplication with random feature matrices using implementations in the GPU's Cupy library (v 10.4.0) and NumPy library (v 1.22.3) (results are shown in Tables 6 and 7). In short, for a moderate number of random features, this implementation is faster than matrix multiplication, and its performance advantage increases as expected with the number of random features, enabling rapid scaling to large numbers of random features. For comparison with other fast transforms, we also compared the CPU Discrete Cosine Transform (DCT) implemented in the SciPy library (v 1.8.0) (results are shown in Table 1).

[0248] Tables 5 and 6 show the results in milliseconds. In Table 5, we evaluate the time required to generate random features when the number of columns in the input matrix is ​​the same as the number of features—an unusual case that significantly demonstrates the speedup of the procedure based on the Fast Hadamard Transform. In Table 6, a more realistic scenario is considered: the number of input features is fixed at 1024, and an increasing number of random features must be generated.

[0249] We used an input matrix with 2000 rows because this block size is common when processing real-world datasets. For the GPU implementation, we used 32-bit (single-precision) floating-point, as single-precision matrix multiplication offers significant performance improvements on most GPUs; while for the CPU tests, we used 64-bit (double-precision) floating-point. However, generally, unless the user chooses otherwise, xGPR will use single-precision for both CPU and GPU random feature generation.

[0250] Table 5: Comparison of the speed of SORF operations implemented in this study (highlighted in the table) with the original random characteristics and the general SciPy fast transform.

[0251]

[0252]

[0253] The results shown here have been rounded to reflect precision, include CPU and GPU times (where applicable), and are calculated as the average time for each operation after 100 repetitions of a specified operation. All studies were conducted on an Intel i5-7500 CPU and a GeForce GTX 1070 GPU. The methods used in this study are shown in bold.

[0254] Table 6: Comparison of the speed of SORF operations implemented in this study (highlighted in the table) with the original random characteristics and the general SciPy fast transform.

[0255]

[0256] The results shown here have been rounded to reflect precision, include CPU and GPU times (where applicable), and are calculated as the average time for each operation after 100 repetitions of a specified operation. All studies were conducted on an Intel i5-7500 CPU and a GeForce GTX 1070 GPU. The methods used in this study are shown in bold.

[0257] S4: Preconditions for conjugate gradients and stochastic gradient descent

[0258] Consider the well-known conjugate gradient algorithm, which solves for the system Aw = b of w. It can be proven that the upper bound of the residual in the nth iteration is:

[0259]

[0260] Where κ is the condition number of A, and ||e0|| is the error of the initial guess. If κ is very large, the algorithm will converge very slowly. If the inverse of the preconditioner is M... -1 If the inverse of A is very close, then the condition number of A can be reduced and the convergence speed can be accelerated. Crucially, there must be a matrix Z that does not require explicit construction. T Z can be used to construct preconditioners because it requires expensive matrix multiplication, and for a large number of random features, the matrix is ​​too large to be stored in memory.

[0261] The kernel matrix of a Gaussian process is often approximated using the Nystrom method, which provides the following low-rank approximation:

[0262] A≈A m，q A q -1 A m，q T

[0263] Where m is the number of rows in A, and q is a randomly selected subset of columns in A. The Nystrom approximation has low construction cost and has successfully solved many problems in practice. Although sampling from a uniform distribution has traditionally been the most popular method, its performance is sensitive to the sampling method used to select columns and rows in A.

[0264] We found that the randomized Nystrom approximation method performs better. In this scheme, for the system:

[0265] (Z T Z+λ 2 I) -1 Z T y = w;

[0266] Z T The low-rank approximation of Z = A is given by the following formula:

[0267]

[0268] For Ω, an independent and identically distributed matrix plotted from a zero-mean, unit-variance normal distribution is used; Ω is of size M×L, where M is the number of random features. This method provides a rank-L approximation of A and offers a numerically stable generation method (implemented in the library) as shown in algorithm S2. Importantly, this method can be constructed as illustrated without explicitly constructing matrix Z. T Instead of using Z, a single loop is used on the dataset to load one data chunk at a time into memory.

[0269] Algorithm S2: Constructing randomized Nystrom preconditioners

[0270] Input: Dataset X, divided into blocks x1...x n The format is stored on disk, and the required preconditioner size L, GP kernel, and number of random features M are specified.

[0271] Output: Factorization preconditioner of form UΛU T , where U is orthogonal and Λ is diagonal.

[0272] 6. Generate an M×L matrix Ω by sampling N(0,1) from it.

[0273] 7. Ω=qr(Ω)QR decomposition

[0274] 8. Generate an M×L matrix Q and initialize it to all zeros.

[0275] 9. For i from 1 to n, load x i

[0276] a. Generate z using algorithm S1 i

[0277] b.Q+=z i T (z i Ω)

[0278] 10.

[0279] 11.Q v =Q+vΩ is used for numerical stability

[0280] 12. C = cholesky (Ω) T Q v )

[0281] 13. CB = Q v Trigonometric solution B

[0282] 14.U, Σ, ~=svd(B) SVD of B

[0283] 15. Λ=max(0,Σ 2 -vI)

[0284] 16. Return to U, Λ

[0285] Here, a novel modification is introduced. The M×L matrix Ω is replaced with a subsampled randomized Hadamard transform, defined as follows:

[0286]

[0287] Where D is a diagonal matrix whose elements are derived from the Lamacher distribution; H is a normalized Hadamard matrix; and S is an L-row subset extracted from the identity matrix with equal probability. In practical applications, the Fast Hadamard Transform is used instead of the Hadamard matrix. For step 9b, it is replaced with:

[0288] Q+=SRHT(Z T )Z

[0289] In step 12, for the product Ω T Q v , replacing sSRHT(Q v Finally, to ensure numerical stability, the shift is replaced, and then the above steps (steps 11 and 15) are subtracted, where SRHT(Q) is not performed in step 12. v Instead of using the Cholesky decomposition of B, we use SVD to find the self-adjoint square root of the matrix, and then use it to solve for B in the next step.

[0290] While previous literature has proposed using SRHT to construct low-rank matrix approximations, it has not been used to construct preconditioners for CG / stochastic gradient descent as suggested in this paper. Since matrix multiplication with O(NML) scaling is replaced by fast transformation operations, it significantly reduces the cost of preconditioner construction; this is particularly advantageous on CPUs. Empirical evidence shows that even for datasets with only 40,000 data points, the SRHT-based preconditioner construction procedure is significantly faster, and this advantage becomes more pronounced as the dataset size and the number of random features increase.

[0291] In various situations encountered, experience has shown that the performance of preconditioners constructed using the SRHT method suggested here is as good as or better than that of preconditioners constructed using the unmodified algorithm S2. Figure 4 As shown, preconditioners were constructed using either SRHT or the unmodified algorithm S2, and the number of iterations required for conjugate gradient convergence was compared to a predetermined threshold using various L values. In all cases, the SRHT-constructed preconditioner achieved the same or greater acceleration, which was indeed the case throughout the study. Currently, both the modified and unmodified preconditioner construction algorithms are provided as options for users to select in xGPR, with the SRHT-based algorithm used by default due to its superior speed and empirical performance.

[0292] Furthermore, the preconditioner can be further modified to improve its performance. In step 9 of algorithm S2, a product Z is generated. T (ZΩ), where Z is a random feature matrix of shape N×M in the training data, containing M random features and N data points. Matrix Q can be obtained from the QR decomposition of this product and used to perform a second traversal of the dataset, thereby generating a second product Z. T (ZQ). This product can then replace matrix Q in step 9 in algorithm S2 (this process can even be repeated a third time, although the improvement to the preconditioner performance is negligible). For example... Figure 5 As shown, the resulting preconditioner reduces the number of iterations required for conjugate gradient convergence in all studies, by 20-25% in most cases. However, it is more expensive to build because it requires two traversals of the dataset. Whether this is advantageous 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 choose this option is provided in the documentation.

[0293] Regardless of which option is used to construct it, the inverse of the preconditioning matrix takes the following form:

[0294] M -1 =(β) L +λ 2)U(Λ+λ 2 I) -1 U T +(I-UU T )

[0295] Where β L λ is the smallest eigenvalue of the preconditioner, and λ is the hyperparameter of the GP kernel. Matrix M -1 It is not explicitly constructed, but since all the necessary values ​​have been generated, the matrix-vector product M is... -1 v is available on demand, which is exactly what is needed for preconditioned conjugate gradients or preconditioned stochastic gradient descent.

[0296] The larger L is, the better the approximation effect; if β1≥β2≥β3… are eigenvalues ​​of A, then the faster the eigenvalue decays, the better the approximation effect; and for β… j≥L The smaller the eigenvalue, the better. Let UΛU T A is of rank L nystrom The eigenvalue decomposition, let the preconditioner P be... And let β L It is A nystrom The smallest eigenvalue of . It can be proven that if κ2 is the condition number of the preconditioned system, then the probability is:

[0297]

[0298] It is ≥1-δ, where δ is a certain fault tolerance capability.

[0299] How to choose L? For a given L, it can be estimated by constructing a preconditioner. The value of L is then doubled until... This continues until the number of iterations required to fit the model is deemed acceptable. For example... Figure 1A As shown, for a given expected tolerance, the logarithm of the number of iterations typically exhibits a roughly linear relationship with the expected tolerance. Therefore, for a given The number of iterations for fitting can be estimated. If the final number of iterations is considered too high, L can be adjusted and the preconditioner reconstructed.

[0300] Down Figure 1B This illustrates the relationship between the value of L and the number of iterations required for fitting using CG.

[0301] Preconditioners can quickly fit models, but they cannot solve the problem of hyperparameter fine-tuning. Calculating the marginal likelihood of a Gaussian process requires calculating the logarithmic determinant of an M×M matrix of M random features, which becomes computationally very expensive if M is large. Therefore, we develop a method for estimating the marginal likelihood of a large number of random features. S5: Logarithmic Determinant Approximation and Conjugate Gradient

[0302] The hyperparameters of a Gaussian process can be fine-tuned using performance on a validation set (just like any other model), a method implemented as an option in xGPR, but they can also be fine-tuned by maximizing the marginal likelihood. This process is more robust to overfitting and does not require a validation set, so this method should be used whenever possible. Nevertheless, fine-tuning using marginal likelihood does present a unique challenge. Recall that for a Gaussian process, the expression for the negative marginal log-likelihood is as follows:

[0303] -

[0304] Where K y =K(X, X)+λ 2 I. In a randomized characteristic scheme, the approximate value Z T Z is replaced by K(X, X) and rearranged to obtain:

[0305]

[0306] As shown in S4 above, it can efficiently calculate (Z) T Z+λ 2 I) -1 Z T y, without explicitly forming the matrix |Z using preconditional conjugate gradients. T Z+λ 2 I|. However, if the number of random features is large, then log|Z is calculated. T Z+λ 2 I| is more challenging. In Part S6, efficient methods were developed to fine-tune hyperparameters with a small amount of randomness (e.g., <4000, which often occurs in practice). However, a method is first introduced that uses preconditioned conjugate gradients to approximate log|Z. T Z+λ 2 I|, no need to explicitly form Z T Z+λ 2 The I matrix. Even with a large number of random features, the same tools can be used for fitting, allowing for efficient fine-tuning of hyperparameters.

[0307] The literature outlines several methods for approximating the logarithmic determinant of a large matrix A, including the stochastic Lanzos quadrature method, the stochastic trace estimation method combined with Chebyshev polynomials, and the subspace iteration method. The subspace iteration method is only applicable when a small subset of the eigenvalues ​​of A is dominant, i.e., there exists an eigenvalue λ. i >>λ i+1 , λ i+2 ...λ N The larger feature gaps exist, but this does not always hold true during fine-tuning.

[0308] The stochastic Lanzos quadrature (SLQ) method is particularly attractive in this regard because, for a large number of stochastic features, the model can be fitted using preconditioned conjugate gradients, and the tridiagonal matrix required by the SLQ method can be automatically generated during the conjugate gradient optimization process with negligible additional cost. It has been shown that using incomplete Cholesky decomposition to construct preconditioners allows for accurate estimation of the logarithmic determinant of an exact (non-approximate) Gaussian process.

[0309] Algorithm S5 below demonstrates how to combine this method with stochastic feature approximation and stochastic Nystrom approximation to accurately estimate logarithmic determinant terms. Simultaneously calculate performance loss (The remaining terms in the marginal likelihood are either constants or do not require data retrieval, making them easy to compute.) The final algorithm exhibits linear scaling with the dataset size, subquadratic scaling with the number of random features, and does not require loading multiple data blocks into memory at once.

[0310] First, rewrite the logarithmic determinant and use Hutchison's trace estimator. 21 :

[0311]

[0312] Where, p i It is a probe vector (usually) extracted from a La Mach or Gaussian distribution.

[0313] Next, when using the preconditioned conjugate gradient method and the preconditioner P, we are actually solving the preconditioned system log|P. -1 / 2 Z T Z+λ 2 IP -1 / 2 |. It can be rewritten using the same method as the exact (non-approximate) GP:

[0314] log|Z T Z+λ 2 I|=log|P -1 / 2 (Z T Z+λ 2 I)P -1 / 2 |+log|P|

[0315]

[0316] Therefore, it is necessary to apply the preconditioner to the probe vectors. This is easily done if the probe vectors follow a Gaussian distribution; they can be derived solely from the normal distribution of the covariance matrix, which is the preconditioner. These preconditioned probe vectors are now called c. i .

[0317] The term is evaluated as follows:

[0318]

[0319] This is implemented using an orthogonal-based procedure but with a preconditional conjugate gradient algorithm; by storing the α and β coefficients generated in each conjugate gradient iteration, a tridiagonal matrix constructed using the Lanczos algorithm is generated. At this point, only the matrix containing Z needs to be processed. T Run conjugate gradients on multiple target vectors of the y vector to compute the performance loss (Z). T Z+λ 2 I) -1 Z T y; and also includes preconditioning probe vectors to generate n v Tridiagonal matrix A simplified version of the complete program is outlined in Algorithm S5 below.

[0320] Algorithm S5: Used for simultaneous estimation (Z) T Z+λ 2 I) -1 Z T y and log|Z T Z+λ 2 Preconditional conjugate gradient method of I|

[0321] Input: Dataset X, where N data points are divided into blocks x1...x n and y1...y n The format is stored on disk, and the initialized kernel has fixed hyperparameters, which can be the input point x. i Generate a dimensional random feature representation k(x) of length M. i The preconditioner P is constructed using algorithm S2 or a suitable variant; the number of probe vectors n v Convergence threshold tol

[0322] 1. Create an empty array Z of shape M. T y

[0323] 2. For i from 1, 2…n:

[0324] aZ T y+=k(x i ) T y i

[0325] 3. For i from 1, 2…n v :

[0326] a. Extract the probe vector c from N(0,P). i

[0327] 4. Create a shape with (M, n) vEmpty array r of (+1) j

[0328] 5. Use to fill r j for the first n v columns. Copy Z T y into the remaining columns.

[0329] 6. Let where ||r jk || is the Euclidean norm of the k-th column in r, so L1 is a row vector of length n j +1. v +1.

[0330] 7. Create an initial weight matrix w of shape (M, n v +1); set all elements of the matrix to zero. j ; Set all elements of the matrix to zero.

[0331] 8. Create an empty array B of shape (M, n v +1). Let 1 M be a vector of 1s of length M. Let err be a vector of 1s of length n v +1.

[0332] 9. p j = z j = P -1 r j ; j = 1

[0333] 10. When max(err) > tol and j < maxiter:

[0334] a. Set all elements of B to zero

[0335] b. For i from 1, 2... n:

[0336] i. B += k(x i )k(x i ) T p k

[0337] c. B += λ 2 p j

[0338] d. Perform element-wise division so that α j is a row vector of length n v +1

[0339] e.

[0340] f. r j+1 = r j-(1 M α j )*B

[0341] g. Let Where ||r jk || is r j The Euclidean norm of the k-th column, therefore err is a column of length n. v +1 row vector, where ||r ji ||is r j The norm of the i-th column.

[0342] h.err = err / L1 (elemental division)

[0343] iz j+1 =P -1 r j

[0344] j. Perform element-wise division so that β j It is of length n v +1 row vector

[0345] k.

[0346] l. Storage β j α j

[0347] m.j+=1

[0348] 11. Let α ik It corresponds to the k-th iteration and the probe vector i (corresponding to p). j α of the (i+1)th column k The elements. Let β ik It corresponds to the k-th iteration and the probe vector i (corresponding to p). j β of the (i+1)th column k Element.

[0349] 12.Γ=0

[0350] 13. For i from 1, 2…n v :

[0351] a. Create an empty matrix T of shape (j, j) and set all elements of the matrix to zero.

[0352] b. For k from 1, 2…j:

[0353] i

[0354] c. For k from 1, 2…j-1:

[0355] i.

[0356] ii.

[0357] d.[Y, Θ] = Eig(T); τ m =e1 T y m for m in 1, 2...j

[0358] e.

[0359] 14.

[0360] 15. Return the estimated logarithmic determinant Γ and vector Z. T y and w j The first column; the last one corresponds to (Z) T Z+λ 2 I) - 1 Z T y is a performance term used to calculate the negative marginal log-likelihood.

[0361] To evaluate the accuracy of this program, eight combinations of hyperparameters were randomly selected, each targeting seven different datasets ranging in size from several hundred to half a million data points. For each combination, the negative marginal log-likelihood of 4,096 random features was calculated using Algorithm S5 or matrix factorization. In all cases, the preconditioner rank of 512, the convergence threshold of 1e-5, and the "srht_2" preconditioner construction algorithm were used. Either 40 or 25 probe vectors were used; the distribution of absolute percentage error is as follows: Figure 7 As shown, 40 probe vectors do not provide a significant improvement over 25 probe vectors. In this study, 25 probe vectors are sufficient to achieve high accuracy.

[0362] To further evaluate the accuracy of this procedure, we combined approximate marginal likelihood with Bayesian optimization to fine-tune the hyperparameters of two tutorial / example experiments in the user documentation. The first experiment involved the UCI protein dataset, and the second involved the QM9 dataset with one-hot encoded inputs. In both cases, the fine-tuning procedure improved the best validation set performance achieved by fine-tuning fewer random features using matrix factorization. We further fine-tuned the hyperparameters of tabular data and small molecule experiments on the aforementioned QM9 dataset using the same method to achieve highly competitive performance. These results are consistent with... Figure 7 The data shown are consistent, demonstrating that given a sufficiently large preconditioner rank, a sufficiently small convergence threshold, and a sufficient number of probe vectors, algorithm S5 can estimate the negative marginal log-likelihood with high accuracy for the dataset considered in this paper.

[0363] S6: Improved strategies for hyperparameter fine-tuning

[0364] S5 above demonstrates how to efficiently compute marginal likelihood using the preconditioned conjugate gradient method when dealing with a large number of random features. However, for noisy data, sufficient performance is often achieved by fine-tuning the hyperparameters using only 1000 to 3000 random features. If it is indeed decided to use fewer random features than the fitted version for fine-tuning, a fast and efficient strategy can be introduced. Recall that the marginal likelihood in the random feature approximation is given by:

[0365]

[0366] Z T Z is a Hermitian matrix, therefore as long as it is full rank, it can generate the eigenvalue decomposition PΛP. T Such that P is orthogonal, i.e., PP T =I. Then, the above expression can be rearranged:

[0367]

[0368] For three and four hyperparameter kernels, this simplification can be used in a variety of ways, such as Algorithm S3 and Algorithm S4.

[0369] Algorithm S3: Efficient Grid Search

[0370] Input: Dataset X, where N data points are divided into blocks x1...x n and y1...y n Stored on disk, the initialized kernel has two, three, or four hyperparameters (λ, β, σ1, σ2 if applicable) with M random characteristics, and selected values ​​for hyperparameters (if applicable) σ1 and σ2.

[0371] 1. Initialize the M×M array W and the M×1 array S, setting them to zero.

[0372] 2.a=0

[0373] 3. For i from 1 to N:

[0374] a. Load x j y j And use algorithm S1 to use δ i Formation Q j =z j / β

[0375] b.W+=Q T Q; S+ = Q T y i

[0376] c.a+=y j T y j

[0377] 4. U, Λ=eigh(W+δI) A small constant δ is added to maintain numerical stability.

[0378] 5.Λ-=δ

[0379] 6. S = U T S

[0380] 7. Define loss_fn(β,λ):

[0381] a. Return

[0382] 8.λ best ,β best =argmin λ，β loss

[0383] 9. Score = loss(λ) best ,β best )

[0384] Return rating λ best ,β best

[0385] Algorithm S4: Bayesian Hyperparameter Optimization (Thompson Sampling)

[0386] Input: Dataset X, where N data points are divided into blocks x1...x n and y1...y n Stored on disk, the initial kernel has three or four hyperparameters (λ, β, σ1, or σ2 if applicable) with M random features, maxiter is the maximum number of iterations, each iteration has n candidates, and the convergence criterion is tol.

[0387] 1. For 5-10 randomly selected kernel-specific hyperparameters ω1...ω 10 The evaluation score λ is calculated using algorithm S3. best ,β best .

[0388] 2. Let Ω be the set of all kernel-specific hyperparameter combinations currently evaluated; let ω best For the current optimal combination 3, for niter from 10…maxiter:

[0389] a. Fit an accurate Gaussian process S using a Matern kernel v = 5 / 2.

[0390] b. From Uniform(σ) max , σ min Extract n kernel-specific hyperparameter combinations ω1...ω from ) n

[0391] c. From the location ω1...ω n m samples are drawn from S at each location; select

[0392] ω new =argmin ω (m)

[0393] d. If ||ω new -ω best ||2<tol:

[0394] i. Termination

[0395] e. Evaluation ω new Loss function of algorithm S3

[0396] f. If loss new <loss best :

[0397] i.ω new =ω best

[0398] 4. Return ω best

[0399] Note that S4 uses Thompson sampling, but other sampling functions can be easily replaced (e.g., as expected improvements). A simple grid search for kernel-specific hyperparameters can also be used, which is also implemented in xGPR. The above algorithm can also be improved by calculating Z. T The eigenvalue decomposition of Z is modified, and this calculation, which uses singular value decomposition of z, is the most stable method.

[0400] The polynomial kernel contains only two hyperparameters (λ, the noise parameter, and β, the amplitude parameter, both of which are common to all kernels). Using the above strategy, fine-tuning can be completed in a single data traversal, thus the hyperparameter fine-tuning speed of the polynomial kernel is extremely fast.

[0401] S7: Dataset used in the experiments described in the main text

[0402] All datasets used in this embodiment are publicly available data, and details of each dataset are shown in the table below.

[0403] Table 7. Details of the dataset used in the study described in this paper

[0404]

[0405]

[0406] S8; Comparison of xGPR with stochastic variational inference and exact Gaussian processes

[0407] The stochastic variational approximation of Gaussian process models and their variants appears frequently in the literature. Theoretically, it possesses many attractive properties, but how does it perform in practice? Table 8 compares the performance of the SVGP implementation in the GPyTorch library with the stochastic feature approximation implemented in xGPR; for sufficiently small datasets, it is also compared with the performance of the exact Gaussian process implementation in GPyTorch (when training the exact Gaussian process on a large dataset, an out-of-memory error occurs on the UCI protein dataset containing approximately 8,000 data points, making it clearly impractical for larger datasets). In all cases, the stochastic variational inference was trained in GPyTorch for 40 training epochs, and the RBF kernel was used in all cases because it did not consider whether sequence data convolution kernels compatible with SVGP were efficient enough to handle datasets exceeding 100,000 sequences in a short time, and no such kernels were implemented. In all cases, training for 20 training epochs resulted in a performance degradation.

[0408] For xGPR, hyperparameters are fine-tuned using the "fast" strategy in Algorithm S3 of S6 above with 3,000 random features. Then, tabular data is fitted with either 16,384 or 32,768 random features, and sequence data is fitted with either 8,192 or 16,384 random features (the number of random features used for sequence data is smaller, corresponding to the number of features used for comparison with CNN deep learning models on sequence data). It's important to note that several datasets used as benchmarks here are protein sequence datasets, which are also used as benchmarks for comparison with CNNs. Although sequence kernels that do not require input alignment are used, RBF kernels are used because no sequence kernels are available in GPyTorch or scikit-learn. Therefore, in this comparison, for protein sequence datasets, each dataset is converted to multiple sequence alignment and one-hot encoded before being used as input to the RBF kernel, which differs from the modeling method used for sequence datasets.

[0409] The data in Table 8 clearly shows that xGPR outperforms SVGP on 8 out of 9 datasets, and the advantage is usually significant. While their performance is roughly equivalent when SVGP is competitive (Song UCI), SVGP training time is 80 times longer. xGPR typically performs better with the same kernel and at a significantly lower cost; furthermore, xGPR performance can often be further improved by using more random features for fitting or by fine-tuning hyperparameters with more random features. For example, on the Kin40K dataset, "fine-tuning" hyperparameters using the marginal likelihood approximation strategy from Section S5 above requires an additional 10 minutes and reduces the mean absolute error by an additional 5-8% (depending on the number of random features used for fitting). In contrast, improving SVGP performance is very difficult because increasing the number of induced points in SVGP makes training on larger datasets unacceptable.

[0410] It's possible to improve SVGP's performance by modifying the settings for stochastic gradient descent in the Adam optimizer, for example, by tuning the hyperparameters of the Adam algorithm for a specific dataset. However, in a sense, this is just another drawback of SVGP: xGPR doesn't require extensive experimentation with optimizer settings to achieve good results. Given these results, the library did not implement stochastic variational inference, but instead opted for the methods described above.

[0411] Table 8. Performance of SVGP (implemented in the GPyTorch library), random features (implemented in xGPR), and exact GP (implemented in scikit-learn) on various benchmark datasets.

[0412]

[0413] S9: Lemma 1.1

[0414] Lemma 1.1: set up for On the positive definite kernel, and set From X to Let X be any mapping of u, w ∈ X, where X is some non-empty set. Then, for u, w ∈ X, K(f(u), f(w)) is a positive definite kernel on X.

[0415] prove: set up From X to Let X be any mapping of X, where X is a non-empty set. yes Positive definite kernel on, that is, for any and for any

[0416]

[0417] Because of this property Any set of points x1, x2, x3...x n Both are true, and the mapping is true. Map the input to Therefore, for any x1, x2, x3...x n For X, we have:

[0418]

[0419] In other words, this is just another way of saying that anything that maps input to... Feature engineering processes or functions can all be "plugged into" a It is a positive definite kernel that retains its positive definiteness. This can be used to build "plug-in" kernels, such as the FastConv1d kernel mentioned above.

[0420] S10: Detailed information on QM9 molecular modeling

[0421] The complete QM9 dataset contains approximately 133,000 molecules. First, molecules whose geometrically optimized structures failed the consistency test were removed. Then, using a random seed of 123, the remaining data was randomly divided into 110,000 training sets, 10,000 validation sets, and 10,831 test sets.

[0422] For SOAP descriptors implemented in the descriptor library, several hyperparameters need to be fine-tuned:

[0423] 1) Used to represent the Gaussian standard deviation (sigma) of atoms in the neighborhood of the central atom.

[0424] 2) rcut, which is the radius by which an atom is considered "outside its neighborhood" and is no longer included in the representation.

[0425] 3) Weighting, which reduces the weight of atoms that are far from the central atom (several different schemes are provided).

[0426] In addition, there are the number of radial basis functions, n_max, and the number of spherical harmonic functions, l_max. Increasing these numbers improves the accuracy of the representation but also significantly increases the size of the descriptor vector. n_max is set to 12, and l_max to 9; both values ​​are the same because a larger descriptor vector size was not used in this study. However, even with this configuration, approximately 18,000 features are generated per atom.

[0427] Willatt et al. (Phys.Chem.Chem.Phys.,20(47),29661-29668,2018) described the determined weighting schemes and weights as optimal, which used power-law weighting, where atoms within a distance rcut are weighted according to their distance using the following expression:

[0428]

[0429] For all SOAP feature hyperparameters, values ​​selected by Willatt et al. based on grid search were used, with two exceptions: r0 from the aforementioned weighted scheme, and rcut (cutoff distance), which were optimized starting from the values ​​of Willatt et al. r0 and rcut were fine-tuned by evaluating the marginal likelihood of a randomly selected subset of training data consisting of 25,000 data points. Initially, it was found that generating features only for heavy atoms (i.e., non-hydrogen) and omitting hydrogen atoms could achieve almost the same accuracy. This is advantageous because otherwise the amount of disk space required to store the training set data would be very large. After adjustment, the following settings were adopted and described:

[0430] N_max:12

[0431] L_max:9

[0432] Sigma: 0.25

[0433] R_cut:3.25

[0434] Weighted: Power Mode

[0435] R0:1.5

[0436] M:9

[0437] M and r0 determine the rate at which the weights decrease from the central atom. SOAP features are generated individually for each atom; since each molecule in this dataset contains at most 9 heavy atoms, each molecule has at most 9 descriptor vectors corresponding to these atoms. These descriptor vectors can be used as input to the graph convolution kernel or graph polynomial kernel described above. The SOAP descriptor of each atom is divided by the norm of the descriptor vector, making each descriptor vector an identity norm; this ensures that if the GraphPoly kernel is used, it will have the same interpretation as the original SOAP kernel.

[0438] To fine-tune the kernel hyperparameters, we used the workflow described in S5 and S6 above. First, a starting point was found using a matrix factorization-based marginal likelihood method with a small number of random features, followed by Bayesian optimization with a large number of random features (16,384). This process took 4 hours on a single A6000 GPU. All hyperparameter adjustments were targeted at an internal energy of 298K. As shown in the figure, hyperparameters that perform well at an internal energy of 298K are also suitable for other targets. Then, the final model was fitted using 16,384, 32,768, or 66,536 random features. The results are shown in Table 4 above.

[0439] Using more random features can yield a small performance improvement, but even 16,384 random features are sufficient to achieve a mean absolute error that is 5 times lower than the accuracy of chemistry. Therefore, the practicality of using more random features is limited in practice; for this problem, 16,384 random features are perfectly adequate. Furthermore, the most important target for further improvement is likely the SOAP features themselves. The large number of SOAP features generated for each atom significantly increases computational cost because many features must be generated, saved, and loaded from disk when processing each mini-batch of data during iterative dataset processing. Moreover, it limits the ability to improve model accuracy by increasing the number of basis functions. Several authors have proposed more compact representation methods to achieve similar or higher accuracy. More compact representation methods can increase the number of basis functions, thereby improving model accuracy, while reducing the length of feature vectors, thus lowering their computational cost.

[0440] It's worth noting that the Conv1d kernel based on the Fast Hadamard Transform (FHTConv1d) also offers advantages. Using the raw random features, generating 66,536 random features from 18,000 input features would require a matrix containing approximately 1.2 billion floating-point numbers. In contrast, the FHTConv-1d kernel only needs to store approximately 288,000 8-bit integers and 72,000 floating-point numbers—significantly reducing memory usage.

[0441] S11: Clustering and kernel PCA using the representations generated by xGPR

[0442] The QM9 dataset was used here. From approximately 130,000 molecules in the dataset, 20,000 molecules were randomly selected as the test set. Each atom and its neighbors were encoded using one-hot encoding, up to the 15th farthest neighbor (an arbitrary cutoff value). The weight of each one-hot encoded neighbor was 1 / distance**6; the 6th power used here corresponds to the weights used in the London matrix, which is sometimes used as a representation for small molecule machine learning. Following a standard workflow, the hyperparameters were first tuned using a matrix factorization-based procedure with a small number of random features (2048), and then “fine-tuned” using an approximate marginal likelihood procedure with a large number of random features (8192). As shown, the model achieved a mean absolute error of <1.1 kcal / mol, depending on the number of random features used for fitting.

[0443] Next, the kernel_xpca tool provided by the visualization_toolkit module in xGPR version 0.0.2.1 is used to retrieve the first 500 principal components. For example... Figure 8 As shown, the kernel PCA plot is generated using only the first two principal components for the test set data.

[0444] As shown, the top principal components correlate well with the predicted true labels (energy 298K). However, this is not always the case, as using a 2D representation instead of the 16,384-dimensional representation generated by xGPR discards a significant amount of information. To further explore the data, we used the scikit-learn library, performing k-means clustering with the first 500 principal components as input, and used the sum of squared errors to evaluate the performance of a given number of clusters. This operation is approximated by performing kernel K-means clustering on the original training data using the GraphConv-1d kernel.

[0445] Elbow diagram as follows Figure 9A As shown, this result indicates that 5 clusters may be appropriate, therefore the final K-means model is also fitted using 5 clusters. Figure 9B The energy distribution of each cluster at 298K is shown. In many drug development scenarios, the ability to cluster data using the same model as the fitted data can be useful in itself, as it can be used to efficiently retrieve similar molecules from the training data, as described above.

[0446] S12: Molecules from the test set and the most similar molecules from the training set.

[0447] Figure 10-14Five random molecules from the QM9 test set and the 50 most similar molecules from the training set are shown, with similarity quantified using a (approximate) GraphRBF kernel function. The molecules are represented using one-hot encoding as described above. It is noteworthy that this molecular representation method has relatively little information content, which is why the model performs so well (MAE approximately 1 kcal / mol) and is able to retrieve molecules that do indeed share many similarities in many cases, which is somewhat surprising.

[0448] in conclusion

[0449] Gaussian processes represent a principle-based Bayesian machine learning approach capable of directly quantifying uncertainty. Their most significant drawbacks are unacceptable scaling and a lack of efficient kernel functions for sequences and graphs. This embodiment describes efficient methods for fitting Gaussian process models using stochastic feature approximation kernel functions, achieving linear scaling with training set size and with molecule size or sequence length. These tools are implemented in the publicly available Python library xGPR, which has proven competitive across a wide range of protein and small molecule property prediction tasks.

[0450] As described in this article, GPs equipped with these cores are able to generate well-calibrated uncertainty estimates for their predictions, unlike any other deep learning model used as a comparison. This property is crucial for protein engineering and drug development; since evaluating predictions is often costly, high-confidence predictions are naturally preferred whenever possible.

[0451] Furthermore, GPs equipped with these kernels are more interpretable than deep learning models used as comparisons. For protein sequences and small molecules, each kernel can explain in detail how the model determines the similarity between any two data points, and thus how it makes predictions. This helps in understanding the model's limitations and how to improve it. Current popular deep learning architectures for sequence or graphical data cannot generate such detailed explanations because they are all "black boxes." Moreover, as with small molecules, the kernels of a trained GP regression model can be used to cluster the input data or perform kPCA, which can provide more insights into the distribution of the input data. This capability is particularly useful for drug design, allowing small molecules to be classified into subgroups.

[0452] This combination of competitive accuracy, uncertainty quantification, and improved interpretability suggests that the xGPR library may be useful for a range of protein and small molecule property prediction tasks.

[0453] The kernels discussed in this embodiment all accept a single protein or small molecule as input. For certain tasks (e.g., protein-small molecule interactions), it is preferable to use “paired” kernels that accept two inputs (e.g., a small molecule and a protein sequence) and compare them with other pairs. The kernels discussed here can be easily modified to generate multiple “paired” kernels that may be suitable for certain problems. For example, the first layer of the FastConv-1d kernel discussed here can be used to generate characterizations of small molecules and proteins; these characterizations are then concatenated and used as input to an RBF kernel. Alternatively, a GraphRBF kernel can be used to generate random features for small molecules, while an FHTConv-1d or FastConv-1d kernel can be used to generate random features for proteins; these random features are then concatenated into a single vector to generate predictions. This last arrangement is equivalent to using the summation of two kernels. It is possible to design other efficient kernels specifically for such pairwise comparisons.

[0454] Example 2

[0455] The Rapid Protein Engineering System (RESP) AI model accelerates the recognition of tightly bound antibodies.

[0456] This embodiment describes an easy-to-implement machine learning-assisted process for identifying high-affinity antibodies that address these challenges. A novel autoencoder model was trained based on over 3 million B-cell receptor sequences, demonstrating that this learned representation method delivers better results than state-of-the-art embedding schemes on the target task. Next, a Bayesian neural network was developed and trained for ordinal regression analysis to model the relationship between sequences and binding affinity or dissociation rate using directed evolutionary data. Importantly, the model is able to estimate the uncertainty in its predictions, which is not achievable with current deep learning methods. To benchmark the method, it achieves competitive accuracy on a literature dataset while providing useful uncertainty information not provided by deep learning models in the original relevant studies. Furthermore, simulated mutagenesis was performed using a simulated annealing strategy to explore sequences not present in the mutation library, and their binding affinity was evaluated for experimental assessment. The autoencoder, the Bayesian network for ordinal regression, and the search strategy together constitute the key computational components of this RESP process.

[0457] To demonstrate the function of RESP, we attempted to engineer atezolizumab (trade name Tecentriq) heavy chain mutants to enhance their affinity for programmed death-ligand 1 (PD-L1), thereby improving the affinity of known antibodies for well-studied antigens. We largely randomized the atezolizumab heavy chain sequences, binned the mutants using yeast display and FACS, and sampled from each bin for sequencing. We used an autoencoder model to convert the mutant atezolizumab sequence data into a low-dimensional characterization, and trained a Bayesian ordinal regression model using the encoded mutant atezolizumab sequences, scoring each sequence based on its probability of being a strong PD-L1 binder. Furthermore, an improved simulated annealing algorithm was used to select sequences for testing. Results showed that this process identified a group of atezolizumab scFv mutants with improved PD-L1 dissociation rates, and one identified mutant showed an approximately 10-fold reduction in the dissociation rate between human PD-L1 and atezolizumab. D This represents a 17-fold increase. The method described here can be used as a general approach to improve antibody-antigen interactions while reducing experimental workload. Mutants discovered using this method could be useful agents for treating PD-L1-positive tumors.

[0458] method

[0459] software

[0460] Analysis and modeling were performed using Python 3.9, with PyTorch version 1.8.1, NumPy version 1.19.5, SciPy version 1.5.4, and scikit-learn version 0.24.2. 0 was used as a random seed for model weight initialization, training / test splitting, etc.

[0461] Antibody-encoded autoencoder model

[0462] The training set used by the autoencoder came from the cAb-Rep database, which contains a library of antibody sequences 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 high-depth antibody library dataset was used. This dataset was further filtered by numbering all sequences using ANARCI software with Chothia numbers. During this process, ANARCI software used HMMer software to align the input sequences to the domain-specific Hidden Markov Model (HMM) profiles of human antibodies. Any sequences with low bit scores resulting from this alignment were rejected, minimizing the possibility of non-antibody proteins in the dataset. After ANARCI numbering and filtering, 2,725,492 sequences remained. These sequences are antibodies because they are derived from antibody sequence databases and have been screened to identify those most likely generated by antibody MSA. It is important to note that only the VH heavy chain is used here to match the sequences used in the experimental work.

[0463] The autoencoder accepts one-hot encoded sequences as input. To number the positions of each antibody, the Chothia numbering scheme is used. In this scheme, a heavy chain may contain up to 132 amino acids, but due to the different CDR sizes of different antibodies, our mutants contain fewer amino acids. Therefore, each one-hot encoded sequence is a matrix with 132 rows and 21 columns. Although only 20 amino acids are present, the 21st position represents a blank because the heavy chain does not contain all 132 amino acids, thus resulting in blanks at some Chothia numbering positions. It is also possible to simply leave blank positions as all zeros, but for clarity, we prefer to explicitly encode blanks, and this method is adopted in this paper. Each row contains a 1 at one position to indicate which amino acid is present (or a blank if no amino acid is present).

[0464] Typically, when selecting and designing sequences, it is crucial to determine which positions are most important for affinity or other desired properties. This consideration dictates that the learned representation generated by the autoencoder should have the same number of rows (positions) as the input. Furthermore, to ensure that a particular row of the encoding contains information related to that position, the decoder should be able to use the information from that position (and possibly neighboring positions) to reconstruct the amino acids present at each position.

[0465] To meet these constraints, the encoder part of the model is designed to compress the input from a 132×21 matrix to a 132×3 matrix. The encoder module consists of a convolutional network with two convolutional layers with appropriate padding and a third linear layer. The first convolutional layer contains 40 kernels, each with a width of 21, while the second convolutional layer contains 20 kernels, each with a width of 11. Each convolutional layer uses a gated convolutional activation function to improve the performance of the convolutional neural network on language modeling tasks. For this activation function, the sigmoid function is applied to the first half of the layer's output, and then multiplied element-wise with the second half of the output. In other words, if the output of a layer for a given sequence is a 132×40 matrix, then for the gated convolutional activation function, the sigmoid function is applied to the first 20 columns, and then multiplied element-wise with the remaining columns to obtain a final output of size 132×20.

[0466] The final output of the linear layer in the autoencoder is a 132×3 matrix, which is the encoding of the original sequence. This encoding serves as the input to both the decoder and prediction modules. The decoder module consists of a single layer of the following form:

[0467] (1) softmax(a·W+b)

[0468] Here, 'a' represents a line of the encoder output. Since the encoder output is a 132×3 matrix, 'a' is a three-dimensional vector. 'b' is the learned bias vector, and 'W' is the learned 3×21 weight matrix. Softmax is the softmax function.

[0469] (2)

[0470] Where z is a 21-element vector obtained by a·W+b. The decoder layer is applied to each line of the encoder output. This design imposes a strict constraint on the autoencoder: the model needs to reconstruct the input using a decoder function that shares relatively few parameters across all positions.

[0471] The encoder's output is also provided to a prediction module, which distinguishes between human B-cell receptor sequences and junk sequences. This prediction module consists of a simple logistic regression model, where the probability that the sequence is an antibody is given by the following equation:

[0472] (3)

[0473] Where 'a' is the complete output of the encoder, flattened from a 132×3 matrix into a vector of length 396; and 'w' is the learned weight vector of length 396 plus the learned bias term. Since both the prediction and decoder modules have very few parameters (compared to typical deep learning models), they cannot learn the complex mapping from the encoded sequence to the input, thus placing the burden of generating the most informative and relevant encodings as possible on the encoder.

[0474] To generate junk sequences, we generated mutant versions of 2,725,492 sequences selected from the cAb-Rep database. Our goal was to force the model to distinguish between true antibody sequences and closely related sequences, thereby compelling the model to integrate information about the relative abundance of specific amino acids at specific positions into its encoding. Random sequences could also be generated simply, but this does not force the model to generate informative encodings; random sequences differ so much from antibodies that logistic regression models used to predict whether a sequence is an antibody can distinguish them without modifying the representation generated by the encoder. Therefore, it is best to generate decoy sequences that are similar enough to be difficult to distinguish, yet distinct enough to exhibit modifications rarely found in true antibody sequences.

[0475] This experiment showed that seven mutations provided a good balance; therefore, copies were replaced with randomly selected amino acids at seven randomly chosen positions. The end result of this process was a library containing approximately 6 million sequences, half of which were human B-cell receptor sequences and the other half were not. Therefore, an autoencoder model could be trained to encode both the input sequence and embed canonical features observed in real antibody sequences. This autoencoder was implemented using the PyTorch library in Python 3.6.9 and trained on the complete dataset of 6 million sequences until convergence. The code for this step, as well as all other steps described in this paper, is available online at (10.5281 / zenodo.7508853).

[0476] Accuracy on the prediction and reconstruction tasks was evaluated using a reserved test set. These metrics are for diagnostic purposes only, as they assess the autoencoder's ability to reconstruct the input. The real test for the autoencoder is the impact of its generated encodings on the prediction accuracy of a model trained using that learned representation as input. To evaluate this more critical metric, we encoded wild-type and mutant atezolizumab libraries (constructed as follows) using the aforementioned autoencoders, one-hot encoding, the ProtVec encoding scheme, the UniRep encoding scheme, the FAIR-ESM encoding scheme, the AbLang sequence encoding scheme, and the AntiBertY embeddings, respectively. For AbLang, the seq-coding option was used. For the AntiBertY embeddings, the average of all residue-specific representations in the sequence was taken. For AntiBertY, initially all residue-specific embeddings were used without averaging, which resulted in poor performance. A Bayesian neural network (structure described below), a fully connected neural network with the same number of layers and weights as a Bayesian network, was trained, and a random forest model using these available encodings was used as a baseline. Five-fold cross-validation on the training set and independent test sets were used to evaluate the classification performance of each model under each encoding scheme. Matthews correlation coefficient (MCC) and AUC-ROC were used as metrics to distinguish RH03 from others. Compared to one-hot encoding in all three model types (especially the Bayesian network, which is the most important in this study), a good encoding should improve or at least not harm performance.

[0477] Generation of WT (wild-type) and mutant atezolizumab scFv libraries

[0478] First, WT atezolizumab scFv was cloned into the pYD1 yeast display vector to test its function on the yeast surface (see [link]). Figures 24A-24B For the first atezolizumab library, PCR amplification was performed using the WT plasmid as a template. The WT light chain was prepared using Q5 hot-start DNA polymerase (NEB), or the mutant heavy chain was prepared using Taq polymerase (Invitrogen) via error-prone PCR (as described above). The two PCR products were assembled into one product by overlap extension PCR and co-electroplated with the linearized pYD1 vector into EBY100 yeast (as described). Based on the colony counts after serial dilutions on selection plates, the library transformation yielded 7.8 × 10⁸ colonies. 7 A transformant.

[0479] Atezolizumab scFv library was screened using yeast surface display.

[0480] The optimal competition time for dissociation rate screening was determined using the previously described method. WT or mutant yeast libraries were thawed, inoculated into selective growth medium, and cultured at 30°C for 22 hours. The library / WT was induced at 20°C in galactose-induced medium for 42 hours. Subsequently, the WT or library was labeled with biotin-PD-L1 (Sino Biological 10084-H08H-B) in TBS-BSA at room temperature for 3 hours, and then incubated with an excess of non-biotinylated PD-L1 (Sino Biological 10084-H08H) for the determined competition time (room temperature, TBS-BSA). For FACS, yeast cells were labeled in TBS-BSA at 4°C with anti-V5 antibody (Thermo Fisher R960-25, formerly known as 46-0705), and then labeled in TBS-BSA on ice for 30 minutes with SA-PE (BD554061) and goat anti-mouse IgG2a AF647 (Thermo Fisher A21241). Cells were sorted according to fast, medium, and slow dissociation rates (see Figure 16A). 1-2 and Figure 34 The target organisms were cultured at 30°C to a high density and then cryopreserved at -80°C.

[0481] Preparation of mutant DNA libraries with faster WT levels and slower dissociation rates compared to PD-L1

[0482] To isolate plasmids from the initial library screening, plasmids were harvested from yeast using the Zymoprep Yeast Plasmid Mini-Extraction II Kit (Zymo Research), eluted with ddH2O, concentrated using the DCC-5 (DNA Clean and Concentrator 5) Kit (Zymo Research), and then subjected to 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 using Q5 Hotstart DNA polymerase (NEB) with a moderate number of PCR cycles (17 cycles, determined to be optimal by qPCR) to avoid library over-amplification; the PCR products were concentrated using the DCC-5 kit and purified by agarose gel extraction. The second round of PCR amplification was performed using different primers (low / medium / high binding NGS F and NGS R) to barcode the binding groups of the MiSeq Nano PE250 (see Table 11), with KAPA HiFi Hotstart Premixed Reagent (KK2601) used, and a cycle count of 5 (determined optimal by qPCR). The PCR products were purified using AmpureXP beads (A63880) and submitted for QC using a TapeStation analysis system (Agilent HSD1000), and then sent to the Institute for Genomic Medicine at the University of California, San Diego for MiSeq sequencing.

[0483] Objective 21 - Generation and screening of mutant libraries (see also S2.4 below)

[0484] Twenty-one gene fragments (IDTs) from 21 mutant heavy chains were fused to the light chains via overlap extension PCR (Q5 hot start), and then co-electroplated into EBY100 yeast cells with the linearized pYD1 vector. The method was similar to that used in the first library (>10). 7 The transformants were identical. The screening method for the obtained library was basically the same as that for the first library, but the final sorting process required a very strict competition time (39 hours at room temperature). The target was obtained by small-scale extraction of yeast plasmid, transformed into GC10 competent cells, and then obtained by small-scale bacterial extraction. Finally, sequencing was performed using the standard Sanger sequencing method.

[0485] WT atezolizumab, mutant 4, durvalumab, and averumab scFv were cloned into the pYD1 vector.

[0486] Gene fragments of durvalumab and avelumab were purchased from IDT. Their codons were optimized for yeast use, and PCR amplification was performed using the following primers: Atez LE F and Durv LE R for durvalumab, and Atez LE F and Avel LE R for avelumab (Table 11). The WT / mutant 4 gene was amplified from its plasmid by PCR using Atez LE F / R primers. The PCR products were double-digested with XhoI / NheI-HF (NEB) and ligated into a linearized / dephosphorylated pYD1 vector (using the same enzymes for digestion and dephosphorylation / ligation using 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 Small Scale Preparation Kit (Zymo Research) and sequenced before transformation into EBY100 yeast. These constructs were used for yeast k off and K D Measurement.

[0487] WT / Mutant 4 / Dvorumab / Averumab on the surface of yeast K D and k off Determination of value

[0488] The method described in Section S2.3 below was used to determine k in a TBS-BSA. off The obtained data were fitted to a single-phase decay model using GraphPad Prism 9.3.0 software, where the following formula was used (Y is the yeast fraction binding to biotin-PDL1, X is time, Y0 = 1 (the yeast fraction binding to biotin-PDL1 at the 0s competition time point), the plateau value is a constant based on non-specific yeast binding, and K is K0. off ):

[0489] (4) Y = (Y0 - plateau value) * exp(-K * X) + plateau value

[0490] K on the surface of yeast D The values ​​are basically as described. 49 The measurements were performed, and the data were fitted to the following equation:

[0491] (5)

[0492] (Bmax is the maximum MFI value, and X is the PD-L1 concentration). Because mutant 4 has a very slow dissociation rate, yeast needs to be incubated with PD-L1 at room temperature for 6 days.

[0493] Cloning and purification of WT / mutant 4 atezolizumab scFv (see also S2.6 below)

[0494] The scFv sequence was amplified by PCR 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 protein was purified primarily according to the described method. The scFv was stored frozen in PBS + 10% glycerol at -80°C. Purity was assessed using SDS-PAGE. Figure 28A and 28B The protein concentration was determined by BCA analysis (Pierce #23227).

[0495] K was determined using the biolayer interferometry (BLI) method. D value

[0496] BLI was performed at 25 °C on an Octet Red96 at the Biophysics and Biochemistry Core Laboratory of The Scripps Research Institute. Binding reactions were carried out in 1X kinetic buffer (Sartorius, 18-1105) containing 20 mM phosphate-buffered saline (pH 7.6), 2 mM KCl, 150 mM NaCl, 0.02% Tween 20, 0.1% BSA, and 0.05% sodium azide. Biotinylated PDL1 (Sartorius 10084-H08H-B) was immobilized on the SA biosensor (Sartorius, 18-5019) by immersing the sensor in 100 nM b-PDL1 solution until signal saturation. A blank loading channel was used as a negative control. Kinetic experiments were performed using wild-type scFv and scFv mutant 4. The experiment used gradient concentrations, with wild-type scFv concentrations ranging from 475 nM to 0.6 nM and scFv mutant concentrations ranging from 160 nM to 0.2 nM. Binding kinetics were evaluated using a 2:2 binding model with Octet Data Analysis HT software version 12.

[0497] Mass spectrometry analysis

[0498] Mass spectrometry analysis was performed at room temperature on a Refeyn TwoMP mass spectrometer at the 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 identification was performed in 1X PBS (Cytiva BR100672 phosphate buffer). For each experiment, the scFv or PD-L1 (or both) was diluted to the final concentration with PBS. 60-second images were generated for each experiment, 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).

[0499] Sequence processing

[0500] The raw paired-end reads from the atezolizumab dataset underwent quality checks (see the "Sequence Quality" section for details on the screening criteria). After sequence processing, the sequences were split into an 80% training set and a 20% test set. All training and cross-validation were performed only on the 80% training set.

[0501] In many cases, sequences appear in multiple categories, but with varying frequencies across each category. Clearly, the confidence level of category assignment decreases when the frequency of a sequence in a given category is not significantly different from its frequencies in other categories. To quantify the confidence level, a weighting method is used: increment the frequency of each sequence in its respective category by 1, divide by the total number of occurrences, and then increment by 3. This is the posterior probability calculated using a multinomial likelihood function and a uniform Dirichlet prior (i.e., a Dirichlet distribution with parameter α = [1,1,1]).

[0502] The ordinal regression model used here employs binary cross-entropy loss:

[0503] (6)

[0504] Here, for data point i, p is the probability assigned by the model, representing the probability that the model classifies the sequence as belonging to a more tightly bound category than category j—involving M possible categories and N data points; y(i,j) is 0 or 1, used to indicate whether the sequence belongs to the more tightly bound category (see below for details on the calculation method of p). The loss function of sequence i is multiplied by the weight of sequence i. Thus, when the confidence of the class assignment is low, the penalty for misclassified sequences by the model will be correspondingly weakened.

[0505] Construction of Ordinal Regression Model

[0506] Ordinal regression analysis based on a Bayesian neural network was performed on the atezolizumab mutant library data. Each atezolizumab mutant sequence was first encoded by an autoencoder to generate a 132×3 matrix, which served as input to the atezolizumab model. The model architecture is shown in Figures 23A-23B. Overall, it is similar to the backpropagation Bayesian architecture described by Blundell et al. (ArXiv150505424Cs Stat, 2015), but the model has been adapted for ordinal regression analysis, as described below.

[0507] In traditional neural networks, each parameter is learned and remains fixed once the model is trained. In contrast, in variational networks, each parameter has an associated Gaussian distribution, described by its mean and standard deviation, whose optimal value is learned during training. To generate predictions, instead of using fixed learned parameters as in traditional networks, we sample N weight distributions from the weight distribution to generate N predictions. The variance of these predictions measures our uncertainty about the final prediction. This approach requires a more complex training process compared to the traditional neural networks discussed below.

[0508] In short, we seek the parameters θ of the weight distribution w, which will minimize the lower bound of evidence or ELBO, as follows:

[0509] (7)argmin θ KL[q(w|θ)|||p(w)]-E q(w|θ) [log(p(D|w))]

[0510] Where q(w|θ) is the distribution of the weights, approximating the true Bayesian posterior probability, p(w) is the prior probability, and p(D|w) is the likelihood of the data. The parameter θ is the normal distribution parameter of all weights in each hidden layer, as shown in Figure 23A. 1-2 As shown in Figure 23B. The second term is the negative log-likelihood, while the first term is the regularization term, which measures the difference between the weight distribution learned by the model and the prior probabilities. We chose Cauchy prior probabilities with a unit scale and zero positions, which adds some extra flexibility since some weights are now expected to be outliers.

[0511] Costs are estimated using a Monte Carlo sampling procedure. Using the current parameter set θ (mean and standard deviation of the weight distribution), n weighted samples are drawn for each mini-batch to derive the following approximate cost function:

[0512] (8)

[0513] Where M is the number of mini-batches, and other terms are the same as above. In other words, for each mini-batch of training data, n sets of sample weights are plotted, and then the average of these weights and the approximate cost function of all mini-batches in the training set is taken. Backpropagation of this cost function can be handled using Kingma and Welling's reparameterization technique. Samples are taken from a standard normal distribution with a mean of zero and a standard deviation of 1, and then the means are summed and multiplied by the standard deviation of the distribution with weights j. To ensure that the standard deviation is always positive, we do not use the standard deviation itself as a parameter, but instead parameterize each distribution with parameter ρ and convert it to the standard deviation using the softplus function:

[0514] (9)σ=log(1+e ρ )

[0515] Using this reparameterization technique, the gradient of the approximate cost function for each weight and bias term in each hidden layer with respect to the mean and ρ of the distribution can be easily computed. This model structure has a dual advantage. First, it imposes strong regularization on the model parameters, which can improve performance on certain tasks. Second, it estimates the uncertainty in the predictions, thus assigning greater weight to the most reliable predictions when selecting test sequences.

[0516] All these characteristics of this model are the same as most Bayesian neural network architectures. Therefore, this model is used to perform ordinal regression. As shown in Figures 23A-23B, the last hidden layer of the Bayesian network outputs a latent score. For the case of M classes, this score is added to M-1 learned thresholds. Then, a sigmoid function is applied to each of the M-1 outputs to generate an output vector containing M-1 probability values. Each element i of this output vector is the model-assigned probability that the input sequence belongs to a more stringent combination class than i. Since there are three combination classes (RH01, RH02, and RH03 or weak, medium, and strong), the output vector is two-dimensional; the first element represents the probability that the sequence belongs to RH02 or RH03, and the second element represents the probability that the sequence belongs to RH03. The model is trained by minimizing the binary cross-entropy loss described above.

[0517] This arrangement ranks the categories hierarchically: sequences with higher potential scores are assigned to higher categories. It's important to note that this does not predict actual dissociation rates or binding affinity. Instead, the score reflects the confidence of a sequence relative to other sequences in the training set that they are strong binders. This method has been used in the literature to rank sequences for protein engineering tasks and select them for experimental evaluation.

[0518] Sequence scoring and selection

[0519] The classic simulated annealing algorithm has been modified and equipped with a trained model for computer-directed evolution, such as... Figure 17 As shown. First, the frequency of each amino acid at each position in the entire dataset is calculated, all values ​​are incremented by 1, and then divided by the total number of sequences in the dataset plus 20 to obtain the marginal probability of each amino acid at each position. By incrementing all frequencies by 1, the marginal probability of unobserved amino acids in the dataset can be obtained; this probability, while not zero, is very small.

[0520] Next, the top 500 sequences with the highest scores are selected from the atezolizumab dataset, and the 10 positions with the highest mutation frequencies among these sequences are determined. The selection of these 10 sites is arbitrary, as more or fewer sites can be used as needed. In each iteration, any one of these top ten sites is selected with equal probability. The selected site is randomly reassigned a new amino acid; the probability of selecting any given new amino acid is determined by the marginal probability calculated as described above. For example, suppose that 80% of all sequences observed in the dataset carry arginine at position 100, another 10% carry lysine, and so on. If position 100 is selected in a given iteration, a mutation will occur at that position, with an 80% probability of converting to arginine, a 10% probability of converting to lysine, and so on.

[0521] Both the current sequence and the mutated proposed sequence are encoded using an autoencoder, and a score is assigned to each sequence using a trained ordinal regression model. The probability of a proposed sequence being accepted is as follows:

[0522] (10)

[0523] Where T is temperature, S best and S proposed These are the best score to date and the score of the proposed sequence, respectively.

[0524] The algorithm is further explained as follows:

[0525] Improved simulated annealing algorithm

[0526]

[0527] This process is a simple method for exploring the sequence space while ensuring it doesn't deviate too far from the training set. To ensure repeatability and avoid random fluctuations, the Bayesian neural network uses the mean of the distribution of each weight as the weight value to generate predictions during this modeling phase.

[0528] S2.1: Sequence quality

[0529] Sequence reads containing one or more base pairs and a Phred quality score <10 were discarded, as these reads might be unreliable in such cases. If the overlapping regions at both ends did not match, both reads were discarded to prevent mismatches between paired ends. Reads meeting these quality criteria were merged and translated to produce a mutant atezolizumab sequence of 118 amino acids. Each sequence appearing more than once in any given category was assigned a frequency for that category, indicating the number of times it was found. Sequences with mutations in the first 30 or last 8 positions were excluded from further consideration because these positions were not targeted for mutation, and mutations at these positions were extremely rare, making it difficult to assess their importance.

[0530] S2.2: Generation of WT and mutant atezolizumab scFv libraries

[0531] To test the function of WT atezolizumab scFv in yeast display form, the gene for WT atezolizumab scFv was purchased from Integrated DNA Technologies in the form of a gene block (using yeast-optimized codons). The gene block was cloned into the modified pYD1 (Addgene #73447) yeast display vector by PCR amplification using AtezF / AtezR primers (Table 11). The digested gene block was then double-digested / ligated into the digested vector using Esp3I (Thermo Scientific) and T4 DNA ligase (New England Biolabs (NEB)) via Golden Gate Assembly. The sequence was validated before conversion to EBY100 yeast by heat shock. For the atezolizumab library, WT scFv from the pYD1 vector was used as a template. High-fidelity PCR was performed on the light chain region using Q5 hot-start DNA polymerase (NEB) (primers IF1F / IF1R), and error-prone PCR was performed on the heavy chain region using Taq DNA polymerase (Invitrogen 18038-018, IF2F / IF2R). 8-oxo-dGTP (TriLink N-2034-1) and dPTP (TriLink N-2037-1) were used in the reaction. To improve the mutation rate, error-prone PCR was performed for 30 cycles, with a final concentration of 20 μM for each dPTP and 8-oxo-dGTP (starting concentration 200 μM). The PCR products of the light and heavy chains were concentrated using the DCC-5 (DNA Clean and Concentrator 5) kit (Zymo Research) and purified using an agarose gel extraction kit (Zymo Research gel extraction kit). Two DNA fragments (WT light chain, mutant heavy chain) were assembled into single DNA molecules via overlap extension PCR using primers OF / OR (all primers in Table 11) and Q5 hot-start DNA polymerase. The resulting products were concentrated again using a DCC-5 kit and purified by gel extraction. TA cloning (Invitrogen TA Cloning Kit, K202020) was performed on the products with 3'A added using Taq polymerase to examine the mutation rate in the heavy chain. Random sequencing of the clones showed 0–8 mutations per gene. In preparation for yeast electroporation, the pYD1 vector was double-digested with XhoI and EcoRI-HF, concentrated using a DCC-5 kit, and the gel was extracted / eluted into ddH2O.

[0532] The scFv library was transformed into yeast via electroporation (3 μg of pYD1 digested, 9 μg of scFv gene per electroporation), and then transformed via homologous recombination. 6Assembly was performed in EBY100 cells (ATCC MYA-4941), and the yield of 7.8 × 10⁸ cells was determined by serial dilutions on plates containing selective growth medium (16.7 g / L BD bacterial agar on solid medium). 7 One transformant. The library was passaged several times in selective medium to ensure that each cell contained only one plasmid, and then frozen in equal portions (3.6 × 10⁻⁶ cells) at -80°C in 85% ddH₂O, 10% glycerol, and 5% DMSO. 8 (cells / tube). As a sorting control, the WT scFv gene / pYD1 was also electroporated into EBY100 to ensure its position in the vector was identical to that of the library.

[0533] S2.3: Screening atezolizumab scFv libraries via yeast surface display

[0534] First, to determine the optimal competition time, we determined the koff (dissociation rate) of WT scFv and PD-L1 using essentially the same method as before. Yeast expressing WT scFv was washed with TBS-BSA and then subjected to a reaction at 16 × 10⁻⁶ ppm. 6 Cells were resuspended in 400 μL. Then, 12 μL of b-PDL1 (final concentration 178 nM, Sinopharm 10084-H08H-B) was added to the yeast, and the cells were labeled at room temperature for 3 hours. Afterward, the yeast was washed with TBS-BSA and resuspended in 800 μL, and 174 nM of unlabeled PDL1 (Sinopharm 10084-H08H) was added. 2 × 10⁻⁶ cells were added to each Eppendorf tube. 6 Cells were collected and rotated at room temperature for x hours (0–21 hours). For each time point, cells were washed three times with TBS-BSA (cold) and stored in 198 μL TBS-BSA on ice until the experiment was completed. After collecting all data points, 1:100 V5 antibody (R960-25, formerly known as 46-0705) was added, and yeast was labeled on ice for 1 hour. Yeast was washed with TBS-BSA (cold) and resuspended in 200 μL 1:100 SA-PE (BD 554061) and goat anti-mouse AF647 (Thermo Fisher Scientific A21241) and labeled on ice for 30 minutes. Subsequently, yeast was washed with TBS-BSA (cold) and resuspended in TBS-BSA (cold) for flow cytometry analysis. The optimal competition time was determined to be 29 hours at room temperature. To screen for libraries targeting PD-L1, thaw the WT or library and add it to selective growth medium (20 × 10⁻⁶). 6One cell was added to 50 mL of growth medium + ampicillin and incubated at 30°C with shaking for 22 hours. Cells were then induced by cell precipitation, washed once with sterile ddH2O, and resuspended in the same type of medium, but with 20 g / L galactose instead of glucose. The pH was adjusted to 5.6 with ampicillin and the OD value was set to [missing value]. 600 =0.2). scFv was induced by shaking at 20°C for 42 hours. After 42 hours, the culture flask was placed on ice, and 40×10⁻⁶ cells were added. 6 Yeast cells (WT or library type) were centrifuged at 17,000x G for 30 seconds to precipitate, resuspended in 1 mL TBS-BSA (25 mM Tris, 150 mM NaCl, pH 7.5, 5 mg / mL BSA), and resuspended again in TBS-BSA buffer at a concentration of 40 million cells / mL. Subsequently, WT or library cells were labeled with biotin-PD-L1 (Sinochem 10084-H08H-B) at a final concentration of 178 nM and incubated at room temperature for 3 hours. The cell pellet / washed with TBS-BSA was then resuspended in non-biotinylated PD-L1 (Sinochem 10084-H08H, 174 nM) and incubated at room temperature for varying times (mostly 29 hours, but one selection was done at 7 hours). To facilitate FACS sorting, cells were precipitated / washed three times with cold TBS-BSA, and then labeled on ice for 60 minutes at a 1:100 dilution with anti-V5 antibody (Thermo Fisher Scientific R960-25, specificity verified by the supplier, lot number 2001339 / 2249078 / 2212258). 7 Cells / mL were precipitated / washed with cold buffer, and then labeled on ice for 30 minutes with various secondary detection reagents at a concentration of 1:100 (PE streptavidin BD 554061, goat anti-mouse IgG2a AF647 (Thermo Fisher Scientific A21241)). 7 Cells / mL. Cells were then precipitated / washed with cold buffer and resuspended in cold TBS-BSA buffer for use on BDFACSDiva v.8.0.1 software at BD FACSAria II (UC San Diego Moore Cancer Center), see [link to BDFACSDiva software]. Figure 34 Example of an FACS plot. To set sorting gating to collect mutants with faster, wild-type, or slower dissociation rates to PD-L1, see [link to FACS plot]. Figure 34 Example of gating in the study. Hits were collected in selective growth medium and cultured at 30°C to high density, then aliquoted into cryopreservation solution at -80°C. For conjugates with faster dissociation rates or wild-type levels, the process was repeated twice; for mutants with slower dissociation rates, the process was repeated four times.

[0535] S2.4: Generation and screening of the target 21 mutant libraries

[0536] Twenty-one gene blocks encoding 21 heavy chain mutants were purchased from IDT. After PCR amplification and overlap extension PCR (Q5 hot-start DNA polymerase), 21 full-length atezolizumab scFv mutant genes were generated. The sequences of the 21 heavy chain mutants are shown in Table 12. These genes (along with linearized pYD1) were electroporated into EBY100 yeast to assemble the target library, and passaged / frozen according to the initial library preparation method. After titration to selective medium plates, transformants >10 were found. 7 .

[0537] To screen for the most improved (slowest dissociation rate) mutants from the library, progressively longer competition times were used in the presence of excess unlabeled PD-L1 (initially saturating the WT or mutant library with 200 nM biotin-PD-L1, followed by cell precipitation / washing to remove unbound biotin-PD-L1, and resuspending in excess (128 nM unlabeled PD-L1) at room temperature for 7.5 h, 17 h, and then 39 h). The most vigorous clones were collected in growth medium by FACS (using WT scFv as a reference) and propagated to the next round of sorting. After three rounds of selection, plasmids from the hits were harvested via yeast miniprep and transformed into GC10 competent cells, single colonies were selected, and bacterial miniprep was performed. Plasmids were sequenced using Eton Bioscience sequencing services.

[0538] S2.5: WT / mutant 4 / dvorumab / avirumab on the surface of yeast K D and k off Determination of value

[0539] The room temperature Koff (dissociation rate) assay for the scFv-PD-L1 complex was performed as described previously, using biotinylated / non-biotinylated PD-L1 (the same type used in library screening), TBS-BSA as buffer, and the same antibody reagents as in library screening. Yeast was cultured / induced in the same manner as in the screening protocol and treated with 200 nM b-PD-L1 at 4 × 10⁻⁶ ppm. 7 Label cells / mL, precipitate / wash, then add 64 nM unlabeled PD-L1 (10⁻⁶ cells / mL). 7Cells / mL) were collected and rotated at room temperature. The mean fluorescence intensity (MFI) of V5-positive yeast was recorded at each time point, and the data were fitted to a single-phase decay model using GraphPad Prism 9.3.0 software (Y = (Y0 – plateau value) * exp(-K * X) + plateau value, where Y is the fraction of yeast bound to biotin-PDL1, X is time, Y0 = 1 (fraction of yeast bound to biotin-PDL1 at the 0s competitive time point), the plateau value is a constant based on non-specifically bound yeast, and K is Koff). K on the yeast surface... D The MFI of the V5+ population was determined using the previously described method, whereby it was plotted relative to antigen concentration and biotin-PD-L1, the same method used for library screening. D Important considerations for the measurement include: using 10 for each data point 5 Cells were resuspended in 100 μL of cold TBS-BSA for 1 hour for each data point, then incubated with 1:200 anti-V5 antibody in V5 mAb. Cells were washed with cold TBS-BSA, then labeled with 1:100 SA-PE / IgG-AF647 on ice for 20 minutes. Cells were washed and resuspended in cold TBS-BSA for flow cytometry. The data were fitted to 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 is subtracted from all values). Because the dissociation rate of mutant 4 is very slow, the sample needs to be incubated at room temperature for 6 days to allow the system to reach equilibrium.

[0540] S2.6: Cloning and purification of WT / mutant 4 atezolizumab scFv

[0541] The gene encoding the scFv sequence was amplified from the pYD1 vector by PCR using Q5 hot-start DNA polymerase (NEB) and primers containing NcoI-HF and XhoI restriction sites (scFv F / R, Table 11). The PCR product was double-digested (along with the pET27b(+) vector (69863-3, MilliporeSigma)) (the vector was also dephosphorylated using Quick CIP (NEB)). The digested products were ligated using T4 DNA ligase (NEB), and the ligated products were transformed into NEB 5-α cells. The plasmids were purified and sequenced using the Zippy plasmid mini-extraction kit.

[0542] The scFv plasmid was transformed into Rosetta (DE3) cells (70954-3, Millipore Sigma) on LB-0.5% glucose plates containing 30 μg / mL kanamycin and 34 μg / mL chloramphenicol, and cultured with shaking at 37°C for 5 hours in 10 mL LB-0.5% glucose medium supplemented with kanamycin / chloramphenicol (30 μg / mL kanamycin and 34 μg / mL chloramphenicol) until OD was reached. 600 =0.8-1. Then, 4 mL of culture medium was added to 400 mL of self-induction medium (71300-4, Millipore Sigma, containing 100 μg / mL kanamycin and 34 μg / mL chloramphenicol), and cultured at 25°C with shaking at 200 RPM for 24 hours, and antibiotics were added. Subsequently, the scFv-His6 fusion protein was purified essentially according to the method described above, but using 5 mM MgSO4 instead of MgSO2. The supernatant was mixed with 0.3 M NaCl (not 0.5 M), and eluted by centrifugation using 5 mL Ni-NTA agarose (Qiagen 30210) filtered through a Vivaspin-20 (10,000 MWCO, Cytiva 28932360) filter. The protein was then exchanged for PBS buffer (140 mM NaCl, 2.7 mM KCl, 10 mM Na2HPO4, 1.8 mM KH2PO4, pH 7.4). Furthermore, the Ni-NTA agarose resin was washed with 50 mM HEPES, 300 mM NaCl, and 50 mM imidazole (pH 7.5), followed by washing with the same buffer (containing 70 mM imidazole), and finally eluted with the same buffer (containing 250 mM imidazole). The concentrated scFv was diluted to a glycerol content of 10%, then flash-frozen and stored at -80°C. The chromatogram was determined by SDS-PAGE. Figure 28A and 28B Protein purity was assessed, and concentration was determined by BCA analysis.

[0543] Statistics and repeatability

[0544] No statistical methods were used to predetermine the sample size. When processing the raw sequence data, unreliable reads (containing one or more bases and a phred quality score <10, or reads with mismatched overlapping regions) were discarded before further analysis or processing. These steps were designed to ensure that only reliable reads were used for analysis. No data was excluded from any subsequent analysis or model training. The test set used to evaluate 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 numpy library version 1.19.5 with a seed value of 0. When evaluating model performance using cross-validation, the cross-validation splits were generated by randomly dividing the dataset into 5 equal-sized partitions using the KFold function in Python scikit-learn library version 0.24.2.

[0545] The final evaluation of the model performance was conducted using a "blind" method, which generates predictions for sequences not included in the current dataset and validates these predictions according to the aforementioned experimental methods.

[0546] Data availability

[0547] The construction of the cAbRep database is described by Guo et al. (Front. Immunol. 10, 2365, 2019). The source data is provided with this article.

[0548] Code availability

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

[0550] result

[0551] RESP Process Overview

[0552] The RESP workflow comprises four key components. First, we developed a simple, novel encoding scheme in which an autoencoder was designed to learn a characterization capable of distinguishing human B-cell receptor (BCR) sequences from closely related sequences. This learned characterization is able to model trends in fluorescence-activated cell sorting (FACS) data more efficiently and accurately compared to characterizations learned from other popular protein and antibody sequences (see Table 9). This part of the workflow is applicable to any antibody sequence and can be reused in any project.

[0553] Table 9 provides a comparison of 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.

[0554] Table 9. Comparison of classification performance of different encoding types and model architectures on the atezolizumab dataset

[0555]

[0556] Secondly, yeast surface display libraries of mutants were constructed against the starting antibody sequence (e.g., atezolizumab here, where residues in the antibody heavy chain are randomly mutated). The mutant libraries were incubated with labeled target antigens (e.g., PD-L1 here), and mutants with slower dissociation kinetics (slower dissociation rate, i.e., potentially tighter binding) were screened by incubation in the presence of excess unlabeled target antigen. Flow cytometry experiments collected mutants with lower, intermediate, and faster dissociation rates compared to the WT antibody, thus dividing the population into a series of sequencing groups (…). Figure 15 And Figure 16A 1-2 and 16B 1-4 Then, sequences are determined by sequencing, and their group identities are distinguished using encoding schemes.

[0557] Third, by constructing a variational Bayesian neural network, experimental data were modeled using ordinal regression to predict the dissociation rate of sequences. This affinity model takes encoding sequences from a yeast surface display library as input and attempts to predict the probability that a given sequence has a slow dissociation rate (i.e., is a tight binding). The variational Bayesian architecture provides powerful regularization capabilities, minimizing the risk of overfitting and estimating the model's uncertainty for each prediction. Using ordinal regression, each sequence is mapped to a one-dimensional latent score, which indicates the model's confidence in the ordination group to which the sequence belongs. This approach naturally considers the ordination order of the groups, and as previously discussed in protein engineering problems... 41 As demonstrated in [the study], it is a useful method. Affinity models should be retrained for new antigens or different starting antibody sequences while their architecture remains unchanged.

[0558] Furthermore, the simulated annealing algorithm was modified to develop a computer-directed evolutionary algorithm. This algorithm utilizes the first three parts of the process to effectively explore the sequence space surrounding the training set, as well as sequences that may exhibit dissociation rates significantly lower than those of the identified parent sequences. The final analysis process is as follows: Figure 17 As shown. This component provides an efficient method for exploring the sequence space, which can be performed using a Bayesian neural network trained on any target antibody-antigen pair.

[0559] This study employed the RESP process to enhance the binding affinity of atezolizumab, an antibody targeting PD-L1. Twenty-one candidate sequences with predicted dissociation rates slower than wild-type sequences were selected. One of these sequences exhibited a dissociation rate approximately 10 times slower than that of wild-type atezolizumab scFv. D It is 17 times more potent (closer) than wild-type atezolizumab scFv.

[0560] Furthermore, we evaluated the core computational components of RESP based on data obtained by Mason et al. (Nat. Biomed. Eng. 5(6), 600–612, 2021). The results show that the model described in this paper has the same accuracy as the CNN model in the original study, while providing more useful uncertainty information. The individual components of the process will be discussed in more detail below.

[0561] Encode antibody sequences using a new autoencoder model.

[0562] First, a novel autoencoder model was developed, trained to characterize antibody sequences and distinguish real antibody sequences from other closely related sequences. The requirement to distinguish antibodies from closely related sequences may force the encoder to embed typical positional bias information into its learned representation. To this end, a convolutional autoencoder consisting of three main modules or components was constructed, such as... Figure 18 As shown, its performance was compared with one-hot encoding, UniRep, ESM-1b, ProtVec, AntiBertY, and AbLang. The components of this process are universal and do not require retraining for neoantigens or wild types.

[0563] The training set for the autoencoder was extracted from the cAb-Rep database, containing 121 sequenced B-cell receptor sequences from human donors, representing real antibody sequences. Specifically, the high-depth sequence library dataset was filtered using the ANARCI package to remove incomplete sequences, and the remaining sequences were numbered using the Chothia numbering scheme, resulting in a dataset containing 2,725,492 sequences. This dataset was then augmented with the same number of decoy sequences, generated by copying each real antibody sequence and randomly modifying it at seven positions. This number was chosen to strike a balance: increasing the number of mutations would reduce the likelihood of any mutant matching a real human B-cell receptor sequence in nature. However, too many mutations would make the classifier easily distinguish decoy sequences from the rest, causing the encoder unit to be less forced to learn information-rich input representations. Experiments with different numbers of mutations showed that seven mutations provided a good balance.

[0564] To ensure that the bait sequences truly differed from typical human antibody sequences, we randomly selected approximately 50,000 original sequence samples and used three different models—BioPhi, AbLSTM, and ANARCI—for scoring and bait processing (see [link to relevant documentation]). Figures 19A-19C Using any of these three models / tools, the decoy sequences exhibited significantly different humanization scores compared to the original sequences (two-sided 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 indicated that the p-values ​​were too small to be reliably calculated due to floating-point errors in double-precision operations and a return value ≈ 0.0. Therefore, these decoy sequences were explicitly rated by the model as having low humanization, indicating their suitability for use as decoys.

[0565] The autoencoder was trained on the full one-hot encoded cAb-Rep decoy dataset, with 200,000 sequences reserved for the test set, for eight training epochs, at which point the training loss converged. Reconstruction accuracy on the test set was >99.99%, while prediction accuracy for the B-cell receptor versus mutant task was 97.4%. These results demonstrate that the autoencoder can compress input sequences while retaining all the information needed for reconstruction or classification. The "Atezolizumab Modeling" section below will evaluate the performance of the autoencoder and use it as input to a model predicting binding affinity.

[0566] Generation and screening of atezolizumab scFv heavy chain libraries with improved dissociation rates

[0567] To develop antibodies against a target antigen using the RESP process, specific training data for that target antigen must first be generated so that the model can be trained to identify the types of sequences that bind to it. Here, the test example is the antigen PD-L1, starting with the single-chain variable region (scFv) heavy chain of the anti-PD-L1 antibody atezolizumab, hereinafter referred to as wild-type (WT).

[0568] First, the binding of wild-type atezolizumab (WT atezolizumab) scFv to human PD-L1 was tested, and its functionality was then assessed in yeast display mode. Figures 24A-24BThe results showed that it had good binding stability with the 6.1 nM antigen. The variable heavy chain (VH) domain of atezolizumab was chosen for mutagenesis because the Fab structure of atezolizumab binding to PD-L1 indicated that the heavy chain was primarily involved in the binding interaction (PDB numbers 5XXY, 5X8L). Furthermore, only the heavy chain was mutated, which is beneficial for deep sequencing because the region required for MiSeq sequencing is relatively short (approximately 316 bp in this case). Error-prone PCR was used to randomly generate gene mutations in the heavy chain region, which were then transformed into EBY100 yeast to produce a library of up to 78 million atezolizumab scFv variants. Due to the very strong / slow binding interaction and dissociation rate between WT scFv and PD-L1 (K0), the mutagenesis was not performed. D =1.75nM, k off =1.56×10 -4 s -1 Under these conditions, yeast exhibits a better fit for reducing the dissociation rate than for directly improving K. D Because maintaining a high molar excess of antigen relative to scFv on the yeast surface requires a large volume, and when the dissociation rate is already very slow, screening for improved K... D A very long incubation time is required. WT scFv is used as a control during the sorting process, sorted from FACS (Figure 16A). 1-2 Variants with faster WT levels and slower dissociation rates were isolated from the sample, and their plasmids were harvested and subjected to multiple rounds of PCR for MiSeq Nano PE250 barcoding of each group. Notably, selecting clones with slower dissociation rates appeared to split into two populations, one stronger than WT and the other weaker than WT (Figure 16B). 1-4 (“HB” sub-image).

[0569] The atezolizumab dataset contains mutant sequences from three different dissociation rate decreasing bins: RH01, RH02, and RH03. RH01 contains weak binders with a relatively fast dissociation rate; RH02 contains moderate binders with a dissociation rate similar to the wild type; and RH03 contains strong binders with a relatively slow dissociation rate. After quality screening, a total of 92,553 unique sequences were identified (out of a total of 550,215 sequences, as a sequence may appear multiple times in one bin or multiple times in different bins), of which 15,004 sequences were discarded because they contained mutations outside the randomized target protein region. RH03, RH02, and RH01 contain 15,070, 26,122, and 34,439 unique sequences, respectively (Table 10).

[0570] Table 10: Read counts by sorting category

[0571]

[0572]

[0573] Each sequence was assigned to the category with the highest frequency and given a weight equal to its frequency in the assigned category plus one, divided by the total frequency plus three. If a sequence had equal frequencies in both categories, it was discarded because it could not be explicitly assigned. After this filtering process, 75,631 unique sequences were retained. The weight of each sequence corresponds to its posterior probability of belonging to a given category, based on multinomial likelihood and uniform Dirichlet prior. This weight is crucial because it allows the model to distinguish which sequences can reliably be classified into a single category and which cannot. For example, a sequence that appears 5 times in both RH02 and RH03 is significantly less likely to be a strong binding than a sequence that appears 5 times only in RH03.

[0574] Atezolizumab modeling

[0575] The next component of this process is the affinity model, which takes a representation of a candidate sequence as input and predicts the ordination category (RH01, RH02, or RH03) to which the sequence should belong. This component requires experimental data for the target antigen, and if a new target antigen is selected, it needs to be retrained using the new experimental data.

[0576] Unlike most protein-engineered classifiers, affinity models are trained for ordinal regression, and the last layer of the network outputs a latent score. Traditional classification models treat categories as nominal (i.e., they have no specific order), so the model cannot learn any way to rank the categories. This is obviously not suitable for ordinally ranking data. Ordinal regression solves this problem by ranking the categories and using latent scores to determine which category a sequence should belong to. This approach provides a direct means of ranking sequences and selecting them for experimental evaluation. Past experiments have shown that for protein data with more than 3 binding categories, ordinal regression performs better in correctly predicting which sequences will appear in more rigorous future rankings, especially the sequence ranking method proposed by Liu et al. (Bioinformatics 36(7), 2126–2133, 2020). Since there are only 3 binding categories, these experiments are difficult to reproduce here, but ordinal regression is used here in light of past results.

[0577] The output score of the affinity model quantifies the model's belief that the sequence is likely to be a strong binder; a higher score indicates that the model is more certain the sequence should belong to a more stringent category. Next, this score is added to M-1 learned thresholds for M categories, and then a sigmoid function is applied to generate an output vector of M-1 probability values. Each element i in this output vector is a probability assigned by the model, indicating that the input sequence belongs to a more stringent binding category than i.

[0578] In short, the affinity model described in this paper employs an architecture similar to the backpropagation Bayesian algorithm, except that it is suitable for performing ordinal regression. Traditional fully connected neural networks learn a specific value for each weight and bias term in each hidden layer during training. In contrast, Bayesian neural networks treat each weight as a (typically Gaussian) distribution and learn the parameters of each weight distribution during training. Instead of using maximum likelihood for fitting as is common in neural networks, the Bayesian architecture approximates the posterior probability distribution using variational methods. By repeatedly sampling the weight distribution, uncertainty in the prediction can be estimated—or, if desired, the average of the Gaussian distribution of each weight can be used to generate point estimates of the prediction. In addition to generating uncertainty estimates, Bayesian neural networks also strongly regularize the model parameters by penalizing the bias from the prior.

[0579] The uncertainty information provided by the model can obviously be used in two ways. Predictions with high correlation uncertainty correspond to sequences located in relatively unexplored regions of the input space. In practice, considering cost and time, only predictions with relatively low correlation uncertainty should be selected for experimental evaluation, which is the strategy adopted in this paper.

[0580] The results of a standard fully connected neural network trained on simple classification were compared with those of a Bayesian network trained on ordinal regression across encoding types (Table 9). It is important to note that both neural networks have the same number of hidden layers and weights. While the non-Bayesian network shows a slight performance improvement, the Bayesian network provides additional information about the uncertainty associated with a given prediction, which is crucial for analyzing ordinal data as it allows us to measure the relative reliability of the model's predictions for candidate sequences.

[0581] Affinity models can use one-hot encoded sequences as input or other representations. To determine whether the representations generated by the autoencoder could be used for encoding affinity models, we encoded atezolizumab sequences using the following methods: 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 datasets were split into training (80%) and test (20%) sets, and 5-fold cross-validation was performed on the training set for each encoding type using a Bayesian neural network trained with ordinal regression (Table 9).

[0582] While different encoding methods may have advantages for specific tasks, autoencoders are the only characterization method that can rival or even outperform one-hot encoding for any model type (Table 9). Compared to one-hot encoding, it offers comparable performance despite significantly reduced model size and computational cost. It's important to note that autoencoders have approximately 20,000 parameters, while the FAIR-ESM model has 750 million parameters and the AntiBERTy model has 26 million parameters.

[0583] Therefore, autoencoders offer at least two specific benefits. First, they significantly reduce model size and complexity compared to simple one-hot encoding, and are much less expensive than language models due to their smaller number of parameters. Second, they provide a real-valued representation for each sequence, which can be used to cluster sequences or determine their similarity in later stages of the process. Given these advantages, autoencoders are preferred over one-hot encoding despite their similar performance.

[0584] Computer-directed evolution is used to select the most promising candidates.

[0585] The atezolizumab dataset was encoded using an autoencoder, and then a Bayesian neural network (affinity model) was trained for 30 training epochs on the dissociation rate based on the complete dataset. At this point, the workflow used the trained affinity model to search the sequence space surrounding the training set to generate new candidate antibody sequences not present in the training set.

[0586] At this stage, it's necessary to reduce the size of the search space. If a 118-amino acid sequence mutates in a computer simulation, the search space that needs to be covered would be incredibly large. Therefore, the focus is on a subset of smaller positions that significantly contribute to binding affinity—those positions that mutate most frequently in the highest-scoring sequences (i.e., the sequences the model considers most promising). For this dataset, we found that using the top 500, 1000, 2000, or 4000 sequences with the highest scores to select the top 10 most frequently mutated positions retrieves the same set of positions used in subsequent search steps. This step in the process is flexible and can use other criteria to select positions (e.g., only selecting positions with frequent mutations present in the CDR).

[0587] In the subsequent search steps, the model randomly selects a position in the wild-type sequence in each iteration and randomly mutates it with any amino acid. It's important to note that all mutants generated through this process have the same length as the wild-type sequence. The probability of selecting any specific amino acid at a given position is equal to the frequency of that amino acid in the training set plus one, divided by the total number of sequences in the training set plus 20. New candidate sequences are encoded using an encoder and scored using a trained affinity model. This score is evaluated using the classic simulated annealing criterion, and the probability of accepting a candidate sequence is determined by its score, the score of the previously accepted candidate sequence, and the temperature. The temperature is initially set to a high value to explore a larger sequence space and decreases in each iteration, thus reducing the probability of accepting candidate sequences that do not improve their scores.

[0588] Since the average sequence in the training set contains only 7 mutations, while the average high-scoring sequence contains even fewer (3 mutations on average), the most common amino acids at most positions are those present in the wild type. Therefore, the algorithm tends to sample a large number of sequences similar to the wild type. However, given the algorithm's randomness, it can explore combinations not present in the training set, thus finding new high-scoring sequences. Figure 20A 1-2 , 20B And 20C shows how the best score obtained so far changes during the optimization process.

[0589] Ten (10) simulated annealing chains were run; all chains converged within 1000 iterations. Candidate sequences with scores above the 90th percentile were selected from each chain, and duplicate sequences were removed. The uncertainty of the scores was evaluated by sampling each sequence 1000 times from the Bayesian neural network. The top 50% of sequences with the largest specified score standard deviations were removed. The 50% threshold was arbitrary, and this process produced a total of 127 sequences.

[0590] Then, median hierarchical clustering was used to cluster these sequences (the generated dendrogram is shown below). Figure 33 (As shown). The results clearly show that the selected sequences can be divided into two main subgroups. Marginal distribution (see Figure 20A) 1-2 The results (20B and 20C) indicate that the two main subgroups differ primarily at position 77, but also include subgroups with other interesting differences. For example, most of the selected sequences show the R98C mutation, but a subset of cluster 1 shows no change at this position compared to the wild type.

[0591] Cutting the tree at a lower height produces 11 subclusters, and the two highest-scoring sequences from each cluster are selected (or one of them if there is only one), resulting in 21 candidate sequences. This threshold is chosen to obtain a manageable number of sequences for experimental evaluation.

[0592] Mutation examination of 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 ).

[0593] Table 12. Amino acid sequences of 21 selected mutants for further evaluation.

[0594]

[0595]

[0596]

[0597] The amino acid sequences of the full-length IgG scFv sequences of the twenty-one mutants are shown below:

[0598] Mutant 1

[0599] DIQMTQSPSSSLSASVGDRVTITCRASQDVSTAVAWYQQKPGKAPKLLIYSASFLYSGVPSRFSGSGSGTDFTLTISSLQPEDFATYYCQQYLYHPATFGQGTKVEIKRGGGGSGGGGSGGGGSGG GGSEVQLVESGGGLVQPGGSLRLSCAASGFTFSDSWIHWVRQAPGQGLEWVAWISPYGGSTYYADSVKGRFTISADTSKNTTYLQMNSLRAEDTAVYYCVCRHWPGGFDYWGQGTLVTVSS(SEQ ID NO:22)

[0600] Mutant 2

[0601] DIQMTQSPSSLSASVGDRVTITCRASQDVSTAVAWYQQKPGKAPKLLIYSASFLYSGVPSRFSGSGSGTDFTLTISSLQPEDFATYYCQQYLYHPATFGQGTKVEIKRGGGGSGGGGSGGGGSGGGGSEVQLVESGGGLVQPGGSLRLSCAASGFTFSDSWIHWVRQAPGQGLEWVAWISPYGGSTYYADSVKGRFTASADTSKNTTYLQMNSLRAEDTAVYYCVCRHWPGGFDYWGQGTLVTVSS(SEQ ID NO:23)

[0602] Mutant 3

[0603] DIQMTQSPSSLSASVGDRVTITCRASQDVSTAVAWYQQKPGKAPKLLIYSASFLYSGVPSRFSGSGSGTDFTLTISSLQPEDFATYYCQQYLYHPATFGQGTKVEIKRGGGGSGGGGSGGGGSGGGGSEVQLVESGGGLVQPGGSLRLSCAASGFTFSDSWIHWVRQAPGRGLEWVAWISPYGGSTYYADSVKGRFTTSADTSKNTTYLQMNSLRAEDTAVYYCVRRHWPGGFDYWGQGTLVTVSS(SEQ ID NO:24)

[0604] Mutant 4

[0605] DIQMTQSPSSLSASVGDRVTITCRASQDVSTAVAWYQQKPGKAPKLLIYSASFLYSGVPSRFSGSGSGTDFTLTISSLQPEDFATYYCQQYLYHPATFGQGTKVEIKRGGGGSGGGGSGGGGSGGGGSEVQLVESGGGLVQPGGSLRLSCAASGFTFSDSWIHWVRQAPGKGLEWVAWISPYGGSTYYADSVKGRFTASADTSKNTTYLQMNSLRAEDTAVYYCVRRHWPGGFDYWGQGTLVTVSS(SEQ ID NO:25)

[0606] Mutant 5

[0607] DIQMTQSPSSLSASVGDRVTITCRASQDVSTAVAWYQQKPGKAPKLLIYSASFLYSGVPSRFSGSGSGTDFTLTISSLQPEDFATYYCQQYLYHPATFGQGTKVEIKRGGGGSGGGGSGGGGSGGGGSEVQLVESGGGLVQPGGSLRLSCAASGFTFSDSWIHWVRQAPGEGLEWVAWISPYGGSAYYADSVKGRFTTSADTSKNTTYLQMNSLRAEDTAVYYCVCRHWPGGFDYWGQGTLVTVSS(SEQ ID NO:26)

[0608] Mutant 6

[0609] DIQMTQSPSSLSASVGDRVTITCRASQDVSTAVAWYQQKPGKAPKLLIYSASFLYSGVPSRFSGSGSGTDFTLTISSLQPEDFATYYCQQYLYHPATFGQGTKVEIKRGGGGSGGGGSGGGGSGGGGSEVQLVESGGGLVQPGGSLRLSCAASGFTFSDSWIHWVRQAPGRGLEWVAWISPYGGSAYYADSVKGRFTTSADTSKNTTYLQMNSLRAEDTAVYYCVCRHWPGGFDYWGQGTLVTVSS(SEQ ID NO:27)

[0610] Mutant 7

[0611] DIQMTQSPSSLSASVGDRVTITCRASQDVSTAVAWYQQKPGKAPKLLIYSASFLYSGVPSRFSGSGSGTDFTLTISSLQPEDFATYYCQQYLYHPATFGQGTKVEIKRGGGGSGGGGSGGGGSGGGGSEVQLVESGGGLVQPGGSLRLSCAASGFTFSDSWIHWVRQAPGKGLEWVAWISPYGGSTYYADSVKGRFTASADTSKNTTYLQMNSLRAEDTAVYYCVCRHWPGGFDYWGQGTLVTVSS(SEQ ID NO:28)

[0612] Mutant 8

[0613] DIQMTQSPSSLSASVGDRVTITCRASQDVSTAVAWYQQKPGKAPKLLIYSASFLYSGVPSRFSGSGSGTDFTLTISSLQPEDFATYYCQQYLYHPATFGQGTKVEIKRGGGGSGGGGSGGGGSGGGGSEVQLVESGGGLVQPGGSLRLSCAASGFTFSDSWIHWVRQAPGKGLEWVAWISPYGGSAYYADSVKGRFTTSADTSKNTTYLQMNSLRAEDTAVYYCVCRHWPGGFDYWGQGTLVTVSS(SEQ ID NO:29)

[0614] Mutant 9

[0615] DIQMTQSPSSLSASVGDRVTITCRASQDVSTAVAWYQQKPGKAPKLLIYSASFLYSGVPSRFSGSGSGTDFTLTISSLQPEDFATYYCQQYLYHPATFGQGTKVEIKRGGGGSGGGGSGGGGSGGGGSEVQLVESGGGLVQPGGSLRLSCAASGFTFSDSWIHWVRQAPGKGLEWVAWISPYGGSAYYADSVKGRFTISADTSKSTVYLQMNSLRAEDTAVYYCVCRHWPGGFDYWGQGTLVTVSS(SEQ ID NO:30)

[0616] Mutant 10

[0617] DIQMTQSPSSLSASVGDRVTITCRASQDVSTAVAWYQQKPGKAPKLLIYSASFLYSGVPSRFSGSGSGTDFTLTISSLQPEDFATYYCQQYLYHPATFGQGTKVEIKRGGGGSGGGGSGGGGSGGGGSEVQLVESGGGLVQPGGSLRLSCAASGFTFSDSWIHWVRQAPGKGLEWVAWISPYGGSAYYADSVKGRFTISADTSKDTVYLQMNSLRAEDTAVYYCVCRHWPGGFDYWGQGTLVTVSS(SEQ ID NO:31)

[0618] Mutant 11

[0619] DIQMTQSPSSLSASVGDRVTITCRASQDVSTAVAWYQQKPGKAPKLLIYSASFLYSGVPSRFSGSGSGTDFTLTISSLQPEDFATYYCQQYLYHPATFGQGTKVEIKRGGGGSGGGGSGGGGSGGGGSEVQLVESGGGLVQPGGSLRLSCAASGFTFSDSWIHWVRQTPGRGLEWVAWISPYGGSAYYADSVKGRFTISADTSKDTTYLQMNSLRAEDTAVYYCACRHWPGGFDYWGQGTLVTVSS(SEQ ID NO:32)

[0620] Mutant 12

[0621] DIQMTQSPSSLSASVGDRVTITCRASQDVSTAVAWYQQKPGKAPKLLIYSASFLYSGVPSRFSGSGSGTDFTLTISSLQPEDFATYYCQQYLYHPATFGQGTKVEIKRGGGGSGGGGSGGGGSGGGGSEVQLVESGGGLVQPGGSLRLSCAASGFTFSDSWIHWVRQTPGEGLEWVAWISPYGGSAYYADSVKGRFTISADTSKDTTYLQMNSLRAEDTAVYYCTCRHWPGGFDYWGQGTLVTVSS(SEQ ID NO:33)

[0622] Mutant 13

[0623] DIQMTQSPSSLSASVGDRVTITCRASQDVSTAVAWYQQKPGKAPKLLIYSASFLYSGVPSRFSGSGSGTDFTLTISSLQPEDFATYYCQQYLYHPATFGQGTKVEIKRGGGGSGGGGSGGGGSGGGGSEVQLVESGGGLVQPGGSLRLSCAASGFTFSDSWIHWVRQTPGRGLEWVAWISPYGGSAYYADSVKGRFTASADTSKDTTYLQMNSLRAEDTAVYYCVCRHWPGGFDYWGQGTLVTVSS(SEQ ID NO:34)

[0624] Mutant 14

[0625] DIQMTQSPSSLSASVGDRVTITCRASQDVSTAVAWYQQKPGKAPKLLIYSASFLYSGVPSRFSGSGSGTDFTLTISSLQPEDFATYYCQQYLYHPATFGQGTKVEIKRGGGGSGGGGSGGGGSGGGGSEVQLVESGGGLVQPGGSLRLSCAASGFTFSDSWIHWVRQVPGRGLEWVAWISPYGGSAYYADSVKGRFTASADTSKDTTYLQMNSLRAEDTAVYYCVCRHWPGGFDYWGQGTLVTVSS(SEQ ID NO:35)

[0626] Mutant 15

[0627] DIQMTQSPSSLSASVGDRVTITCRASQDVSTAVAWYQQKPGKAPKLLIYSASFLYSGVPSRFSGSGSGTDFTLTISSLQPEDFATYYCQQYLYHPATFGQGTKVEIKRGGGGSGGGGSGGGGSGGGGSEVQLVESGGGLVQPGGSLRLSCAASGFTFSDSWIHWVRQVPGKGLEWVAWISPYGGSAYYADSVKGRFTASADTSKDTTYLQMNTLRAEDTAVYYCVCRHWPGGFDYWGQGTLVTVSS(SEQ ID NO:36)

[0628] Mutant 16

[0629] DIQMTQSPSSLSASVGDRVTITCRASQDVSTAVAWYQQKPGKAPKLLIYSASFLYSGVPSRFSGSGSGTDFTLTISSLQPEDFATYYCQQYLYHPATFGQGTKVEIKRGGGGSGGGGSGGGGSGGGGSEVQLVESGGGLVQPGGSLRLSCAASGFTFSDSWIHWVRQVPGKGLEWVAWISPYGGSAYYADSVKGRFTASADTSKDTIYLQMNSLRAEDTAVYYCVCRHWPGGFDYWGQGTLVTVSS(SEQ ID NO:37)

[0630] Mutant 17

[0631] DIQMTQSPSSLSASVGDRVTITCRASQDVSTAVAWYQQKPGKAPKLLIYSASFLYSGVPSRFSGSGSGTDFTLTISSLQPEDFATYYCQQYLYHPATFGQGTKVEIKRGGGGSGGGGSGGGGSGGGGSEVQLVESGGGLVQPGGSLRLSCAASGFTFSDSWIHWVRQTPGKGLEWVAWISPYGGSAYYADSVKGRFTTSADTSKDTTYLQMNSLRAEDTAVYYCVCRHWPGGFDYWGQGTLVTVSS(SEQ ID NO:38)

[0632] Mutant 18

[0633] DIQMTQSPSSLSASVGDRVTITCRASQDVSTAVAWYQQKPGKAPKLLIYSASFLYSGVPSRFSGSGSGTDFTLTISSLQPEDFATYYCQQYLYHPATFGQGTKVEIKRGGGGSGGGGSGGGGSGGGGSEVQLVESGGGLVQPGGSLRLSCAASGFTFSDSWIHWVRQVPGKGLEWVAWISPYGGSAYYADSVKGRFTTSADTSKDTTYLQMNSLRAEDTAVYYCVCRHWPGGFDYWGQGTLVTVSS(SEQ ID NO:39)

[0634] Mutant 19

[0635] DIQMTQSPSSLSASVGDRVTITCRASQDVSTAVAWYQQKPGKAPKLLIYSASFLYSGVPSRFSGSGSGTDFTLTISSLQPEDFATYYCQQYLYHPATFGQGTKVEIKRGGGGSGGGGSGGGGSGGGGSEVQLVESGGGLVQPGGSLRLSCAASGFTFSDSWIHWVRQAPGKGLEWVAWISPYGGSTYYADSVKGRFTVSADTSKDTTYLQMNSLRAEDTAVYYCVRRHWPGGFDYWGQGTLVTVSS(SEQ ID NO:40)

[0636] Mutant 20

[0637] DIQMTQSPSSSLSASVGDRVTITCRASQDVSTAVAWYQQKPGKAPKLLIYSASFLYSGVPSRFSGSGSGTDFTLTISSLQPEDFATYYCQQYLYHPATFGQGTKVEIKRGGGGSGGGGSGGGGSGG GGSEVQLVESGGGLVQPGGSLRLSCAASGFTFSDSWIHWVRQAPGKGLEWVAWISPYGGSTYYADSVKGRFTASADTSKDTTYLQMNSLRAEDTAVYYCVRRHWPGGFDYWGQGTLVTVSS(SEQ ID NO:41)

[0638] Mutant 21

[0639] DIQMTQSPSSSLSASVGDRVTITCRASQDVSTAVAWYQQKPGKAPKLLIYSASFLYSGVPSRFSGSGSGTDFTLTISSLQPEDFATYYCQQYLYHPATFGQGTKVEIKRGGGGSGGGGSGGGGSGG GGSEVQLVESGGGLVQPGGSLRLSCAASGFTFSDSWIHWVRQAPGKGLEWVAWISPYGGSAYYADSVKGRFTTSADTSKDTAYLQMNFLRAEDTAVYYCACRHWPGGFDYWGQGTLVTVSS(SEQ ID NO:42)

[0640] Analysis of the existing structures of atezolizumab Fab and PD-L1 (PDB numbers 5XXY and 5X8L) revealed that none of the mutated residues were in contact with PD-L1 in the WT sequence. Figure 21 This suggests that these mutations enhance binding affinity by altering the heavy chain conformation, rather than directly improving the binding contact. Another possibility is that these mutations stabilize the conformation, rather than altering it. For example, it has been shown that stabilizing the conformation of a particular antibody can be beneficial for enhancing affinity in certain situations.

[0641] To explore these possibilities, we used antigen-free IgFold technology to predict the structure of the highest-scoring mutants from twenty-one (21) sequences. Mutations in the high-scoring sequences altered the conformation of the protein’s two antigen-contact regions, particularly CDR H3. Figures 30A-30B and Figure 31The conformational changes of R98C (a mutation present in many of the highest-scoring sequences) and its surrounding residues are shown, and a distinct hydrophobic pocket is formed in the highest-scoring mutants compared to the wild type. It is important to note that the mutations described here do not directly contact the antigen. These predicted structures suggest that the mutations selected by this process function by altering the conformation of the contact region, rather than by directly forming new contacts with the antigen.

[0642] The trained affinity model was also used to calculate the importance of individual mutations: A) the contribution of a single mutation to binding when introduced into the WT alone without any other mutations; and B) the change in binding affinity when a single mutation is removed from a beneficial mutation combination. Full results are shown in Tables 13 and 14.

[0643] Table 13. Changes in affinity model allocation scores for each mutation found in the final 21 selected sequences if the mutation is introduced into the wild type without other modifications.

[0644] mutation Fractional changes (positive values ​​= more likely to form a strong combination) K

[43] Q 0.164 A

[79] T 2.941 A

[97] V 1.442 R

[98] C 5.013 I

[70] A 3.351 K

[43] R 0.013 I

[70] T 1.677 K

[43] E -0.050 T

[58] A 2.125 N

[77] S 0.074 A

[79] V 0.215 N

[77] D 0.839 A

[40] T -0.004 A

[97] T -0.047 A

[40] V -0.007 S

[85] T 0.026 A

[79] I 5.071 I

[70] V 0.142 S

[85] F -0.074

[0645] Table 14. Changes in affinity model assignment scores (positive = more likely to form a strong binder) if a specified mutation is removed from the selected sequence used for experimental evaluation.

[0646]

[0647]

[0648]

[0649]

[0650]

[0651]

[0652] The results do indeed show that some mutations may be more important than others, and also show that of the 19 introduced mutations, 12 are predicted to be beneficial even when they exist alone, i.e., increasing the likelihood that the related sequences become tightly bound sequences. The remaining mutations are of little harm or do not cause any change when they exist alone, although they are predicted to have beneficial effects in some cases. The model predicts that mutations such as A

[79] I or A

[79] T, R

[98] C, I

[70] A, A

[97] V and T

[58] A will have significant beneficial effects even when they exist alone.

[0653] Interestingly, the model is able to predict synergistic effects between individual mutations. For example, it predicts that if K43Q is introduced alone into the wild type, its beneficial effect is almost negligible; however, if it is removed from the mutants K

[43] Q, A

[79] T, A

[97] V, and R

[98] C, the negative impact on the score will be much greater. A

[97] T is minimally harmful when introduced alone, but it is predicted that A

[40] T, K

[43] E, T

[58] A, N

[77] D, A

[79] T, A

[97] T, and R

[98] C will benefit. For example, while it predicts that A

[79] T will be beneficial regardless of its introduction, the magnitude of its effect on the score will vary depending on the other mutations present, and may be up to twice as significant.

[0654] Validation of predicted tight-binding antibodies

[0655] The model generated 21 (21) mutants (all located on the heavy chain), which were purchased as gene blocks and fused with wild-type light chain sequences by PCR, then transformed into EBY100 yeast. Small libraries of the 21 mutants were screened in a manner similar to that of the natural libraries, but with a longer competition time (39 hours) used in the final screening to attempt to isolate the mutants with the slowest dissociation rates from the remaining 21 mutants. Figures 25A-25D The process failed to yield clearly isolated mutants. Sequencing of random clones after the most rigorous sorting revealed 12 independent mutants (with similar frequencies) among 17 random sequences. This indicates that a significant portion of the 21 mutant members had a significantly slower dissociation rate than the WT (WT exhibited significantly reduced fluorescence compared to the final, most rigorously sorted mutant library (Figures 25C-25D, bottom)). Upon examination of the original / unsorted 21 members after 8 hours of dissociation rate competition at room temperature, it appeared that virtually all 21 mutants selected by the model dissociated slower than WT (Figures 25A-25B, top subplots).

[0656] From the various mutants selected from the 21 member libraries after the most rigorous sorting, the I70A / A79T / A97V mutant (named mutant 4, sequence shown in Table 11) was selected for further characterization.

[0657] Table 11. WT atezolizumab, mutant 4, durvalumab and averumab scFv gene blocks, cloning primers, primers for generating mutant libraries and DNA libraries for next-generation sequencing

[0658]

[0659]

[0660]

[0661] The amino acid sequence of the full-length IgG heavy chain / constant region of mutant 4 is shown below:

[0662] MGWSCIILFLVATATGVHSEVQLVESGGGLVQPGGSLRLSCAASGFTFSDSWIHWVRQAPGKGLEWVAWISPYGGSTYYADSVKGRFTASADTSKNTTYLQMNSLRAEDTAVYYCVR RHWPGGFDYWGQGTLVTVSSASTKGPSVFPLAPSSKSTSGGTAALGCLVKDYFPEPVTVSWNSGALTSGVHTFPAVLQSSGLYSLSSVVTVPSSSLGTQTYICNVNHKPSNTKVDKKV EPKSCDKTHTCPPCPAPELLGGPSVFLFPPKPKDTLMISRTPEVTCVVVDVSHEDPEVKFNWYVDGVEVHNAKTKPREEQYASTYRVVSVLTVLHQDWLNGKEYKCKVSNKALPAPIE KTISKAKGQPREPQVYTLPPSREEMTKNQVSLTCLVKGFYPSDIAVEWESNGQPENNYKTTPPVLDSDGSFFLYSKLTVDKSRWQQGNVFSCSVMHEALHNHYTQKSLSLSPGK(SEQ ID NO:67)

[0663] The amino acid sequence of the full-length IgG light chain / constant region of mutant 4 is shown below:

[0664] MGWSCIILFLVATATGVHSDIQMTQSPSSSLSASVGDRVTITCRASQDVSTAVAWYQQKPGKAPKLLIYSASFLYSGVPSRFSGSGSGTDFTLTISSLQPEDFATYYCQQYLYHPATFG QGTKVEIKRTVAAPSVFIFPPSDEQLKSGTASVVCLLNNFYPREAKVQWKVDNALQSGNSQESVTEQDSKDSTYSLSSTLTLSKADYEKVYACEVTHQGLSSPVTKSFNRGEC(SEQ ID NO:68)

[0665] This mutant differs from other mutants that underwent rigorous sequencing; it lacks the introduction of cysteine ​​residues in its sequence, making it more attractive because its mutation does not form additional disulfide bonds. To characterize the improved dissociation / binding rate of the isolated mutant to PD-L1, the K-values ​​of WT and mutant 4scFv to PD-L1 were determined using a yeast display model. D and k off Value. First, the apparent dissociation rate of the yeast surface at room temperature was determined (Figure 22A), and it was found that mutant 4 had a higher k value. off It is about 10 times slower than WT atezolizumab (WT is 6.3×10). -5 s -1 The half-life is 3.04 hours; mutant 4 has a half-life of 6.5 × 10⁻⁶ hours. -6 s -1 (Half-life is 29.8 hours). Then, the k of mutant 4 was... off The values ​​were compared with two other FDA-approved anti-PD-L1 mAbs (dvorumab and avirumab, which are converted to scFv form on the yeast surface) (Figure 22B). The dissociation rate of mutant 4 was significantly slower than that of wild-type atezolizumab and duvorumab, and significantly slower than that of avirumab. The dissociation rate order of the FDA-approved mAb PD-L1 scFv (dvorumab < atezolizumab < avirumab) is consistent with previous results of these scFv values ​​obtained by surface plasmon resonance (SPR) binding assays.

[0666] Then, the binding affinity (Kb) of WT and mutant 4scFv was determined by titrating PD-L1 onto each scFv on the yeast surface. D K value. K value of mutant 4. D It was approximately 17 times denser than WT (on average) (5.3 pM vs. 92 pM, Figure 22C). It is noteworthy that the K+ on the yeast surface was determined at such a low antigen concentration. D It is difficult because maintaining a molar excess of yeast surface antigen to scFv at low antigen concentrations requires a large volume. 49 This also makes it difficult to repeat the data points. However, three independent K... D Measurements showed a significant (>10-fold) improvement in binding affinity (Figure 22C).

[0667] To assess whether the improved binding affinity and dissociation rate observed on the yeast surface could be replicated after purifying each scFv and reversing the direction of the binding assay (instead of fusing scFv to the yeast surface and binding it to soluble PD-L1, but immobilizing PD-L1 and binding it to soluble scFv), a BLI (biolayer interference) binding assay was performed. The assay revealed that mutant 4 exhibited a significantly slower dissociation rate than the wild-type, consistent with the results observed on the yeast surface. In this assay, the dissociation of the mutant from PD-L1 was minimal even after 2 hours at room temperature, while the wild-type showed more pronounced dissociation within approximately 60 minutes. However, BLI yielded two K values ​​for the wild-type. D The values ​​(219 pM for major species and 5.93 nM for minor species) were different from those of mutant 4. D The value is unmeasurable due to its very slow dissociation rate, making it impossible to fit / determine the kinetic constant value (see Figure 26A). 1-2 (and 26B). Minor species may be caused by a small fraction of unfolded scFv in the protein formulation. This observation prompted us to investigate whether scFv and PD-L1 are monomers or formed oligomers, so we detected them by mass spectrometry (a method that can image the distribution of protein oligomers in solution). Imaging analysis of PD-L1 (residue Met1-T239, expected molecular weight 35-38 kDa, molecular weight range determined by glycosylation) showed that at the tested concentrations (7.4 nM, 18.5 nM, and 74 nM), the protein was mainly monomeric ( Figure 27A-27G This is consistent with previous findings on human PD-L1, which identified PD-L1 as a monomer. These data support the K levels measured on the yeast surface. D The values ​​ranged between those of monomeric PD-L1 and scFv immobilized on the yeast surface. Imaging of WT or mutant 4scFv in solution showed that, at 100 nM for each scFv, both were mixed populations of monomers, dimers, trimers, and possibly higher molecular weight oligomers. WT and mutant 4scFv were purified twice (using the exact same method) and detected twice by mass spectrometry. WT scFv was a mixture of monomers and dimers in both preparation methods, while the mutant, for unknown reasons, was either predominantly monomeric (first batch) or a mixture of monomers / dimers / trimers (second batch). Figure 27A-27G Furthermore, the mutant scFv appeared to aggregate more readily than the WT scFv. Taken together, the oligomerization and aggregation of scFv in solution renders binding assays such as SPR (surface plasmon resonance) or BLI (where scFv is in solution and PD-L1 is immobilized) unsuitable for K-type assays. D Measurement is difficult because it is hard to know the concentration of soluble monomers in scFv.

[0668] Evaluation of the dataset by Mason et al. (Nat. Biomed. Eng., 5(6), 600-612, 2021)

[0669] To ensure the universality of RESP, it was then tested based on data from Mason et al. They started with the non-binding mutant of trastuzumab and successfully recovered an antibody with similar affinity to the original wild-type trastuzumab. Using the same training and test sets, an autoencoder was used to encode the sequence data, and a variational Bayesian network with the same data structure was trained to predict binding and non-binding states.

[0670] The model achieved the same accuracy as their counterparts on the test set (Matthews correlation coefficient 0.68, AUC-ROC curve 0.91) and correctly predicted sequences that were experimentally tested as binding compounds (see Figures 32A-32C). The model also assigned high scores to most of these sequences. For example, trastuzumab was assigned a higher score than 94% of the sequences in the training set, indicating a high confidence level in its role as a tight binding compound. The score distributions for the training set and experimentally evaluated sequences are shown in Figures 32A-32C.

[0671] However, unlike the process described by Mason et al., the process discussed in this embodiment provides additional uncertainty information not available in their model. The uncertainty of the affinity model predicting incorrectly on the test set is significantly higher than the uncertainty of predicting correctly on the test set (p<1e-23, two-sided Mann-Whitney U test). Therefore, uncertainty can be used to help determine whether the prediction is reliable, something that the process developed by Mason et al. cannot do, and this advantage is obtained without sacrificing accuracy.

[0672] discuss

[0673] This study provides an efficient workflow for antibody development, with the following innovations. First, a simple learning encoding for antibodies was developed. This encoding not only includes easily recoverable information from the original sequence but also contains additional encoded information describing key sequence features that distinguish human antibodies from the surrounding sequence space. Results show that when training models using antigen-specific experimental data to capture binding affinity trends, encoding the input sequence using the characterization generated by the autoencoder is more accurate than encoding the input sequence using popular advanced encodings such as UniRep, ESM-1b, AntiBertY, and AbLang. Notably, for this specific task, UniRep, ESM-1b, AntiBertY, and AbLang all outperformed one-hot encoding among all the models considered. 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 the performance of antibody affinity prediction compared to simple one-hot encoding.

[0674] Next, the training set was fitted using a model designed to provide intuitive sequence ordering and uncertainty quantification. The results show that the cross-class sequence read distribution can be used as data point weights, incorporated into the model fitting in a principled and intuitive way. The Bayesian ordinal regression model generates estimated posterior probabilities, thus providing additional information that traditional deep learning classifiers cannot offer, as the cross-class probability distributions predicted by traditional deep learning classifiers do not indicate the reliability of a given prediction.

[0675] Furthermore, an algorithm was designed to explore the sequence space covered by the training set. By evaluating the reliability of each prediction and limiting the search to the space covered by the training set, the time and cost wasted in evaluating poor candidate sequences were minimized.

[0676] Furthermore, the robust functionality of the disclosed process in this embodiment was validated. By training the model on a single large library, mutants with 10-fold / 17-fold increases in dissociation rate / binding affinity, respectively, were screened. This is consistent with previous results from CBX1 protein engineering, indicating that similar strategies produce binding affinity enhancements comparable to longer directed evolution processes. The mutant 4 identified and screened in this paper can be used as a scFv for cancer therapy, similar to previously reported high-affinity anti-PD-L1 protein (a mutant of PD-1 with a binding KD of 110 pM to PD-L1). Like PD-1, mutant 4 scFv is significantly smaller than the monoclonal antibody (30 kDa vs. 150 kDa), and therefore may be more effective in tumor penetration.

[0677] The procedure disclosed in this embodiment may have several important advantages over purely experimental methods. Phage and yeast display technologies only allow screening of a small number of clones, not individual clones, thus requiring additional experiments (e.g., ELISA or yeast K). D (Assays) are used to evaluate clones with the tightest binding affinity. These techniques can only identify strong binders present in the original library, thus often requiring very large libraries and / or multiple libraries to maximize sequence space coverage. Typically, the best binder identified through this process still does not possess sufficiently strong affinity, therefore, so-called affinity maturation is usually necessary to achieve the desired affinity. In this process, random mutations are introduced at selected sites, and the resulting target library is further screened and experimentally evaluated. It is important to note that antibodies produced through this process do not guarantee other desirable properties, such as solubility or stability.

[0678] Taking the initial discovery of atezolizumab as an example (as described in US Patent 8217149B2), firstly, phage display libraries screened for the PDL1 target underwent four rounds of panning, yielding 96 enriched clones. Then, two further libraries were constructed at two selected sites for affinity maturation to improve binding efficiency. These libraries were then used for a single plate sort, followed by five to six rounds of solution sorting. Finally, the enriched clones obtained from the final sorting were screened using high-throughput ELISA to identify the best candidates. It is worth noting that while this procedure is reliable, it is time-consuming and costly. Antibody engineering typically requires the construction of multiple libraries and the high-throughput ELISA detection of thousands of clones; for other examples of antibodies currently in clinical trials, see US Patents 20180086848A1 and 8313746B2, etc.

[0679] In contrast, the method disclosed in this paper only requires constructing a large library, followed by dissociation rate sorting and binning based on FACS to screen for tightly binding antibodies. This method does not require any high-throughput K... D Determination or subsequent focused library screening. Therefore, the significant time required to construct multiple libraries and the high throughput of ELISA screening / K D The time and cost required for assays are eliminated. Instead, the method requires only a few computational steps that are easy to implement and can be run quickly on a single computer equipped with a GPU. Bayesian neural networks need to be retrained each time a new antigen is selected based on data specific to that antigen (as described in this embodiment). In contrast, autoencoders can be reused without retraining. Furthermore, the method disclosed in this embodiment identifies tight bindings not present in the original library, unlike traditional phage and yeast display methods, which can only screen for sequences present in the library.

[0680] It is noteworthy that this method can be easily modified to introduce computer-assisted screening for stability, solubility, and other desired properties, unlike purely experimental techniques which cannot easily optimize these other properties simultaneously. Adding additional filters in the search step of the process is straightforward: rejecting candidates suggested by the improved simulated annealing algorithm if the predicted solubility or immunogenicity is poor. In contrast, purely experimental methods must optimize these properties individually, incurring considerable additional costs. The ability to achieve higher affinity while easily adding additional filters as needed is a key advantage of machine learning-assisted methods.

[0681] Recent literature describes other computer-aided antibody design strategies; however, the method disclosed in this paper offers several advantages. For example, Mason et al. achieved a 3-fold increase in antibody affinity compared to wild-type trastuzumab, while most of the mutants selected by their algorithm as promising actually showed weak affinity for the target, whereas the method disclosed in this paper increases affinity by 17-fold. Mason et al. performed multiple rounds of mutagenesis and library generation, including using rational design to select further mutagenesis sites based on data from previous selection steps. They required additional effort to limit the search space. In contrast, the improved simulated annealing algorithm disclosed in this paper enables computer simulation-based constraints. Unlike their proposed CNN, the Bayesian neural network disclosed in this paper provides an estimate of uncertainty related to prediction accuracy, thereby eliminating potentially unreliable predictions.

[0682] They describe a plausible design approach (as opposed to machine learning-based methods). However, the accuracy of their plausible design components is only comparable to the predictions of Rosetta modeling software, resulting in low prediction accuracy and an inability to determine which predictions are most reliable. For example, in an attempt to improve the binding of an antibody called G6 to its target VEGF, they used a computational search program that required approximately 250 CPU days, several orders of magnitude more expensive than the computational program used in this paper and published elsewhere. Of the mutants screened for experimental evaluation using this AbLift program, 60% were worse than the wild-type, and the vast majority of the remaining mutants were only comparable to the wild-type. Only one design in their approach truly improved K... D This improved it by five times, which is smaller than the improvement effect disclosed in this paper. It is worth noting that most or all of the 21 (21) mutants screened by the model disclosed in this paper had a slower dissociation rate than the wild type (and therefore the binding affinity may be higher).

[0683] Gaussian processes are used to model the trend of binding affinity as a function of the input sequence. Instead of experimentally validating their method, the authors used affinity predictions from the Absolut! software to determine whether a sequence is a strong binder; however, the Absolut! software itself is based on docking-generated structures and affinity, so its accuracy in tracking experimental data is unclear. Importantly, their method is affected by some well-known limitations of exact Gaussian processes. These models suffer from limitations in the number of training points (or O(N) in some more efficient modern implementations). 2 O(N) is shown on the surface. 3 The scaling factor is limited, making it completely infeasible for datasets larger than approximately 5,000-10,000 sequences. In contrast, the variational Bayesian network-based method presented in this paper is not subject to any such limitations.

[0684] While the method disclosed in this paper has been validated using atezolizumab as a starting point and PDL1 as a target, it is not target-specific and can therefore be easily adapted to other targets and challenges. The described workflow can be readily modified to optimize only a single CDR or a subset of available positions and to combine with other models to rank candidate sequences for desired properties beyond affinity. Therefore, we anticipate that this workflow and its modified versions could become useful tools for accelerating antibody discovery and development.

Claims

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

2. The method as described in claim 1, wherein, The directed evolution, sequencing, and machine learning analysis steps include encoding the directed evolution sequences using a trained autoencoder computer model.

3. The method as described in claim 2, wherein, The directed evolution, sequencing, and machine learning analysis steps also include using a trained variational Bayesian neural network to perform ordinal regression on the encoded directed evolution sequences to quantify the desired features of interest.

4. The method of claim 2, wherein, The directed evolution, sequencing, and machine learning analysis steps also include using an approximate Gaussian process (GP) computer model, which approximates the kernel of the encoded directed evolution sequence using random features for ordinal regression.

5. The method of claim 4, wherein, The GP computer model reduces the number of computations required when approximating the kernel using the random feature by utilizing the preconditional conjugate gradient method.

6. The method of claim 3, wherein, Computer mutagenesis was performed using a simulated annealing strategy to explore directed evolutionary sequences that were not present in the mutant library and to evaluate the properties of interest.

7. The method of claim 1, wherein, The proteins are selected from the group consisting of: therapeutic proteins, protein domains, enzymes, and proteins used to enhance CRISPR efficiency and enzyme activity.

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

9. The method of claim 7, wherein, The protein domains are used as imaging probes.

10. A method for designing and engineering an antibody comprising the desired mutation and having improved binding affinity or dissociation rate to a target antigen, comprising: a) Establish an autoencoder model based on the coding mutation sequence data of known antibodies that bind to the target antigen; b) Develop a Bayesian ordinal regression model that scores the likelihood of each mutated sequence forming a strong binding compound with the target antigen; and c) An improved simulated annealing algorithm, tested for binding affinity or dissociation rate to the target antigen, is used to select antibodies containing the desired mutations.

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 includes the scFv amino acid sequence as shown in SEQ ID NO:1-21.

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

25.

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

16. An antibody or a composition thereof comprising a mutated sequence of a known antibody against a target antigen, wherein the antibody is engineered by the method of any one of claims 1-14.

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

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

25.

19. The antibody as claimed in claim 17 or 18, wherein, The antibody contains a mutated sequence of atezolizumab.

20. The antibody according to any one of claims 16-19, wherein, The antibody is used to treat PD-L1 positive tumors.

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

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

Images