Dynamic-Backbone Protein-Ligand Structure Prediction with Multiscale Generative Diffusion Models

Molecular complexes formed by proteins and small-molecule ligands are ubiquitous, and predicting their 3D structures can facilitate both biological discoveries and the design of novel enzymes or drug molecules. Here we propose NeuralPLexer, a deep generative model framework to rapidly predict protein-ligand complex structures and their fluctuations using protein backbone template and molecular graph inputs. NeuralPLexer jointly samples protein and small-molecule 3D coordinates at an atomistic resolution through a generative model that incorporates biophysical constraints and inferred proximity information into a time-truncated diffusion process. The reverse-time generative diffusion process is learned by a novel stereochemistry-aware equivariant graph transformer that enables efficient, concurrent gradient field prediction for all heavy atoms in the protein-ligand complex. NeuralPLexer outperforms existing physics-based and learning-based methods on benchmarking problems including fixed-backbone blind protein-ligand docking and ligand-coupled binding site repacking. Moreover, we identify preliminary evidence that NeuralPLexer enriches bound-state-like protein structures when applied to systems where protein folding landscapes are significantly altered by the presence of ligands. Our results reveal that a data-driven approach can capture the structural cooperativity among protein and small-molecule entities, showing promise for the computational identification of novel drug targets and the end-to-end differentiable design of functional small-molecules and ligand-binding proteins.


Introduction
Protein structures are dynamically modulated by their interactions with small-molecule ligands, triggering downstream responses that are crucial to the regulation of biological functions [1][2][3]. Proposing ligands that selectively target protein conformations has become an increasingly important strategy in small-molecule-based therapeutics [4][5][6]. However, computational prediction of protein-ligand structures that are coupled to receptor conformational responses is still hampered by the prohibitive cost of physically simulating slow protein state transitions [7,8], as well as the static nature of existing protein folding prediction algorithms [9,10]. While several schemes have been proposed to remedy these issues [11][12][13][14][15][16][17][18][19][20], such methods often require case-specific expert interventions and lack a unified framework to predict 3D structures in a systematic and cooperative fashion.
Here we propose NeuralPLexer, a Neural framework for Protein-Ligand complex structure prediction. NeuralPLexer leverages diffusion-based generative modeling [21,22] to sample 3D structures from a learned statistical distribution. We demonstrate that the multi-scale inductive bias in biomolecular complexes can be feasibly integrated with diffusion models by designing a finite-time stochastic differential equation (SDE) with structured drift terms. Owing to this formulation, NeuralPLexer can generalize to ligand-unbound or predicted protein structure inputs once trained solely on experimental protein-ligand complex structures that are not paired to alternative protein conformations. When applied to blind protein-ligand docking, NeuralPLexer improves both the geometrical accuracy and structure quality compared to baseline methods; when applied to ligand binding site design, an inpainting version of NeuralPLexer can Sampling from NeuralPLexer. The protein (colored as red-blue from N-to C-terminus) and ligand (colored as grey) 3D structures are jointly generated from a learned SDE, with a partially-diffused initial state q T * approximated by the protein backbone template and predicted interface contact maps. (c-e) Key elements of the NeuralPLexer technical design. (c) Ligand molecules and monomeric entities are encoded as the collection of atoms, local coordinate frames (depicted as semi-transparent triangles), and stereospecific pairwise embeddings (depicted as dashed lines) representing their interactions. (d) The forward-time SDE introduces relative drift terms among protein Cα atoms, non-Cα atoms and ligand atoms, such that the SDE erases local-scale details at t = T * to enable resampling from a noise distribution. (e) Information flow in the equivariant structure diffusion module (ESDM). ESDM operates on a heterogeneous graph formed by protein atoms (P), ligand atoms (L), protein backbone frames (B) and ligand local frames (F) to predict clean atomic coordinatesx 0 ,ŷ 0 using the coordinates at a finite diffusion time t > 0.
accurately repack 44% of failed AlphaFold2 [9] binding sites with up to 60% success rate improvements compared to the method in Rosetta [23]. Furthermore, NeuralPLexer only requires molecular graphs as ligand inputs, therefore can enable end-to-end gradient-based design for functional small-molecules and ligand-binding proteins when coupled to recently-proposed differentiable protein sequence [24][25][26] and molecular graph generators [27,28].

