Learning Harmonic Molecular Representations on Riemannian Manifold
Authors: Yiqun Wang
Yuning Shen
Shi Chen
Lihao Wang
Fei YE
Hao Zhou
Venue: ICLR 2023 poster (2023)
Keywords: Riemannian manifold
molecular surface
harmonic analysis
functional map
binding site prediction
rigid protein docking

Abstract: Molecular representation learning plays a crucial role in AI-assisted drug discovery research. Encoding 3D molecular structures through Euclidean neural networks has become the prevailing method in the geometric deep learning community. However, the equivariance constraints and message passing in Euclidean space may limit the network expressive power. In this work, we propose a Harmonic Molecular Representation learning (HMR) framework, which represents a molecule using the Laplace-Beltrami eigenfunctions of the molecular surface. HMR offers a multi-resolution representation of molecular geometric and chemical properties on 2D Riemannian manifold. We also introduce a harmonic message passing method to realize efficient spectral message passing over the surface manifold for better molecular encoding. Our proposed method shows comparable predictive power to current models in small molecule property prediction, and outperforms the state-of-the-art deep learning models for the rigid protein docking challenge, demonstrating its versatility in molecular representation learning.

Published as a conference paper at ICLR 2023
LEARNING HARMONIC MOLECULAR REPRESENTA- TIONS ON RIEMANNIAN MANIFOLD
Yiqun Wang', Yuning Shen', Shi Chen"; Lihao Wang', Fei Ye', Hao Zhou
BByteDance Research, UUniversity of Wisconsin-Madison, Institute for AI Industry Research (AIR), Tsinghua University
(yiqun. wang, yuning. shen, wangl 1217, yefei, joyce))byttderrce.com schen636@wise edu, zhouhao@air tsinguuu en
ABSTRACT
Molecular representation learning plays a crucial role in Al-assisted drug discov- ery research. Encoding 3D molecular structures through Euclidean neural net- works has become the prevailing method in the geometric deep learning commu- nity. However, the equivariance constraints and message passing in Euclidean space may limit the network expressive power. In this work, we propose a Har- monic Molecular Representation learning (HMR) framework, which represents a molecule using the Saplace-Beltramy eigenfunctions of its molecular surface. HMR offers a multi-resolution representation of molecular geometric and chem- ical features on 2D Riemannian manifold. We also introduce harmonic mes- sage passing method realize efficient spectral message passing over the surface manifold for better molecular encoding. Our proposed method shows compara- ble predictive power to current models small molecule property prediction, and outperforms the state-of-the the-art deep learning models for ligand-binding protein pocket classification and the rigid protein docking challenge, demonstrating its versatility in molecular representation learning.
INTRODUCTION
Molecular representation learning is a fundamental step in Al-assisted drug discovery. Obtaining good molecular representations is crucial for the success of downstream applications including pro- tein function prediction (Gligorijevćć et al.) 202T) and molecular matching, e.g., protein-protein docking (Ganea al. 2021). In general, ideal molecular representation should well integrate both geometric (e.g.. 3D conformation) and chemical information (c.g., electrostatic potential). Ad- ditionally, such representation should capture features in various resolutions accommodate dif- ferent tasks, e.g., high-level holistic features for molecular property prediction, and fine- fieeggained features for describing whether two proteins can bind together at certain interfaces.
Recently, geometric deep learning (GDL) (Bronstein et al. 2017 2021 Monti et al. 2017 has been widely used in learning molecular representations (Atz et al. captures necessary information by performing neural message passing (NMP) on common struc- tures such as 2D/3D molecular graph (Klicpera et al. 2020 Stokes et al. 2020). 3D voxel (Liu et al. 2021 Land point cloud (Unke et 2001). Specifically, GDL encodes: a) geometric features by modeling atomic positions in the Euclidean space, and b) chemical features by feeding atomic information into the message passing networks. High-level features could then be obtained by ag- gregating these atom-level features, which has shown promising results empirically.
2021 Townshend et al. (02T). GDL
However, we argue that current molecular representations via NMP in the Euclidean space is not necessarily the optimal solution, which suffers from several drawbacks. First, current GDL ap- proaches need employ equivariant networks (Thomas et al. (2018) to guarantee that the molecular representations transform accordingly upon rotation and translation (Fuchs et al. 2020), which could undermine the network expressive power (Cohen et al. 2018 Li et al. (021). Therefore, developing a representation that could properly encode 3D molecular structure while bypassing the equivari- ance requirement is desirable. Second, current molecular representations in GDL are learned
"This work was conducted during internship at ByteDance Research.

---

Published as a conference paper at ICLR 2023
bottom-up manner, which are hardly able to provide features in different resolutions for different tasks. Specifically, NMP in Euclidean space typically achieves long-range communication between distant atoms by stacking deep layers or increasing the neighborhood radius. This would hinder the effective representation of macromolecules with tens of thousands of atoms (Battiston et al.) 2020, Boguna et al. 202T). To remedy this, residue-level graph representations are commonly used for large molecules (Jumper et al., 2021] Gligorijević al., 2021). Hence designing efficient multi- resolution message passing mechanisms would be ideal for encoding molecules with distinct sizes.
On the other hand, the molecular surface a high-level representation of molecule's shape, which has been widely used to study intermmolecula interactions (Richards, 1977; Shulman-Peleg et al.) 2004). Intuitively, the interaction between molecules commonly described as a "key-lock pair", where both shape complementarity (Li et al. (013) and chemical interactions (e.g., hydrogen bond) determine whether the key matches the lock molecule. has been shown that the molecular surface holds key information about inter-molecula interactions (Gainza et al. 2020), which makes it an ideal candidate for molecular representation learning (Sverrisson et al.
2021 Somnath et al. 2021
Inspired by the idea of Shape- DNA (Reuter et al.] (2066), hereby we propose Harmonic Molecular Represen- tation learning (HMR) by utilizing the Laplace-Belteer Beltrami cigenfunctions on the molecular surface (a 2D manifold). Our representation has the following advan- tages: a) HMR works on 2D Riemannian manifold instead of in the 3D Euclidean space, thus the resulting molecular rep- resentation is by design roto-translation invariant; b) HMR represents a molecule in top-down manner, and is capable of offering multi-resolution features that accommodate various target molecules (i.c.. from small molecules to large proteins), thanks the smooth nature of the Laplace-Beltrami eigenfunctions (Fig. I): c) HMR naturally integrates geometric and chemical features the molecular shape defines the Riemannian manifold (i.e., the underlying domain equipped with a metric), and the atomic configurations determine the associated functions distributed on the manifold (e.g., electrostatics).
reduced spatial resolution
truncated Laplace- Beltrami eigenfunctions
Figure 1: Multi-resolution molecular surface represen- tation. Showing the electrostatic potential (blue regions being negatively charged, PDB ID: 3V6F) at different resolutions under our representation. See (Appendix B for technical details about tuning resolution.
To demonstrate that HMR is generally applicable to different downstream tasks including molecu- lar property prediction and molecular matching. we propose two specific techniques: (1) manifold harmonic message passing for realizing holistic molecular representations, and (2) learning regional functional correspondence for molecular surface matching. Without loss of generality, we apply the proposed techniques to solve three drug discovery- problems: QM9 small molecule property regression, bogand-binding protein pocket classification, and rigid protein docking pose prediction. Our proposed method shows comparable performance for small molecule property prediction to NMP-based models, while outperforming the state-of-the-art deep learning models in protein pocket classification and the rigid protein docking challenge.
RELATED WORK
Molecular Surface Representation The molecular surface representation is commonly adopted for tasks involving molecular interfaces (Duhovny et al.] (002). where non-covalent interactions (e.g., hydrophobic interactions) play decisive role (Sharp) 1994p. Non-Euclidean convolutional neural networks (Monti et al. 2017) and point cloud-based learning models (Sverrisson et al.) 2022) have been applied to encode the molecular surface for downstream applications, e.g., protein bunding site prediction (Mylonas et al.) 2021). However, existing methods apply filters with fixed sizes and highly dependent the surface mesh quality, which limit the expressive power for molecular shape representation across different spatial scales (Somnath et al. 2021 Isert et al. 2022).

---

Published as a conference paper at ICLR 2023
Geometry Processing Our work was inspired by some pioneering work for shape encoding and shape matching in the geometry processing research community (Litman & Bronstein 2013 Bi- asotti et al. 2016). The surface of a 3D object is typically discretized into a polygon mesh with vertices and faces. Intrinsic properties of the surface manifold used encode the shape (Sun eet al. 2009) Broostein & Kokkinos] 2010). Functional maps (Ovsjanikov et al. [012 Litany et al.)
2017b) have been proposed establish spectral-space functional correspondence between two man-
ifolds. Recently, deep learning has been applied to learn representative features to facilitate shape
recognition and matching (Litany al. 2017a] Donati al. 2020 Attaiki al. 202T).
Spectral Message Passing Our proposed method decomposes surface functions/features as the linear combination of some basis functions (hence the name "harmonic") and realizes message pass- ing by applying various spectral filters. Graph convolutional network (GCN) is closely related to our proposed method, which operates the graph Laplacian eigenspace (Kipf Welling 2016) [Shen] et al. (02T). The major difference is that graph is discrete by construction, while our method works with continuous Riemannian manifold (Coifman & Lafon 2006). other words, the underlying manifold and its spectrum remain the same with different surface discretizations, hence is robust representation of the surface shape (Coifman et al. 2005). See "Shape-DNo" in Sec. 3
3 PRELIMINARIES
The goal of this work is to propose representation learning method using the molecular surface, which could properly encode both geometric and chemical features. Intuitively, the molecular sur- face defines the shape of molecule in 3D) Fuclidean space (i.e., geometry), while chemical features (e.g., hydrophobicity) can be treated as functions distributed on the molecular surface. The surface and these associated functions co-determine the underlying molecular properties, e.g., whether an antibody could bind with an antigen. Therefore, we propose to represent a molecule as a set of (learndd) functions/features defined its surface.
Before moving on, we first explain some basic concepts behind our proposed representation learning introduce molecular "Shape-DNA" and set of basis functions inherent
framework. In [.e..... to the surface manifold. In Seec. 3.2] we illustrate how to apply harmonic analysis to decompose function (defined on the molecular surface) into the linear combination of its "Shape-DNA" ba- sis functions, as shown in Fig. 2] These concepts enable us to represent a molecular surface as 2D Riemannian manifold with associated functions, which form the cornerstone of our proposed framework.
Co1
+ Colo
+ Cotoaa
electrostatic potential
Φ1
010/
1000
Figure 2: Illustrating manifold harmonic analysis. Left-hand side: the simulated electrostatic poten- tial function on the protein surface. Right-hand side: the linear combination of Saplace-Beltrame eigenfunctions (ø1...., .10....) with corresponding coefficients (Co,.... Co Only few selected terms are explicitly shown from the infinite sum. Note that these eigenfunctions exhibit different spatial frequencies (resolutions).
3.1 THE SHAPE-DNA
The molecular surface can be viewed as 2D Riemannian manifold (M), which adopts a discrete set of eigenfunctions (x) that solves
Ao. Xoo i= 0,
(I)
Here, A is the Laplace-Beltrami (LB) operator acting surface scalar fields, defined as Aʃ div V f)).. #o.... are a set of orthonormal eigenfunctions (i.c., (Øx, ;,) = ōtz). And 0 = dn S \1 ... are the corresponding cigenvalues. These set of cigenvalues {A:] } of the LB operator
For instance, an isosurface of its electron density field, or the solvem-aaceessileleeccsssss surface, etc.

