Unsupervised Protein-Ligand Binding Energy Prediction via Neural Euler's Rotation Equation
Authors: Wengong Jin
Siranush Sarkizova
Xun Chen
Nir Hacohen
Caroline Uhler
Venue: NeurIPS 2023 poster (2023)
Keywords: Energy-based Models
Denoising Score Matching
Equivariant Neural Networks
TL;DR: We develop a SE(3) denoising score matching method for unsupervised antibody-antigen binding affinity prediction.

Abstract: Protein-ligand binding prediction is a fundamental problem in AI-driven drug discovery. Previous work focused on supervised learning methods for small molecules where binding affinity data is abundant, but it is hard to apply the same strategy to other ligand classes like antibodies where labelled data is limited. In this paper, we explore unsupervised approaches and reformulate binding energy prediction as a generative modeling task. Specifically, we train an energy-based model on a set of unlabelled protein-ligand complexes using SE(3) denoising score matching (DSM) and interpret its log-likelihood as binding affinity. Our key contribution is a new equivariant rotation prediction network called Neural Euler's Rotation Equations (NERE) for SE(3) DSM. It predicts a rotation by modeling the force and torque between protein and ligand atoms, where the force is defined as the gradient of an energy function with respect to atom coordinates. Using two protein-ligand and antibody-antigen binding affinity prediction benchmarks, we show that NERE outperforms all unsupervised baselines (physics-based potentials and protein language models) in both cases and surpasses supervised baselines in the antibody case.