Method
We assume the model inputs are a receptor protein backbone template containing the amino acid sequence s and (N, Cα, C) atomic coordinatesx ∈ R nres×3×3 , and a set of ligand molecular graphs {G k } K k=1 containing atom/bond types and stereochemistry labels (e.g., tetrahedral or E/Z isomerism [29]). We aim to sample (x, y) ∼ q φ (·|s,x, {G}) from a generative model q φ with predicted 3D heavy-atom coordinates of the protein x ∈ R n×3 and that of the ligands y ∈ R m×3 . It can be understood as a conditional generative modeling problem for partially-observed systems.
NeuralPLexer adopts a two-stage architecture for protein-ligand structure prediction (Figure1a). The input protein backbone template and molecule graphs are first encoded and passed into a contact predictor that iteratively samples binding interface spatial proximity distributions for each ligand in {G}; the output contact map parameterizes the geometry prior, a finite-time marginal of a designed SDE that progressively injects structured noise into the data distribution. An equivariant structure diffusion module (ESDM) then jointly generates 3D protein and ligand structures by denoising the atomic coordinates sampled from the geometry prior through a learned reverse-time SDE (Figure1b).

Protein-ligand structure generation with biophysics-informed diffusion processes
Diffusion models [22] introduce a forward SDE that diffuses data into a noised distribution and a neural-networkparameterized reverse-time SDE that generate data by reverting the noising process. To motivate the design principles for our biomolecular structure generator, we first consider a general class of linear SDEs known as the multivariate Ornstein-Uhlenbeck (OU) process [30] for point cloud Z ∈ R N ×3 : where Θ ∈ R N ×N is an invertible matrix of affine drift coefficients and W t is a standard 3N -dimensional Wiener process. The forward noising SDEs used in standard diffusion models [31,32] can be recovered by setting Θ = θI, converging to an isotropic Gaussian prior distribution at the t → ∞ (often expressed as t → 1 with reparameterized t [33]) limit. In contrast, we design a multivariate SDE with data-dependent drift matrix Θ(Z 0 ) and truncate the SDE at t = T * < ∞ such that the final state of forward noising process is a partially-diffused, structured distribution q T * that can be well approximated by a coarse-scale model. We propose a set of SDEs depicted by Figure1d and detailed in Table A1, with separated lengthscale parameters σ 1 , σ 2 such that the forward diffusion process erases residue-scale local details but retains global information about protein domain packing and ligand binding interfaces, yielding the following time-dependent transition kernels: where we use an exponential scheduleτ = (σ 2 min /σ 2 1 )e t with truncation T * = 2 log(σ 2 /σ min ). c is a softmaxtransformed contact map as detailed in Sec. 2.2, which attracts the diffused ligand coordinates y(t) towards binding interface Cα atoms while preserving SE(3)-equivariance. We choose σ 1 = 2.0 Å to match the average radius of standard amino acids with task-specific σ 2 > σ 1 such that at t = T * : (a) the terms involving x nonCα (0) and y(0) approximately vanishes thus are set to zeros to initialize the reverse-time SDE, and (b) the Cα-atom coordinate marginal q T * x Cα (t)|x(0) is sufficiently close to which approximated by the backbone template q T * x Cα (t)|x , guided by the theoretical result proposed in [34]. Proofs regarding SE(3)-equivariance are stated in the Appendix A.1.2.

Contact map prediction and sampling from the truncated reverse-time SDE
Given protein-ligand coordinates (x, y), we define the contact map L ∈ R nres×m with matrix elements where j runs over all protein atoms in amino acid residue A and α = 0 To sample from the reverse-time SDE, we use the contact predictor to generate inferred contact mapsL and parameterize the geometry prior q T * (·|x,L) -the initial condition of reverse-time SDE -by replacing x(0) in q T * with the backbone templatex and the ligand-Cα relative drift coefficient c with the predicted c(L). Note that in the general multivariate OU formulation, this corresponds to replacing the clean-datadependent drift coefficients Θ(Z 0 ) by a model estimationΘ. To account for the multimodal nature of protein-ligand contact distributions, the contact predictor models L as the logits of a categorical posterior distribution over a sequence of one-hot observations {l} K k=1 sampled for individual molecules in {G}. The forward pass of contact predictor ψ takes an iterative form: where k ∈ {1, · · · , K} and we setL . . =L K . All results reported in this study are obtained with K = 1 due to the curation scheme of standard annotated protein-ligand datasets, but we note that the model can be readily trained on more diverse structural databases with multi-ligand samples.