---

Published as a conference paper at ICLR 2023
are called the "Shape-DNa and their corresponding eigenfunetions { are unique to each shape. In other words, different molecules adopt different LB eigenfunctions. Notably, two different 3D conformations of the same molecule may not have the same LB eigenfunctions, but a pair of chiral molecules (i.e., mirror image of each other) share the same LB spectrum.
Fig. shows a few selected eigenfunctions (01, 610. 1000) of a protein surface manifold on the right-hand side. These eigenfunctions are intrinsic properties of the surface manifold, which remain invariant under rigid transformations to the molecule. The eigenvalues X; reflect the surface Dirich- let energy (defined (Af, M),, which measures the smoothness of eigenfunction &i over M. The eigenvalues increase linearly, whose slope is roughly inversely proportional to the surface area (known as Weyl's asymptotic law). See Seeenndp Afro more details about the "Shape-DNA"..
3.2 BASICS OF MANIFOLD HARMONIC ANALYSIS
Now, we introduce basic manifold harmonic analysis, merit of using the Riemannian manifold representation. "Harmonic analysis" refers to the representation of functions as the superposition of some basic waves. Specifically in our case, given the molecular surface manifold M and its LB eigenfunctions {or];-o, any scalar-valued function f that square-integrable on M can be decomposed into generalized Fourier series:
f(x) - (f.))
(2)
In other words, I can be represented as the linear combination of the LB eigenfunctions. The linear combination coefficient (J. M can be interpreted as the "projection" of onto the eigenfunction 6,, which reflects the contribution of this particular cigenfunction to synthesizing function f.
Interestingly. it is easy to notice that these eigenfunctions display different spatial frequencies Fig. 2). i.c., di varies slowly over the surface, while 1000 oscillates at a much higher frequency. This is analogous to the set of sin(k:)) functions in the ID case, where high frequency waves with larger & values exhibit more oscillations within the period of length 2x. Hereafter we refer to LB eigenfunctions with (relatively) small/large eigenvalues as the low/high frequency components.
Therefore, by manipulating the lincar combination coefficients, we could control the contribution of different frequency components to synthesizing function f. For instance, with only low frequency components, the synthesized function will have a lower spatial resolution on the surface (i.e., more smoothed out, Fig. riight panel, see/Appe Boro technical details.), and vice versa. Synthesizing a new function f with some selected frequency components is commonly known as wave filtering in signal processing, which will be applied in our representation learning framework to realize multi- resolution encoding of the molecular surface. We refer the readers to some excellent review papers for more details about geometry processing (Bronstein al., 2017 Rosenberg 1997).
METHODOLOGY
We realize that the aforementioned geometry processing methods suitable for molecular surface representation. However, a clear distinction between the 3D objects used in geometry processing research and the molecular systems is that molecules are not simply shapes their underlying atomic structures beneath the surface govern the molecular functionality. In other words, both geometry (i.e., shape) and chemistry have to be considered for molecular encoding. It has been shown that the molecular surface displays chemical and geometric patterns which fingerprint protein's 's mode of interactions with other biomolecules (Gainza al.) (2020). Therefore, we formulate both geometric and chemical properties of molecule as functions distributed on its surface manifold.
Then comes the question: how do we properly learn these geometric and chemical features on the molecular surface? One viable solution is to emulate the message passing framework commonly used GNNs, whose goal is to propagate information between distant surface regions to encode surface features at different scales. To that end, we present two methods under the HMR frame- work. In Sec. 4.1, we introduce manifold harmonic message passing, which makes use of harmonic
Two shapes may share the same I. eigendecomposition (isometries). The uniqueness eigendecompo- sition up to isometries is associated with the question "Can One Hear the Shape of a Drum?" (Kac) 1966)

---

