DatasetsSKEMPIWe used the SKEMPI v219 dataset for training and validation. The dataset contains 7085 ΔΔGbind measurements on 348 complexes. We performed preprocessing following13,18, discarding data with ambiguous ΔΔGbind values (e.g., the mutant is a nonbinder without exact KD measurements) or high ΔΔGbind variance across multiple measurements (>1 kcal/mol), and arrived at 5729 distinct mutations, with their ΔΔGbind measurements, on 340 complexes. See Supplementary Table S6 for a list of discarded, high-variance SKEMPI mutations. For each mutation, we sampled the mutant structure with FoldX 49 based on the wild-type crystal structure. We used PDBFixer v1.837 to fix the PDB structures beforehand if the raw structure could not be processed by FoldX. These FoldX-derived SKEMPI structures are used to train deep-learning models, including Bind-ddG, GearBind, and GearBind+P. The resulting dataset was split into five subsets with roughly the same size using a split-by-PDB strategy, in order to perform cross-validation.CATHFor pretraining, we use a non-redundant subset of CATH v4.3.0 domains, which contains 30,948 experimental protein structures with less than 40% sequence identity. We also remove proteins that exceed 2000 AAs in length for efficiency. During pretraining, we randomly truncate long sequences into 150 AAs for efficiency. It is important to note that, currently, our pretraining exclusively utilizes single-chain proteins. The information learned by single-chain pretraining can be transferred to downstream tasks on protein complexes and we have found that this approach alone is sufficient to yield improvement.HER2 bindersThe HER2 binders test set was collected from ref. 27. The raw data include SPR data for 758 binders and 1097 non-binders. As all benchmarked methods only support amino acid substitutions, we filter out the binders that have different lengths compared to the wild-type antibody (Trastuzumab), leaving 419 Trastuzumab mutants. ΔΔGbind values are calculated by \(\Delta \Delta {G}_{{{{\rm{bind}}}}}=-RT\ln \left({K}_{{{{\rm{D}}}}}^{({{{\rm{mt}}}})}/{K}_{{{{\rm{D}}}}}^{({{{\rm{wt}}}})}\right)\) based on the SPR-measured binding affinity. Note that we only use this dataset as a test set to evaluate physics-based models (FoldX, Flex-ddG) and deep-learning models (Bind-ddG, GearBind, and GearBind+P) trained on SKEMPI.GearBind implementationGiven a pair of wild-type and mutant structures, GearBind predicts the binding free energy change ΔΔGbind by building a geometric encoder on a multi-relational graph, which is further enhanced by self-supervised pretraining. Note that the key feature that makes the neural network geometric is that it considers the spatial relationship between entities, i.e., nodes in a graph. In the following sections, we will discuss the construction of multi-relational graphs, multi-level message passing, and pretraining methods.Constructing relational graphs for protein complex structuresGiven a protein–protein complex, we construct a multi-relational graph for its interface and discard all other atoms. Here a residue is considered on the interface if its Euclidean distance to the nearest residue from the binding partner is no more than 6 Å. Each atom on the interface is regarded as a node in the graph. We add three types of edges to represent different interactions between these atoms. For two atoms with a sequential distance lower than 3, we add a sequential edge between them, the type of which is determined by their relative position in the protein sequence. For two atoms with spatial distance lower than 5 Å, we add a radial edge between them. Besides, each atom is also linked to its 10-nearest neighbors to guarantee the connectivity of the graph. Spatial edges that connect two atoms adjacent in the protein sequence are not interesting and thus discarded. The relational graph is denoted as \(({{{\mathcal{V}}}},\, {{{\mathcal{E}}}},\, {{{\mathcal{R}}}})\) with \({{{\mathcal{V}}}},\, {{{\mathcal{E}}}},\, {{{\mathcal{R}}}}\) denoting the sets of nodes, edges, and relation types, respectively. We use the tuple (i, j, r) to denote the edge between atom i and j with type r. We use one-hot vectors of residue types and atom types as node features for each atom and further include sequential and spatial distances in edge features for each edge.Building geometric encoder by multi-level message passingOn top of the constructed interface graph, we now perform multi-level message passing to model interactions between connected atoms, edges and residues. We use \({{{{\bf{a}}}}}_{i}^{(l)}\) and \({{{{\bf{e}}}}}_{(j,i,r)}^{(l)}\) to denote the representations of node i and edge (j, i, r) at the l-th layer. Specially, we use \({{{{\bf{a}}}}}_{i}^{(0)}\) to denote the node feature for atom i and \({{{{\bf{e}}}}}_{(j,i,r)}^{(0)}\) to denote the edge feature for edge (j, i, r). Then, the representations are updated by the following procedures:$${{{{\bf{a}}}}}_{i}^{(l)}\leftarrow {{{\rm{AtomMP}}}}\left({{{{\bf{a}}}}}_{i}^{(l-1)}\right),$$
(1)
$${{{{\bf{e}}}}}_{(j,i,r)}^{(l)}\leftarrow {{{\rm{EdgeMP}}}}\left({{{{\bf{e}}}}}_{(j,i,r)}^{(l-1)}\right),$$
(2)
$${{{{\bf{a}}}}}_{i}^{(l)}\leftarrow {{{{\bf{a}}}}}_{i}^{(l)}+{{{\rm{AGGR}}}}\left({{{{\bf{e}}}}}_{(j,i,r)}^{(l)}\right),$$
(3)
$${{{{\bf{a}}}}}_{{{{\rm{C}}}}\alpha (i)}^{(l)}\leftarrow {{{{\bf{a}}}}}_{{{{\rm{C}}}}\alpha (i)}^{(l)}+{{{\rm{ResAttn}}}}\left({{{{\bf{a}}}}}_{{{{\rm{C}}}}\alpha (i)}^{(l)}\right).$$
(4)
First, we perform atom-level message passing (AtomMP) on the atom graph. Then, a line graph is constructed for the message passing between edges (EdgeMP) so as to learn effective representations between atom pairs. The edge representations are used to update atom representations via an aggregation function (AGGR). Finally, we take the representations \({{{{\bf{a}}}}}_{{{{\rm{C}}}}\alpha (i)}^{(l)}\) of the alpha carbon as residue representation and perform a residue-level attention mechanism (ResAttn), which can be seen as a special kind of message passing on a fully connected graph. In the following paragraphs, we will discuss these components in detail.Atom-level message passingFollowing GearNet22, we use a relational graph neural network (RGCN)38 to pass messages between atoms. In a message-passing step, each node aggregates messages from its neighbors to update its own representation. The message is computed as the output of a relation (edge type)-specific linear layer when applied to the neighbor representation. Formally, the message-passing step is defined as:$${{{\rm{AtomMP}}}}\left({{{{\bf{a}}}}}_{i}^{(l-1)}\right)={{{{\bf{a}}}}}_{i}^{(l-1)}+\sigma \left({{{\rm{BN}}}}\left({\sum}_{r\in {{{\mathcal{R}}}}}{{{{\bf{W}}}}}_{r}^{(a)}{\sum}_{(j,i,r)\in {{{\mathcal{E}}}}}{{{{\bf{a}}}}}_{j}^{(l-1)}\right)\right),$$where BN( ⋅ ) denotes batch norm and σ( ⋅ ) is the ReLU activation function.Edge-level message passing and aggregationModeling sequential proximity or spatial distance alone is not enough for capturing the complex protein–protein interactions (PPI) contributing to the binding. Multiple works have demonstrated the benefits of incorporating angular information using edge-level message passing20,22,39. Here we construct a line graph40, i.e., a relational graph among all edges of the above atom-level graph. Two edges are connected if and only if they share a common end node. The relations, or edge types, are defined as the angle between the atom-level edge pair, discretized into eight bins. We use \(({{{{\mathcal{V}}}}}^{{\prime} },\, {{{{\mathcal{E}}}}}^{{\prime} },\, {{{{\mathcal{R}}}}}^{{\prime} })\) to denote the constructed line graph. Then, relational message passing is used on the line graph:$${{{\rm{EdgeMP}}}}\left({{{{\bf{e}}}}}_{x}^{(l-1)}\right)=\sigma \left({{{\rm{BN}}}}\left({\sum}_{{r}^{{\prime} }\in {{{{\mathcal{R}}}}}^{{\prime} }}{{{{\bf{W}}}}}_{{r}^{{\prime} }}^{(e)}{\sum}_{(y,x,{r}^{{\prime} })\in {{{{\mathcal{E}}}}}^{{\prime} }}{{{{\bf{e}}}}}_{y}^{(l-1)}\right)\right),$$where x and y denote edge tuples in the original graph for abbreviation.Once we updated the edge representations, we aggregate them into its end nodes. These representations are fed into a linear layer and multiplied with the edge type-specific kernel matrix \({{{{\bf{W}}}}}_{r}^{(a)}\) in AtomMP:$${{{\rm{AGGR}}}}\left({{{{\bf{e}}}}}_{(j,i,r)}^{(l)}\right)=\sigma \left({{{\rm{BN}}}}\left({\sum}_{r\in {{{\mathcal{R}}}}}{{{{\bf{W}}}}}_{r}^{(a)}{\sum}_{(j,i,r)\in {{{\mathcal{E}}}}}{{{\rm{Linear}}}}\left({{{{\bf{e}}}}}_{(j,i,r)}^{(l)}\right)\right)\right),$$which will be used to update the representation for atom i as in Eq. (3).Residue-level message passingConstrained by the computational complexity, atom and edge-level message passing only consider sparse interactions while ignoring global interactions between all pairs of residues. By modeling a coarse-grained view of the interface at the residue level, we are able to perform message passing between all pairs of residues. To do this, we design a geometric graph attention mechanism, which takes the representations of the alpha carbon of residues as input and updates their representations with the output as in Eq. (4). Here we follow the typical definition of self-attention to calculate attention logits with query and key vectors and apply the probability on the value vectors to get residue representations ri:$${\alpha }_{ij}= {{{{\rm{Softmax}}}}}_{j}\left(\frac{1}{d}\cdot {{{{\rm{Linear}}}}}_{q}\left({{{{\bf{a}}}}}_{{{{\rm{C}}}}\alpha (i)}^{(l)}\right)\cdot {{{{\rm{Linear}}}}}_{k}\left({{{{\bf{a}}}}}_{{{{\rm{C}}}}\alpha (j)}^{(l)}\right)\right),\\ {{{{\boldsymbol{r}}}}}_{i}^{(l)}= {\sum}_{j}{\alpha }_{ij}\cdot {{{{\rm{Linear}}}}}_{v}\left({{{{\bf{a}}}}}_{{{{\rm{C}}}}\alpha (i)}^{(l)}\right),$$where d is the hidden dimension of the representation \({{{{\bf{a}}}}}_{{{{\rm{C}}}}\alpha (i)}^{(l)}\) and the Softmax function is taken over all j.Besides traditional self-attention, we also include geometric information in the attention mechanism, which should be invariant to roto-translational transformations on the global complex structure. Therefore, we construct a local frame for each residue with coordinates of its Nitrogen, Carbon and alpha Carbon atoms:$$\begin{array}{rcl}&&{{{{\bf{v}}}}}_{i1}={{{{\bf{x}}}}}_{{{{\rm{N}}}}(i)}-{{{{\bf{x}}}}}_{{{{\rm{C}}}}\alpha (i)},\\ &&{{{{\bf{v}}}}}_{i2}={{{{\bf{x}}}}}_{{{{\rm{C}}}}(i)}-{{{{\bf{x}}}}}_{{{{\rm{C}}}}\alpha (i)},\\ &&{{{{\bf{u}}}}}_{i1},{{{{\bf{u}}}}}_{i2}={{{\rm{GramSchmidt}}}}({{{{\bf{v}}}}}_{i1},{{{{\bf{v}}}}}_{i2})\\ &&{{{{\bf{R}}}}}_{i}=\left[\frac{{{{{\bf{u}}}}}_{i1}}{\parallel {{{{\bf{u}}}}}_{i1}\parallel },\frac{{{{{\bf{u}}}}}_{i2}}{\parallel {{{{\bf{u}}}}}_{i2}\parallel },\frac{{{{{\bf{u}}}}}_{i1}}{\parallel {{{{\bf{u}}}}}_{i1}\parallel }\times \frac{{{{{\bf{u}}}}}_{i2}}{\parallel {{{{\bf{u}}}}}_{i2}\parallel }\right],\end{array}$$where we use x to denote the coordinate of an atom and GramSchmidt( ⋅ ) refers to the Gram-Schmidt process for orthogonalization. Then, the geometric attention is designed to model the relative position of beta carbons of all residues j in the local frame of residue i:$${{{{\bf{p}}}}}_{i}^{(l)}={\sum}_{j}{\alpha }_{ij}{{{{\bf{R}}}}}_{i}^{\top }\left({{{{\bf{x}}}}}_{{{{\rm{C}}}}\alpha (i)}-{{{{\bf{x}}}}}_{{{{\rm{C}}}}\beta (j)}\right),$$where \({{{{\bf{p}}}}}_{i}^{(l)}\) is the spatial representations for the residue i. When the complex structure is rotated, the frame Ri and relative position xCα(i) − xCβ(j) are rotated accordingly and the effects will be counteracted, which guarantees the rotation invariance of our model.The final output is the concatenation of residue representations \({{{{\boldsymbol{r}}}}}_{i}^{(l)}\) and spatial representations \({{{{\bf{p}}}}}_{i}^{(l)}\):$${{{\rm{ResAttn}}}}\left({{{{\bf{a}}}}}_{{{{\rm{C}}}}\alpha (i)}^{(l)}\right)={{{\rm{Concat}}}}\left({{{{\boldsymbol{r}}}}}_{i}^{(l)},\, {{{{\bf{p}}}}}_{i}^{(l)}\right).$$After obtaining representations for each atom, we apply a mean pooling layer over representations of all alpha carbons aCα(i) to get protein representations h. An antisymmetric prediction head is then applied to guarantee that back mutations would have the exact opposite predicted ΔΔGbind values:$$\widetilde{\Delta \Delta {G}_{{{{\rm{bind}}}}}}={{{\rm{MLP}}}}\left.\right({{{{\bf{h}}}}}^{({{{\rm{wt}}}})},\, {{{{\bf{h}}}}}^{\left.\right({{{\rm{mt}}}}}\left.\right)-{{{\rm{MLP}}}}({{{{\bf{h}}}}}^{({{{\rm{mt}}}})},\, {{{{\bf{h}}}}}^{({{{\rm{wt}}}})}),$$
(5)
where h(wt) and h(mt) denote the representations for wild-type and mutant complexes and \(\widetilde{\Delta \Delta {G}_{{{{\rm{bind}}}}}}\) is the predicted ΔΔGbind from our GearBind model.Modeling energy landscape of proteins via noise contrastive estimationAs paired binding free energy change data is of relatively small size, it would be beneficial to pretrain GearBind with massive protein structural data. The high-level idea of our pretraining method is to model the distribution of native protein structures, which helps identify harmful mutations yielding unnatural structures. Denoting a protein structure as x, its distribution can be modeled with Boltzmann distribution as:$$p({{{\bf{x}}}};{{{\boldsymbol{\theta }}}})=\frac{\exp (-E({{{\bf{x}}}};{{{\boldsymbol{\theta }}}}))}{A({{{\boldsymbol{\theta }}}})},\, A({{{\boldsymbol{\theta }}}})=\int\exp (-E({{{\bf{x}}}};{{{\boldsymbol{\theta }}}}))d{{{\bf{x}}}},$$
(6)
where θ denotes learnable parameters in our encoder, E(x; θ) denotes the energy function for the protein x and A(θ) is the partition function to normalize the distribution. The energy function is predicted by applying a linear layer on the GearBind representations h(x) of protein x:$$E({{{\bf{x}}}};{{{\boldsymbol{\theta }}}})={{{\rm{Linear}}}}({{{\bf{h}}}}({{{\bf{x}}}})).$$
(7)
Given the observed dataset {x1, . . . , xT} from PDB, our objective is to maximize the probability of these samples:$${{{\rm{maximize}}}}\,\,\frac{1}{2T}{\sum}_{t}\log p({{{{\bf{x}}}}}_{t};{{{\boldsymbol{\theta }}}}).$$
(8)
However, direct optimization of this objective is intractable, since calculating the partition function requires integration over the whole protein structure space. To address this issue, we adopt a popular method for learning energy-based models called noise contrastive estimation24. For each observed structure xt, we sample a negative structure yt and then the problem can be transformed to a binary classification task, i.e., whether a sample is observed in the dataset or not.$${{{\rm{minimize}}}}\,\,\frac{1}{2T}{\sum}_{t}\log \left[\sigma (E({{{{\bf{x}}}}}_{t};{{{\boldsymbol{\theta }}}})-E({{{{\bf{y}}}}}_{t};{{{\boldsymbol{\theta }}}}))\right],$$
(9)
where σ( ⋅ ) denotes the sigmoid function for calculating the probability for a sample xt belonging to the positive class. We could see that the above training objective tries to push down the energy of the positive examples (i.e., the observed structures) while pushing up the energy of the negative samples (i.e., the mutant structures.For negative sampling, we perform random single-point mutations on the corresponding positive samples and then generate its conformation by keeping the backbone unchanged and sampling side-chain torsional angles at the mutation site from a backbone-dependent rotamer library25. Besides, to further enhance the model’s capability to distinguish structural noises, we randomly choose 30% residues to randomly rotate torsional angles when generating negative samples.After pretraining on the CATH database, we finetune the GearBind encoder on downstream tasks for prediction to avoid overfitting.Cross-validation on SKEMPIDuring cross-validation, a model is trained and tested five times, each time using a different subset as the test set and the remaining four subsets as the training set. Results are calculated for each test set, and their mean and standard error of mean are reported as the final cross-validation performance. During the process of cross-validation, each individual data point is incorporated into the test set precisely once. This ensures that a comprehensive “test result table” is compiled, which includes predictive values for each data point when it is part of the test set. Subsequent performance analyses are done by splitting this table by various criteria and evaluate performance on each subset.After cross-validation on SKEMPI, we obtain five sets of model parameters. During inference, we use the mean of the predicted values of these five checkpoints as the model prediction result.Baseline implementation detailsFoldXIn this work, we use FoldX 49 for mutant structure generation. First, each PDB file is processed with the RepairPDB command for structural corrections. Then, the wild-type, mutant structure pair is built using the BuildModel command. We use the AnalyseComplex command to get the FoldX ΔΔGbind prediction based on the wild-type and mutant structures.Flex-ddGWe run Flex-ddG with its official implementation at https://github.com/Kortemme-Lab/flex_ddG_tutorial. Each PDB file is processed with PDBFixer v1.837. Using the default parameters, we sample 35 structure models for each mutation, with the number of backrub trails set to 35,000 and the energy function set to fa_talaris2014. The final ΔΔGbind values are predicted with a generalized additive model that reweights the score terms.GearBind(+P)We implement GearBind with the TorchDrug library41. For message passing, we employed a 4-layer GearBind model with a hidden dimension of 128. Regarding edge message passing, the connections between edges are categorized into 8 bins according to the angles between them. To predict the ΔΔGbind value from graph representations, we utilized a 2-layer MLP.The model was trained using the Adam optimizer with a learning rate of 1e-4 and a batch size of 8. The training process is performed on 1 A100 GPU for 40 epochs. For pretraining, we use the same architecture with 4-layer GearBind model with a hidden dimension of 128. The pretraining was conducted using the Adam optimizer with a learning rate of 5e-4 and a batch size of 8, employing 4 A100 GPUs for 10 epochs.Bind-ddGTo ensure a fair comparison, we re-implement and re-train the Bind-ddG model on our SKEMPI splits. We follow the configuration of the original implementation to set the dimensions of hidden and pair representations at 128 and 64, respectively. Our validation performance indicates that the optimal configuration for our setup includes a two-layer geometric attention mechanism and a four-layer MLP predictor. We trained the model using an Adam optimizer with a learning rate of 1e-4 and a batch size of 8, on a single A400 GPU, for a total of 40 epochs.In silico affinity maturation of CR3022 and anti-5T4 UdAbPDB 6XC342, in which chains H and L comprise antibody CR3022 and chain C is the SARS-CoV-2 RBD, was chosen as the starting complex for CR3022 affinity maturation. To better simulate the CR3022 interaction with Omicron RBD, we constructed the complex structures for BA.4 and BA.1.1 mutants with SWISS-MODEL35. We then performed saturation mutagenesis on the CDRs of CR3022 and generated mutant structures using FoldX9 and Flex-ddG10. Specifically, residues 26–35, 50–66, 99–108 from the heavy chain H and residues 24–40, 56–62, 95–103 from the light chain L are mutated. This totals 1400 single-point mutations (if we count the self-mutations). We use our ensemble model to rank the mutations and select the top-ranked mutants for synthesis and subsequent experimental validation. Mutations are ranked by the modified z-score (where values are subtracted by the median rather than the mean to be less sensitive to outliers) averaged across multiple ΔΔGbind prediction methods.An unpublished complex structure was used to optimize anti-5T4 UdAb. As the single-domain antibody binding two distinct epiopes on 5T4 (Fig. 5b), anti-5T4 UdAb has a larger interface region compared to traditional antibodies. After analyzing its interface with 5T4, we decided to perform saturation mutagenesis on residues 1, 3, 25, 27–30, 31–33, 39–45, 52–57, 59, 91–93, 95, 99, 100–102, 103, 105, 110, 112, 115–117. This totals 780 single-point mutations (if we count the self-mutations) that go through the same ranking and selection strategies as described above.Antigen preparationThe gene encoding SARS-CoV RBD was synthesized by Genscript (Nanjing, China) and subcloned into pSectag 2B vector with C-terminal human IgG1 Fc fragment and AviTag. The recombinant vector was transfected to Expi 293 cells and cultured at 37 °C for 5 days, followed by centrifugation at 2200×g for 20 min. The supernatant was harvested and filtered through a 0.22-μm vacuum filter. The protein G resin (Genscript) was loaded into the column, washed by PBS, and flow the supernatant through to fully combine the resin. Then the targeted protein was eluted with 0.1 M glycine (pH 3.0) and neutralized with 1 M Tris-HCL (pH 9.0), followed by buffer-exchanged and concentrated with phosphate buffered saline (PBS) using an Amicon ultra centrifugal concentrator (Millipore) with a molecular weight cut-off of 3 kDa. Protein concentration was measured using the NanoDrop 2000 spectrophotometer (Thermo Fisher), and protein purity was examined by sodium dodecyl sulfate-polyacrylamide gel electrophoresis (SDS-PAGE). The Delta RBD protein was purchased from Vazyme (Nanjing, China) and Omicron S protein was purchased from ACROBiosystems (Beijing, China). The biotinylated human TPBG/5T4 and human TPBG/5T4-Fc antigen were purchased from ACROBiosystems (Beijing, China).Preparation for mutant and wild-type CR3022 antibodiesThe heavy chain and light chain genes of different CR3022 antibodies were synthesized and subcloned into expression vector pcDNA 3.4 in IgG1 format. These constructed vectors were transfected into CHO cells and purified by Protein A. All antibodies were produced by Biointron Biological Inc. (Shanghai, China).Generation of mutant and wild-type anti-5T4 UdAbsThe pComb3x vector encoding the gene of wild-type anti-5T4 UdAb was constructed in previous work and preserved in our laboratory. All anti-5T4 UdAb mutants with single-point mutation were constructed with QuickMutation™ Site-directed Gene Mutagenesis Kit (Beyotime, Shanghai, China) following the manufacturer’s protocol. The expression of different mutant and wild-type anti-5T4 UdAbs were performed in E. coli HB2151 bacterial culture at 30 °C for 16 h accompanied by 1 mM isopropyl b-d-1-thiogalactopyranoside (IPTG). The cells were harvested and lysed by polymyxin B at 30 °C for 0.5 h, followed by centrifugation at 8800×g for 10 min. The supernatant was collected, filtered through 0.8 μm polyethersulphone membranes by sterile syringes and purified by Ni-NTA (Smart Lifesciences) following the manufacturer’s instructions. Briefly, the filtered supernatant was loaded over the column with Ni-NTA. The resin was washed by washing buffer (10 mM Na2HPO4, 10 mM NaH2PO4 [pH 7.4], 500 mM NaCl, and 20 mM imidazole), and proteins were eluted in elution buffer (10 mM Na2HPO4, 10 mM NaH2PO4 [pH 7.4], 500 mM NaCl, and 250 mM imidazole). The collected pure fractions were immediately buffer-exchanged into PBS and concentrated with Amicon ultra centrifugal concentrators (Millipore). Protein concentration was measured using the NanoDrop 2000 spectrophotometer (Thermo Fisher), and protein purity was examined by SDS-PAGE.Enzyme-linked immunosorbent assay (ELISA)For comparison of different CR3022 mutants, the RBD of Delta (B.1.617.2) strain and S protein of Omicron (B.1.1.529) strain at 100 ng per well was coated in 96 wells half area microplates (Corning #3690) overnight at 4 °C. The antigen-coated plate was washed by three times with PBST (PBS with 0.05% Tween-20) and blocked with 3% MPBS (PBS with 3% skim milk) at 37 °C for 1 h. Following three times washing with PBST, 50 μL of threefold serially diluted antibody in 1% MPBS was added and incubated at 37 °C for 1.5 h. The HRP-conjugated anti-Fab and anti-Fc (Sigma-Aldrich) secondary antibodies were used for the detection of different tested antibodies. After washing with PBST for five times, the enzyme activity was measured after the addition of ABTS substrate (Invitrogen) for 15 min. The data were acquired by measuring the absorbance at 405 nm using a Microplate Spectrophotometer (Biotek) and the EC50 (concentration for 50% of maximal effect) was calculated by GraphPad Prism8.0 software. To verify different UdAb mutants, the same experimental protocol as mentioned above was adopted. Briefly, the human TPBG/5T4-Fc antigen was coated in 96 wells half area microplates, then blocked with 3% MPBS and added serial diluted antibodies. The HRP-conjugated anti-Flag (Sigma-Aldrich) secondary antibody was used, followed by adding ABTS substrate and detected at 405 nm. The reported EC50 value is the mean value from three duplicates on a single experiment.Bio-layer interferometry (BLI) binding assayThe binding kinetics of different antibodies to SARS-CoV-2 Omicron S and 5T4 antigens were measured by BLI on an Octet-RED96 (ForteBio). Briefly, the his-tagged Omicron S protein at 8 μg/ml and biotinylated 5T4 protein at 5 μg/ml were separately loaded onto Ni-NTA and streptavidin-coated (SA) biosensors. The antigen immobilized sensors were incubated with threefold serially diluted CR3022 candidates or two-fold serially diluted anti-5T4 UdAbs starting at 300 nM in 0.02% PBST for 300 s for association, and then immersed into 0.02% PBST for another 300 s at 37 °C for dissociation. All the curves were fitted by 1:1 binding model using the Data Analysis software 11.1. All KD values were determined with R2 values of greater than 98% confidence level.Protein structure and ΔΔG
bind contribution analysisProtein structure analysis is conducted by python scripts. The antibody–antigen complex structure after mutation was obtained from Rosetta Flex-ddG relaxation10. The relaxed protein structure can provide more accurate side-chain conformations, which are critical for accurate contact and conformational analysis. The improved accuracy of such analyses enables a deeper understanding of the underlying binding mechanisms and can facilitate the identification of key characteristics involved in protein–protein interactions. The contribution scores are derived by using Integrated Gradients (IG)43, a model-agnostic attribution method, on GearBind to obtain residue-level interpretation following22. All protein structure figures are created with PyMOL v3.0.Molecular dynamics simulationFor antibody mutation structural analysis, we conducted molecular dynamics simulations of the wild-type and mutant antibody–antigen complex. Initial structures were taken from the Rosetta Flex-ddG relaxed structures used by GearBind. The LEaP module in the AMBER 22 suite was used for building starting structures and adding ions and solvent for the simulation44. The protonation states of the molecules were kept at the default settings as assigned by LEaP during the initial structure preparation. All systems were simulated with the ff19SB protein force field45 and solvated in boxes of water with the OPC346 solvent model. Simulated systems were solvated using a 10 Angstrom solvent box. All bonds involving hydrogen atoms were constrained with the SHAKE algorithm47. The particle mesh Ewald (PME) algorithm was used to calculate long-range electrostatic interactions48. Initial structures were relaxed with maximum 10,000 steps of minimization before convergence, then subjected to heating for 20 ps and equilibrating for 10 ps in the NPT ensemble with PMEMD. The CUDA version of PMEMD was used to accelerate the simulations49. The simulation temperature was set to room temperature 298K. All systems were simulated for 1 μs production MD with one replica, and samples in the first 50 ns were not used for analysis. CPPTRAJ module in AMBERTools was used to analyze the simulation results, including calculating the root mean square deviation (RMSD), room-mean square fluctuation (RMSF), hydrogen bonding, and dihedral angle50.Reporting summaryFurther information on research design is available in the Nature Portfolio Reporting Summary linked to this article.
Pretrainable geometric graph neural network for antibody affinity maturation