Architecture overview
Here we outline the key neural network design ideas and defer the featurization, architecture, and training details to the Appendix. To enable stereospecific molecular geometry generation and explicit reasoning about long-range geometrical correlations, NeuralPLexer hybridizes two types of elementary molecular representations (Figure1c): (a) atomic nodes and (b) rigid-body nodes representing coordinate frames formed by two adjacent chemical bonds. For small-molecule ligand encoding, we introduce a graph transformer with learnable chirality-aware pairwise embeddings that are constructed through graph-diffusion-kernel-like transformations [35]; such pairwise embeddings are pretrained to align with the intra-molecular 3D coordinate distributions from experimental and computed molecular conformers. The protein backbone template encoding module and the contact predictor are built upon a sparsified version of invariant point attention (IPA) adapted from AlphaFold2 [9] and are combined with standard graph attention layers [7,36] and edge update blocks.
The architecture of ESDM (Figure1e) is inspired by prior works on 3D graph and attentional neural networks for point clouds [37,38], rigid-body simulations [39] and biopolymer representation learning [9,[40][41][42]. In ESDM, each node is associated with a stack of standard scalar features f s ∈ R c and cartesian vector features f v ∈ R 3×c representing the displacements of a virtual point set relative to the node's Euclidean coordinate t ∈ R 3 . A rotation matrix R ∈ SO(3) is additionally attached to each rigid-body node. Geometry-aware messages are synchronously propagated among all nodes by encoding the pairwise distances among virtual point sets into graph transformer blocks. where NeuralPLexer accurately inpaints the binding site protein-ligand structure, while directly aligning AlphaFold2 prediction to the ground-truth complex resulted in steric clashes between the ligand and binding site residues. (f) Summary of binding site accuracy (measured by the all-atom lDDT-BS score) and ligand clash rate over the test dataset. 32 conformations are sampled for each protein-ligand pair; dots indicates the median value and errorbars indicates 25% and 75% percentiles. (g) Success rates compared to baseline methods. A success is defined as: lDDT-BS > 0.7, ligand RMSD < 2.0 Å, and clash rate = 0.0. The pink "true contact map" curves are obtained by initializing the geometry prior q T * using the true protein-ligand contact map, while the gold curves are obtained by generating both protein and ligand conformations end-to-end.
Explicit non-linear transformation on vector features f v is solely performed on rigid-body nodes through a coordinateframe-inversion mechanism, such that the node update blocks are sufficiently expressive without sacrificing equivariance or computational efficiency. On the contrary, 3D coordinates are solely updated for atomic nodes while the rigid-body frames (t, R) are passively reconstructed according to the updated atomic coordinates, circumventing numerical issues regarding fitting quaterion or axis-angle variables when manipulating rigid-body objects. The nontrivial actions of a parity inversion operation on rigid-body nodes ensure that ESDM can capture the correct chiral-symmetry-breaking behavior that adheres to the molecular stereochemistry constraints.

Results
Fixed-backbone protein-ligand docking. In this setting the ground-truth receptor protein backbone is given as input x, and both ligand coordinates and protein sidechain coordinates are predicted using the input protein backbone and ligand graphs. Results are compared to a recent learning-based method EquiBind [43]; for reference, we also include results from a physics-based blind docking method CB-Dock [44] obtained with ground-truth all-atom receptor inputs and using a computing budget similar to learning-based methods. Models are trained and tested on the PDBBind-2020 [45] dataset split used in [43], with additional test dataset processing to ensure a reasonable comparison to docking-based methods (see Appendix A.8.1). As shown in Figure 2a-c, NeuralPLexer achieves both improved geometrical accuracy (reported as the ligand heavy atom root-mean-squre-deviation (RMSD)) and lower steric clash rate (the fraction of ligand heavy atoms with a Lennard-Jones energy > 100 kcal/mol, using UFF [46] parameters). We found that good ligand structure quality and geometrical accuracy can be achieved using as few as 10 integrator steps (0.2 second per conformation on a single V100 GPU).
Ligand-coupled binding site repacking. Here we apply a diffusion-based inpainting strategy to jointly sample ligand and protein structure for a cropped region within 6.0 Å of the ligand conditioning on the uncropped parts of the protein. Protein binding site accuracy is measured by the lDDT-BS metric [47] with cutoff parameters consistent with CAMEO [48]. Input backbones are obtained using template-free AlphaFold2 (AF2) predictions of 154 selected chains whose TM-score [49]>0.8 and lDDT-BS<0.9 out of the abovementioned PDBBind test set, a subset representing cases where AF2 correctly predicts the global protein folding but unable to reproduce the exact bound-state binding site structure. We found 82% of structures contain steric clash with the ligand when directly aligned to reference complex structure in PDB, while NeuralPLexer is able to rescue 44% of these AF2 binding sites with joint protein-ligand inpainting (Figure 2e-g). Comparing to an energy-based flexible ligand-receptor modeling method RosettaLigand [23], NeuralPLexer increases success rate by up to 60% on the combined metric for ligand accuracy, binding site accuracy and physical plausibility.
Cryptic pockets and binding-induced protein conformation transitions. Lastly, we assessed NeuralPLexersampled structures for 31 systems from the PocketMiner dataset [50] which represents proteins with substantial ligand-binding-induced conformation changes. As a preliminary examination, we use the ligand-unbound (apo) crystal structure from PDB as the input backbone template and fix the ligand conformation to ground-truth coordinates along sampling. We found NeuralPLexer shifts the sampled ensemble toward bound-state (holo) structures when performing joint protein-ligand generation, compared to unconditioned protein-only sampling results (Figure 3a). Human evaluations reveal that NeuralPLexer correctly predicts biologically-relevant motions as illustrated by examples in Figure 3b-c, but a more systematic examination is currently hampered by the sensitivity of TM-Score and lDDT-BS to binding-irrelevant fluctuations. We note that native contact analysis algorithms [51] may provide improved metrics for interpreting protein generative models and consider that a future direction.

Discussion
We have presented a learning-based method for dynamic-backbone protein-ligand structure prediction, establishing an accuracy and sampling efficiency advantage relative to baseline approaches. We anticipate the incorporation of state-of-the-art protein representation learning techniques such as the use of sequence evolutionary signals, pretrained language models or higher-level attention mechanisms [9,24,25] and training on large-scale structure datasets to further improve the methodology and facilitate applications in various downstream molecular design problems.
The forward-time and reverse-time SDEs The forward-time SDEs in NeuralPLexer are summarized in Table A1. For generality, we introduce an effective time stamp τ such that the drift and diffusion coefficients are constant θ(t) = θ, σ(τ ) = σ. The symbolic conventions are as following: • x Cα ∈ R nres×3 denotes the collection of alpha-carbon coordinates in the protein, following the standard nomenclature for amino acid atom types: • x nonCα ∈ R (n−nres)×3 denotes the set of coordinates for all non-alpha-carbon protein atoms (backbone N, C, O, and all side-chain heavy atoms); • y ∈ R m×3 denotes all ligand heavy atom coordinates. Note that m . . = K k=1 m k with m k being the number of heavy atoms in each ligand molecule G k .

A.1.1 Transition kernel densities and sampling
Following the general result for Ornstein-Uhlenbeck processes [52] given the effective time-homogeneous diffusion process described in Table A1, for internal coordinates x nonCα − x Cα : since the Brownian motions w 1 , w 2 are independent, we obtain the transition kernel for the a finite time interval s: Similarly, for the ligand degrees of freedom the transition kernel is = N e −θs (y(τ ) − c T x Cα (τ )); (1 − e −2θs ) σ 2 2θ 2 (I + c T c) The transition kernel for alpha-carbon atoms is a standard Gaussian Defining σ 2 1 = σ 2 2θ , σ 2 2 = σ 2 · τ (T * ), andτ = 2θτ , we recover (2)(3)(4). For model training in practice, we use an exponential noise schedule defined by τ = τ 0 e t and τ 0 = σ 2 min σ 2 with σ min being a minimum perturbation scale as commonly adopted in variance-exploding (VE) [22] SDEs. For completeness, the SDEs defined in the transformed time horizon t ∈ [0, T * ] is given by replacing the drift coefficient θ and the diffusion coefficient σ with the following time-dependent counterparts: and To sample from the marginal distribution q t . . = p data * q 0:t derived from the forward SDEs: x where α(t) = e −2θτ (t) .
For the reverse-time SDE the ESDM φ predicts the denoised observationsx(0),ŷ(0) usingx(t),ŷ(t) which is formally equivalent to estimating the score function ∇ Z log q t (Z) [53]. Given a time discretization schedule with interval s, we obtain the expression for the predicted observation meanZ(φ, t − s) in one denoising step Z(t) → Z(t − s): Cα (0)) standard ODE-based or SDE-based integrators can then be adapted to sample from (19).