Published as a conference paper at ICLR 2023
Propagated Functions
Feature Extraction
frequency filters
Aut
heat operator
(A.0)
Agt
Eigenfunction Computation
FmP
+
Agt
00ncet
Dxg
Forward
Feed
Aut
Manifold Harmonic Message Passing
Figure 3: HMR workflow. Given molecular surface mesh with vertices, we compute the first Laplace-Beltrami eigenfunctions (column-wiss stacked into an array h, and क denotes its Moore- Penrose pseudo-inverse, see Appendix BJ for discrete calculations) with ascending eigenvalues, and extract initial surface features FmuP through MLP. Then, we apply neural network-learned spectral filters to propagate the features over the surface to achieve message passing. Note that each feature channel has unique Gaussian frequency filter and propagation time t. Relevant tensor sizes indicated in parentheses. Multiple message passing blocks can be stacked for better representations.
analysis techniques to allow efficient multi-range communication between regions on the molecu- lar surface regardless of its size. This enables HMR to encode information within a molecule and be applied to molecule-level prediction tasks. In [Sec. 4.2 we propose HMR-powered rigid pro- tein docking pipeline. We use this docking challenge to demonstrate the potential of our proposed representation learning method for modeling interactions between large biomolecules.
4.1 LEARNING HARMONIC MOLECULAR REPRESENTATIONS
Given the 3D atomic structure molecule, HMR returns (1) a discretized molecular surface man- ifold (i.e., triangle mesh) with N vertices (XI,...,XX) .,) C R3 and the corresponding faces; (2) set of n per-vertex features, E R N.Nx which can viewed as n learned surface functions. The discretized surface manifold well these features represent the geometry and chemistry of the underlying molecule, and can be used for various downstream prediction tasks.
Surface Preparation We use MSMS (Ewing & Hermisson, 2010) to compute the molecular solvent-excluded surface as a triangle mesh with N vertices. Then, we compute the first k LB eigenfunctions { k-1 with ascending eigenvalues as described in Reuter et al. (2009), and stack i20 them into an array @ E IR Nxk where each column stores eigenfunction.
Geometric features can be readily calculated given the surface mesh. We compute the per-vertex mean curvature, Gaussian curvature, and the Heat Kernel Signatures as described in Sun et al. (2099. Local chemical environment captured using a simple multilayer perceptron (MLP). For each vertex, we encode its neighboring atoms within predefined radius (e.g.. 6 σ) through MLP, then sum over the neighbors to obtain its chemical embedding. We use another MLP to combine the per-vertex initial features Fiup + MLP (coneat (Fgeom Fchem)) )) Fimp Nx These n. features reflect the local geometric and chemical environment of each surface vertex, which will be used as input to the harmonic message passing module. See more implementation details in Appenddix C
The output of the surface preparation module includes (1) the molecular surface triangle mesh, (2) the surface Caplace-Beltrante eigenfunctions 4 and (3) the per-vertex features Fimp.
Harmonic Message Passing Our proposed harmonic message passing mechanism is closely re- lated the heat diffusion process on an arbitrary surface. Joseph Fourier developed spectral analysis methods to solve the heat equation 0f + A = 0, where f is some heat distributed on the sur- face. This concise partial differential equation describes how heat distribution evolves over time, for
whose solution can be expressed using the heat operator exp(-A1), A/), i.e., s(() initial heat distribution fo at f 0. Intuitively, heat will flow from hot regions to cool regions on the surface. As time approaches infinity, the heat distribution f will converge to constant value (i.e., the global average temperature on the surface), assuming that total energy is conserved.
In fact, heat diffusion be thought of message passing process, where surface regions with different temperatures communicate with each other and propagate the initial heat distribution de- terministically. The heat exchange rate is dependent the difference in temperature (determined by the LB operator), while the message passing distance is determined by the heat diffusion time I.

---

Published as a conference paper at ICLR 2023
Following this idea, we generalize the heat diffusion process by proposing a function propagation operator P with neural network-learne frequency filter (a()):
pj -==Eon (f(A))ffoo)mm)
(3)
11)= exp(-Xt) where (,, 0²
(A(A) = exp
(4)
As shown in Eq.. 3] the input function f is first expanded as the linear combination of the LB eigen- functions with coefficients (f. M- We then learn spectral-space frequency filter Fo (A), which is a function of the corresponding eigenvalue A. Here we abuse the usage of frequency. which actu- ally refers the LB eigenvalues. The filter (A) consists of components (Eq. 4) a Gaussian frequency filter (with parameters and σ), and the heat operator part e (with parameter !). The matrix representation of the function propagation operator plotted Fig.
For each input function f within Finp, the neural net- work learns unique set of parameters (H,, a. t) through backpropagation. In Fig. 4] we showcase a few exam- ples of surface functions obtained by applying different frequency filters an initial delta function.
All-frequency
The Gaussian frequency filter allows the network to prop- agate the input function along some selected eigenfunc- tions with eigenvalues close to μa, while the number of selected eigenfunctions is determined by the filter width a. Owing to the multi-resolution nature of the LB eigen- functions with different spatial frequencies, this filter will help capture surface functions at different resolutions.
High-frequency
Low-frequency
The heat operator part governs the communication dis- tance. With longer propagation time, function ʃ will be more averaged towards the global mean, leading smoothed function. In addition, the heat operator is by definition low-pass filter, where components with larger eigenvalues decay faster. This makes eigenfunctions with large eigenvalues contribute less significantly during mes-
Figure 4: Examples of versatile mes- sage passing outcomes using different frequency filter settings.
sage passing. Therefore, the combination of the Gaussian frequency filter and heat operator could help the network focus on some higher frequency components (Aubry al. 2011).
As demonstrated, HMR is able to represent variety of surface functions through harmonic message passing. The output of this module has the same size as the input, Fmp E NR, where each channel respectively contains the propagated version of its input function. These features represent the neighboring geometric and chemical environment across multiple spatial scales, and can be used for surface property- related tasks, e.g., protein binding site prediction (see [Sec. 4.2) In addition, molecule-Ieve representations could be obtained through global pooling (see Sec. 5.1).
4.2 LEARNING SURFACE CORRESPONDENCE FOR RIGID PROTEIN DOCKING
In this section, we demonstrate how to learn surface correspondence for molecular matching. Specif- ically, we introduce a surface-based rigid protein docking workflow Fig. powered by HMR. Rigid protein docking a significant problem in structural biology, whose goal is to predict the pose of the protein complex based on the structure of the ligand and receptor proteins. It has been shown that protein-protein interfaces exhibit similar geometric and chemical patterns (Gainza et al., 2020 In other words, two proteins may interact if part of their surfaces display similar shapes and chemical functions (i.c., functional correspondence). This is similar to solving a puzzle problem, where both shape and pattern of the missing piece have match in order to complete the puzzle.
Following this idea, we propose to realize rigid protein docking in two consecutive steps: (1) given two protein surfaces, prediet the region where binding might occur (i.e., binding site prediction, lo- cating the missing piece); (2) establish functional correspondence between the ligand/receptor bind- ing surfaces, and convert it to real-space parte-too-verte correspondence (shape/pattern matching). Rigid docking could then be achieved by aligning the corresponding binding site surface vertices.

---

Published as a conference paper at ICLR 2023
Binding Site Prediction Module
Rigid Docking Module
rocestor interface
Surface Preparation
HMR
Cross Attention
Binding Site Classifier
Ma
Tell
Kabsch
rotation translation
Mi.
Iiaan
Surface
HMR
Cross Attantion
Binding Site
Functional Correspondence
Interface
Preparation
Leyer
Clasifior
Corsespendence
Predicted Protein Complex Structure
Ipand
Figure 5: The rigid protein docking pipeline. Given the surface of the ligand and receptor proteins, apply multiple HMR and cross attention layers to allow communication within and between the surfaces. The learned surface representations are then used to predict the binding interfaces and establish functional correspondence using functional maps. Rigid docking is achieved by converting the functional correspondence to rigid transformations which aligns the predicted interfaces.
Binding Site Prediction Given the ligand and receptor protein surface meshes, we first predict the regions where they interact, which is per-vertex binary classification problem. We iteratively apply HMR and cross-attention layers (Fig. 5) to encode the surfaces with intra- and inter-surface communications, Next, use the learned features on each vertex classify whether it belongs to the binding interface. Detailed descriptions of this module are available in Appendix D
The output of the binding site prediction module includes: a) the Receptor (Ligand) surface mesh of the binding interface, which is submanifold of the entire protein surface, represented as MA (Mz) with NR (NL) vertices, and neural network-leande features or so-called surface functions FMR and F.M distributed on the receptor and ligand binding interfaces, respectively.
Rigid Docking with Functional Maps We know that protein interfaces MR and Mr exhibit similar geometry and also adopt set of corresponding functions F.R E VR and F.M e Ne Fig. 5). Intuitively, if we could somehow align these of corresponding functions, then we have found way to align the protein binding interfaces. To that end, we employ functional maps to "align" F.M. and F.M. in spectral domain.
Specifically, given the truncated LB eigenfunctions HM E Na and G- E N2 of the re- ceptor and ligand interface manifolds, we respectively compute the spectral representation of learned functions as A (G and B Ma M +M1, where + denotes the Koore-Penrose pseudo-inverse, and A, E Rkxn Functional correspondence (C) can be obtained by minimizing:
= arg min ||CA -BI(F,
CER
where ||. ||r denotes the Frobenius norm. Once the functional mapping (i.e., the C. matrix) is recovered through numerical optimization, vertex-to-vertex correspondence between the receptor and ligand surfaces can be established by mapping indicator functions of vertices on MR those on Mr. In practice, we adopted a slightly more complicated functional maps approach, which is illustrated (mpppendix Finally, we perform rigid docking by aligning the proteins according to the vertex-to- to-exxt-t--etex surface correspondence using the Kabsch algorithm (Kabsch) 1976).
EXPERIMENTS
Table QM9 Mean Absolute Error
5.1 QM9 MOLECULAR PROPERTY REGRESSION
We employ HMR to perform property regression tasks on the QM9 dataset (Ramakrishnan al.) 2014). We use the same data split and atomic features as Satorras et al. (2021). We compare our results with both invariant and equivarient networks shown in Table I including SchNet (Schütt et al. 2018). NMP (Gilmer et al. (017). Cormorant (Anderson et al. 2019), SE(3)-Transformer (Fuchs et al.) 2020).
Task Unit
a
Ag
9номо meV
C.
bohr
meV
D
cal/(mol K)
SchNet
.235
63
41
.033
.033
NMP
.092
69
43
.030
.040
Cormorant
.085
61
34
.038
.026
SE()--Tr.
.142
35
.051
.054
SEGNN
.060
42
24
.023
.031
HMR
.102
59
.037
.040

---

Published as conference paper at ICLR 2023
and SEGNN (Brandstetter et al. 2021). See Appendix F for the complete table and experimental setup. Interestingly, despite HMR completely discards the bonding information and only performs massage passing over the molecular surface, it still shows comparable performance in predicting these molecular properties.
5.2
HIGAND-BINDING POCKET CLASSIFICATION
HMR
MaSiF-ligand
Next, showcase the representation power of HMR in predicting the ligand-binding preference protein pock- ets. The preference of protein binding specific small molecule (ligand) is associated with the geometric and chemical environment at the binding region (i.e., bind- ing pocket). As defined in Gainza et al. (202). we clas- sify a protein pocket as one seven ligand- binding types (ADP, COA, FAD, HEM, NAD, NAP, and SAM) using its surface information. The dataset contains 1,634 training. 202 validation, and 418 test cases. We compare the per- formance of HMR to MaSIF-ligand (Gainza et al. 2020), which is geodesic convolutional neural network-based model. set of features similar to MaSIF-ligand are used in our model: hydrophobicity score and partial charge chemical features; mean/Gaussian curvature and the Heat Kernel Signatures geometric features.
1.0
Balanced Accuracy
0.81
0.8
0.71
0.66
0.85
0.6
0.4
0.2
0.0
G+C
Geom
Chem
Figure 6: Balanced accuracy of ligand- binding protein pocket classification. Models with different input features are compared: "G+C" includes both ge- ometric and chemical features, while "Geom"/" """Cemm only contains geomet- ric/chemical features.
As shown Fig. 6 the balanced accuracy of our model consistently better than that reported in Gainza al.) (2020), suggesting that the proposed HMR framework can more effectively encode protein surface information through harmonic message passing over the molecular surface Riemannian manifold. In addition, we draw similar conclusion that both geometric and chemical information of the binding pockets are important in predicting the type of its binding molecules. See Appendix for implementation details and more result analysis.
5.3 RIGID PROTEIN DOCKING
Finally, evaluate HMR on more challenging task: rigid protein docking.
Experimental setup HMR is trained on a modified version Database Interacting Protein Structures (DIPS) (Twnnseend et al. 2019) and evaluated on the gold-standard Docking Bench- mark 5.5 (DB5.5) (Guest et al. 2021). We compare our model with the state-of-the-an the-art GNN-based deep learning model EQUIDOCK (Ganea et al., 202I) and two traditional docking methods, AT- 2015) and HDOCK (Yan et al. 2020). To evaluate docking performance,
TRACT (de Vries et al. compute the Complex and Interface noot-mea-sqqaal deviation (RMSD) following Ganea al.) (202I). and calculate DockQ following Basu & Wallner (2016). We also report the success rate in- dicating whether the result achieves "Acceptable" higher according Lensink & Wodak (203))
Results Model performance are summarized in Table 2 and Fig. 7a. HMR (Top 1) outperforms GNN-based EQUIDOCK model under all metrics. Notably, HMR achieves much higher success rate, indicating more test cases have results close to the ground truth complex structure. We note that traditional methods still exceed in terms of docking performance but at greater computational
Table 2: Rigid Prediction Results on Docking Benchmark 5.5
Complex RMSD |
Interface RMSD
DockQ f
Success rate 1
Inference
Model
Median
Mean
Median
Mean
Median
Mean
(2 Acceptable)
time (sec)
HMR (Top
12.01
12.09
10,90
11.43
0.06
0.22
0.28
5.8
HMR (Top 3)
8.44
9.68
7.50
8,70
0.10
0.26
0.33
6,6
EQUIDOcK*
15.90
17.18
14.45
14.85
0.02
0.04
0.00
3.9
ATTRACT
8.99
11.11
10.31
11.67
0.07
0.42
0.44
882.7
HDOCK
0.39
5.97
0.32
5.68
0.97
0.71
0.73
884.9
*EQUIDoCk model is trained on DIPS, provided by the authors

---

Published as a conference paper at ICLR 2023
Conolty RM53
Dookk
Pesroon'
20
0.6
Accepadee
High
GUlftinaniigg
correspondence
PDB ID: 1RB5
(MSD:
Deckq:
RMSD:
CAPRI gi: Hig!
Receptor
Predicted Interfaces
Figure 7: Rigid docking performance analysis. a. Distribution of Complex RMSD and DockQ scores for poses predicted using HMR (ours, Top 1), EQUIDoCK, HDOCK, and ATTRACT. b. Cor- relations between learned functions on the predicted ligand-receple interfaces. For each test case, calculate the correlation the Pearson's between the learned function values on ligand inter- face vertices and their nearest receptor vertices (corresponding vertices), averaged over 128 hidden channels. A showcase of functional correspondence between interfaces for well-docked case.
cost (Table H.T). Further experiments show that the harmonic message passing mechanism learns to
propagate information different scales (Fig. H.I and is critical the effectiveness of the model (Table H.2). Our feature ablation studies show that chemical information is particularly important in the rigid protein docking task Table H.3). Compared EQUIDocK, the novel framework of HMR and optimized training dataset collaboratively contribute to the higher performance (Talle H.4
We further examine the learned features the predicted binding interfaces. As shown in Fig. 7pc, surface functions on the two binding sites are highly correlated, confirming the good alignment achieved using functional maps. Cases with higher docking quality show stronger interface feature correlations, suggesting our HMR is able to learn complex surface interaction patterns, which also supports the claims in Gainza et al. (2020) about the significant role of protein surface.
HMR predicts multiple binding sites for some proteins (either due to certain protein symmetries or model uncertainty). Therefore, we also assess the model performance by including candidate poses from top 3 binding site pairs, ranked by the mean probability predicted by the binding site classifer. The best scores from top 3 poses show that the performance of HMR is competitive to ATTRACT. See Appendix for more detailed analysis the rigid protein docking experiment.
CONCLUSIONS AND OUTLOOK
We presented HMR, powerful surface manifold-based molecular representation learning frame- work. By integrating geometric and chemical properties as functions distributed on the molecular surface manifold, and applying harmonic message passing in spectral domain, we achieve multi- resolution molecular representations. HMR shows promising performance in molecular property prediction, protein pocket classification, and molecular matching tasks. Our work highlights im- portant aspect of molecular "structure-aaciivity" relationship, that is "shape-activity" relationship. We think this is particularly significant for large biomolecules, where the surface shape and chemical patterns determine some fundamental biological activities, such protein-protein interactions.
The HMR framework could serve as complementary method GNN/NMP- models in solv- ing challenges for complex biological systems, which exhibits its unique advantages and shortcom- ings. One of the foreseeable challenges in further developing the HMR framework is the efficient computation of the "Shape-DNA" (i.e., solving second-order partial differential equation) in der to incorporate protein dynamics. Since either a change of the protein backbone or some side chains near the surface may alter the entire "Shape-DNA", which needs to be recomputed upon pro- tein conformational change. This is particularly important for large biomolecules where the surface shape undergoes significant changes (c.g., the papplementtrtt--Deeenniie body upon binding with an antigen, or the allosteric site of some enzymes). To that end, we call for more research attention surface manifold-bar based molecular representation learning.
Regions of an anti-

---

Published as a conference paper at ICLR 2023
ACKNOWLEDGEMENTS
The authors thank Dr. Hang Li and Dr. Quanquan Gu for their insightful comments. Hao Zhou is supported by Vanke Special Fund for Public Health and Health Discipline Development, Tsinghua University N10.20221880055 Guoqiang Research Institute General Project, Tsinghua University (No. 2021GQG1012).
REPRODUCIBILITY STATEMENT
The code and data are available https//ghhho.ccm/ooooooollllglgr QM9 raw
dataset is provided at https: / http://sppinnngenurr..e... figshare.................. hompaare.cm//nee hwtpm//noownodder...
files/3195389 The dataset for the ligand-binding pocket classification is provided at and the split used by MaSIF is at https:
DIPS dataset can be downloaded from the following website https: /tttp:://githbe com/ EQUIDOCK model and checkpoints can be ATTRACT can be downloaded from http:///taaaaaaaaaaaaaaaaaaaaaaa public HDOCK is implemented using its local version HDOCKlite, which be downloaded from
downloaded from bttps://gitnn:
http:a/ww.....oo...............r be downloaded http:///ttaaeeggt.t.oo..o///o//a
DockQ can
REFERENCES
Brandon Anderson, Truong Son Hy, and Risi Kondor. Cormorant: Covariant molecular neural networks. Advances in neural information processing systems, 32, 2019.
Souhaib Attaiki, Gautam Pai, and Maks Ovsjanikov. DPFM: Deep partial functional maps. In 2021 International Conference on 3D Vision (3DV), pp. 175-185, IEEE, 2021.
Kenneth Atz, Francesca Grisoni, and Gisbert Schneider. Geometric deep learning on molecular representations. Nature Machine Intelligence, 3(12): 1111000210222 2021.
Mathieu Aubry, Ulrich Schlickewei, and Daniel Cremers. The wave kemel signature: A quantum mechanical approach to shape analysis. In 2011 IEEE international conference on computer vision workshops (ICCV workshops), pp. 1626-1633 1633. IEEE, 2011.
Sankar Basu and Björn Wallner. DockQ: A Quality Measure for Protein-Protein Docking Models.
PLOS ONE, (1):6111,, 2016.
Federico Battiston, Giulia Cencetti, Iacopo lacopini, Vito Latora, Maxime Lucas, Alice Patania, Jean-Gabriel Young, and Giovanni Petri Networks beyond pairwise interactions: structure and dynamics. Physics Reports, 874:1-92, 2020.
Silvia Biasotti, Andrea Cerri, Alex Bronstein, and Michael Bronstein. Recent trends, applications, and perspectives in 3d shape similarity assessment, In Computer graphics forum, volume 35, pp- 87-119. Wiley Online Library, 2016.
Marian Boguna, Ivan Bonamassa, Manlio De Domenico, Shlomo Havlin, Dmitri Krioukov, and 2021.
M Serrano. Network geometry. Nature Reviews Physics,
Johannes Brandstetter, Rob Hesselink, Elise van der Pol, Erik Bekkers, and Max Welling, Ge- ometric and physical quantities improve e (3) equivariant message passing. arXiv preprint arXiv:2110.02900 2021.
Michael M Bronstein and lasonas Kokkinos. Scale-invariant heat kernel signatures for non-rigid shape recognition. In 2010 IEEE computer society conference on computer vision and pattern recognition, pp. 1704-1711. IEEE, 2010.
Michael M Bronstein, Joan Bruna, Yann LeCun, Arthur Szlam, and Pierr Vandergheynst. Geomet- ric deep learning: going beyond Euclidean data. IEEE Signal Processing Magazine, 34(4):18 42, 2017.
10

---

Published as a conference paper at ICLR 2023
Michael Bronstein, Joan Bruna, Taco Cohen, and Petar Veličković. Geometric deep learning: Grids, groups, graphs, geodesies, and gauges. arXiv preprint arXfv:2104.1347,, 2021.
Taco Cohen, Mario Geiger, Jonas Köhler, and Max Welling. Spherical CNNs. arXiv preprint arXivv 2018.
Ronald R Coifman and Stéphane Lafon. Diffusion maps. Applied and computational harmonic analysis, 21(1):5-30, 2006.
Ronald R Coifman, Stephane Lafon, Ann B Lee, Mauro Maggioni, Boaz Nadler, Frederick Warner, and Steven W Zucker. Geometric diffusions as a tool for harmonic analysis and structure defini- tion of data: Diffusion maps. Proceedings of the national academy of sciences, 102(21):7426- 1431,2005,
Sjoerd J. de Vries, Christina E. M. Schindler, Isaure Chauvot de Beauchêne, and Martin Zacharias. Web Interface for Easy Flexible Protein-Protein Docking with ATTRACT. Biophysical Journal, 108(3):462-4653 February 2015.
Todd J Dolinsky, Jens E Nielsen, J Andrew McCammon, and Nathan A Baker. Pdb2pqr: au- tomated pipeline for the setup of poisson-boltzmaan electrostatics calculations. Nucleic acids research,
2004.
Nicolas Donati, Abhishek Sharma, and Maks Ovsjanikov. Deep geometric functional maps: Ro- bust feature learning for shape correspondence. Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 8592-8601, 2020.
Dina Duhovny, Ruth Nussinov, and Haim J Wolfson. Efficient unbound docking of rigid molecules. In International workshop on algorithms in bioinformatics, pp. 185-200. Springer, 2002.
Gregory Ewing and Joachim Hermisson. Msms: coalescent simulation program including recom- bination, demographic structure and selection at single locus. Bioinformatics, 26(16):2064- 2065, 2010.
Fabian Fuchs, Daniel Worrall, Volker Fischer, and Max Welling. SE(3)-Transfong 3D roto- translation equivariant attention networks. Advances in Neural Information Processing Systems, 33: 33:1970-1981, 2020.
Pablo Gainza, Freyr Sverrisson, Frederico Monti, Emanuele Rodola, D Boscaini, MM Bronstein, and BE Correia. Deciphering interaction fingerprints from protein molecular surfaces using geo- metric deep learning. Nature Methods, 17(2):14-1122 2020.
Octavian-Eugen Ganea, Xinyuan Huang, Charlotte Bunne, Yatao Bian, Regina Barzilay, Tommi Jaakkola, and Andreas Krause. Independent -Equivaqui models for end-to-end rigid protein docking. arXiv preprint arXiv:2/11.07780
2021.
Justin Gilmer, Samuel Schoenholz, Patrick F Riley, Oriol Vinyals, and George E Dahl. Neural message passing for quantum chemistry. In International conference on machine learning, pp. 1263-1272. PMLR, 2017.
Vladimir Gligorijević, P Douglas Renfrew, Tomasz Kosciolek, Julia Kochler Leman, Daniel Beren- berg, Tommi Vatanen, Chris Chandler, Bryn C Taylor, Ian M Fisk, Hera Vlamakis, et al. Structure- based protein function prediction using graph convolutional networks. Nature communications, 12(1):1-14, 2021.
Johnathan D. Guest, Thom Vreven, Jing Zhou, lain Moal, Jeliazko R. Jeliazkov, Jeffrey J. Gray, Zhiping Weng. and Brian G. Pierce. An expanded benchmark for antibody-antigen docking and affinity prediction reveals insights into antibody recognition determinants. Structure, 29(6):606~ 621.e5, June 2021.
Clemens Isert, Kenneth Atz, and Gisbert Schneider. Structure-based drug design with geometric deep learning. arXiv preprint arXiv:22 10. /250, 2022.
Alec Jacobson and Daniele Panozzo. Libigl: Prototyping geometry processing research in c++. In SIGGRAPH Asia 2017 courses, pp. 1-172. 2017."

---

Published as a conference paper at ICLR 2023
John Jumper, Richard Evans, Alexander Pritzel, Tim Green, Michael Figurnov, Olaf Ronneberger, Kathryn Tunyasuvunakool, Russ Bates, Augustin Zidek, Anna Potapenko, al. Highly accurate protein structure prediction with alphafold. Nature, 596(7873):58355888 2021.
Wolfgang Kabsch. A solution for the best rotation to relate two sets of vectors, Acta Crystallo- graphica Section A: Crystal Physics, Diffraction, Theoretical and General Crystallography, 32 (5):922-923, 1976.
Mark Kac. Can one hear the shape of a drum? The american mathematical monthly, 73(4P2):1-23, 1966.
Thomas N Kipf and Max Welling. Semi-supervised classification with graph convolutional net- works. arXiv preprint arXiv: (00900297, 2016.
Johannes Klicpera, Janek GroB, and Stephan Günnemann. Directional message passing for molec- ular graphs. arXiv preprint arXiv:2003.03122 2020.
John M Lee. Smooth manifolds. In Introduction to smooth manifolds, pp- 1-31. Springer, 2013.
Mare Lensink and Shoshana J. Wodak. Docking, scoring, and affinity prediction in CAPRI. Proteins: Structure, Function, and Bioinformatics, 81(12):0082-2095, 2013.
Xianzhi Li, Ruihui Li, Guangyong Chen, Chi-Wine Chi-Wing Fu, Daniel Cohen-On and Pheng-Amn Heng. A invariant framework for deep point cloud analysis. IEEE Transactions on Visualization and Computer Graphics, 2021.
Ye Li, Xianren Zhang, and Dapeng Cao. The role of shape complementarity in the protein-protee protein interactions. Scientific reports, 311:1177,011. 2013.
Or Litany, Tal Remez, Emanuele Rodola, Alex Bronstein, and Michael Bronstein. Deep functional maps: Structured prediction for dense shape correspondence. In Proceedings of the IEEE inter- national conference computer vision, pp. 5659-5667, 2017a.
Or Litany, Emanuele Rodolà, Alexander M Bronstein, and Michael M Bronstein. Fully spectral partial shape matching. In Computer Graphics Forum, volume 36, pp. 247-258. Wiley Online Library, 2017b.
Roee Litman and Alexander M Bronstein. Learning spectral descriptors for deformable shape cor- respondence. IEEE transactions on pattern analysis and machine intelligence, 36(1):171-180, 2013.
Qinqing Liu, Peng-Shuai Wang, Chunjiang Zhu, Blake Blumenfeld Gaines, Tan Zhu, Jinbo Bi, and Minghu Song. OctSurf: Efficient hierarchical voxel-hased molecular surface representation for protein-ligand affinity prediction. Journal of Molecular Graphics and Modelling, 105:107865, 2021.
Federico Monti, Davide Boscaini, Jonathan Masci, Emanuele Rodola, Jan Svoboda, and Michael M Bronstein. Geometric deep learning graphs and manifolds using mixture model CNNs. In Proceedings of the IEEE conference on computer vision and pattern recognition, pp. 5115-5124, 2017.
Stelios Mylonas, Apostolos Axenopoulos, and Petros Daras. DeepSurf: surface-based deep learning approach for the prediction of ligand binding sites proteins. Bioinformatics, 37(12): 1681-1690, 2021.
Raúl Méndez, Raphaël Leplae, Leonardo De Maria, and Shoshana J. Wodak. Assessment of blind protein-protein
predictions of protein protein interactions: Current status of docking methods. Proteins: Struc- ture, Function, and Bioinformatics, (2(1):51-61 2003.
Raúl Méndez, Raphaël Leplae, Mare F. Lensink, and Shoshana J. Wodak. Assessment CAPRI predictions in rounds 3-5 shows progress in docking procedures. Proteins: Structure, Function, and Bioinformatics, 60(2):150-169, 2005.

---

Published as a conference paper at ICLR 2023
Maks Ovsjanikov, Mirela Ben-Chen, Justin Solomon, Adrian Butscher, and Leonidas Guibas. Func- tional maps: flexible representation of maps between shapes. ACM Transactions Graphics (ToG), 31(4):1-11, 2012.
Peter Petersen, Riemannian geometry, volume 171. Springer, 2006.
Raghunathan Ramakrishnan, Pavlo Dral, Matthias Rupp, and Anatole Von Lilienfeld. Quantum chemistry structures and properties of 134 kilo molecules. Scientific data, 1(1):1-7,2014
Jing Ren, Adrien Poulenard, Peter Wonka, and Maks Ovsjanikov. Continuous and orientation- preserving correspondences via functional maps. ACM Transactions on Graphics (ToG), 37(6): 1-16, 2018.
Martin Reuter, Franz-Erich Wolter, and Niklas Peinecke, Laplace-beltrami spectra 'shape-dna's of surfaces and solids. Computer-Aided Design, 38(4):342-34) -36, 2006.
Martin Reuter, Silvia Biasotti, Danicla Giorgi, Giuseppe Patanè, and Michela Spagnuolo. Discrete laplace-beltrami operators for shape analysis and segmentation. Computers & Graphics, 33(3): 381-390, 2009.
Frederic M Richards. Areas, volumes, packing, and protein structure. Annual review of biophysics and bioengineering, 6(1):151-176, 1977.
Steven Rosenberg. The Laplacian Riemannian manifold: an introduction to analysis mani- folds. Number 31. Cambridge University Press, 1997.
Victor Garcia Satorras, Emiel Hoogeboom, and Max Welling. E(n) equivariant graph neural net- works. In International conference machine learning, pp. 9323-9332. PMLR, 2021.
Kristof T Schütt, Huziel E Sauceda, P-J Kindermans, Alexandre Tkatchenko, and K-R Müller. Schnet-a deep learning architecture for molecules and materials. The Journal of Chemical Physics, 148(24)244122, 2018.
Kim A Sharp. Electrostatic interactions in macromolecules. Current Opinion in Structural Biology, 4(2):234-22) 239, 1994.
Yangmei Shen, Wenrui Dai, Chenglin Li, Junni Zou, and Hongkai Xiong. Multi-scale graph convo- lutional network with spectral graph wavelet frame. IEEE Transactions Signal and Information Processing over Networks, 210, 2021.
Alexandra Shulman-Peleg Peleg, Ruth Nussinov, and Haim Wolfson. Recognition of functional sites in protein structures. Journal of molecular biology, 339(3)x077-633, 2004.
Vignesh Ram Somnath, Charlotte Bunne, and Andreas Krause. Multi-scale representation learning on proteins. Advances in Neural Information Processing Systems, 34:25244-25255, 2021.
Jonathan Stokes, Kevin Yang, Kyle Swanson, Wengong Jin, Andres Cubillos-Ruiz, Nina M Donghia, Craig R MacNair, Shawn French, Lindsey A Carfrae. Zohar Bloom-Ackermann, et al. A deep learning approach to antibiotic discovery. Cell, 180(4):688-700 2020.
Jian Sun, Maks Ovsjanikov, and Leonidas Guibas. A concise and provably informative multi-scale signature based on heat diffusion. In Computer graphics forum, volume 28, pp. 1383-1392. Wiley Online Library, 2009.
Freyr Sverrisson, Jean Feydy, Bruno E Corrcia, and Michael M Bronstein. Fast end-to-end learning protein surfaces. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 15272-15281, 2021.
Freyr Sverrisson, Jean Feydy, Joshua Southern, Michael M Bronstein, and Bruno Correia. Physics- informed deep neural network for rigid-body protein docking. In ICLR2022 Machine Learning for Drug Discovery, 2022.
Nathaniel Thomas, Tess Smidt, Steven Kearnes, Lusann Yang, Li Li, Kai Kohlhoff, and Patrick Riley. Tensor field networks: Rotation-and translation-equivrriant neural networks for 3d point 2018.
clouds. arXiv preprint
13

---

Published as a conference paper at ICLR 2023
Raphael J. L. Townshend, Rishi Bedi, Patricia A. Suriana, and Ron O. Dror. End-to-End Learn- ing 3D Protein Structure for Interface Prediction. Technical report, arXiv, December 2019. arXiv:1807.01299 [cs, q-bio, stat] type: article.
Raphael JL Townshend, Stephan Eismann, Andrew M Watkins, Ramya Rangan, Maria Karelina, Rhiju Das, and Ron O Dror. Geometric deep learning of structure. Science, 373(6558):
1047-1051, 2021.
Oliver Unke, Mihail Bogojeski, Michael Gastegger, Mario Geiger, Tess Smidt, and Klaus-Robert Müller. SE())-quivarriat prediction molecular wavefunctions and electronic densities. Ad- vances in Neural Information Processing Systems, 34:14434-14447, 2021.
Petar Veličković, Guillem Cucurull, Arantxa Casanova, Adriana Romero, Pietro Liò, and Yoshua Bengio. Graph Attention Networks, February 2018. htt.ttpp://axxli.orr/bbbbb
10903 arXiv:17 109003 [cs, stat].
J Michael Word, Simon C Lovell, Jane S Richardson, and David Richardson. Asparagine and glutamine: using hydrogen atom contacts in the choice of side-chain amide orientation. Journal of molecular biology, 283(4):1735-1777, 1999.
Saining Xie, Jiatao Gu, Demi Guo, Charles R Qi. Leonidas Guibas, and Or Litany. Pointcontrast: Unsupervised pre-training for 3d point cloud understanding. In European conference computer vision, pp. 574-591. Springer, 2020.
Yumeng Yan, Huanyu Tao, Jiahua He, and Sheng- You Huang. The HDOCK server for inte- protein-protein grated protein protein docking. Nature Protocols, 18(5:1182-11822 May 2020. doi: 10.1038/ s4 1596-020-03112 Number: 5 Publisher: Nature Publishing Group.
Zhou. Pymesh- processing library for python. Software available for downloud at hhtps//giiaahccoom///eeee//eeehe 2019.

---

Published as a conference paper at ICLR 2023
Appendix
A THE MOLECULAR SHAPE-DNA
Riemannian Manifold A manifold is space that is locally flat but not necessarily globally flat. Formally speaking, a d-dimensional manifold M is topological space where each point p E M has a neighborhood that is homeomorphic to a d-dimensional Euclidean space (Lee) 2013), which equivalent to the tangent space at p and is denoted by T.MM.
We can further assign a positive definite inner product g TpM x T,M + R on every tangent space, and this inner product is called a Riemannian metric. A manifold equipped with Riemannian metric called Riemannian manifold. Intuitively, the Riemannian metric provides a measurement of the velocity when a particle moves on the manifold, and many other quantities defined. For example, for any tangent vector Xp Tp.M, the quantity [Xpl v (XXp. Xp) can be interpreted as the traveling speed of particle when passing through p. Hence, the traveling distance along curve (i.e., the length of curve) can be defined the integral of |Xpl along the curve. On Riemannian manifold, quantities that can be expressed in terms of the Riemannian metric are called intrinsic (e.g., geodesic distance).
therefore be
When a manifold is realized in the Euclidean space, a natural Riemannian metric can be induced from the ambient Euclidean space. We always refer to this induced metric when talking about Riemannian manifold in the following appendices.
The Lalpace- Salpace-Beltrami Beltrami Operator We denote a real-valued scalar function on the manifold M by f. Given two functions f1, f2 on the manifold, we can define the inner product (f1, f2) = [M (=)))))))) where the area element du induced by the Riemannian metric. We denote by (M) = (ʃ: M - IR (S: J' < ]o) the space of square-integrable functions on M. On Riemannian manifold M, we can generalize the usual Euclidean gradient V f and the positive semi- definite Laplace operator Af = -div f) to the intrinsic gradient vf and the Laplace-Beline Beltrami (LB) operator Amf, respectively (Petersen 2006).
The LB operator admits a discrete set of eigenfunctions that solves
(m) = (x))(())
EM
(A.1)
with homogeneous Neumann boundary conditions if . has boundary. Here, 0 A1 S X2 S ... are eigenvalues and Φ1, d22,.. ... are the corresponding eigenfunctions. The LB eigenfunctions form an orthonormal basis of T,² (MM) i.e., (6., 0;)M d.j.
Molecular Shape DNA Intuitively, the LB eigenfunctions are the smoothest functions on the surface, where smoothness measured by the Dirichlet energy ((Afff) Eigenfunctions with smaller eigenvalues are more smooth (i.e., with lower spatial frequency). However, these eigenfunc- tions are not directly comparable across different manifolds. Fig. A.I. shows the first 300 eigen-
b
ligand
receptor
0.5
$10
Eigenvalues
0.4
0.3
Aιo 0.017
Aso 0.008
0.2
2200
AVR = 0.344
0
50
100
150
200
250
300
Aana 0.161
Number of eigenfunctions
Figure A.1: Weyl's asymptotic law, where the protein with larger surface area exhibits slower cigenvalue growth rate. The 10th and th eigenfunction of the ligand and receptor surface manifold, respectively. Eigenfunctions with similar eigenvalues exhibit similar spatial frequencies regardless of the surface area, since they have similar smoothness in the sense of Dirichlet energy.
15

---

Published as a conference paper at ICLR 2023
values of a pair of ligand and receptor protein surfaces (PDB ID: 3V62). The protein with smaller surface area exhibits a larger slope, known Weyl's asymptotic law. Fig. A.Ip presents few eigen- functions with correspondending eigenvalues for the ligand and receptor molecules. Since different eigenfunctions have their unique spatial resolutions, in experiments tions with eigenvalues smaller than predefined value (determined empirically) to guarantee that molecules of different sizes have the same highest spatial resolution in their eigenfunctions.
compute all eigenfunc-
B IMPLEMENTING MANIFOLD HARMONIC ANALYSIS
Laplace-Beltrami Eigendecomposition Given Riemannian manifold M and its Laplace- Beltrami operator A,, the Laplacian eigenvalue problem we consider states as
Amf Af,
with homogeneous Neumann boundary condition.
To realize discrete calculations, we approximate the manifold with triangle mesh consisting of N vertices } C M, and the corresponding faces. We then solve the discretized eigenvalue problem using linear finite element method (FEM) (Reuter et al., 2009). The discretized equation we obtain is the following generalized eigenvalue problem
A f = ABf, f (f(x)))
(B.1)
where Acot is the stiffness matrix with cotagent weights,
2 (cot + cot Biy)
if (i,j) is an edge
Aco (i, = EkEN(a) (1, ifi=j
otherwise,
and B an N N sparse mass matrix which is associated with the weight each surface vertex,
T2 ((4 (12))
if (i,j) is an edge
B(i,))
E
/41 if
0
otherwise
Here ai and By are the two angles opposite to the edge (i,j), and M(i) denotes the vertices that are adjacent to vertex The set T(i) contains all the triangles that have as its vertex, and [t&] the area of the triangle t,. We also use & and to to denote the triangles that share the edge (i. j).
Such generalized symmetric eigenvalue problem can be solved with commonly used numerical sim- ulaiton packages (e.g., scipy). We can find non-negative eigenvalues A and eigenvectors Z such that
Z A
ZZ 1.
I is the identity matrix, A == diag( (oo AI....) is the diagonal matrix of the eigenvalues, and the matrix Z E IR NxN has the eigenvectors of Eq. B. as its column vectors. Note that, unlike other conventional orthonormal basis (i.c., Q Q here Z forms orthonormal basis w.r.t. the mass matrix B, that is, Z BZ I.
Practically, do not need to store the entire N x N eigenvector matrix Z. Just like in Fourier series expansion, a truncated Fourier basis with finite terms can be used to approximate the original signal. The number of basis that we keep determines the resolution of the synthesized signal, where typi- cally high-frequency components are truncated. In our case, we empirically determine the number of eigenvectors to keep for different molecular systems, which also task-dependent.

---

Published as a conference paper at ICLR 2023
Resolution Tuning with Harmonic Analysis We now explain how to realize surface resolution tuning under our representation framework (e.g., how to make Fig. 1.
Given a molecular surface triangle mesh, we first compute its Laplace- Belttramm cigendecomposition as described above. Under the discrete setting, we refer to the cigenfunctions eigenvectors. We
obtain a set of truncated Laplace- Laplace-Beltrami eigenvectors Z E TR. x (k is the number of eigenvectors we keep), the associated eigenvalues {A, 1 and the sparse mass matrix B E RNN..
Let f TR- be the initial function of interest (stored as an N-dimensional array), which is distributed on the underlying molecular surface (e.g.. electrostatic potential, scalar value at each surface ver- tex). The spectral representation of function can be calculated as
fapec = Z Bf,
papec R²
which is similar to a discrete Fourier transform. To project the function back to real space (inverse Fourier transform), we simply do
= f' PERN.
Resolution tuning achieved by controlling the number of basis we use (i.e., tuning k) in spectral space. As shown in Fig. 1, the electrostatic potential function resolution can be tuned by manipulat- ing the number of Laplace-Beltrant eigenfunctions.
However, what about the resolution of the shape itself? In Fig. we see that the surface smoothness can also be tuned. It is important to realize that manifold is an abstract concept, which does not necessarily have a particular realization in the Enclidean space. The surfaces that we visualize in Fig. are realizations of the underlying manifold in the Euclidean space, where each surface vertex associated with some Cartesian coordinates. These coordinates are extrinsic properties of the manifold, thus can be treated in the same way as other surface functions (e.g.. the electrostatic potential). Therefore, in order to reconstruct the molecular surface with lower spatial resolution, we can simply calculate the smoothed Cartesian coordinates:
x = Zx ZZ Bx,
and do the same for and z to obtain the smoothed surface coordinates (x, y').). In short, we realize resolution tuning with a series of (sparse) matrix multiplications bringing the surface functions back- and-forth between the real space and the generalized Fourier space.
C SURFACE PREPARATION
The raw input to the HMR framework is simply the 3D atomic structures (e.g., xyz files for small molecules, or PDB files for proteins). In other words, we only need to know where these atoms and their atomic species. For proteins with only heavy atoms (since hydrogen atoms are almost invisible under X-ray diffraction detectors), we use reduce (Word et al. 1999) (or alternatively PDB2POR (Dolinsky al. 2004)) software to add hydrogen atoms.
Next, we employ MSMS (Ewing & Hermisson 2010) to calculate the solvent-excluded surface of the molecule (with probe radius A. sampling density 3.0 for small molecules and 1.0 for proteins) as a triangle mesh. We use PyMesh (Zhou) 2019) to further refine the surface mesh in order to reduce the number of vertices and fix poorly meshed areas. Degenerate vertices or disconnected surfaces would lead to numerical issues for solving the generalized eigenvalue problem in the next step, thus should be fixed beforehand. We then compute the truncated Caplace-Beltrant eigenvectors, eigenvalues, and the mass matrix as described in Appendix B
Initial geometric features can be directly calculated given the surface triangle mesh, where we use the libigl (Jacobson & Panozzo 2017) package to calculate the mean and Gaussian curvatures, and compute the Heat Kernel Signatures as described in Sun et al. (2009). These geometric fea- tures capture shape-related properties of the molecular surface, and stored scalar-type array Freom E R N Nxp, where is the number of initial geometric features (a user defined variable).
Chemical features projected from atoms to their neighboring surface vertices. We first obtain an initial descriptor vector for each atom (e.g., atomic number, charge, etc.). Then, for each surface
17

---

Published as a conference paper at ICLR 2023
vertex X;, we compute its v nearest neighbor atoms centered at fai. a with features {ui.... }. We apply multilayer perceptron (MLP) the vector u | a ]I] for each neighboring atom, then compute the average over the neighbors to obtain the chemical feature vector Frhem E RNx,, where q is user defined variable indicating the dimension of initial chemical features. In short, the initial chemical features of each surface vertex are determined by its neighboring atomic species and their distance, which learned by MLP.
We use another MLP to combine the per-vertex initial features
FiuP MLP(concat (Facmm Fchem)), FiuP
Nx
These n. features reflect the local geometric and chemical environment of each surface vertex, which will be used as the input the harmonic message passing module.
The output of the surface preparation module includes (1) the truncated Daplace-Beltrane eigenvec- tors Z R N the corresponding eigenvalues {Ax [==) and the sparse mass matrix B E NRN, (2) per-vertex scalar features Finp E R Nxn
k-1
D THE BINDING SITE PREDICTION MODULE
First, we define the binding site as the protein surface region which is within 3 A to its counterpart surface, and obtain the of all corresponding surface points P the nearest neighbor vertices between the ground truth protein interfaces.
The binding site prediction module stacks feature propagation blocks, each consists of three HMR layers (introduced in Secc. 4.1) and a cross attention layer. Given the propagated receptor features F and ligand features G, the cross attention layer enables communication between proteins:
(FWQ) Vdz
F softmax
(GWv);
(GWa) (FWv),
G' softmax
Vdz
where F NRxd, G Nr xdv denotes the propagated features the receptor/ligand protein surface, dg denotes the dimension of features, and Wq. WK and Wy are the parameter matrices for the query, key, and value in attention computation, respectively.
The loss function consists of two components. The first is a binary cross entropy loss, which en- courages the model correctly predict the binding site:
Lace (i) = ly, log x; + (1 - v.) log(1 x;)],
where and v, are the label and predicted probability of whether vertex belongs to the binding site.
The second term is PointInfoNCE loss (Xie et al. 2020). contrastive matching loss that mini- mizes the distance between the features of corresponding surface points and maximizes the distance between non-corresponding point features:
exp(f, g,/T)
Σ log El,neP exp(f, g+/T)
mmttch
(KJ))
where the set all corresponding surface points and is the temperature factor (a hyperparam- eter). Here f; and gy the neural network-learned feature vectors at point and j, which belong the the receptor and ligand surface, respectively.
The total loss is the weighted sum of the two loss terms:
[ Lhee + Xmmatch
where empirically we set to 0.1 our docking experiments.

---

Published as a conference paper at ICLR 2023
E DETAILS ON FUNCTIONAL MAPS
In this section, we present the details of functional maps used in Sec. 4.2
Let us be given two manifolds M and N. The aim of functional maps is to find bijective mapping T: M + N to align these two manifolds. Unlike traditional methods that try to recover the point- to-point correspondence directly, functional maps lift the mapping T to correspondence between the functional spaces on the two manifolds. Formally, let (-) be the functional space of square integrable functions manifold, infer the functional correspondence Tp (2M (N) induced by the mapping T', and is defined by Tpf = f T for any f E (MM..
To compute the functional correspondence, we need to utilize the Laplace-Beltne Beeltce- basis on the manifolds. Actually, such functional correspondence has a concise expression in spectral domain: k2
M given the respective truncated LB eigenfunctions { 1320 on M and { 1i20 on N, the func- tional correspondence Tp can be (approximately) represented as a change of basis matrix:
C (Ci))aaxk = ((d) Trd)
( (
Now given a of q corresponding functions ffi...., , fg) C (MM and f91.....9. C L2(N), we denote their spectral representations by coefficients A (ai))ki where ai (0r fz) and B (bey (2p where bız (0, .9)N.. The matrix can be obtained by solving the following quadratic minimization problem:
C argminceR |[C A B + aalc 88 8N CII-- + + E IC -iM Ay ClI (E.1)
where ||- || denotes the Frobenius norm, and the first term on the right hand side is the change of basis constraint. Two more regularization terms are introduced into the formula with tunable weights a, 3 > 0. |[C- 6. - ōN C/|r enforces the isometry of the two manifolds where the matrices dM and da are the spectral representation the LB operators. The matrices M and Age- are called the orientation operator (Ren et al., 2018) and are defined by the functions f; and Ji. respectively. The commutator with the orientation operator incorporates extrinsic properties into the formulation and enforces the orientation of functional maps (i.e., which side of the 2D surface is pointing "outward").
Once the functional mapping (i.e., the C matrix) is recovered, point-to-point correspondence T between manifold M and A can be obtained by mapping indicator functions of vertices M to the corresponding functions on A because Trom ÓT(m) for any vertex m E MM (Ovsjanikov et al. 2012).
F DETAILS ON QM9 PROPERTY REGRESSION
Dataset For the QM9 molecular property regression task, we align input data with EGNN (Satorras et al. 2021) (with a total 130,831 instances). 181 molecules failed during molec- ular surface extraction (MSMS failed generate reasonable surface mesh for some molecules). We follow the same data split as EGNN, leading to 130,650 instances (99,862 for training. 17,719 for validation, and 13,069 for test). Only the atomic number and atomic positions are used as the ini- tial molecular information, no extra handcrafted features are involved. We compute the molecular surface as described in Appendix C with an average of 439 vertices and 47 eigenfunctions per molecular surface.
Model Architecture and Performance Different from other graph-based models, we do not ex- plicitly encode inter-atomic distance or chemical bonding information. Instead, we feed the molec- ular surface manifold as the model input, which contains local chemical information (Appendix C). Surface geometric features (i.e., curvature and the Heat Kernel Signatures) are not included, since we find these features have no contribution to predicting molecular properties. We stack layers of HMR, followed by global average pooling layer to aggregate information from all surface points to make a final property prediction. Batch normalization is applied to all MLPs. The mean absolute prediction error of all 12 properties shown in Table F.I results averaged over three independent runs with different random seeds.

---

Published as a conference paper at ICLR 2023
Table F.1: Model performance on the QM9 dataset, reporting the Mean Absolute Error (MAE).
Task
a:
Á
EhOMO
ELUMO
C.
G
H
y²
U
Uo
ZPVE
Units
bohr
meV
meV
meV
D
cal/(mol K)
me V
meV
bohr²
me
meV
meV
NMP
.092
69
43
38
.030
.040
19
17
.180
20
20
1.50
SchNet
235
63
41
34
.033
.033
14
14
.073
19
14
1.70
Cormorant
.085
61
34
38
.038
.26
20
21
.961
21
22
2.02
SE(3)-Tr.
.122
53
33
.051
.054
EGNN
.071
48
29
25
.029
.031
12
.106
12
1.55
SEGNN
.060
42
24
21
.023
.031
15
16
.660
13
15
1.62
HMR
.102
59
37
31
.037
.040
40
41
.653
44
42
2.97
Resource Consumption Computing the molecular surface and its Laplace- Beltrami eigenfune- tions is computationally intensive. We performed data preprocessing in parallel on 64 CPUs, where molecular surface computation took 198 seconds, surface mesh refinement took 59 minutes, and solving eigenfunctions took 4.5 hours. With average of 439 mesh vertices and 47 eigenfunc- tions, the entire QM9 dataset consumes 12.1 GB disk space (approximately kB per molecule). The model for production contains 356.609 learnable parameters, most of which are linear trans- formation coefficients in MLPs. However, manifold harmonic message passing involves matrix multiplications bringing features back-and-forth between the real space and the generalized Fourier space, which is performed in serial manner (instead of batch processing). This is because the num- ber of surface vertices and Laplace-Beltrami Therefore, we only perform batch operations on the MLPs, but not message passing layers. We trained our model on NVIDIA A100 GPUs with 80 GB memory with a batch size of 32, which on average takes 240 seconds train single epoch (99,862 molecules), and 16 seconds for inference on the test set (13,069 molecules). For each molecule we keep its N x n feature matrix (N ver- tices and features), and the N x eigenvectors (Z matrix) and its inverse matrix (Z matrix, Appendix B), which does not consume much GPU memory.
Beltrami eigenfunctions are different across different molecules.
G DETAILS ON LIGAND-BINDING POCKET CLASSIFICATION
Dataset The dataset obtained from Gainza (2020) with 1,438 non-redundant protein struc- tures that each corresponds to list of bound ligands with their atom coordinates. We first generate protein surface meshes as described in Appendix C extract ligand binding pockets, we identify pocket vertices on the surface mesh that are within 4 A to any atom of the ligand and extract the largest connected component of pocket vertices as the binding pocket. Binding pockets that contain < 100 vertices are removed and the LB eigenfunctions for the remaining pocket surfaces are calcu- lated. protein complexes failed in protein surfaces generation and 63 protein complexes failed in binding pocket extraction due disconnected surface or too few vertices after surface refinement. In total, 2,254 binding pockets are extracted that bind to one of seven ligands: ADP (634), HEM (396), FAD (338), COA (263), NAD (232), SAM (222), and NAP (169). On average, 1,199 ver- tices and 190 eigenfunctions per pocket surface are obtained. The same training, validation, and test split as Gainza al. (2020) are used and obtained 1,634 training pockets (in 986 protein complexes), 202 validation pockets (in 112 protein complexes), and 418 test pockets (in 274 pro- tein complexes). Similar Gainza al. (2020). we provide chemical information (hydrophobicity score, partial charge) as well as local geometric information (mean/ Gaussian curvature and the Heat Kernel Signatures) as input features.
Model Architecture and Performance The HMR-based classification model contains 6 layers of HMR propagation layers followed by global average pooling aggregate information from all surface points of pocket. A simple 2-layer MLP is used classify pockets into seven ligand- binding classes. HMR trained to minimize the cross-entropy loss for 400 epochs (approx. 20,000 iterations with batch size of 32) and the with the best balanced accuracy score the validation set is selected. To compare with MaSlF-ligand, we calculated the balanced accuracy for multi-class classification (Fgg. 6). Per class performance of our full model (geometric chemical features) is shown as a confusion matrix Fig. G.T.
20

---

Published as a conference paper at ICLR 2023
ADP
0.92
0.00
0.01
0.02
0.03
0.00
0.02
COA
0.07
0.62
0.30
0.12
0.10
0.07
0.03
0.8
Ground truth
FAD
0.00
0.04
0.87
0.03
0.04
0.01
0.03
0.8
HEM
0.90
NAD
0.00
0.00
0.4
SAM NAP
0.05
3.05
0.30
0.00
0.18
0.73
0.00
0.24
3.00
0.00
0.00
0.07
0.03
0.79
ADP
CDA
FAD
HEM Predicted
NAD
Figure G.1: Confusion matrix of ligand specificity predicted by HMR on the test set.
Resource Consumption Similar to QM9, we performed data preprocessing in parallel on 64 CPUs. For all 2,254 protein pocket cases: the protein surface computation took 306 seconds, surface pocket extraction and refinement took 330 seconds, and solving the eigenfunctions for the extracted pockets took 26 minutes. The processed dataset (N=2,254) takes 2.4 GB disk space (about 1.1 MB protein pocket). The full model contains 407,047 trainable parameters. We trained our model on NVIDIA A100 GPUs with 80 GB memory with a batch size 32, which on average takes 33 seconds to train single epoch (1,634 protein pockets), and 9 seconds for inference on the test set (418 protein pockets). The model structure is almost identical to QM9 that shares the same batch processing design and has similar GPU memory consumption.
H DETAILS ON RIGID PROTEIN DOCKING
H.1 DATASET
The original DIPS dataset (Townshend et al. 2019) contains interacting protein chains extracted from experimentally determined complex structures in the RCSB PDB database (https: www. resb.org// While DIPS is designed represent the interactions in protein complexes including homomeric proteins (i.e., protein complexes composed of identical proteins), we are particularly interested in interactions between different proteins in biological processes (e-g., antibody-antigen, enzyme- -niiiie-n enzyme-substrates, etc.). To better resemble such interactions, we apply a modi- fied pipeline based on DIPS, referred as DIPS-Het.
Specifically, we include the PDB entries in DIPS as well as newly deposited entries from the RCSB PDB database that (1) are determined using diffraction-bart methods electron microscopy, (2) have < 3.5 A resolution, (3) are not hybrid protein complexes (e.g., protein-DN. complexes are excluded), and (4) have
protein chains. Different from the original DIPS pipeline, we further remove the PDB entries with only homomeric interfaces (i.e., interactions between identical pro- teins), and b) adopt a one-vs-rest strategy that uses each one of the protein chains in complex the "ligand protein" and the rest protein chains as the "receptor protein" form ligand-recepenn pairs.
Similar to Ganea et al. 2021) 1, we remove the cases with any protein chain sharing the same protein sequence clusters (at 30% sequence similarity) with DB5.5. Finally, we cluster proteins based on protein chain sequence similarity and separate these clusters into the training and validation set. This sequence-based split helps selecting a model with better generalizability. After removing proteins failed in surface mesh generation, DePS-Het contains 11,373 training cases and 508 validation cases.
H.2 EVALUATION METRICS
We evaluate the rigid protein docking results using four metrics: Complex RMSD, Interface RMSD, DockQ, and Success rate. Complex RMSD and Interface RMSD are used in Ganea et al. (2021) let
21

---

Published as a conference paper at ICLR 2023
Z' E (3((t))) and Z E (3x(m+m) be the a-carbon coordinates of the ground truth and predicted protein complexes, respectively, were and are the number of a-carbons in the receptor and ligand protein. After superimposing the complex structures using the Kabsch algorithm, Complex
I 2 Similarly, Interface RMSD is calculated using the same RMSD calculated as m+
procedure but with the a-carbon coordinates of interface residues (< 8 Á to the other protein's residues). Smaller RMSD value means the predicted structure is closer the ground truth structure.
In addition, we evaluate the overall quality of docking using DockQ (Basu & Wallner 2016) and
the success rate of achieving "Acceptable" or higher according to Méndez al. (2003 2005).
Both metrics are based on three standardized criteria used by Critical Assessment of PRedicted Interactions (CAPRI): Lens the ligand (the smaller protein) RMSD calculated based on backbone atoms, after superimposing receptor's 's backbone atoms; I is the backbone RMSD of interface residues, after superimposing the interface residues (residues with any atom is 10 Á to atoms in the other protein); and Jnat is the recall in recovering residue-residue contacts between the proteins, where two residues are "in contact" if any pair of atoms from two residues has distance < Å. DockQ is continuous score between O and I (the higher the better), derived from Lyms, Luss and fnat- A docking result is considered as a "success" if ranked "Acceptable" or higher according to CAPRI's criteria, that is
fma 2 0.1 (Lems S 10.0 V Amss I 4.0)
OR
fmn 0.3
Both DockQ score and CAPRI ranking is calculated using the DockQ package (https:/) ctchub..com/boonoo wwwerwawlleer/Dockkeee
H.3 RESOURCE CONSUMPTION
Dataset We report the rigid protein docking resource consumption using the Docking Benchmark 5.5 (DB5.5) dataset (Guest et al. (021), which contains 253 pairs of representative protein complex structures. This is the gold-standard test set to evaluate protein docking model performance.
We preprocessed the DB5.5 dataset in parallel 64 CPUs. Adding hydrogen atoms the PDB data took 34 seconds (using reduce), computing the solvent-excluded molecular surface with MSMS took 6 seconds, triangle mesh refinement by PyMesh took 3.4 minutes, and computing the Laplace- Beltrami eigenfunctions took 18 minutes with s scipy (eigsh The average number of vertices per protein (i.e., ligand receptor) 3,380, with 215 calculated Laplace-Beltrami eigenfunctions (maximum eigenvalue capped to 0.3, which is determined empirically). The total storage space for the processed dataset is GB (about 8 MB per protein, much larger than QM9 molecules). For our training set with 11 11,781 protein complexes, the storage space is 101 GB.
Inference The HMR inference time in predicting DB5.5 proteins complexes is shown in Table H. in comparison to EQUIDOCK, HDOCK, and ATTRACT. For HMR, data preprocessing took 80% of inference time, functional maps (docking pose prediction) took 19%, while binding site prediction took only 1% of time.
Table H.I: Rigid protein docking inference time on DB5.5 dataset (per protein complex docking time averaged over 253 cases).
Method
Time
HDock
884.9 sec
ATTRACT
882.7 sec
EquiDock
3.9
HMR (Top 1)
5.8 sec
HMR (Top 3)
6.6 sec
SSolving the generalized eigenvalue problem should be done with less parallel processes to allocate more CPU resources to each solver for better efficiency.
22

---

Published as a conference paper at ICLR 2023
H.4 DISTRIBUTIONS OF LEARNED PROPAGATION PARAMETER
We visualize the distribution of propagation parameters: propagation time frequency filter mean and variance o, learned from data. The histogram of selected layers are shown in Fig. H.I
Propagation Layer
60
Channels
Chennia
20
40
10
10
8.0
Propagation
Fihor
Propagation Layer 3
75
Charels
Cheradd
60
28
20
0.0
0.2
0.4
8.0
0.1
0.2
0.3
0.4
Propagation timw
Filer
Filtor d
Propagation Layer
100
Clames
90
Chemaas
Chomes
20
80
110
0.0
0.2
8.0
Propagation
Filer
Figure H.I: Distribution of learned propagation time / and frequency filter (f and ()) at selected propagation layers.
H.5 ANALYSIS OF DIFFERENT MESSAGE PASSING MECHANISMS
We study the effectiveness of proposed Harmonic Message Passing mechanism by comparing with spatial graph message passing, a "Graph Mcssage Passing" model where Harmonic Massage Passing is replaced by Graph Attention Network (Veličković al. 2018), and with "No Message Passing" model where only a multilayer perceptron (MLP) is used to encode per-vertex features. Experimental results show that the Harmonic Message Passing model outperforms Graph Message Passing and No Message Passing models, suggesting its higher efficiency in propagating information on protein surface meshes.
Table H.2: Rigid Protein Docking Performance with Different Message Passing Models
Complex RMSD
Interface RMSD
DockQ
Model
Median
Mean
Median
Mean
Median
Mean
Harmonic Message Passing
12.01
12.09
10.90
11.43
0.06
0.22
Graph Message Passing
15.41
15.67
13.38
14.37
0.03
0.08
No Message Passing
16.44
17.23
14.11
16.11
0.02
0.04
H.6 ABLATION STUDIES ON INPUT FEATURES
To assess the importance of input chemical and geometric features for rigid protein docking, we compare the performance of HMR ("Full model") to models trained with only chemical or geo- metric features: the "Chem only" model only includes chemical features (i.e., atom types, residue types, residue's hydrophobicity, and whether the atom is an a-carbon); the "Geom only" model only contains features related to manifold geometry, including Gaussian curvature, mean curvature, and Heat Kernel Signatures. Table H.3 shows that the model performance drops upon removing either chemical geometric features. And the chemical features is more critical for HMR in the rigid protein docking task.
23

---

Published as a conference paper at ICLR 2023
Table H.3: Ablation studies on input features
Complex RMSD
Interface RMSD
DockQ
Model
Median
Mean
Median
Mean
Median
Mean
Full model
12.01
12.09
10.90
11.43
0.06
0.22
Chem only
13.96
13.95
12.73
13.33
0.04
0.17
Geom only
17.87
19.15
16.58
18.46
0.01
0.03
H.7 CROSS EXPERIMENTS WITH EQUIDOCK
Compared to EQUIDoCK, HMR differs both in the model and the dataset used for training. To un- derstand which aspect contributes more the performance of HMR, conduct to train HMR on the original DIPS (used in EQUIDOCK) and EQUIDOCK trained on DIPS-Het (pro- posed by us, used in HMR). The same training and validation sets are used, and all models evaluated on DB5.5. As summarized in Table H.4, model structure and training dataset are both important: HMR achieved higher performance than EQUIDocK when both models are trained on DIPS dataset; the performance of HMR further improves when it is trained on DIPS-Het dataset.
experiments
Table H.4: Compare EQUIDoCK and HMR trained DIPS or DIPS-Het datasets
Complex RMSD
Interface RMSDD
DockQ
Model
Training data Median
Mean
Median
Mean
Median
Mean
DIPS
15.68
16.04
14.81
15.38
0.02
0.07
HMR
DIPS-Het
12.01
12.09
10.90
11.43
0.06
0.22
DIPS
17.27
18.09
15.04
16.11
0,02
0.03
EQUIDocK
DIPS-Het
17.22
17.65
14.27
14.76
0.02
0.04
H.8 DB5.5 TEST CASES FROM EQUIDOCK
For fair comparison, further compare HMR with EQUTDocK on the 25 cases in DB5.5 selected as a test set in Ganea et al. (202T). We report the EQUIDocK model fine-tuned on DB5.5, as orig- inally reported Ganea et (2021). As shown in Table H.S HMR still outperforms EQUIDOCK, despite not being fine-tuned on DB5.5.
Table H.5: Rigid docking performance DB5.5 test set
Complex RMSD Interface RMSD
DockQ
Success rate
Model
Median
Mean
Median
Mean
Median
Mean (> Acceptable)
HMR (Top 1)
13.10
12,42
11.05
11.68
0.07
0.15
0.20
HMR (Top 3)
11.35
11.18
10.02
10.03
0.08
0.16
0.20
EQUIDOCK
14.13
14.72
11.97
13.23
0.04
0.05
0.00
H.9 MODEL PERFORMANCE DIFFERENT RESOLUTIONS
We perform the rigid docking experiment using HMR at different resolutions. Specifically, we cap the largest eigenvalue the ligand and receptor protein surface manifolds used the model. Since eigenfunctions with larger cigenvalues exhibit higher spatial resolutions, restricting the eigenvalues effectively constrains the model resolution. The results are shown in Table H.6 In general, we observe worse model performance with fewer eigenfunctions for harmonic message passing and functional maps. Interestingly, the metrics of the binding site prediction module (i.e., AUC and AP)
24

---

Published as a conference paper at ICLR 2023
are less sensitive to resolution than the rigid docking module (powered by functional maps). This means the model could still infer where the binding site is with lower resolution, yet the quality of learned functional correspondence decreases, leading to worse docking power.
Table H.6: HMR rigid protein docking performance at different spatial resolutions (quantified by the largest eigenvalue of surface manifold). AUC (area under the receiver operating characteristic curve) and AP (average precision) are metrics of the binding site prediction module.
max eigenvalue
AUC
AP
CRMSD (Top 1)
CRMSD (Top 3)
0.30
0.84
0.54
12.11
8.36
0.25
0.84
0.53
13.14
8.9
0.20
0.84
0.54
11.91
9.53
0.15
0.83
0.53
12.19
9.91
0.10
0.82
0.52
12,90
10.02
0.05
0.51
13.84
11.32
H.10
TRAINING DETAILS
Table H.7: Hyperparameter choices of HMR and the training phase settings
Hyperparameters
Values
Binding Site Prediction Module
Number of Feature Propagation Blocks
2
Number of HMR Layers
Dimension of Propagated Features d;
128
Number of Attention Heads
Dropout Rate
0.1
NCELoss Temperature
10
NCELoss Number of Point Samples
50
NCELoss Weight
0.1
Training
Batch Size
4
Epoch Learning Rate
80
5x10-4
Optimizer
Adam
Learning Rate Scheduler
Cosine Annealing
F-score
0.4
0.8
Valu
0.2
0.1
0.5
Salit
AUC
vold
0.95
0.8
0.85
0.6
Fooch
Epach
Figure H.2: HMR learning curve for the rigid protein docking task. Model with the best validation average precision (AP) score in binding site prediction selected for testing on DB5.5.
25


Reviews:
[{"review_id": "t2QSn3ccUx", "reviewer": "Reviewer_tWAd", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "The authors study the molecule representation problem from a novel aspect based on the surface Riemannian manifold. Specifically, the authors adopt the Shape-DNA techniques to obtain surface representation and propose the information aggregation for the surface. The authors further design the frameworks for molecular property prediction and binding prediction.", "summary_of_the_review": "Overall I like the idea in this paper. Using representation and message passing of surface is a pretty natural approach to the binding problem. The proposed idea is novel. However, there are concerns that some important claims are not fully justified and experimental results are not strong enough.", "strength_and_weaknesses": "\n### Strengths\n\n(+) The studied problem of the “shape-activity” relationship as an alternative to the “structure-activity” relationship is interesting and can provide some insights into understanding the system.\n\n(+) The proposed method itself, including surface representation and information aggregation, is interesting and novel for geometric graph learning. I’m sure it can also be applied to many other scenarios than molecules.\n\n(+) The proposed rigid protein docking framework has additional technical contributions.\n\n### Weaknesses\n\n(-) The use of surface only is not well-motivated nor justified. Although the authors claim that “molecular surface displays chemical and geometric patterns which fingerprint a protein’s mode of interactions with other biomolecules”, some analysis or empirical results are needed to justify the claim and make the proposed surface representation more convincing. Moreover, I agree that the surface may be important for binding problems to some degree, but the claim “this possibly explains the success of predicting these properties using only molecular surface information, since the surface is approximately an isosurface of the molecular electron density field” is not convincing enough and may need additional justifications.\n\n(-) Current experimental results are not convincing enough to show the significance and effectiveness of the proposed methods. The experimental results on QM9 do not show a significant improvement to the baselines (some are worse). It indicates the practical application of the proposed method may be limited to only a few specific tasks. For the binding task, the EquiDock can be a relatively weak baseline. The authors may want to include stronger baselines such as EquiBind [1].\n\n(-) There is a lack of details or comparisons on the method's computation complexity in terms of runtime and memory cost. Also, it would be better if the authors include results on how increasing the dimensions (surface resolution) can improve the performance and increase the computation cost.\n\n[1] Stark et al. EquiBind: Geometric Deep Learning for Drug Binding Structure Prediction. ICML 2022.\n", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "2: Several of the paper’s claims are incorrect or not well-supported.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "I8YTCjfaL5o", "reviewer": "Reviewer_tNJq", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "This paper proposes propose a Harmonic Molecular Representation learning (HMR) framework,which offers a multi-resolution representation of molecular geometric and chemical properties on 2D Riemannian manifold. In order to realize efficient spectral message passing over the surface manifold for better molecular encoding, a harmonic message passing method is proposed. The proposed method outperforms the state-of-the-art deep learning models for the rigid protein docking challenge, and shows comparable predictive power to current models in small molecule property prediction tasks.", "summary_of_the_review": "The paper is well written in general, and the proposed method was examined in comparison with state-of-the-art in various experimental analyses.\n\nThe recommended updates can improve presentation of the paper and elucidate dynamics of the proposed learning methods and networks.", "strength_and_weaknesses": "The paper addresses an important problem of encoding 3D molecular structures in deep neural networks. The paper is well written and the proposed methods are examined in various tasks.\n\nSome parts of the presentation can be improved and additional results can be provided as suggested below.", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "9Zgh9ZOPKv", "reviewer": "Reviewer_GJXL", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "This paper develops a harmonic representation of a molecular surface using the Laplace-Beltrami eigenfunctions. The approach is applied to predict the properties of small molecules and predict the docking pose of proteins.", "summary_of_the_review": "In summary, this paper develops a novel harmonic representation of the molecular surface and applies the approach to two tasks: QM9 and protein rigid docking. Despite the potential improvement in the rigid docking task, the HMR approach may benefit from the additional geometric and chemical features in the surface representation. Comparing it to another surface representation baseline which encodes the same information is needed to demonstrate the advantage of the harmonic representation from the reviewer’s perspective.", "strength_and_weaknesses": "Strengths:\n- The Laplace-Beltrami operator as a harmonic filter provides a principled approach to encode interactions at different resolutions of molecular surfaces.\n- The model outperforms EquiDock in the protein rigid docking task. \n\nWeaknesses:\n- The comparison with GNN baselines is not very informative. It is well-known that molecular surface representation captures very different information compared with GNNs. The HMR also uses surface geometry features and chemical features which are not present in the GNN baselines. It is not surprising that HMR can outperform the GNN based EquiDock baseline due to the additional information. I think comparing it with another protein-surface representation baseline like MaSIF [Gainza et al., Nature Methods 17.2 (2020): 184-192.] or Multi-Scale Protein Representation [Somnath et al., NeurIPS 2021] is needed to demonstrate the usefulness of this new molecular surface representation approach.\n- In Table 2, the runtime of the HMR approach is not reported. The key advantage of EquiDock is its significantly faster inference compared conventional docking software. The surface construction might be slow for the HMR method. How does the inference time compare with EquiDock, Attract, and Hdock?\n", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "-wTEY7Yhm4", "reviewer": "Reviewer_R5jV", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "This paper proposes a framework to learn Riemannian manifold representation using harmonic message passing. To this end, the authors (1) pre-process molecules (or macromolecules like proteins) as a surface function using triangulation, (2) decompose the surface function into components with LB eigenfunctions and coefficients, and (3) perform spectral message passing between the component-wise features. Importantly, the proposed harmonic message passing is motivated by the interpretation of heat diffusion as a message passing between surface regions (LB eigenfunctions and coefficient) with different temperatures (eigenvalues). The proposed methodology is applied for molecular property prediction, binding site prediction, and rigid docking simulation. The empirical results are promising and competitive with respect to the recently proposed molecular representation learning algorithms. ", "summary_of_the_review": "This paper proposes a very novel idea for molecular representation learning. While the experiments could be slightly improved, I think this paper is potentially very impactful and interesting for molecular representation learning community ", "strength_and_weaknesses": "Strength. \n\nThis paper proposes a novel molecular representation learning framework inspired by molecular surface representation and harmonic analysis. To my knowledge, various geometric deep learning works have focused on message passing between mesh vertices rather than the surface regions (LB eigenfunctions). If my knowledge is correct, I believe this work to be quite significant. \n\nI also appreciate how the experiments are done on both molecules & proteins with various tasks (molecular property prediction, binding site prediction, rigid docking simulation). This makes the empirical benefit of the proposed method convincing.\n\nThe proposed method looks quite scalable, so I think this method might be potentially useful in large-scale scenario where computation time is important.\n\nWeakness. \n\nIn the experiments, it is not clear whether if the complexity of preprocessing molecular surface is considered. I find this confusing since the EquiDock authors mentioned how they avoided molecular surface preprocessing due to high complexity. \n\nThe experimental results for molecular property prediction is not exactly state-of-the-art, but I think this is fine given novelty of the work. Nevertheless, I suggest comparing with more recent works like \"Geometric and physical quantities improve e(3) equivariant message passing, ICLR 2022\" for molecular property prediction (QM9) for completeness of the paper.\n\nI think clarity the description of harmonic message passing could be improved. For example, the input and the output of harmonic message passing is expressed as a function. However, to my understanding, the features being actually updated correspond to different coefficients for the LB eigenfunctions. I think the matrix representation of the function operator (Fig 3) should be described in the text. Also, the symbols in Fig 3 has no corresponding explanation, e.g., what is $\\Phi$? \n\nThe authors claim this method to be scalable to large scale molecules and proteins. However, although the logic is convincing, the authors do not design any experiements to support this fact. When HMP is compared to Equidock in term of preprocessing + inference time, is it faster?\n\n(minor) I find that the experimental setting differs from Equidock. Could you explain the reasoning behind this?\n\n\n", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}]

Author rebuttal:

Decision: Accept: poster