Unsupervised Protein-Ligand Binding Energy Prediction via Neural Euler's Rotation Equation
Wengong Jin, Siranush Sarzikova, Xun Chen, Nir Hacohen, Caroline Uhler
Broad Institute of MIT and Harvard
[wj in, sarkizov xun, nwwcohah, hute org
Abstract
Protein-ligand binding prediction is fundamental problem in Al-driven drug dis- covery. Previous work focused on supervised learning methods for small molecules where binding affinity data abundant, but it is hard to apply the same strategy to other ligand classes like antibodies where labelled data is limited. In this paper, we explore unsupervised approaches and reformulate binding energy prediction as generative modeling task. Specifically, we train an energy energy-based model on a complexes using SE(3) denoising score matching
set of unlabelled (DSM) and interpret its log- ligliiol as binding affinity. Our key contribution is a new equivariant rotation prediction network for SE(3) DSM called Neural Euler's Rotation Equations (NERE). It predicts a rotation by modeling the force and torque between protein and ligand atoms, where the force is defined as the gradient of an energy function with respect to atom coordinates. Using two protein-ligan antibody-antigen binding affinity prediction benchmarks, we show that NERE out- performs all unsupervised baselines (physics-basee potentials and protein language models) in both cases and surpasses supervised baselines the antibody case.
and
1
Introduction
One of the challenges in drug discovery is to design ligands (small molecules antibodies) with high binding affinity to target protein. In recent years, many protein-ligand binding prediction models have been developed for small molecules [13] 19 43 These models typically trained large set of crystal structures labeled with experimental binding data. However, it is difficult to apply this supervised learning approach to other ligand classes like antibodies where binding affinity data is limited. For example, the largest binding affinity dataset for antibodies (31) has only 566 data points. Transfer learning from small molecules antibodies is also hard their structures very different.
Given this data scarcity challenge, we aim to develop general unsupervised binding energy prediction framework for small molecules and antibodies. The basic idea is to learn energy -aaedd model (EBM) [6] pootein-ligande complexes by maximizing the log-likelihood of crystal structures in the Protein Data Bank (PDB). This generative modeling approach is motivated by recent advances of protein language models (PLMs) which show that the likelihood of protein sequences correlated with protein mutation effects 220 [1]. However, PLMs are not applicable small molecules and they model the likelihood of protein sequences rather than structures. We suppose the likelihood of protein complex structures is better predictor of binding affinity, which is determined by the relative orientation between a protein and a ligand. Indeed, our method shares the same spirit as traditional physics-basee based energy functions ( designed by experts so that crystal structures have low binding energy (i.c., high likelihood). Our key departure is to use neural networks to learn more expressive energy function in data-driven manner. Importantly, our model is applicable to large-scale virtual screening because it does not require expensive molecular dynamic simulations 21]
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

Specifically, we train EBMs using SE(3) denoising score matching (DSM) B3 15 as maximum likelihood estimation is hard. For each training example, we create a perturbed protein-ligant complex ligand and ask our model to reconstruct the rotation/translation noise. The main challenge of SE(3) DSM in our context how to predict rotation noise from the score of an EBM in an equivariant manner. Motivated by physics, we develop a new equivariant rotation prediction network called Neural Euler's Rotation Equation (NERE). The key observation is that the score of EBM equals the force of each atom, which further defines the torque between protein and a ligand. Using Buler's rotation equation, we convert this torque into the angular velocity of ligand and its corresponding rotation matrix calculate SE(3) DSM loss. Importantly, NERE is compatible with any (E((3)-invariin) networks and guarantees equivariance without any further requirement network architectures.
by randomly rotating/rranslating
We evaluate NERE on two protein-ligams and antibody-antigen binding affinity benchmarks from PDBBind 3d and Structural Antibody Database (SAbDab) [31. We compare NERE with unsu-
pervised thsiics-bsed models like MM/GBSA [21], protein language models (ESM-1v 1v (0 and ESM-IF (TI,, and variety of supervised models regressed on experimental binding affinity data. To simulate real-world virtual screening scenarios, we consider two settings where input complexes are crystallized or predicted by docking software. NERE outperforms all unsupervised baselines across all settings and surpasses supervised models in the antibody setting, which highlights the benefit of unsupervised learning when binding affinity data is limited.
2
Related Work
Protein-ligand binding models fall into two categories: supervised and unsupervised learning.
Supervised models are trained on binding affinity data from PDBBind [36] 29 3 32, 40, [9, 43]. They typically represent protein-ligand 3D complex represented as a geometric graph and encode into a vector representation using neural network for affinity prediction. Unsupervised models are either based on physics or statistical potentials. For example, molecular mechanics generalized Born surface area (MM/GBSA) [2]] calculates binding energy based on expensive molecular dynamics. Our work is closely related to statistical potentials like DrugScore201s [6]. It learns an energy score for each atom pair independently based on their atom types and distance alone, which ignores the overall molccular context. Our key departure is to use neural networks to learn a more expressive energy function powered by context-aware atom representations.
Antibody-antigen binding models are traditionally based on scoring functions in protein docking programs [25] I [0 (39 []. They are typically weighted combinations of physical interaction terms whose weights are tuned so that crystal structures have lower energy than docked poses. Recently, protein language models like ESM-1v (20 and ESM-IF [1] have been successful in predicting mutation effects for general proteins and we evaluate their performance for antibodies this paper. Due to scarcity of binding affinity data, very few work has explored supervised learning for antibody binding affinity prediction. Existing models [23] are based on simple hand mandtfefttdd features using binding affinity data from SAbDab. For fair comparison, we implement own supervised baselines using more advanced neural architectures 38 28).
Denoising score matching (DSM) is powerful technique for training score-based generative models 33 and has been applied to protein structure prediction 41 and protein-ligand docking [5]. Our departure from standard DSM is two fold. First, standard DSM is based on Gaussian noise while our method uses rigid transformation noise based on an isotropic Gaussian distribution for SO(3) rotation group [15] Second, we parameterize the score as the gradient of an -ased model rather than the direct output of a neural network. We choose this parametcrization scheme because we need to predict energy and score at the same time.
Equivariant rotation prediction. Motivated by molecular docking, there has been growing interest in predicting rotations with equivariant neural networks. For instance, Ganea et al. [9 Stärk et al. (35 developed a graph matching network that first predicts a set of key points between protein and
a ligand, and then uses the Kabsch algorithm [14] to construct rigid transformation. Corso et al. [S] predicted rotations based on tensor field networks (37 and spherical harmonics. Our method (NERE) has two key differences with prior work. First, it does not require key point alignment and directly outputs an rotation matrix. Second, NERE does not require spherical harmonics and can be plugged into any SE(3) 3E(3)-invariant encoder architectures for equivariant rotation prediction.

---

Pocket
b
Jeposusn wwmmuuttis
Antigen
One-
ESM2
Energy
Energy B(A,X)
E(A.X)
ESM2
Antibody
X
Force 0E/oX
Force dETOX
Figure 1: The architecture of (AA,)) for small molecules (left) and antibodies (right). For small molecules, the input a binding pocket and a molecular graph encoded by a message passing network (MPN). For antibodies, the input contains antibody CDR residues and epitope residues, where residue features A come from an ESM2 protein language model.
Unsupervised Binding Energy Prediction
protein prooten-liigand (or an antibody-antigen) complex is geometric object that involves protein and ligand (a small molecule or an antibody). It is denoted as tuple (A, X), with atom features A == Aa ... an and atom coordinates X I1..... In. (column-wise concatenation). The key property of complex is its binding energy F,A, X). A lower binding energy means a ligand binds more strongly to a protein. In this section, we deseribe how to parameterize (AA, X) and design proper training objectives to infer the true binding energy function from a list of crystal structures without binding affinity labels.
3.1 EBM Architecture
An energy function E(A, is composed of protein encoder and output layer. The encoder is frame averaging neural network [28] that learns SE(3)-invariant representation h; for each atom. We choose this architecture because of its simplicity (detailed in the appendix). The output layer do a ferw-fowward neural network with one hidden layer. predicts the interaction energy do (h., h2) for each pair of atoms. Finally, we define E( X) as the sum of pairwise interaction energies:
E(A,X)= = Σ
go h2)
i,j:d <d
Since atomic interaction vanishes beyond certain distance, we only consider atom pairs with distance dı < d. So far, we have described our method in generic terms. We now specify the input features and preprocessing steps tailored small molecules and antibodies.
Small molecules. When the ligand a small molecule, the input to our model a looein--ligand complex where the protein is cropped to its binding pocket (residues within 10Å from the ligand). The model architecture is illustrated in Figure [.. On the protein side, each atom in the binding pocket is represented by one-hot encoding of its atom name (Ca, Cg, N, O, etc.). We consider all backbone and side-chain atoms. On the ligand side, the atom features A are learned passing network (MPN) 42 based on a ligand molecular graph. The MPN and energy function are optimized jointly during training.
message
Antibodies. When the ligand an antibody, the input to our model is an antibody-antigen binding interface. The binding interface is composed of residues in the antibody complementariy--deeermin region (CDR) of both heavy and light chains and top antigen residues (epitope) closest to CDR. The model architecture is depicted in Figure [p.. Different from the small molecule case, we only consider C. atoms of each antibody/antigen residue. Each residue is represented by a 2560 dimensional pre-trained ESM-2 embedding [18]. For computational efficiency, we freeze ESM-2 weights during training.
3.2 Training FBMs with Denoising Score Matching
Given that our training set does not have binding affinity labels, we need to design an unsupervised training objective different from a supervised regression loss. Our key hypothesis that we can infer the true binding energy function (up affine equivalence) by maximizing the likelihood of crystal
3

---

Algorithm Training procedure (single data point)
Require: A training complex X)
1: Sample a noise level Ø.
2: Sample rotation vector w NSso() with variance o²
3: Sample translation vector ~ (0,
4: Perturb the coordinates X by applying rigid transformation (w, t) to the original complex.
5: Compute the score of energy function (w, t) based on its gradient (force) -0E/OX.
6: Minimize DSM objective ddsm-
structures in our training set. The motivation of our hypothesis is that a crystal structure the lowest energy state of a -rigain-ligand complex. The maximum likelihood objective seeks to minimize the energy of crystal structures since the likelihood
complex is p( A, X) x exp(-E(A, X)).
While maximum likelihood estimation (MLE) is difficult for EBMs due to marginalization, recent works 33 34 has successfully trained EBMs using denoising score matching (DSM) and proved
that DSM is a good approximation of MLE. In standard DSM, we create a perturbed complex by adding Gaussian noise to ligand atom coordinates, i.c., X = X + e, where é ~ p(E) = (O, o² DSM objective tries to match the score of our model -ÔE xx and the score of the noise distribution Ve p(e) -c/0²:
4 = E[ -
Intuitively, lg forces the gradient to be zero when the input complex is a crystal structure (c = 0). As result, a crystal structure pose will be the local minima of an EBM under the DSM objective.
3.3 Perturbing Complex with Rigid Transformation Noise
Nevertheless, adding Gaussian noise is not ideal for -igeen binding because it may create nonsensical conformations that violate physical constraints (e.g., an aromatic ring must be planar). A better solution is to create perturbed complex (A, X) via random ligand rotation and translation, similar molecular docking. To construct a random rotation, we sample rotation vector w from N Nso(3). an isotropic Gaussian distribution over SO(3) rotation group [15) with variance o². Each S30(3) has the form w Aúo, where W a vector sampled uniformly from a unit sphere and 0 e [0, m] rotation angle with density
(())
cos
Σ
20
(21 + 1)e 11+1)00 sin ((1 + 1/2)0)
π
sin(0/2)
Likewise, we sample a random translation vector t from normal distribution o(0,021) Finally, we apply this rigid transformation the ligand and compute its perturbed coordinates X R. X +t, where R., is the rotation matrix given by the rotation vector w (W., Wy: wz).
w
Wy
R (x)) Ww
W= Wy
0
Ww
(4)
Wr
Here exp means matrix exponentiation and W., is an infinitesimal rotation matrix. Since W., is skew symmetric matrix, its matrix exponential has the following closed form
R = exp(pppp + 6 W + 02 W
(5)
sin ||م|||
1 cos |||
C1
C2
|
Moreover, we do not need to explicitly compute the matrix exponential R since W., is the linear mapping of cross product, i.e. cross product operations that are very efficient:
War. Therefore, applying rotation matrix only involves
Rw x, + Ciw n. c2w X (w x.)

---

Input complex
Residue features
Residue coords
Random rigid
(4
transformation
Energy- based Model
Energy Ex()
Rotation p()
Translation p(
Euler's rotation equation
Aggregate forces
Gradient
V ExE
DSM Loss
Translation
Rotation à
Torque
B.)
Forces
Figure 2: NERE predicts rotation in three steps. It first calculates pairwise residue force based on the gradient f: = (-0E/Oz)". then predicts the corresponding torque and angular velocity by solving Euler's rotation equation. Lastly, it converts the angular velocity vector to a rotation via matrix exponential map.
3.4 Training EBMs with SE(3) Denoising Score Matching
Under this perturbation scheme, DSM aims to match the model score -0E/aX and the score of rotation and translation noise V. log p(w), V. log p(t). Our SE(3) DSM objective losses: fasm = l + l,, where la and l, correspond to a translation and rotation DSM loss. The translation part is straightforward since follows normal distribution and V. log p(t) t//22
of two
i i-- E, 0E/02:
4 = E [[] V Voggg(())))
The rotation part of DSM is more complicated. As w is sampled from a uniform distribution over sphere (whose density is constant), the density and score of p(w) is
p(w) (()), V logr(w) Vo log s(0)
In practice, we calculate the density and score by precomputing truncated infinite series in f(0). However, the main challenge is that the model score -0E/OX E R=x3 is defined over atom coordi- nates, which is not directly comparable with V log)(w) R$ they have different dimensions. To address this issue, we need map -0E/0X a rotation vector G via rotation predictor F and perform score matching in the rotation space:
E E ت
V log)(w)))?)
== F(-E///)
(10)
The rotation predictor F is based Neural Fuler's Rotation Equation, which will be described in the next section. The overall training procedure is summarized Algorithm and Figure
Neural Euler's Rotation Equation
In this section, present Neural Fuler's Rotation Equation (NERE), which be used in SE(3) DSM to infer rotational scores from the gradient -0E/OX. To motivate method, let us first review some basic concepts related to Euler's Rotation Equation.
4.1 Euler's Rotation Equations
In classical mechanics, Euler's rotation equation is first-order ordinary differential equation that describes the rotation rigid body. Suppose a ligand rotates around its center mass f with angular velocity w. Euler's rotation equation in an inertial reference frame is defined as
di dt
=T.
Σ, (2. 11) f:
(11)
where In R3x3 is the inertia matrix of a ligand, T the torque received, and fi the force applied to ligand atom The inertia matrix describes the mass distribution of a ligand and the
5

---

torque needed for desired angular acceleration, which is defined as
IN xx 1 11 I For short period of time At, approximate the angular velocity wi=Ai by
(.. p)) - 11)
(12)
du Wi=at di - W IN'T
(13)
Al
Since we assume the system is in an inertial reference frame (Wi=D = 0), we have Wi=At IN AAt (we set AL 0.1). We note that calculating the inverse I I is cheap because is a x matrix.
4.2 Applying NERE to SE(3) DSM
In SE(3) DSM, our goal is to map the gradient -0E/Ox; to a rotation vector w. In physics, the gradient -0E / 0o: is the force fi of atom i. Therefore, the rotation predictor F in Eq. 10 is simple application of the Euler's s Rotation Equation
w F(-0E//0 IN T&t,
T T= -=E, (z. (1) x (-OE/Ox,)
(14)
A nice property of NERE is that it is equivariant under SO(3) rotation group because it is derived from physics. We formally state this proposition as follows and its proof in the appendix. Proposition 1. Suppose we rotate a ligand so that its new coordinates become x² Ray The new force f', torque T', iinertia matrix N, and angular velocity w' for the rotated complex are
RR: = RT, -Rx w' Rww In other words, NERE is equivariant under SO(3) rotation group.
Once we establish SO(3) equivariance, is easy to satisfy SE(3) equivariance by first placing the rotation center at the origin (æ; 2; p1), applying the predicted rotation via NERE, and then adding y back to each atom.
Experiments
We evaluate our model on two drug discovery applications: protein-ligand binding and antibody- antigen binding affinity prediction. The experimental setup described follows.
5.1
Protein-Ligand Binding
Data. Our training data has 5237 ligand complexes from the refined subset of PDBbind v2020 database [36]. Our test set has 285 complexes from the PDBbind
with binding affinity labels converted into log scale. Our validation set has 357 complexes randomly sampled from PDBbind by Stärk al. (35] after excluding all test cases. The final training set has 4806 complexes (without binding affinity labels) after removing all ligands overlapping with the test set.
Metric. We report the Pearson correlation coefficient between true binding affinity and predicted energy E(A, X). We do not report root mean square error (RMSE) because our model does not predict absolute affinity values. In fact, shifting [AA, X)) by any constant will be equally optimal under the DSM objective. We run our model with five different random seeds and report their average.
Baselines. We consider three sets of baselines for comparison:
Physics-based potentials calculate binding affinity based energy functions designed by ex- perts. We consider four popular methods: Glide (8) AutoDock,jm. (7]. DrugScoreznss [6], and MM/GBSA 21 Among these methods, MM/GBSA is the most accurate but computationally expensive. It takes one hour to calculate energy for just one complex on a 64-core CPU server. Unsupervised models. Since unsupervised learning is relatively underexplored in this area, we implement two unsupervised EBMs using the same encoder architecture NERE but different training objectives. The first baseline is the standard Gaussian DSM (see Section 3.2) The second baseline is contrastive learning 4. For each crystal structure (A, X),
, we apply K random rigid transformations to obtain K perturbed protein-ligand complexes X i..... X k as negative samples. Suppose - F, X;) the predicted energy for X,, we train our EBM to maximize the likelihood of the crystal structure - (epp-E(AAX)) X))/ E. exp(-E( E(A,X)),

---

Sup
Crystal
Docked
Sup
Crystal
Docked
IGN
0.837
0.780
FANN
0.325
0.326,019
TankBind
n/a
0.824
FANN transfer
0.350 033
0.359 .099
PLANET
0.811
0.799
APSPISA
0.323
0.144
Glide
0.467
0.478
ZRANK
0.318 X
0.163
Autodockvina
0.604
0.578
PYDOCK
0,248 X
0.164
DrugScorez018
0.602
0.540
ESM-Iv
0.038
0,038
MM/GBSA
0.647
0.629
ESM-IF
0.024
0.025
Contrastive
0.625 .012
0.623,002
Contrastive
X
0.308 .37
0.304 ,21
Gauss DSM
0.638.017
0.632.016
Gauss DSM
0.335 .038
0.330
NERE DSM
0.656.a12
0.651.m1 .011
NERE DSM
X
0.361 .05
0.360
Table 1: Pearson correlation on PDBBind (left) and SAbDab test sets (right). We report standard deviations for all models implemented in this paper. "Sup"" indicates whether a method is supervised or unsupervised. "Crystal"""Docked" means a model evaluated on crystal or docked structures. The performance of TankBind crystal structures is not available because it performs docking and affinity prediction simultancously in one model.
Supervised models. Most of the existing deep learning models for binding affinity prediction belong to this category. They are typically trained on the entire PDBBind database with over 19000 binding affinity data points. The purpose of this comparison is to understand the gap between supervised and unsupervised models when there is abundant labeled data. For this purpose, we include three top- peppeffforing methods (IGN [3]. TankBind recent survey by /Zhang al. (2023).
19] and PLANET (3)) taken from a
Results (crystal structures). In our first experiment, we evaluate all methods on the crystal structure of protein-ligand complexes, As shown in Table (eft,, crystal column), NERE outperforms all physics-based and unsupervised models. We note that NERE is orders of magnitude faster than the second best method MM/GBSA (100ms v.s. Ihr per complex). As expected, the three supervised models perform better than NERE because they are trained on 19000 labeled affinity data points. Nonetheless, NERE recovers almost 80% of their performance without using any labeled data, which implies that binding affinity is closely related to the geometric structure of an input complex.
Results (docked structures). The previous evaluation considers an idealistic setting because crystal structures are not available real-world virtual screening projects. To this end, we use AutoDock Vina (7) to dock each ligand in the test set its binding pocket and evaluate all methods docked protein-ligand complexes. As shown in Table docked column), NERE consistently outperforms the baselines this setting. Its performance is quite close the crystal structure setting because docking error is relatively low (median test RMSD=3.55), In summary, these results suggest that NERE has learned accurate binding energy function useful for structure-based virtual screening.
5.2
Antibody-Antigen Binding
Data. Our training and test data come from the Structural Antibody Database (SAbDab) 31 which contains 4883 non-redundant antibody-antigen complexes. Our test set is subset SAbDab (566 complexes) that has binding affinity labels. Our training set has 3416 complexes (without binding affinity labels) after removing antigen/antibody sequences that appear in the test set. Our validation set comes from Myung al. [23]. which has 116 complexes after removing antibodies antigens overlapping with the test set.
Baselines. Similar to the previous section, we consider three sets of baselines for comparison:
Physics-based models. We consider eight biophysical potentials implemented in the CCharPPI webserver 22 and report the top three models (full results are listed the appendix).
Unsupervised models. Besides contrastive learning and Gaussian DSM baselines used in Sec- tion 5.1] we consider two state-of-the-art protein language models (ESM-Iv (0 and ESM-IF [T]. These two models have been successful in predicting mutation effects or binding affinities for
7

---

general proteins and we seek to evaluate their performance for antibodies. The input to ESM-1v is the concatenation of an antibody and an antigen sequence and we take the pseudo log- likellihood of pair. The input to ESM-IF is
antibody CDR residues as the binding affinity of an antibody-antiges complex and we take the conditional log-likelihood of
the crystal structure of an antibody CDR residues given the backbone structure as its binding affinity.
Supervised models. We also compare NERE with supervised models trained on external binding affinity data. The model is based on the same frame- frame-averaging averaging neural network (FANN) (28) and pre-trained ESM-2 residue embeddings used by NERE Since all labeled data from SAbDab are in the validation and test sets, we draw additional data from an protein protein-protein binding mutation database (SKEMPI) [2]. We obtain 5427 binding affinity data points after removing all complexes appeared in the test set, from which 273 instances are antibody-antige
complexes. We explore two training strategies: 1) training on 273 antibody-antigen data points only (FANN ): 2) training on 5427 data points first and then finetuning on 273 antibody-antigent data points (FANNmandee).
Results (Crystal structure). We first evaluate all models on crystal structures and report the Pearson correlation between true and predicted binding energy. As reported in Table (right, crystal column), NERE significantly outperforms physics- based methods and protein language models. ESM-Iv and ESM-IF have almost zero correlation because they only model the likelihood of the antibody sequence. ESM-IF models the conditional likelihood of a sequence given its structure while NERE models the likelihood of the entire protein complex.It It is crucial to model the likelihood of structures because binding affinity depends on the relative orientation between an antibody and an antigen. We also find that NERE outperforms the supervised models, even with a small gain over FANN trained on 5427 labeled data points. Overall, these results highlight the advantage of unsupervised learning in low-data regimes.
Results (Docked structure). In our second experiment, we evaluate all methods on docked complexes to emulate a more realistic scenario. We use the ZDOCK program (26] to predict the structure of all antibody-antigen complexes in the test set (its median RMSD is 19.4). As reported in Table []] (right, docked column), our model still outperforms all the baselines in this challenging setting. The performance of EBMs are quite robust docking errors regardless of the training algorithm (contrastive learning, Gaussian DSM, or NERE DSM). In summary, our model is capable to predict muiibody-antigen binding even when crystal structure is not available.
5.3 Ablation Studies and Visualization
Visualizing energy landscape. First, we study how the learned energy changes with respect to ligand orientations. Given an input complex, we perform a grid search of ligand rotation angles = =- W wa] and plot the predicted energy for each pose. As 3D contour plot hard to visualize, we decompose it into three 2D contour plots by fixing one of the three axis (W), W2, Wa) to zero. Ideally, the crystal structure (w- 10, O, 0)) should be the local minimum because it is physically the most stable conformation. Figure Bpb show contour plots for small molecules and antibody-antige complexes. We find that their crystal structures are located relatively near the local minimum.
Visualizing interaction energy. Since the predicted energy is a summation of all pairwise interactions Eij.d,, Φo (h,,h,)) h,). , seek visualize the contribution of different residues to predicted binding energy. Figure 3 is one example. Each row and column in this heatmap represent an epitope residue and an antibody CDR residue, respectively. Each entry in the heat map is the interaction energy between two residue do (h;, h2). An entry is left blank if the distance di,j > 20A. Interestingly, we find that the model pays the most attention to CDR-H3 and CDR-L3 residues. In most test cases, their energy is much lower than other CDR residues. This agrees with the domain knowledge that CDR-H3 and CDR-L3 residues are the major determinant of binding.
Ablation studies. Lastly, we perform four ablation studies to understand the importance of different model components. We first replace our SE(3)- invariint protein encoder with a invaiiatttnn neural network where atom 3D coordinates are directly concatenated to atom features in the input layer. As shown in Figure 31, the model performance drops significantly when the encoder is not SE(3) invariant. In addition, we run NERE DSM by removing the rotation é, or translation DSM term & from the objective fSM 4 e, We find that removing either of them substantially hurts the model
'We have also tried other neural network architectures like 3D convolutional neural networks and graph convolutional networks, but they did not achieve better performance. Their results are shown in the appendix.
8

---

Light chain
Heary chain
Antizen
460
Full Model
Non-invaraint Encoder
Rotation DSM only
Translation
No
DSM only
ESM
0.8
-15
0.656
0.622
Pearson Correlztion
0.6
0.382
0:30 0.355
0.303 0.312
COTFSSTALIP IIIGTAM TLLGSOYESI IESIIIM THONAEOOTTNYEE
POBBind
SAbDab
Figure 3: a-b) Visualizing learned energy landscape for small molecules and antibodies. We perform a grid search of ligand rotation angles w = [W1, W2, w3] and plot the predicted energy as 2D contour plots with one of the axis w; fixed. The crystal structure to be the local minimum of the energy landscape. c) The heat map of our learned energy function. Each entry represent the binding energy between residues (i,j) (the darker the stronger). An entry left blank (grey) if their distance d, 20A. Our model correctly puts more attention CDR-H3/L3 residues. d) Ablation studies on different model components.
the origin (marked with $) and supposed
performance. Therefore, it is crucial to consider both rotation and translation degree of freedom in unsupervised binding energy prediction. Lastly, we train NERE DSM by replacing ESM-2 amino acid embedding with one-hot encoding (this change only applies to the antibody model), We find that ESM-2 helpful but NERE DSM is still able infer binding affinity with reasonable accuracy after removing language model features.
6 Discussion
In this paper, we developed energy-based model for unsupervised binding affinity prediction. The energy-based model was trained under SE(3) denoising score matching where the rotation score was predicted by Neural Euler's Rotation Equation. Our results show that the learned energy correlates with experimental binding affinity and outperforms supervised models in the antibody case.
Limitations. Indeed, there are many ways to improve our model. Our current model for antibodies only considers their backbone structure (Ca atoms) and ignores all side-chain atoms. While the model for small molecules include all atoms, we have not incorporated the flexibility of side chains protein or rotatable bonds in a ligand. Considering their interaction is crucial for protein-ligant binding prediction. Our future work is to extend our encoder with all-atom structures and our NERE DSM algorithm with side-chaie torsional rotations.
Broader Impacts. Our method is applicable to wide range biological areas like drug discovery, immunology, and structural biology. We believe our work does not have any potential negative societal impacts since our aim is to accelerate scientific discovery.
Code and Data. Our code and data are available liggttub..... com http///nnoogg--oo jin LinnDSMBind

---

Acknowledgements
We would like to thank Divya Nori, Chenyu Wang, and anonymous reviewers for their valuable feedback the manuscript. Wengong Jin is supported by the BroadIgnite Award and the Eric and Wendy Schmidt Center at the Broad Institute of MIT and Harvard. Caroline Uhler was partially supported by NCCIH/NIH (IDP2AT012345), ONR (N0001--2--1-21166,, AstraZeneca, and Simons Investigator Award.
References
[1] R. F. Alford, A. Leaver-Fay, J. R. Jeliazkov, M. J. OeMear, F. P. DiMaio, H. Park, M. V. Shapovalov, D. Renfrew, V. Mulligan, K. Kappel, et al. The rosetta all-atom energy function for macromolecular modeling and design. Journal of chemical theory and computation, 13(6))3031-3048, 2017.
[2] N. Andrusier, R. Nussinov, and H. J. Wolfson. Firedock: fast interaction refinement in molecular docking. Proteins:
Structure, Function, and Bioinformatics, 69(1):139-155)) 2007.
[3] L. Chen, X. Tan, D. Wang. F. Zhong, X. Liu, T. Yang. X. Luo, K. Chen, H. Jiang, and M. Zheng. Transformerepi: improving compound proteen interaction prediction by sequence-based deep learning with self-attention mechanism and label reversal experiments. Bioinformatics, 36(16); 4406 4414, 2020.
[4] T. Chen, S. Korbliihh, M. Norouzi, and G. Hinton. A simple framework for contrastive learning of visual representations. In International conference on machine learning, pages 1597-1607. PMLR, 2020.
[5] G. Corso, H. Stärk, B. Jing, Barzilay, and T. Jaakkola. Diffdock: Diffusion steps, twists, and turns for molecular docking. arXiv preprint arXiv:2210.01776, 2022.
[6] Dittrich, D. Schmidt, C. Pfleger, and H. Gohlke. Converging knowledge-based scoring function: Drugscore2018. Journal of chemical information and modeling. 59(1):509-55) 521, 2018.
[7] Eberhardt, D. Santos-Martins, A. F. Tillack, and Forli. Autodock vina 1.2.0: New docking methods, expanded force field, and python bindings. Journal of Chemical Information and Modeling, 6188):3891-338 -38988, 2021.
[8] R. A. Friesner, R. B. Murphy, M. P. Repasky, L. L. Frye, J. R. Greenwood, T. A. Halgren, P. C. Sanschagrin, and Mainz. Extra precision glide: Docking and scoring incorporating model of hydrophobic enclosure for protein- ligand complexes. Journal of medicinal chemistry, 49(21):6177-61) 6196, 2006.
[9] O.-E. Ganea, X. Huang, C. Bunne, Y. Bian, R. Barzilay, T. Jaakkola, and A. Krause. In- dependent equiiaaqui models for end-to-end rigid protein docking. arXiv preprint arXiv:2111.07780 2021.
[10] Grosdidier, C. Pons, Solernou, and J. FeerFeeenn-eeeeeee Recio. Prediction and scoring of docking poses with pydock. Proteins: Structure, Function, and Bioinformatics, 694))8852-888 2007.
[11] C. Hsu, Verkuil, Liu, Z. Lin, Hie, T. Sercu, Lerer, and Rives. Learning inverse folding from millions of predicted structures. In International Conference on Machine Learning. pages 8946-8970. PMLR, 2022.
[12] J. Jankauskaité, B. Jiménez--araee J. Dapkūnas, J. Fernrndeen-eeeeeeeeeee -ecino, and I. H. Moal. Skempi 2.0: an updated benchmark of changes in protein-protein binding energy, kineties and thermo- dynamics upon mutation. Bioinformatics, 35(3):462-469) 2019.
[13] D. Jiang, C.-Y. Hsieh, Z. Wu, Kang, Wang. Wang, B. Liao, C. Shen, L. Xu, J. Wu, et al. Interactiongraphhnt:: A novel and efficient deep graph representation learning framework for accurate protein-ligand interaction predictions. Journal of medicinal chemistry, 64(24)) 18209-18232, 2021.
10

---

[14] W. Kabsch. A solution for the best rotation to relate two sets of vectors. Acta Crystallograhhica Section A: Crystal Physics, Diffraction, Theoretical and General Crystallography, 32(5):922-
923, 1976.
[15] A. Leach, S. M. Schmon, M. T. Degiacomi, and C. G. Willcocks. Denoising diffusion proba- bilistic models on so (3) for rotational alignment. In ICLR 2022 Workshop on Geometrical and Topological Representation Learning, 2022.
[16] Y. LeCun, S. Chopra, R. Hadsell, M. Ranzato, and F. Huang. A tutorial on energy-based learning. Predicting structured data, 1(0), 2006.
[17] T. Lei. When attention meets fast recurrence: Training language models with reduced compute. arXiv preprint arXiv: 102. 12459, 2021.
[18] Z. Lin, H. Akin, R. Rao, B. Hie, Z. Zhu, W. Lu, A. dos Santos Costa, M. Fazel-Zarandi, T. Sercu, S. Candido, al. Language models of protein sequences the scale of evolution enable accurate structure prediction. bioRxiv, 2022.
[19] W. Lu, Q. Wu, J. Zhang, J. Rao, C. Li, and S. Zheng. Tankbind: Trigonometry-awee neural networks for drug-protedd binding structure prediction. bioRviv, 2022.
[20] J. Meier, R. Rao, R. Verkuil, J. Liu, T. Sercu, and A. Rives. Language models enable zero-shot prediction of the effects mutations on protein function. Advances in Neural Information Processing Systems, 34:29287-29300 2021.
[21] B. R. Miller III, T. D. McGee Jr, J. M. Swails, N. Homeyer, H. Gohlke, and A. E. Roitberg. Mmpbsa.py: an efficient program for end-state free energy calculations. Journal of chemical theory and computation, 3(9):3314-3321,
[22] I. H. Moal, B. Jiménez-Gareía, and J. Fernández-Rees Recio. Ccharppi web server: computational characterization of protein-protein interactions from structure. Bioinformatics, 1231):123-125, 2015.
[23] Y. Myung, D. E. Pires, and D. B. Ascher. Csm-ab: graph-based antibody- binding affinity prediction and docking scoring function. Bioinformatics, 38(4): 141-, 2022.
[24] B. Pierce and Weng. Zrank: reranking protein docking predictions with an optimized energy function. Proteins: Structure, Function, and Bioinformatics, 67(4): 2007.
[25] B. Pierce and Z. Weng. A combination of rescoring and refinement significantly improves protein docking performance. Proteins: Structure, Function, and Bioinformatics, 72(1):270-279 2008.
[26] B. G. Pierce, K. Wiehe, H. Hwang, B.-H. Kim, T. Vreven, and Z. Weng. Zdock server: interactive docking prediction protein-protein complexes and symmetric multimers. Bioinformatics, 30 112):1771-1773, 2014.
[27] C. Pons, D. Talavera, X. De La Cruz, M. Orozco, and J. Fernandez-Recio. Scoring by in- termolecular pairwise propensities of exposed residues (sipper): new efficient potential for protein- protein docking. Journal of chemical information and modeling, 51(():370-377, 2011.
[28] O. Puny, M. Atzmon, H. Ben- Gen-Hamu, E. J. Smith, I. Misra, A. Grover, and Y. Lipman. Frame averaging for invariant and equivariant network design. arXiv preprint arXiv:2110.03336, 2021.
[29] M. Ragoza, J. Hochuli, E. Idrobo, J. Sunseri, and D. R. Koes. Protein-ligand scoring with convolutional neural networks. Journal of chemical information and modeling, 57(4)9942-957, 2017.
[30] D. Ravikant and R. Elber. Pic-efficient filters and coarse grained potentials for unbound docking. Proteins: Structure, Function, and Bioinformatics, 78(2):400-419,
2010.
[31] C. Schneider, M. I. Raybould, and C. M. Deane. Sabdab in the age of biotherapeutics: updates including sabdab-nano, the nanobody structure tracker. Nucleic acids research, 50(DI):D1368. D1372, 2022.
11

---

[32] V. R. Somnath, C. Bunne, and A. Krause. Multi-scale representation learning on proteins. Advances Neural Information Processing Systems, 34:25244-25255, 2021.
[33] Y. Song and S. Ermon. Generative modeling by estimating gradients of the data distribution. Advances in Neural Information Processing Systems, 32, 2019.
[34] Y. Song, C. Durkan, I. Murray, and S. Ermon. Maximum likelihood training of score-basee based diffusion models. Advances in Neural Information Processing Systems, 34:1415-1428, 2021.
[35] H. Stärk, O. Ganea, I.. Pattanaik, R. Barzilay, and T. Jaakkola. Equibind: Geometric deep learning for drug binding structure prediction. In International Conference Machine Learning. pages 20503-20521. PMLR, 2022.
[36] M. Su, Q. Yang, Y. Du. G. Feng, Z. Liu, Y. Li, and R. Wang. Comparative assessment of scoring functions: the casf-2016 update. Journal of chemical information and modeling, 59(2):895-915, 2018.
[37] N. Thomas, T. Smidt, S. Kearnes, L. Yang, L. Li, K. Kohlhoff, and P. Riley. Tensor field networks: Rotation-and conslation-equivarree neural networks for 3d point clouds. arXiv preprint 1tt.:0880000299 2018.
[38] R. J. Townshend, M. Vögele, P. Suriana, A. Derry, A. Powers, Y. Laloudakis, S. Balachandar, B. Jing, B. Anderson, S. Eismann, et al. Atom3d: Tasks on molecules in three dimensions. arXiv preprini arXiv:2012.04033, 2020.
[39] S. Viswanath, D. Ravikant, and R. Elber. Improving ranking of models for protein com- plexes with side chain modeling and atomic potentials. Proteins: Structure, Function, and Bioinformatics, 81(4):592-66 -0-699996666., 2013.
[40] Wang. Zheng, Y. Jiang, C. Li, J. Liu, C. Wen, A. Patronov, D. Qian, H. Chen, and Yang. Structure-aware multimodal deep learning for trug-protein interaction prediction. Journal of chemical information and modeling,
2022.
[41] J. Wu, S. Luo, T. Shen, H. Lan, S. Wang, and J. Huang. Ebm-fold: fully-differenniable protein folding powered by energy- based models. arXiv preprint arXiv:2105.04771, 2021.
[42] K. Yang, K. Swanson, W. Jin, C. Coley, P. Eiden, H. Gao, A. Guzman-Perez, T. Hopper, B. Kelley, M. Mathea, et al. Analyzing learned molecular representations for property prediction. Journal of chemical information and modeling, 5978):3370-3388, 2019.
[43] X. Zhang, H. Gao, H. Wang, Z. Chen, Z. Zhang, X. Chen, Li, Y. Qi, and R. Wang. Planet: A multi-objective graph neural network model for protein-ligand binding affinity prediction. bioRviv, pages 2023-02, 2023.
12

---

A NERE Rotational Equivariance
Proposition Suppose we rotate protein-ligand complex so that new coordinates become x! Ræ, The new force f', torque T', inertia matrix I',, and angular velocity w' under the rotated complex are
= Rf RT, I' RINR W Rw;
Proof. After rotating the whole complex, the energy function E(A, = E'A, X) since the encoder SE(3)-invariant. Given that æ; R a/ and Pax/0xx/ = R the new force becomes
(
T
(
(
dz
Ozz dxy
f =
J: = RRf:
Oa
dz.
dxy
Based on the definition of torque and the fact that cross products satisfy Ræ x Ry = R(x x y), we have
T'
Σ (at k²) x f Σ RRe; Ru) x Rf:
Σ a 11) x f: = RT
Likewise, using the fact that RR I, the new inertia matrix becomes
NN= =
Σ ||,| 12'1211 (x) - 1a' )²)T
Σ | 1ll21 (Rx; Rp)( (Ræ; R4))
Σ R|x IR R(æ;; (x) 11) R
R Σ | uul21 (x: 12) 11) R
RIyR R
For angular velocity, we have
w = CI' CRI CRI Rw
Therefore, NERE layer is equivariant under rotation.
B Experimental Details
Protein encoder architecture. The energy function needs to be SE(3)-invariant and differentiable with respect to X. Thus. we adopt the frame averaging neural network [28) so that E directly takes coordinates X as input rather than distance matrix. To be specific, our energy function is parameterized as follows
(h hn
Σ dn (A.
(15)
SLEG
The encoder da modified transformer network [17] that learns atom representations hy,... h based on atom features A and coordinates X. Specifically, the model first projects the coordinates onto set of eight frames {9x (æx)) defined Puny et (28). concatenate the projected coordinates with atom features a,, encode the vector sequence [a gu (zi)], ... [an, gu (2,)) to their hidden
representations h ... hh, and then average the frame representations h; Σ. h /8 to maintain SE(() invariance.
Model hyperparameters. the small molecule case, our model has two components: molecular graph encoder (MPN) and frame-averaging encoder. For the MPN encoder, use the default
13

---

Crystal
Docked
ZRANK
0.318
0.163
RosettaDOCK
0.064
0.025
PYDOCK
0.248
0.164
SIPPER
-0.138
0.003
AP AP_PISA
0.323
0.144
FIREDOCK
0.101
-0.052
FIREDOCK AB
0.199
0.042
CP PIE
0.234
0.120
3D CNN
0.286.021
0.162.017
GNN
0.244.075
0.249.075
NERE (ours)
0.361.051
0.360.os1
Table 2: Additional baselines on the SAbDab test set.
hyperparameter from Yang et al. [42]. For the protein encoder, we set hidden layer dimension to be 256 and try encoder depth L 3. and distance threshold d (5.0, 10.0}. In the antibody case, we try encoder depth from E {1, 2, 3) and distance threshold E {10.0. 20.03. both cases, we select the hyperparameter that gives the best Pearson correlation the validation set.
Docking protocol. For Autodock Vina, we use default docking parameters with docking grid dimension of 20A, grid interval of 0.375A, and exhaustiveness of 32. For ZDOCK, mark antibody CDR residues as ligand binding site and generate 2000 poses for each antibody-antigen pair. We re-score those 2000 poses by ZRANK2 and select the best candidate.
Additional results. We include additional results for the antibody-antigen binding task. We first compare our method with additional physic-based potentials implemented in the CCharPPI web server [22), including ZRANK (24), RosettaDock (1) PyDock [10], SIPPER (27), AP_ AP_PISA (39]. CP PIE 30 FIREDOCK, and FIREDOCK AB As shown Table the performance of these
models are much lower than NERE.
We have also explored more protein encoder architecture as additional supervised learning baselines. We consider 3D convolutional neural network (3D CNN) and graph convolutional network (GCN) implemented in the Atom3D package (38). For the GCN model, extend its implementation to include ESM-2 residue embedding. We find that the performance of these two models are much worse than FANN when trained 273 antibody-antigen complexes from SKEMPI. Therefore, choose the FANN architecture as our default supervised model architecture for the rest of our analysis.
14

---

Experimental binding affinity (log scale)
CASF (crystal)
CASF (docked)
-8
-10
100000-00000 Predicted binding energy
2300001000000- 100000-75000
Predicted binding energy
Experimental binding affinitty (log scale)
Sabdab (crystal)
Sabdab (docked)
-0
-10
-12
-14000
-6000 Predicted binding energy
Predicted binding energy
Figure 4: Correlation between predicted and experimental binding energy.
15


Reviews:
[{"review_id": "yhb8w2QQd3", "reviewer": "Reviewer_X8hw", "summary": "This paper introduces NERE, an unsupervised method for predicting protein-ligand binding affinity. The authors propose a generative modeling approach that utilizes an energy-based model (EBM) to capture the characteristics of protein-ligand complexes. The EBM is trained by maximizing the log-likelihood of crystal structures from the PDB database. NERE incorporates random rotations and translations of the separate parts of the molecule complexes and predicts the optimal rotation and translation as the EBM's output. To train the EBM, SE(3) denoising score matching is employed. The effectiveness of the proposed approach is evaluated through experiments on both protein-ligand and protein-protein binding. The results reveal a positive correlation between the predicted energy and binding affinity. The authors demonstrate that their method surpasses all unsupervised baselines and even outperforms supervised baselines in the case of antibodies.", "questions": "Please see weaknesses.\n", "limitations": "None.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "Strengths of the review paper:\n\n- The paper introduces a novel approach for unsupervised protein-ligand binding energy prediction by utilizing an equivariant rotation prediction network. This innovative approach contributes to the existing body of knowledge in the field.\n\n- The paper proposes a unique method to predict rotation and translation as part of an unsupervised data augmentation technique. This approach enhances the understanding and modeling of protein-ligand interactions, as well as the correlation with molecule interaction energies.\n\n- The experimental results presented in the paper demonstrate the effectiveness of the proposed approach. The approach outperforms all unsupervised baselines, showcasing its superiority in predicting protein-ligand binding energies. Notably, the method even surpasses supervised baselines in the case of antibodies, indicating its potential for practical applications.\n\n- The paper exhibits clear and well-written content. The authors effectively convey technical concepts and methods, ensuring the readers can grasp the details of the proposed approach and its implementation. The organization of the paper enhances readability and comprehension.\n\n- The authors demonstrate their awareness of the limitations of their work. ", "weaknesses": "- The paper lacks visualizations of the correlation between the predicted energies and binding affinities. Visual representations, such as scatter plots or correlation matrices, would enhance the understanding of the relationship and provide a clearer picture of the results. Additionally, providing information about the distribution of outliers would give insights into the robustness of the proposed approach.\n- The authors could consider spearman correlation and p-value as additional evaluation metrics, as the ranking between different ligands might be more useful for drug discovery and virtual screening instead of absolute binding affinity values.\n- I am curious about the potential impact of the proposed unsupervised method on downstream tasks when serving as a pre-training strategy. It would be valuable to investigate whether the learned representations from the unsupervised approach can be beneficial for supervised downstream tasks. \n- The paper lacks an explicit discussion on the significance and relevance of an unsupervised method in the context of protein-ligand binding energy prediction. It would be beneficial to design experiments or provide theoretical arguments that illustrate how the utilization of unlabeled data, in comparison to traditional supervised methods, contributes to the understanding and performance of the model. Demonstrating the utility of unsupervised learning in this specific area would provide a stronger rationale for the proposed approach.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your insightful comments and suggestions. Please read our response below and let us know if you have more questions.\n\n**Q1**: The paper lacks visualizations of the correlation between the predicted energies and binding affinities.\n\nIn rebuttal Figure 1a, we visualize the correlation between the predicted energies and binding affinities on the small molecule test set (CASF) and the antibody test set (SAbDab). We did not find significant outliers in either crystal or docked settings.\n\n**Q2**: The authors could consider spearman correlation and p-value as additional evaluation metrics\n\nIn the rebuttal Table 1, we report the spearman correlation and p-value of all methods on both test sets. Our model (NERE DSM) remains the best under both metrics.\n\n**Q3**: I am curious about the potential impact of the proposed unsupervised method on downstream tasks when serving as a pre-training strategy. It would be valuable to investigate whether the learned representations from the unsupervised approach can be beneficial for supervised downstream tasks.\n\nWe adopt an HER2 antibody affinity maturation dataset [1] as an additional benchmark to investigate the utility of NERE DSM as a pre-training strategy. This dataset has 422 antibody sequences with experimental binding affinities against HER2. The pre-trained NERE DSM model achieved an average spearman correlation of 0.230 on this dataset, while the best supervised baseline ($\\mathrm{FANN_{transfer}}$) only yields 0.090. To further improve the performance, we selected a subset of the HER2 data to fine-tune the pre-trained NERE DSM model and the supervised $\\mathrm{FANN_{transfer}}$ model. Specifically, we split the HER2 dataset into a training set of $N$ data points and use the rest of the dataset for testing. As shown in rebuttal Figure 1f, NERE DSM consistently outperforms $\\mathrm{FANN_{transfer}}$ when fine-tuned by $N=50, 100, 150, 200$ data points. These results suggest that NERE DSM is a more effective pre-training strategy than supervised pre-training.\n\n**Q4**: The paper lacks an explicit discussion on the significance and relevance of an unsupervised method in the context of protein-ligand binding energy prediction. It would be beneficial to design experiments that illustrate how the utilization of unlabeled data, in comparison to traditional supervised methods, contributes to the performance of the model.\n\nIndeed, the performance of unsupervised methods is much lower than supervised models on the CASF test set (PDBBind core set). However, the supervised models are trained and evaluated on the same PDBBind database and it is still unclear whether they generalize well to new molecules with different structures. To this end, we adopt the free energy perturbation (FEP) benchmark [2] (developed by Merck) as an additional independent test set. It has eight protein targets (cdk8, cmet, eg5, hif2a, pfkfb3, shp2, syk, and tnks2) and 264 ligands in total. As shown in Figure 1b, molecules in this test set have a different distribution from the PDBBind training set. Therefore, a model needs to generalize to a different chemical space to perform well on this test set.\n\nWe applied NERE DSM to this FEP test set, with a small modification in the input layer to include protein side-chain atoms (see our global response above). We calculate the Spearman correlation ($R_S$) between its predicted binding energy and experimental binding affinity for each of the eight targets and report the average correlation as our evaluation metric. Importantly, all hyper-parameters are selected using only the PDBBind validation set to ensure this experiment is a blind, unbiased evaluation. As shown in Figure 1c, our model outperforms not only the unsupervised models (contrastive learning and Gaussian DSM) but also the supervised models. The performance is much better than PLANET & IGN (0.380 vs 0.314 and 0.252) and slightly better than TankBind (0.380 vs 0.376). The Spearman correlation for each target is shown in Figure 1d. Overall, these results suggest that unsupervised models are more robust than supervised models when there is distributional shift (a common problem when models are deployed in the wild). Therefore, unsupervised models can be also useful for protein-ligand binding, even though there is a fair amount of labeled affinity data.\n\nReferences\n1. Shanehsazzadeh et al., Unlocking de novo antibody design with generative artificial intelligence, biorxiv 2023\n2. Schindler et al., Large-scale assessment of binding free energy calculations in active drug discovery projects, 2020\n"}, {"review_id": "9aWATCLFRP", "reviewer": "Reviewer_j9Y6", "summary": "The paper introduces an unsupervised learning approach for predicting protein-ligand binding energy, called NERE (Neural Equivariant Rotation Estimation). The authors employ SE(3) denoising score matching to train an energy-based model on a dataset of unlabeled protein-ligand complexes. The model's effectiveness is demonstrated through evaluations on protein-ligand and antibody-antigen binding affinity benchmarks from PDBBind and the Structural Antibody Database (SAbDab).", "questions": "1. Why was the SE(3) diffusion process chosen over the standard elucidating diffusion process? While the standard diffusion process may generate nonsensical conformations, the proposed diffusion process might not explore the entire conformation space.\n2. In section 4.1, is the center of mass, denoted as $\\mu$, computed based on atomic weight?\n3. How does the model perform in the \"Docked structure\" scenario when the docking error is high?\n4. Why does contrastive learning underperform compared to using score matching to train the energy-based model? Further elaboration on this would be beneficial.\n5. Are there any experiments conducted on the docking benchmark of the CASF test set? It seems that the proposed method is more suitable for identifying the correct binding pose (ranking among multiple poses of a specific molecule) rather than ranking between binding poses of multiple molecules.\n6. Different noise levels appear to lead to different energy functions, while only the minimal noise level can approximate the energy function of the real data distribution. How does the proposed method handle multiple noise scales?", "limitations": "The proposed method appears to be sensitive to docking error. However, in real-world scenarios, finding the true binding pose is a challenging problem. This sensitivity could limit the method's practical applicability.", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. The paper innovatively applies SE(3) score matching to train an energy-based model, positioning it as an unsupervised energy ranking model.\n2. The experimental results highlight the model's effectiveness, particularly when compared with other unsupervised models or physics-based models.\n3. The paper is well-structured and easy to comprehend, making the proposed method and its implications clear to the reader.", "weaknesses": "1. The proposed method only considers rigid transformations (translations & rotations) in the diffusion process, neglecting the degrees of freedom in small molecule ligands, such as additional rotatable bonds.\n2. The method employs a residue-level representation for proteins, which might omit fine-grained information that could be crucial for accurate binding prediction.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your insightful comments and suggestions. Please read our response below and let us know if you have more questions.\n\n**Q1**: The proposed method only considers rigid transformations (translations & rotations) in the diffusion process, neglecting the degrees of freedom in small molecule ligands, such as rotatable bonds. The method employs a residue-level representation for proteins, which might omit fine-grained information that could be crucial for accurate binding prediction.\n\nWe agree that it is beneficial to model the flexibility of ligands and incorporate the side-chain information of proteins. Given the limited time, we partially address this issue by including all protein atoms (including side-chains) as input instead of only $C_\\alpha$ atoms. This modification only changes the input of our model while the model architecture and NERE DSM objective remains the same. Our results suggest that including side-chain atoms substantially improves the performance of our method on some of the datasets. Please read our global response above and the rebuttal Figure 1e for details. In the future, we plan to model the flexibility of rotatable bonds and protein side-chain atoms by introducing a side-chain DSM loss, which adds rotation noise on the rotatable bonds and side-chain atoms and use NERE to predict the rotation noise.\n\n**Q2**: Why was the SE(3) diffusion process chosen over the standard Euclidean diffusion process? While the standard diffusion process may generate nonsensical conformations, the proposed diffusion process might not explore the entire conformation space.\n\nWe have already included results of the standard Euclidean diffusion process in our original manuscript (named as Gauss DSM in Table 1). As shown in the table below, our results suggest that SE(3) diffusion is better than Euclidean diffusion on both test sets.\n\n|   | CASF (crystal) | CASF (docked) | SAbDab (crystal) | SAbDab (docked) |\n|----------------------|----------------|---------------|------------------|-----------------|\n| Eucliean diffusion | 0.638 ± 0.017  | 0.632 ± 0.016 | 0.335 ± 0.038  | 0.330 ± 0.022   |\n| SE(3) diffusion | 0.656 ± 0.012  | 0.651 ± 0.011 | 0.361 ± 0.051 | 0.360 ± 0.051   |\n\n**Q3**: In section 4.1, is the center of mass computed based on atomic weight?\n\nFor small molecule ligands, the center of mass is set as a constant for all atom types because the atomic weight of the most common atoms are quite close (e.g., C=12, N=14, O=16). For antibodies, the center of mass is computed based on amino acid weight because they are more diverse (e.g., Glycine = 57, tryptophan = 204). To validate this modeling choice, we run our model with either constant weight or true atomic weight. As shown in the table below, the effect of using atomic weight is negligible for small molecules but quite helpful for antibodies.\n\n|  | CASF (crystal) | SAbDab (crystal) |\n|---------------|----------------|------------------|\n| Constant weight  | 0.656  | 0.340  |\n| True atomic/residue weight | 0.657  | 0.352 |\n\n**Q4**: How does the model perform in the \"Docked structure\" scenario when the docking error is high?\n\nOur model performance is quite robust to docking error. In the antibody case, the RMSD of ZDOCK is around 19.4 but the performance of NERE DSM under docked structures is quite close to crystal structures (Pearson correlation: 0.361 vs 0.360; Spearman correlation: 0.385 vs 0.363).\n\n**Q5**: Why does contrastive learning underperform denoising score matching (DSM)? Further elaboration on this would be beneficial.\n\nContrastive learning is a simpler objective than DSM because it only teaches the model to assign lower energy to crystal structures than perturbed structures. In contrast, the DSM objective not only teaches the model to assign lower energy to crystal structure, but also to reconstruct the original ligand/antibody pose from perturbed structures. It requires the model to leverage more geometric information in order to succeed in this more difficult task.\n\n**Q6**: Are there any experiments conducted on the docking benchmark of the CASF test set? It seems that the proposed method is more suitable for identifying the correct binding pose (ranking among multiple poses of a specific molecule) rather than ranking between binding poses of multiple molecules.\n\nWe agree that the docking experiment is also relevant to our approach. Our current focus is protein-ligand binding affinity prediction but we will include this experiment in our future work. Thank you for your suggestions.\n\n**Q7**: Different noise levels appear to lead to different energy functions, while only the minimal noise level can approximate the energy function of the real data distribution. How does the proposed method handle multiple noise scales?\n\nIn our implementation, we use multiple noise scales to train our model. In each training step, we first randomly sample a noise scale $\\sigma \\sim [0, 10]$ for the rotation noise and then sample a rotation vector $\\omega \\sim \\mathcal{N}_\\mathrm{SO(3)}(\\sigma)$. This dynamic sampling strategy allows us to sample noise of both small and large scales. The small-scale noise helps the model to approximate the real data distribution, while the large-scale noise allows the model to be trained on a diverse range of perturbed structures.\n"}, {"review_id": "qNAGID5gXr", "reviewer": "Reviewer_3x52", "summary": "The authors address the problem of protein-ligand binding. Specifically, the authors reformulate binding energy prediction as a generative modeling task: they train an energy-based model on a set of unlabelled protein-ligand complexes using denoising score matching and interpret its log-likelihood as binding affinity. The key contribution of this work is a  new equivariant rotation prediction network for SE(3) DSM called Neural Euler’s Rotation Equations (NERE).", "questions": "### Questions\n* I’m curious, why were only the top 50 antigen residues chosen instead of all residues within a certain distance, as done for the protein in the small molecule case?\n\n* Could you provide more details on the nature of the test splits? Are they IID w.r.t. the training set, or are they split by sequence similarity/some other criteria?\n\n###  Questions regarding protein-ligand binding:\n* Why did the authors choose to compare their model to supervised models trained on 19K labeled affinity datapoints? Did I understand correctly that the 4806 datapoints used to train NERE DSM and the 19K are datapoints used to train the supervised models are disjoint? Is there no cases where affinity data is available together with cocrystal data so both models could be trained on the same dataset?\n\n* How long does AutoDock Vina take to dock a protein-ligand pair on average (in your case)?\n\n\n###  Questions regarding Ab-Ag binding:\n* It is not entirely clear to me, are the supervised pre-trained? \n\n* I have the same question here as for protein-ligand binding, is there no cases where affinity data is available together with cocrystal data so both (unsupervised and supervised) models could be trained on the same dataset?\n\n* How long does ZDOCK take to dock a protein-ligand pair on average (in your case)?\n\n###  Questions regarding ablation studies:\n* In lines 292-294, the authors say “Interestingly, we find that the model pays the most attention to CDR-H3 and CDR-L3 residues. In most test cases, their energy is much lower than other CDR residues. This agrees with the domain knowledge that CDR-H3 and CDR-L3 residues are the major determinant of binding “. To support this claim, the authors only show one Ab-Ag example which does not guarantee that the same trend is observed across all samples in the test set. Could you report a global metric, for example, the average difference in the test set between energies in CDRs vs energies in remaining Ab regions (plus statistical significance)?\n\n### Minor comments\n* In lines 72-73, “Due to scarcity of binding affinity data, very few work has explored supervised learning for antibody binding affinity prediction.” Should be “Due to scarcity of binding affinity data, very few works have explored supervised learning for antibody binding affinity prediction.”\n*  Experiments were repeated with 5 random seeds but only the average values were reported in table 1, could you also report the standard deviation with ±?", "limitations": "In my opinion, the authors have not entirely addressed the limitations of their work, specifically, in relation to my comment in the \"Weaknesses\" section\n\n\nI set my assessment to \"borderline reject\" but my opinion can change based on the author's answers and comments", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "* The work is easy to follow and ideas are presented in a clear way\n* The proposed approach (with certain particularities) generalizes to both small and large molecules", "weaknesses": "NERE DSM relies on cocrystal structures for training and cocrystal/docked structures for binding prediction of a new protein-ligand / Ab-Ag pair, which can be limiting in terms of experimental/computational resources. In addition, NERE DSM does not predict binding affinity directly, hence the model is evaluated using correlation instead of RMSE. The utility of NERE DSM in real-world scenarios where we seek to predict affinities for a large set of protein-ligand / Ab-Ag pairs seems limited. Moreover, the extremely slow rate at which new cocrystal data becomes available is an additional limiting factor. Could the authors comment on this?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your insightful comments and suggestions. Please read our response below and let us know if you have more questions.\n\n**Q1**: NERE DSM relies on co-crystal structures for training and co-crystal/docked structures for binding prediction of a new protein-ligand / Ab-Ag pair, which can be limiting in terms of experimental/computational resources.\n\nIn terms of computational resources, modern molecular docking tools (based on neural networks) are much faster than traditional docking tools like AutoDock Vina and ZDOCK. For example, PLANET [1] is a neural network model for protein-ligand docking. On our CASF test set, it takes only about 0.1 second for PLANET to dock one ligand (on average). GeoDock [2] is another neural network model for protein-protein docking. On our SAbDab test set, it takes about 0.1 second to dock one antibody to its antigen. Using these neural network based docking models, we can screen over millions of compounds/antibodies per day. We also evaluated the performance of our model when using the structures docked by PLANET and GeoDock. As shown in the table below, there is only a small decrease in Pearson correlation but the model remains accurate. Therefore, we argue that computational resources are not limiting factors to our unsupervised binding energy prediction framework.\n\n|        | Slow Docking (Vina/ZDOCK) | Ultrafast Docking (PLANET/GeoDock) |\n|----|---|----|\n| CASF | 0.656 ± 0.012 | 0.616 ± 0.009 |\n| SAbDab | 0.360 ± 0.051 | 0.350 ± 0.069  |\n\nIn terms of experimental resources, our method requires co-crystal structures for training. We agree that it can be a limiting factor but we can alleviate this issue by utilizing a larger set of protein-protein complexes to pre-train our model. There are 208,347 co-crystal structures in the Protein Data Bank in total and its growth is about 10,000 per year. We will explore this direction in our future work.\n\n**Q2**: Docking speed of AutoDock Vina, ZDOCK, and other docking models.\n\nThe docking speed of AutoDock Vina and ZDOCK is 25 second and 20 second per complex (using a 64-CPU server). The docking speed is much faster for neural network docking tools like GeoDock and PLANET, with approximately 0.1 second per complex.\n\n**Q3**: Are there no cases where affinity data is available together with co-crystal data so both supervised and unsupervised models could be trained on the same dataset?\n\nFor antibodies, there are only 566 co-crystal structures that have affinity labels. They are all included in the test set and could not be used to train our model. For small molecules, the co-crystal structures in our training set actually come with affinity labels, though we did not use the affinity labels to train our model. For a fair comparison, we trained a supervised model on the same training set and model architecture as NERE DSM, but using affinity labels and a regression training objective. This supervised model achieves an average Pearson correlation of 0.734 ± 0.012 on the CASF test set, which is lower than TankBind or IGN because it uses less training data. This supervised model is better than NERE DSM but the margin is relatively small (0.734 vs 0.656). In this case, our unsupervised method is able to recover nearly 90% of the supervised model performance.\n\n**Q4**: Why were only the top 50 antigen residues chosen as epitopes instead of all residues within a certain distance, as done for the protein in the small molecule case?\n\nWe found that the latter strategy yielded a lower performance (Pearson R = 0.317) because the number of epitope residues are quite uneven when using a distance threshold. The performance becomes better when using a fixed number of epitope residues.\n\n**Q5**: Details on the nature of the test splits\n\nThe current train/test split is based on IID, but we filtered the training set so that the same ligand/antibody/antigen does not appear in both training and test sets. To further study the impact of train/test split strategy, we also train our model with sequence similarity split. In this case, we remove all training instances whose proteins/antibodies/antigens have 40% similarity to the test set. As shown in the following table, we find that the average performance of NERE is not impacted by similarity split (0.508 vs 0.503), which confirms that our model is not making predictions simply based on similarity.\n\n|        | Original split | Similarity split |\n|----|----|-----|\n| CASF (Crystal) | 0.656 ± 0.012 | 0.634 ± 0.009 |\n| SAbDab (Crystal) | 0.360 ± 0.051 | 0.372 ± 0.038 |\n| Average | 0.508 | 0.503 |\n\n**Q6**: Could you report a global metric, for example, the average difference in the test set between the energy of CDR3 residues versus non-CDR3 residues (plus statistical significance)?\n\nThe average difference between the energy of a CDR3 residue and a non-CDR3 residue is -2.18, with a p-value of 5.8e-60 under student t-test. From the rebuttal Figure 1g, we can clearly see that the average difference is shifted towards negative range. Therefore, we conclude that the model pays more attention to CDR3 residues than non-CDR3 residues.\n\n**Q7**: Experiments were repeated with 5 random seeds but only the average values were reported in table 1, could you also report the standard deviation with ±?\n\nThe standard deviation is already reported in table 1 as subscripts (e.g., $0.656_{.011}$ means 0.656 ± 0.011).\n\n**Q8**: In the antibody-antigen case, are the supervised model pre-trained?\n\nYes, the supervised $\\mathrm{FANN_{transfer}}$ model is pre-trained on the SKEMPI protein-protein binding database, which contains approximately 6000 binding affinity data points. This model also uses the pre-trained ESM-2 protein language model embedding for residue features.\n\nReference\n1. Zhang et al. Planet: A multi-objective graph neural network model for protein-ligand binding affinity prediction. bioRxiv 2023\n2. Chu et al. Flexible Protein-Protein Docking with a Multi-Track Iterative Transformer. bioRxiv 2023"}, {"review_id": "wHI980Hefv", "reviewer": "Reviewer_KYrm", "summary": "The authors propose an energy-based model for unsupervised binding affinity estimation, which is trained by SE(3) denoising score matching (DSM). Different from standard DSM, they add noise by random rotations and translations on the ligand. Utilizing Euler's rotation equations, the rotation matrix can be derived from predicted forces on each atom, making it possible to directly supervise on the predicted forces instead of the rotation matrix. Experiments on protein-ligand binding and antibody-antigen binding demonstrate the superiority of the proposed method over existing unsupervised method.", "questions": "1. Under the docked setting, the RMSD is quite large for the generated antibody-antigen complexes (median 19.4), however, the correlation is similar to the crystal setting. Does that mean the model actually did not learn useful geometric interaction information? Otherwise, with such badly docked structures, the model should have much worse performance because the structures can not provide accurate geometric information.\n2. Will the correlation be improved if you first denoise the input structure (regardless whether it is crystal structure or docked structure) with the model, then output the energy from the denoised structure?", "limitations": "The authors have mentioned that the limitation might be only considering the alpha carbon in the antibody-antigen setting. Perturbations on the side chains can be designed with a similar strategy in the future.", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "1. Using Euler's rotation equations to formulate the supervision on the rotations as supervision on the predicted atomic forces is novel and clever.\n2. Unsupervised learning on binding affinity is important and meaningful as precise affinity labels are hard to obtain in practical.\n3. The docked setting is interesting and reals the potential applications of the proposed method in real scenarios where crystal structure is often hard to obtain.\n4. It's very interesting to see that pair-wise energies are mostly lower at CDR-H3 and CDR-L3, which perfectly aligns with the domain knowledge.", "weaknesses": "1. The learnt probability distribution $p(A, X)$ might not be well aligned with the actual $p(A, X)$ because the actual one is more complicated than the prior distribution of rotation and translation noises. In this sense, for the learnt $p(A, X)$, $p(A, X) \\propto \\exp(-E(A, X))$ might not hold.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your insightful comments and suggestions. Please read our response below and let us know if you have more questions.\n\n**Q1**: Under the docked setting, the RMSD is quite large for the generated antibody-antigen complexes (median 19.4), however, the correlation is similar to the crystal setting. Does that mean the model actually did not learn useful geometric interaction information?\n\nThe reason why our model is robust to docking error is that we use a very loose threshold for residue contact. In NERE DSM, the predicted binding energy is a sum of pairwise residue interaction energy: $E(A, X) = \\sum_{i,j: d_{ij} < 20} f_o(h_i, h_j)$. In this formula, we consider two residues to be interacting (a.k.a., forming a contact) if their $C_\\alpha$ distance $d_{ij}<20$. This threshold is much larger than standard practice (5-8Å). In fact, we find that the model becomes a lot more sensitive to docking error if we set the contact threshold to 10Å or lower. With a tighter contact threshold, the true contacts are likely to be lost in the docked setting and the model cannot extract useful geometric information for binding energy prediction.\n\n**Q2**: Will the correlation be improved if you first denoise the input structure (regardless whether it is crystal structure or docked structure) with the model, then output the energy from the denoised structure?\n\nWe found that the performance is roughly the same with or without this denoising step at inference time (see the table below).\n\n|   | CASF (crystal) | CASF (docked) | SAbDab (crystal) | SAbDab (docked) |\n|----------------------|----------------|---------------|------------------|-----------------|\n| without denoising | 0.656 ± 0.012 | 0.651 ± 0.011 | 0.361 ± 0.051 | 0.360 ± 0.051 |\n| with denoising | 0.655 ± 0.013 | 0.649 ± 0.012 | 0.367 ± 0.029  | 0.359 ± 0.043 |\n"}, {"review_id": "eovh7jbh89", "reviewer": "Reviewer_UbeA", "summary": "In this paper, the authors developed an energy-based model for unsupervised binding affinity prediction. The energy-based model was trained under SE(3) denoising score matching where the rotation score was predicted by Neural Euler’s Rotation Equation. Experiments on protein-ligand binding and antibody-antigen binding are conducted. ", "questions": "please see the weakness.", "limitations": "The limitations are well discussed.", "rating": 4, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "1. The paper is well-written and easy to follow.\n2. The proposed method outperforms all unsupervised baselines and supervised baselines in the antibody case.\n3. The code and data are provided.\n4. Systematic ablation studies are performed to show the effectiveness of modules.\n\n", "weaknesses": "1. The idea of using unsupervised generative models for binding affinity prediction is not new. The authors fail to discuss or compare with related works [1,2]\n\n[1] Luo et al., Rotamer Density Estimator is an Unsupervised Learner of the Effect of Mutations on Protein-Protein Interaction, ICLR 23  \n[2] Guan et al., 3D Equivariant Diffusion for Target-Aware Molecule Generation and Affinity Prediction, ICLR 23  \n\n2. The authors only calculate the Pearson correlation coefficient. For the binding affinity prediction, other metrics such as RMSE, MAE, AUROC, and Pearson correlation coefficient.\n\n3. The authors fail to consider the flexibility of proteins and ligands. The sidechains of proteins are also ignored. However, these factors are quite important for protein-ligand/protein binding.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your insightful comments and suggestions. Please read our response below and let us know if you have more questions.\n\n**Q1**: The authors fail to consider the flexibility of proteins and ligands. The side-chains of proteins are also ignored. However, these factors are quite important for protein-ligand and protein binding.\n\nWe agree that it is beneficial to model the flexibility of ligands and incorporate the side-chain information of proteins. Given the limited time, we partially address this issue by including all protein atoms (including side-chains) as input instead of only $C_\\alpha$ atoms. This modification only changes the input of our model while the model architecture and NERE DSM objective remains the same. We refer to this extended model as NERE DSM (all-atom). Our results (rebuttal Figure 1e) suggest that including side-chain atoms substantially improves the performance of our method on some of the datasets. Please read our global response above  for details. In the future, we plan to model the flexibility of rotatable bonds and protein side-chain atoms by introducing a side-chain DSM loss, which adds rotation noise on the rotatable bonds and side-chain atoms and use NERE to predict the rotation noise.\n\n**Q2**: The authors only calculate the Pearson correlation coefficient. For the binding affinity prediction, including other metrics such as RMSE, MAE, and Spearman correlation coefficient would be helpful.\n\nWe report Spearman correlation and p-values as additional metrics in rebuttal Table 1. Our model (NERE DSM) remains the best under all three metrics. Unfortunately, it was not possible to report root mean square error (RMSE) or mean absolute error (MAE) because our model does not predict absolute affinity values. Shifting or scaling $E(A, X)$ by any constant will be equally optimal under the DSM objective.\n\n**Q3**: The idea of using unsupervised generative models for binding affinity prediction is not new. The authors fail to discuss or compare with related works [1,2].\n\nThank you for suggesting these papers and they are quite related. Luo et al. [1] proposed a flow-based generative model to estimate the probability distribution of protein side-chain conformations. Their unsupervised model (named RDE-Linear) uses the learned entropy of side-chains to predict $\\Delta\\Delta G$ of protein mutations. RDE-Linear and NERE DSM are complementary to each other because they consider different degrees of freedom (DOF). RDE-Linear considers the DOF of side-chains and the backbone structure is fixed. NERE DSM considers the DOF of backbone structure (i.e., different docking angles) but not the DOF of side-chains.\n\nTo compare NERE DSM with RDE-Linear, we train our model (the all-atom version) on approximately 27000 non-redundant protein-protein complexes (downloaded from PDB) and evaluate it on the entire SKEMPI mutation effect prediction benchmark. For each downloaded complex, we decompose it into pairs of two chains and remove any pairs of chains with buried surface area less than 500. We then clustered these proteins and took one representative complex from each cluster. In each training step, we randomly rotate one of the proteins in a complex and update the model with our NERE DSM objective. The input to our model includes side-chain atoms because the downstream task is mutation effect prediction.\n\nAt test time, given a pair of wild type and mutated protein complexes, our model first predicts their binding energy $E_\\mathrm{wildtype}$ and $E_\\mathrm{mutated}$ and calculates $\\Delta\\Delta G=E_\\mathrm{mutated} - E_\\mathrm{wildtype}$. Following Luo et al. [1], we report per-structure and overall Pearson and Spearman correlation between predicted and experimental $\\Delta\\Delta G$. To compute per-structure correlation, we group mutations by structure, discard groups with less than 10 mutation data points, and calculate correlations for each structure separately. As shown in the table below, we find that NERE DSM achieves a much higher performance than RDE-Linear on three out of four metrics. We visualize the overall correlation between predicted and experimental $\\Delta\\Delta G$ in rebuttal Figure 1h.\n\n| SKEMPI | Pearson (per structure) | Spearman (per structure) | Pearson (overall) | Spearman (overall) |\n|----|----|---|----|---|\n| RDE-Linear | 0.290 | 0.263 | **0.419** | 0.3514 | \n| NERE DSM (all-atom) | **0.388±0.020**  | **0.402±0.035** | 0.399±0.039 | **0.402±0.013** |\n\nGuan et al. [2] proposed TargetDiff, a diffusion model for 3D small molecule design. Similar to our method, TargetDiff is trained on cocrystal structures of protein-ligand complexes in an unsupervised manner. They demonstrated that the unsupervised features learned by TargetDiff could improve supervised affinity prediction. Specifically, they augmented EGNN [3] with the features learned by TargetDiff and trained the augmented model (EGNN + TargetDiff) on the PDBBind training set. To compare our method with TargetDiff, we followed their training/validation/test split and fine-tuned the pre-trained NERE DSM model on the PDBBind training set. As shown in the table below, our method outperforms the EGNN + TargetDiff baseline, which demonstrates the advantage of our approach. The baseline results are copied from Guan et al. 2023.\n\n| TargetDiff test set  | Pearson | Spearman |\n|----------------------|----------------|---------------|\n| TransCPI | 0.576 | 0.540 |\n| MONN | 0.624 | 0.589 |\n| IGN | 0.698 | 0.641 |\n| EGNN | 0.648 | 0.598 |\n| EGNN + TargetDiff  | 0.680 | 0.637 |\n| NERE DSM (pretrained) | 0.500±0.004 | 0.516±0.013 |\n| NERE DSM (fine-tuned) | **0.703±0.007**  | **0.656±0.017** |\n\nReferences\n1. Luo et al., Rotamer Density Estimator is an Unsupervised Learner of the Effect of Mutations on Protein-Protein Interaction, ICLR 23\n2. Guan et al., 3D Equivariant Diffusion for Target-Aware Molecule Generation and Affinity Prediction, ICLR 23\n3. Satorras et al., E(n) equivariant graph neural networks. ICML 2021"}]

Author rebuttal:
We want to thank all reviewers for their valuable comments and suggestions. We would like to summarize three main results here in response to some common questions/suggestions. The results are included in the attached PDF file (rebuttal Table 1 and Figure 1a-h).

### Additional metrics

As suggested by Reviewer X8hw and UbeA, we report Spearman correlation and p-values as additional metrics in rebuttal Table 1. We could not report root mean square error (RMSE) or mean absolute error (MAE) because our model does not predict absolute affinity values. Shifting or scaling $E(A, X)$ by any constant will be equally optimal under the DSM objective.

### Additional benchmark datasets

We include two additional benchmark datasets to incorporate reviewer X8hw’s suggestions. First, we adopt an HER2 antibody affinity maturation dataset [1] to investigate the utility of NERE DSM as a pre-training strategy. This dataset has 422 antibody sequences with experimental binding affinities against HER2. We report the Spearman correlation between the predicted binding energy and experimental binding affinity on this dataset. To investigate whether the learned representations from the unsupervised approach can be beneficial for supervised downstream tasks, we randomly sample a small number of HER2 data to fine-tune the pre-trained NERE model. Results on this dataset are shown in rebuttal Figure 1f. When fine-tuned on the same number of HER2 data, we find that representations learned by NERE outperform representations learned under supervised pre-training ($\mathrm{FANN_{transfer}}$).

Second, we adopt a free energy perturbation (FEP) benchmark [2] (developed by Merck) to investigate the utility of unsupervised learning in the context of protein-ligand binding. It has eight protein targets (cdk8, cmet, eg5, hif2a, pfkfb3, shp2, syk, and tnks2) and 264 ligands in total. We calculate the Spearman correlation ($R_S$) between the predicted binding energy and experimental binding affinity for each of the eight targets and report the average correlation for each method. Results on this dataset are shown in rebuttal Figure 1b-d. Our main finding is that NERE outperforms supervised baselines on this FEP test set. These results suggest that unsupervised models can be also useful for protein-ligand binding even though there is a fair amount of labeled affinity data. The experimental setup for these two experiments are detailed in our response to reviewer X8hw.

### Including side-chain atoms into the model

Several reviewers have raised concerns that the current model does not model the flexibility of proteins/ligands and the side chains of proteins. Given the limited time, we partially address these concerns by including all protein atoms (backbone + side chains) as input instead of only $C_\alpha$ atoms. This modification only changes the input of our model while the model architecture and NERE DSM objective remains the same. We refer to this extended model as NERE DSM (all-atom) and the original model as NERE DSM ($C_\alpha$-only).

As shown in rebuttal Figure 1e, we found that adding side-chain information gave a similar performance on the original CASF and SAbDab test sets, but substantially improved model performance on the new FEP and HER2 datasets. On the FEP test set, the all-atom model achieved $R_S=0.380$ while the $C_\alpha$-only model had $R_S=0.192$. On the HER2 test set, the all-atom model achieved $R_S=0.230$ while the $C_\alpha$-only model had only $R_S=0.043$. Overall, these results suggest that incorporating side-chain information can be beneficial for protein-ligand/antibody-antigen binding energy prediction. In the future, we plan to model the flexibility of rotatable bonds and protein side-chain atoms by introducing a side-chain DSM loss, which adds rotation noise on the rotatable bonds and side-chain atoms and use NERE to predict the rotation noise.

### References
1. Shanehsazzadeh et al., Unlocking de novo antibody design with generative artificial intelligence, biorxiv 2023
2. Schindler et al., Large-scale assessment of binding free energy calculations in active drug discovery projects, 2020

Meta-review: The paper proposes an unsupervised method for binding classification using denoising score matching. I like this paper. The problem considered is a pertinent one within drug and protein design. The paper is well written, the results are strong and the approach is novel. The additional numerical results presented during the rebuttal period further corroborate the claims of good performance.
Decision: Accept (poster)