A.1.2 Euclidean equivariance
Given group G, a function f : X → Y is said to be equivariant if for all g ∈ G and x ∈ X, f (ϕ X (g)·x) = ϕ Y (g)·f (x). Specifically f is said to be invariant if f (ϕ X (g)·x) = f (x). We are interested in the special Euclidean group G = SE (3) consists of all global rigid translation and rotation operations g · Z . . = t + Z · R where t ∈ R 3 and R ∈ SO(3). To adhere to the physical constraint that p data is always SE(3)-invariant, the transition kernels of forward-time SDE should satisfy SE(3)-equivariance q(Z t+s |Z t ) = q(g · Z t+s |g · Z t ) such that the marginals are invariant q t (Z t ) = q t (g · Z t ) for any time t. The proofs are straightforward: For receptor Cα degrees of freedom For receptor non-Cα degrees of freedom

For ligand degrees of freedom
where we have used c T t = t up to a column-wise broadcasting operation based on the row-wise normalization property of the softmax-transformed contact map c.
Since all transition kernels are SE(3)-equivariant, it then follows that the score ∇ Z log q t (Z) is also SE(3)-equivariant: ∇ Z log q t (Z ) = ∇ Z log q t (Z) · R where Z = t + Z · R and thus the reverse-time SDE is equivariant. While the forward SDE is also E(3)-equivariant as the noising process satisfies q(−Z(τ + s)| − Z(τ )) = q(Z(τ + s)|Z(τ )), it is worth noting that the reverse SDE is only SE(3)-equivariant as parity-inversion transformations i : Z → −Z on the data distribution p data is physically forbidden and thus the score ∇ Z log q t (Z) is of broken chiral symmetry in general:

A.2 Small-molecule featurization and encoding
We consider two types of nodes to construct a graph-based molecular representation: (a) heavy-atoms i ∈ {1, 2, · · · , N atom } and (b) local coordinate frames u ∈ {1, 2, · · · , N frame }, u . . = u(ijk) formed by atom triplets (i, j, k) that are connected by bonds (ij) and (jk). We introduce Path-integral Graph Transformer (PiFormer), an attentional neural network with edge-level operations inspired by the path-integral formulation of quantum mechanics, to infer the long-range inter-atomic geometrical correlations for small molecules based on their graph-topological properties. PiFormer operates on the collection of following classes of embeddings: • Atom representations H ∈ R Natom × c. The input atom representations is a concatenation of one-hot encodings of element group index and period index for the given atom, which is embedded by a linear projection layer R 18+7 → R c ; • Frame representations F ∈ R N frame × c. For a given frame u, F u is initialized by a 2-layer MLP R 4 * 2+18+7 → R c that embed the bond type encodings (defined as [is_single, is_double, is_triple, is_aromatic]) of the "incoming" bond (i(u), j(u)), "outgoing" bond (j(u), k(u)), and the atom type encoding of the center atom j(u); • Stereochemistry encodings S ∈ R N frame ×N frame ×cs . S is a sparse tensor where an element S uv is nonzero only if the pair of frames (u, v) is adjacent, i.e., u and v sharing a common incoming or outgoing bond; • Pair representations G ∈ R N frame ×Natom×cp . G is initialized by an outer sum of H and F which is added to linear-projected S and passed to a 2-layer MLP.
Elements of the stereochemistry encoding tensor S are defined as S uv,9 . . = is_double_or_aromatic(common_bond(u, v)) ∨ is_same_side(u, v) S uv,10 . . = is_double_or_aromatic(common_bond(u, v)) ∨ not_same_side(u, v) is defined as one of the three atoms in frame v is above the plane formed by frame u with normal ) is defined as the two bonds not shared between u, v being on the same side of the common bond, equivalent to v u · v v > 0, vice versa. Our current technical implementations for is_above_plane and is_same_side are based on computing the normal vectors and dot-products using the coordinates from an auxiliary conformer, but we note that in principle all stereochemistry encodings can be generated based on cheminformatic rules without explicit coordinate generations. We additionally denote MASK s as a N frame × N frame logical matrix defined as the adjacency matrix of frame pairs (u, v).
The notion of "frames" in a coordinate-free topological molecular graph is justified by the inductive bias that most bending and stretching modes in molecular vibrations are of high frequency, i.e., most bond lengths and bond angles fall into a small range as predicted by valence bond theory, such that the local frames forms a consistent molecular representation without prior knowledge on 3D coordinates. PiFormer operates solely on the molecular representation defined by the input graph, and the frame coordinates (t, R) are initialized right before the ESDM blocks.
The forward pass of single PiFormer block is expressed as: where K denotes the propagation length truncation for the learnable graph kernel exp(U l ) ≈ (1 + 1 K U l ) K in a single PiFormer block, MLP denotes a 3-layer multilayer perceptron combined with layer normalization [54]. W K , W Q , W S , W G are trainable linear weight matrices. MHAwithEdgeBias(X 1 , X 2 , X edge ) denotes a multi-head cross-attention layer between source node embeddings X 1 and target node embeddings X 1 , with edge embeddings X edge entering attention computation as a relative positional encoding term as in the relation-aware transformer introduced in [7]. For all models descibed in this study, we set l max = 6 and K = 8.

A.2.1 PiFormer model pretraining
In Table A2 we summarize the small-molecule datasets used for training the PiFormer encoder used in the reported NeuralPLexer model. The loss function used in PiFormer pretraining is the following: We use a mixture density network head to encourage alignment between the learned last-layer pair representations G and the intra-molecular 3D coordinate marginals. For a single training sample with 3D coordinate observation y: where where rigidFrom3Points is the Gram-Schmidt-based frame construction operation described in Ref. [9], Alg. 21; we additionally add a numerical stability factor of 0.01 Å to the vector-norm calculations to handle edge cases when computing the rotation matrices from perturbed coordinates. Each component the 3D distance-angle distribution q 3D is parameterized by where PowerSpherical is a power spherical distribution introduced in [62]; m iul . . = (µ, σ, v) iul , and whre 3DMixtureDensityHead is a 3-layer MLP.
Using an equivariant graph transformer similar to ESDM (see Sec. A.6) but with all receptor nodes dropped, we construct a geometry prediction head to perform global molecular 3D structure denoising. We sample noised coordinates y(t) from a VPSDE [22] and introduce a SE(3)-invariant denoising score matching loss based on the Frame Aligned Point Error (FAPE) [9]: L CC−regression is a mean squared loss for fitting the "level 1" chemical checker (CC) [61] embeddings which represents harmonized and integrated bioactivity data, and L CC−ismask is an auxiliary binary cross entropy loss for classifying whether a specific CC entry is available for any molecule in the chemical checker dataset. Model is trained for 20 epochs with 15% masking ratio for atom and bond encodings, 40% masking ratio for stereochemistry encodings, and dropout=0.1; L MLM is a standard cross-entropy loss for predicting the masked tokens which is added to encourage learning on molecular graph topology distributions, but empirically we found L MLM converged within the first two epochs and did not find it to influence the learning dynamics of other tasks.

A.3 Protein sequence and backbone encoding
The inputs to the protein encoder are (i) the one-hot amino-acid type (20 standard residues + 1 "unknown" token) encoding of the 1D sequence s, (ii) the backbone (N, Cα, C) coordinates of a perturbed protein structure x(t) sampled from the forward SDEs described in Table A1, and (iii) a random Fourier encoding of the diffusion time step t. To reduce memory cost, the protein backbone is represented as a sparse graph with each node mapped to each amino acid residue and randomized edges according to the inclusion probability p(add_edge(i, j)) = exp(− x i (t) − x j (t) /10.0 Å) for all residue pairs (i, j). The edge representations are initialized as a random Fourier encoding of the signed sequence distance between two residues (i, j) if i and j are located on the same chain, and are initialized as zeros if (i, j) are located on different chains.
The protein encoder is composed of 4 stacks of invariant point attention (IPA) [9] blocks with two technical modifications: • The attention scores are computed on the sparsified protein graph, instead of the densely-connected graph as in standard self-attention layers; • Each node i is associated with n head replicas of coordinate frames {T } i , instead of a single frame as in a static structure representation. {T } i is initialized as n head copies of the backbone frames constructed by rigidFrom3Points(x N,i , x Cα,i , x C,i ). The layer output is n head × 7 scalars representing the translation vector and the quaternion variable to update the frame associated with each attention head.
the multi-replica design is found to moderately improve model convergence at a fixed network size. For conciseness, we refer to the modified invariant point attention as GraphIPA.

A.4 Contact predictor
As illustrated in Figure 1, the embeddings from the protein and small-molecule ligand graph encoders are passed to the contact predictor to estimated the contact maps L. A protein-ligand graph is created before the contact predictor forward pass, with pairwise intermolecular edges connecting all protein residues and ligand atoms. The contact predictor is composed of 4 modules each comprises of an intra-protein GraphIPA block, a bidirectional intra-ligand-graph self-attention layer, a bidirectional self-attention layer on the protein-ligand intermolecular edges, and a MLP to update protein-ligand edge representations using the attention maps and previous-layer edge representations. The final edge representations are used to predict L as described by Equation 5. The contact predictor weights are shared across all one-hot contact matrix sampling iterations.

A.5 All-atom graph featurization
All protein heavy-atoms nodes (features and 3D coordinates) and the ligand 3D coordinates sampled from the geometry prior q T * are added to the network inputs right before the ESDM block forward pass. Each protein atom representation is initialized as the concatenation of: • The residue-wise representation from the protein backbone encoder; • An one-hot encoding of its atom type as defined by the 37 standard amino acid heavy atom symbols in the PDB format [63]; • A random Fourier encoding of the diffusion time step t.
A random Fourier encoding of the diffusion time step t is also concatenated to the ligand atom representations from the ligand graph encoder and are transformed by a 2-layer MLP.
Given the noised all-atom protein coordinates at diffusion time t, the following edges are added to the protein-ligand graph: • Edges connecting a protein atom node and the residue node that the protein atom belongs to; • Edges connecting two protein atom nodes that are within the same residue; • Edges connecting two protein atom nodes that are within 6.0 Å distance; • Edges connecting a protein atom node and a ligand atom node that are within 8.0 Å distance; The protein-atom-involving edges are initialized as a concatenation of the following features: • A boolean code indicating whether the source node and target node belong to the same residue or the same ligand molecule; • A boolean code indicating whether there is a covalent bond between the source and target nodes. The covalent bonding information for protein-ligand edges are resolved based on the reference protein-ligand complex structure, where an atom pair (i, j) is considered as a covalent bond if the distance satisfies d ij < 1.2σ ij where σ ij = 1 2 (σ i + σ j ) is the average Van der Waals (VdW) radius for the atom pair.
To focus the learning problem on binding-site parts of the protein-ligand complex structure, the following native contact encoding features are added to the protein sub-graph that do not involve residues that are within 6.0 Å of any ligand heavy atom; given two amino acid residues, we define the native contact encoding as the concatenation of clean-protein-structure N − N, Cα − Cα, and C − C distances discretized into  Since the third row of R is a pseudovector as described in rigidFrom3Points, the determinant of the rotation matrix R is unchanged under parity inversion transformations i : x → −x and thus the intermediate quantity f vloc is SE(3)invariant but in general not invariant under parity inversion i. This property ensures that ESDM can learn the correct chiral symmetry breaking behaviors in molecular 3D conformation distributions.
The predicted drift vectors ∆t are added to the input node coordinates; the final coordinate outputs are taken as the predicted denoised observationsx(0),ŷ(0).

A.7 Model training and hyperparameters
The loss function for NeuralPLexer training is: We train the contact predictor ψ to match the posterior distribution defined by the observed contact map q L . . = Categorical nres×m (L) where L . . = k L k with intermediate ligand-wise one-hot matrices l k sampled from q L k : where KL denotes a Kullback-Leibler divergence and JS denotes a Jensen-Shannon divergence. An auxiliary loss is added to the mean term in the predicted geometry prior: The denoising score matching (DSM) loss expressions are given by L DSM−site is defined analogously but averaged for residues that are within 6.0 Å of the ligand in the ground-truth structure. Lastly For the ligand graph encoder, we use 6 PiFormer blocks with a embedding dimension of 512 for atom representation and frame representations, and a dimension of 128 for pair representations. For the protein encoder, we use 4 GraphIPA blocks with a node embedding dimension of 256 and edge embedding dimension of 64. For the contact predictor we use 4 blocks with the same embeddings sizes (256, 64) as in the protein encoder; linear layers are added to project the ligand representations to the length of protein representations before they are passed to the contact predictor. For ESDM, we use a stack of 4 blocks with a embedding dimension of 64 for both node and edge representations, that is, each node i is associated with scalar representations f s,i of size 64 and vector representations f v,i of size [3,64].
The pretrained small-molecule encoder weights are frozen during training. Model is trained with batch size of 8 for 40 epochs, using dropout=0.05, an inital learning rate of 3E-4 with 1000 warmup steps followed by a cosine annealing learning rate decay schedule. On the PDBBind 2020 training set (170k samples), the training run took 20 hours a single NVIDIA-Tesla-V100-SXM2-32GB GPU.

A.7.1 Task-specific fine-tuning
The model used for fixed-backbone protein-ligand docking is fine-tuned on the original PDBBind training dataset, while all backbone atoms (N, Cα, C, O) and Cβ atoms are set to the ground-truth coordinates. Fine-tuning is performed for 20 epochs with a batch size of 8 without teacher forcing for the geometry prior (i.e., sampling the one-hot matrix l from the observed contact map q L = Categorical nres×m (L), using the predicted contact map ψ(l, s,x, G) to parameterize the finite-time transition kernels q t (Z(t)|Z(0)) during model forward pass, and then backpropagating the model end-to-end) using a cosine annealing schedule with a initial learning rate of 1E − 4.
The model used for binding-site inpainting is fine-tuned on all split-chain samples from the original PDBBind training dataset. A protein-chain/ligand pair is included in the fine-tuning dataset if any heavy atom of the ligand is within 10 Å of any heavy atom of the protein chain. All receptor residues that are not within 6.0 Å of the ligand are set to the ground-truth coordinates with the residue-wise and protein-atom-wise time-step encoding set to zeros. Fine-tuning is performed for 40 epochs with a batch size of 10 without teacher forcing for the geometry prior using a cosine annealing schedule with a initial learning rate of 1E − 4.

A.8 Computational details
A.8.1 Test datasets and post-processing While the time-split-based PDBBind 2020 dataset has been used in previous works for studying model generalization to novel protein-ligand pairs, we noticed that the 363-sample test set curated by [43] contains samples with improperly removed alternative ligand conformation ground truths or deleted adjacent chains that strongly interact with the ligand molecule in the full structure (e.g., binding sites near protein-protein interfaces). To ensure a reasonable comparison to docking-based methods, for the test dataset used fixed-backbone ligand conformation prediction experiments we keep all protein chains that are within 10 Å of the ligand from the original PDB file instead of using the receptor PDB files curated by PDBBind; we further removed all covalent ligands and pipetide binders from the test set as such cases are usually tackled by specialized algorithms [64,65], resulting in 275 test samples in total to produce the results presented in Figure 2a-d.
The AlphaFold2 structures used in the ligand-coupled binding site repacking task are predicted using ColabFold [66] with default MSA, recycling, and AMBER relaxation settings, and without using templates in order to best reflect the prediction fidelity of AlphaFold2 on novel targets (since all PDBBind test set samples are deposited before year 2021). The input sequences for all protein chains are obtained from https://www.ebi.ac.uk/pdbe/api/pdb/ entry/molecules/ to avoid issues related to unresolved residues and to represent a realistic testing scenario where the protein backbone models are obtained from the full sequence.

A.8.2 Baseline method configurations
We run CB-Dock [44] with a heuristic low-sampling-intensity configuration (exhaustiveness=1, number of clustered binding sites to start local docking = 1) such that the execution time (43 seconds per ligand on average on single core of an Intel(R) Xeon(R) CPU E5-2698 v4 @ 2.20GHz CPU) is comparable to deep-learning-based methods that were proposed to perform docking at a low computing budget. The top-scored ligand conformations collected for each protein-ligand pair as ranked by Autodock Vina [67] are used to obtain the success rate results in Figure 2a.
EquiBind [43] are launched with the default configuration file, and for each protein-ligand pair 64 ligand conformations are generated using different random RDKit [68] input conformers. We note that the incorporation of side-chain flexibility as provided by AutoDock Vina and the systematic tuning of sampling intensity in docking-based methods may offer a more accurate comparision regarding the accuracy/computational time relationship among physics-based and learning-based methods.
RosettaLigand [23] runs are launched with a configuration modified from the standard protocol. We set the receptor Calpha constraint parameter to 100.0 to enable a fully flexible receptor; the ligand coordinates are initialized using the aligned-ground-truth conformation as obtained by TM-Align [49], with randomized torsion angles using the BCL [69] library as described in the standard protocol. We set the docking box width to 4.0 Å and remove the ligand center perturbation step to ensure the ligand search space during the low-resolution docking stage is constrained to the binding site location. While high-fideltiy physics-based methods such as IFD-MD [18] have been proposed for flexible-receptor ligand docking, such algorithms often incurs orders-of-magnitude higher computational cost thus are not included within the scope of this study.

A.9 Evaluation metric details
All protein structure alignments and TM-Score calculations are performed using TMAlign [49]. All reported TM-Scores are normalized by the chain length of the reference PDB structure. The per-residue all-atom lDDT score is computed using OpenStructure [70]; the lDDT-BS score is then computed by averaging the per-residue scores for ligand binding site residues with a cutoff of 4.0 Å as used in CAMEO [48]. The symmetry-corrected heavy-atom RMSD for ligand structure comparison is computed using the obrms function in OpenBabel [71]. A standard 6-12 Lennard-Jones energy functional form is used for computing the clash rate statistics; the L-J energy and VdW radius parameters are obtained from the UFF parameter file retrieved from https://github.com/kbsezginel/lammps-data-file/ blob/master/uff-parameters.csv.