Harmonic Torsional Diffusion for Protein-Ligand Flexible Docking

arXiv:2608.20366 · q-bio.BM, cs.LG · Submitted 2026-06-18 · Read on arXiv

Listen

Transcript

Introduction to the show: ident: Genomics Radio. Generated commentary on the latest computational biology and genomics papers.

Ines: We're looking at the paper: "Harmonic Torsional Diffusion for Protein-Ligand Flexible Docking".

Marcus: Here's the gist of the paper. The gist The introduction proposes Harmony, a harmonic torsional diffusion framework for flexible protein-ligand docking, which parameterizes ligand and side-chain torsional score fields as derivatives of learned harmonic potentials on the circle to explicitly make periodicity explicit and improve modeling of flexible binding sites<ref:2608.20366#pg2> How it works Harmony introduces a harmonic torsional diffusion framework where ligand and side-chain torsional score fields are parameterized as derivatives of learned harmonic potentials on the circle, whose noise-level dependence is supplied analytically by the heat semigroup of variance-exploding diffusion on the torus<ref:2608.20366#pg5> This construction makes periodicity explicit and gives the model a frequencyaware inductive bias over rotameric motion<ref:2608.20366#pg2> The framework is built upon a Variance-Exploding (VE) diffusion SDE, where the forward process perturbs data through pure stochastic diffusion without any deterministic drift term<ref:2608.20366#pg5> This VE-SDE acts as a heat flow on the torus T, which suggests that the appropriate basis for torsional variables is given by the Fourier modes<ref:2608.20366#pg5> Harmonic Parametrization The model parametrizes torsional angles of ligands and protein side chains in a Fourier basis, which comes with a clean analytical score function formulation and efficient noising structure which gradually suppresses frequency motion components coming from high to low harmonic variations<ref:2608.20366#pg2> For each rotatable bond, the messagepassing backbone produces a local bond representation ue, from which Harmony predicts harmonic coefficients via

ae,1, be,1,..., ae,K, be,K: = gθ(ue), where K is the number of retained Fourier modes<ref:2608.20366#pg5> These coefficients define a latent harmonic function on the circle ψe(ϑ) = X K k=1

ae,k cos(kϑ) + be,k sin(kϑ): <ref:2608.20366#pg5> By Propositions 4.1-4.3, the VE forward process on each torsional coordinate acts as heat flow on T, and Proposition 4.2 diagonalizes this action in the Fourier basis<ref:2608.20366#pg2> The smoothed harmonic function at noise level λ is therefore ψe,λ(ϑ) = X K k=1 exp-1/2 k 2σ 2(λ)

ae,k cos(kϑ) + be,k sin(kϑ): <ref:2608.20366#pg5> Harmony predicts the torsional score field analytically by differentiating the heat-smoothed harmonic expansion with respect to the angle: se(ϑ, λ) = ∂ϑψe,λ(ϑ) = X K k=1 k exp-1/2 k 2σ 2(λ)

-ae,k sin(kϑ) + be,k cos(kϑ): <ref:2608.20366#pg5> Model Architecture and Training Harmony uses an e3nn (Geiger & Smidt, 2022) heterogeneous SE(3)-equivariant graph neural network over three coupled graphs: the ligand atom graph, the receptor Cα graph, and the receptor atom graph in the binding pocket<ref:2608.20366#pg5> The model is trained with the standard diffusion score-matching objective as in (Plainer et al., 2023)<ref:2608.20366#pg2> Training uses OpenEquivariance (Bharadwaj et al., 2025) kernels for Clebsch-Gordan tensor products, reducing runtime from roughly 4.5 days to 2 days on 8 H100 80GB GPUs<ref:2608.20366#pg5> The model utilizes a lightweight confidence head attached to the shared docking encoder which predicts a scalar confidence score for each sampled pose, supervised with proteinligand interface lDDT (PLI-lDDT) (Lu et al., 2024)<ref:2608.20366#pg2> Experimental Results and Ablations Harmony improves over recent generative methods for flexible docking on PDBBind benchmark by ∼24.5% points while following their training setup and model architecture<ref:2608.20366#pg2> On PoseBusters set (Buttenschoen et al., 2024), Harmony shows improvements in chemical validity of generated poses over the large portion of baselines<ref:2608.20366#pg2> Ablation Study on Harmonic Modes indicates that Harmony performs the best with K = 4, suggesting intrinsic low-frequency nature<ref:2608.20366#pg2> Ablation Study on Frequency Damping concludes that in such setting performance degrades and noising of high frequency components is the central part of the proposed model<ref:2608.20366#pg2> Case studies on EBNA1 and KRAS G12D illustrate the method’s behavior on a polar and a shallow binding site, respectively<ref:2608.20366#pg2> We introduced Harmony, a harmonic torsional diffusion module for flexible protein-ligand docking<ref:2608.20366#pg2> By parameterizing ligand and side-chain torsional scores in Fourier space and utilizing variance-exploding diffusion properties on the torus, Harmony makes angular periodicity explicit and improves the modeling of flexible binding sites<ref:2608.20366#pg2> Across standard benchmarks and case studies, this simple geometric change consistently improves pose recovery and all-atom reconstruction<ref:2608.20366#pg2> More broadly, our results suggest that respecting the intrinsic manifold structure of torsional motion is not a minor modeling detail, but a key ingredient for accurate diffusion-based flexible docking<ref:2608.20366#pg2> How it works The pocket-centered docking approach involves identifying nearby residues by computing the minimum distance from each receptor residue to the ligand and selecting residues within a pocket cutoff<ref:2608.20366#pg5> The docking state is parameterized as z = (t, R, τ, χ) ∈ M = R cubed × SO(3) × T ml × T ms, where t ∈ R cubed and R ∈ SO(3) denote ligand translation and rotation, τ ∈ T ml denotes ligand torsions, and χ ∈ T ms denotes flexible pocket side-chain torsions<ref:2608.20366#pg5> The network jointly predicts ligand pose and pocket side-chain arrangement from Gpocket<ref:2608.20366#pg5> Data Preparation For ligands, we generate a fresh conformer using RDKit with a force field refinement, identify the rotatable bonds, and then optimize the ligand torsion angles by differential evolution so that the generated conformer best matches the reference ligand pose<ref:2608.20366#pg5> For receptor side chains, after apo-holo alignment, we optimize each flexible residue independently over its torsional degrees of freedom using differential evolution, with the objective of minimizing the RMSD between the moved side-chain atoms and the corresponding holo coordinates<ref:2608.20366#pg5> Inference Algorithm 2 involves a reverse VE update where t n-1 = t n + gtr(λn) 2∆λˆstr + gtr(λn)√∆λ ξtr, Rn-1 = Exp grot(λn) 2∆λˆsrot + grot(λn)√∆λ ξrot, τ n-1 = wrap τ n + gtor(λn) 2∆ λˆstor + gtor(λn)√Δ λ ξtor, χ n-1 = wrap χ n + gsc(λn) 2∆ λˆssc + gsc(λn)√Δ λ ξsc<ref:2608.20366#pg5> This process iteratively denoises the components to recover the bound ligand pose together with the binding-pocket side-chain arrangement according to backward process (Anderson, 1982)<ref:2608.20366#pg5> The alignment of apo side chains into local backbone frames is performed by computing the rotation that maps the apo frame to the holo frame and applying that rotation to the corresponding local side-chain coordinates<ref:2608.20366#pg5> Evaluation Metrics Our primary ligand metric is the ligand root-mean-square deviation (L-RMSD) between the predicted ligand pose and the reference crystal pose after optimal rigid alignment of the pockets<ref:2608.20366#pg2> To assess whether the full predicted complex is geometrically accurate, including the flexible binding-site side chains, we also measure all-atom RMSD (AA-RMSD) over the receptor atoms retained in the pocket representation<ref:2608.

Paper summary: Ines: So, we're looking at this paper, "Harmonic Torsional Diffusion for Protein-Ligand Flexible Docking." Essentially, the problem they're tackling is that most current diffusion docking models don't handle protein flexibility and ligand twisting well because they use these generic Euclidean heads that just ignore the periodic nature of angles.

Marcus: That mismatch is particularly bad when you're doing flexible docking, where both the ligand and those pocket side chains have to adapt together to form the final complex. This paper introduces something called Harmony, which is a harmonic torsional diffusion framework for this kind of problem.

Yuki: What does this mean for us biologically? It suggests that just treating torsions as simple continuous variables isn't enough; you need a way to make the model aware of how those angles repeat on a circle.

Ines: Exactly, and Harmony does that by parameterizing the ligand and side-chain torsional score fields as derivatives of learned harmonic potentials on the circle. It makes periodicity explicit in the model.

Marcus: And they tie this noise level dependence to the heat semigroup of variance-exploding diffusion on the torus, which is a very specific mathematical structure for how diffusion behaves there. They say this construction gives the model a frequency-aware inductive bias over rotameric motion.

Yuki: So, instead of just predicting one update for a torsion angle, it learns harmonic coefficients that define a Fourier potential. That sounds like it’s giving the model a built-in understanding of the different vibrational modes in those torsional motions.

Ines: Right, and they use this to predict the torsional score field by differentiating that heat-smoothed harmonic expansion with respect to the angle, which gives them an analytical score function formulation. They are basically creating a clean way to describe how those torsions change based on the diffusion process.

Marcus: The architecture itself uses an e3nn heterogeneous SE(three)-equivariant graph neural network over three coupled graphs: one for the ligand atoms, one for the receptor alpha carbons, and one for all the receptor atoms in the binding pocket. This lets it jointly predict both poses and those side-chain arrangements from a single pocket representation.

Yuki: That joint prediction capability is interesting because it means you're not just getting a good ligand pose; you're also getting a better idea of how those flexible side chains fit around the protein, which is crucial for physical validity.

Ines: They tested this on the PDBBind benchmark and found that Harmony improves ligand pose accuracy and pocket all-atom reconstruction compared to recent flexible docking methods. It also shows improvements on PoseBusters regarding the physical validity of generated complexes.

Marcus: I saw their ablation study, and it indicated that Harmony performed best when they kept K equal to four harmonic modes, suggesting an intrinsic low-frequency nature for the motion they are modeling <ref:2608.20366#pg2>. But they also noted that performance degrades if you try to model high frequency components because the noising of those high frequencies is central to the proposed method.

Yuki: That suggests that focusing on a lower frequency representation captures the most important physical interactions in these binding sites, which makes sense when thinking about how flexible molecules interact with pockets.

Ines: So, for someone just listening to this, what this means is that we need to move past treating torsions as simple continuous movements and start respecting the underlying geometry of their motion—the periodicity on a circle.

Marcus: It’s a simple change in how you parameterize the score function that seems to yield tangible improvements on benchmarks like PDBBind. But they did admit a limitation, which is that they aren't explicitly solving for full backbone motion yet; they are focusing specifically on ligand and flexible side-chain torsions.

Yuki: That’s fair; tackling full protein backbone motion is a much harder problem, but solving the torsional part with this harmonic diffusion framework seems like a very effective step toward modeling complex binding dynamics accurately.

Ines: The authors of "Harmonic Torsional Diffusion for Protein-Ligand Flexible Docking" have put forward this specific harmonic torsional diffusion framework to address the limitations of existing models in flexible docking.

Marcus: It really boils down to aligning the score parameterization with the geometry of the diffusion process, which they show is a simple and effective lever for improving flexible docking accuracy across different benchmarks.

Yuki: For us, this suggests that respecting the intrinsic manifold structure of torsional motion is not just a minor modeling detail; it’s actually a key ingredient for getting accurate results in diffusion-based flexible docking.

Conclusion: Ines: So, we're wrapping up this look at Harmony, that harmonic torsional diffusion framework for flexible protein-ligand docking.

Marcus: Yeah, yeah, so it’s about taking these complex rotational motions and making them behave more like they should on a circle because the diffusion process itself is happening on a torus.

Yuki: From my angle, it means we're not just guessing where the ligand sits; we're modeling how those side chains physically twist in response to the binding pocket geometry.

Ines: Exactly, and by using these learned harmonic potentials, they’re giving the model an explicit way to understand rotation that goes beyond just standard Euclidean heads.

Marcus: It’s a simple geometric shift in how you parameterize the score function, but it seems to actually improve pose recovery on benchmarks like PDBBind by about twenty-five percent points.

Yuki: That improvement matters because for population genetics, if we can model those binding sites more accurately, we might get better ideas about how evolutionary pressures shape protein interactions across different species.

Ines: Right, and the authors themselves are showing that this structure respects the intrinsic manifold of torsional motion, which is a pretty fundamental thing for any simulation to get right.

Marcus: They also flagged that while it does a great job with ligand torsions and pocket side chains, it doesn't fully tackle the whole protein backbone flexibility yet.

Yuki: So what they’re showing us now is that focusing on those specific torsional degrees of freedom is a really effective first step before trying to model the entire protein structure at once.

Ines: It’s a solid piece of work for making flexible docking more physically sound, and I think it points toward a way to build better diffusion models.

Maksim Zhdanov, Pavel Strashnov, Vladislav Kurenkov

q-bio.BM, cs.LG

Submitted: 2026-06-18

Updated: 2026-10-18

Comments: Accepted at NeurIPS 2026

Code: https://github.com/dunnolab/harmony

License: http://creativecommons.org/licenses/by/4.0/

The gist: The gist The introduction proposes Harmony, a harmonic torsional diffusion framework for flexible protein-ligand docking, which parameterizes ligand and side-chain torsional score fields as

Key concepts

Harmonic Torsional Diffusion: This is a framework that uses diffusion processes on a circle (the torus) to model torsional angles. It parameterizes the score fields for ligand and side-chain torsions as derivatives of learned harmonic potentials, which explicitly incorporates the periodic nature of rotational motion in molecular structures.
Variance-Exploding (VE) Diffusion: This is a specific type of stochastic differential equation used in Harmony. It perturbs data through pure noise without any deterministic drift term. This process acts as a heat flow on the torus, which naturally suggests that the Fourier basis is the appropriate way to represent these torsional variables.
Fourier Basis Parametrization: The model parameterizes torsional angles using a Fourier basis. This allows for a clean analytical score function formulation and an efficient noising structure. It gradually suppresses high-frequency motion components, which is crucial for accurately capturing the low-frequency, physically relevant rotameric motions.
Frequency-Aware Inductive Bias: By aligning the score parameterization with the geometry of the VE diffusion process on the torus, Harmony gains a frequency-aware inductive bias. This means the model inherently understands and prioritizes different frequencies of motion in torsional variables, leading to more physically valid predictions for flexible binding sites.

Summary

The gist The introduction proposes Harmony, a harmonic torsional diffusion framework for flexible protein-ligand docking, which parameterizes ligand and side-chain torsional score fields as derivatives of learned harmonic potentials on the circle to explicitly make periodicity explicit and improve modeling of flexible binding sites<ref:2608.20366#pg2>

How it works

Harmony introduces a harmonic torsional diffusion framework where ligand and side-chain torsional score fields are parameterized as derivatives of learned harmonic potentials on the circle, whose noise-level dependence is supplied analytically by the heat semigroup of variance-exploding diffusion on the torus<ref:2608.20366#pg5> This construction makes periodicity explicit and gives the model a frequencyaware inductive bias over rotameric motion<ref:2608.20366#pg2> The framework is built upon a Variance-Exploding (VE) diffusion SDE, where the forward process perturbs data through pure stochastic diffusion without any deterministic drift term<ref:2608.20366#pg5> This VE-SDE acts as a heat flow on the torus T, which suggests that the appropriate basis for torsional variables is given by the Fourier modes<ref:2608.20366#pg5>

Harmonic Parametrization

The model parametrizes torsional angles of ligands and protein side chains in a Fourier basis, which comes with a clean analytical score function formulation and efficient noising structure which gradually suppresses frequency motion components coming from high to low harmonic variations<ref:2608.20366#pg2> For each rotatable bond, the messagepassing backbone produces a local bond representation ue, from which Harmony predicts harmonic coefficients via [ae,1, be,1,..., ae,K, be,K] = gθ(ue), where K is the number of retained Fourier modes<ref:2608.20366#pg5> These coefficients define a latent harmonic function on the circle ψe(ϑ) = X K k=1 [ae,k cos(kϑ) + be,k sin(kϑ)]<ref:2608.20366#pg5> By Propositions 4.1-4.3, the VE forward process on each torsional coordinate acts as heat flow on T, and Proposition 4.2 diagonalizes this action in the Fourier basis<ref:2608.20366#pg2> The smoothed harmonic function at noise level λ is therefore ψe,λ(ϑ) = X K k=1 exp-1/2 k 2σ 2(λ) [ae,k cos(kϑ) + be,k sin(kϑ)]<ref:2608.20366#pg5> Harmony predicts the torsional score field analytically by differentiating the heat-smoothed harmonic expansion with respect to the angle: se(ϑ, λ) = ∂ϑψe,λ(ϑ) = X K k=1 k exp-1/2 k 2σ 2(λ) [-ae,k sin(kϑ) + be,k cos(kϑ)]<ref:2608.20366#pg5>

Model Architecture and Training

Harmony uses an e3nn (Geiger & Smidt, 2022) heterogeneous SE(3)-equivariant graph neural network over three coupled graphs: the ligand atom graph, the receptor Cα graph, and the receptor atom graph in the binding pocket<ref:2608.20366#pg5> The model is trained with the standard diffusion score-matching objective as in (Plainer et al., 2023)<ref:2608.20366#pg2> Training uses OpenEquivariance (Bharadwaj et al., 2025) kernels for Clebsch-Gordan tensor products, reducing runtime from roughly 4.5 days to 2 days on 8 H100 80GB GPUs<ref:2608.20366#pg5> The model utilizes a lightweight confidence head attached to the shared docking encoder which predicts a scalar confidence score for each sampled pose, supervised with proteinligand interface lDDT (PLI-lDDT) (Lu et al., 2024)<ref:2608.20366#pg2>

Experimental Results and Ablations

Harmony improves over recent generative methods for flexible docking on PDBBind benchmark by ∼24.5% points while following their training setup and model architecture<ref:2608.20366#pg2> On PoseBusters set (Buttenschoen et al., 2024), Harmony shows improvements in chemical validity of generated poses over the large portion of baselines<ref:2608.20366#pg2> Ablation Study on Harmonic Modes indicates that Harmony performs the best with K = 4, suggesting intrinsic low-frequency nature<ref:2608.20366#pg2> Ablation Study on Frequency Damping concludes that in such setting performance degrades and noising of high frequency components is the central part of the proposed model<ref:2608.20366#pg2> Case studies on EBNA1 and KRAS G12D illustrate the method’s behavior on a polar and a shallow binding site, respectively<ref:2608.20366#pg2>

Conclusion

We introduced Harmony, a harmonic torsional diffusion module for flexible protein-ligand docking<ref:2608.20366#pg2> By parameterizing ligand and side-chain torsional scores in Fourier space and utilizing variance-exploding diffusion properties on the torus, Harmony makes angular periodicity explicit and improves the modeling of flexible binding sites<ref:2608.20366#pg2> Across standard benchmarks and case studies, this simple geometric change consistently improves pose recovery and all-atom reconstruction<ref:2608.20366#pg2> More broadly, our results suggest that respecting the intrinsic manifold structure of torsional motion is not a minor modeling detail, but a key ingredient for accurate diffusion-based flexible docking<ref:2608.20366#pg2>

How it works

The pocket-centered docking approach involves identifying nearby residues by computing the minimum distance from each receptor residue to the ligand and selecting residues within a pocket cutoff<ref:2608.20366#pg5> The docking state is parameterized as z = (t, R, τ, χ) ∈ M = R cubed × SO(3) × T ml × T ms, where t ∈ R cubed and R ∈ SO(3) denote ligand translation and rotation, τ ∈ T ml denotes ligand torsions, and χ ∈ T ms denotes flexible pocket side-chain torsions<ref:2608.20366#pg5> The network jointly predicts ligand pose and pocket side-chain arrangement from Gpocket<ref:2608.20366#pg5>

Data Preparation

For ligands, we generate a fresh conformer using RDKit with a force field refinement, identify the rotatable bonds, and then optimize the ligand torsion angles by differential evolution so that the generated conformer best matches the reference ligand pose<ref:2608.20366#pg5> For receptor side chains, after apo-holo alignment, we optimize each flexible residue independently over its torsional degrees of freedom using differential evolution, with the objective of minimizing the RMSD between the moved side-chain atoms and the corresponding holo coordinates<ref:2608.20366#pg5>

Inference

Algorithm 2 involves a reverse VE update where t n-1 = t n + gtr(λn) 2∆λˆstr + gtr(λn)√∆λ ξtr, Rn-1 = Exp grot(λn) 2∆λˆsrot + grot(λn)√∆λ ξrot, τ n-1 = wrap τ n + gtor(λn) 2∆ λˆstor + gtor(λn)√Δ λ ξtor, χ n-1 = wrap χ n + gsc(λn) 2∆ λˆssc + gsc(λn)√Δ λ ξsc<ref:2608.20366#pg5> This process iteratively denoises the components to recover the bound ligand pose together with the binding-pocket side-chain arrangement according to backward process (Anderson, 1982)<ref:2608.20366#pg5> The alignment of apo side chains into local backbone frames is performed by computing the rotation that maps the apo frame to the holo frame and applying that rotation to the corresponding local side-chain coordinates<ref:2608.20366#pg5>

Evaluation Metrics

Our primary ligand metric is the ligand root-mean-square deviation (L-RMSD) between the predicted ligand pose and the reference crystal pose after optimal rigid alignment of the pockets<ref:2608.20366#pg2> To assess whether the full predicted complex is geometrically accurate, including the flexible binding-site side chains, we also measure all-atom RMSD (AA-RMSD) over the receptor atoms retained in the pocket representation<ref:2608.

Improvements for AI systems

Improved modeling of torsional motion by parameterizing ligand and side-chain torsional score fields as derivatives of learned harmonic potentials on a circle, which makes periodicity explicit and gives the model a frequencyaware inductive bias over rotameric motion. This allows the AI system to naturally represent difficult torsional configurations without learning an arbitrary score field from scratch.
Enhanced accuracy in flexible docking by improving ligand pose accuracy and pocket all-atom reconstruction, as demonstrated by Harmony improving ligand pose accuracy and pocket all-atom reconstruction over recent flexible docking methods on PDBBind benchmark. This leads to better physical validity of generated complexes, as shown by improvements on PoseBusters set.
Increased robustness in handling complex structural variability by explicitly modeling the product space diffusion process involving translation, rotation, ligand torsions, and pocket side-chain torsions independently at noise level lambda. This allows the AI system to jointly predict ligand pose and pocket side-chain arrangement from the graph representation while respecting the geometry of the diffusion process.
Optimized confidence estimation by using a time-weighted regression loss that upweights low-noise samples, emphasizing near-final geometries and adding only minimal computational overhead to predict a scalar confidence score for each sampled pose. This ensures that the system focuses on high-quality predicted structures during sampling.
Reduced training time for equivariant backbones by utilizing OpenEquivariance CUDA kernels specialized to the irreducible representation layout of the model, which substantially improves throughput from roughly 4.5 days to 2 days on 8 H100 GPUs while maintaining the architecture and training objective. This allows for faster iteration and training of large models.

Abstract

Molecular docking requires reasoning jointly about ligand pose and protein flexibility. Most diffusion-based docking models predict torsional updates with generic Euclidean heads that ignore the periodic geometry of angular variables. This mismatch is especially limiting in flexible docking, where ligand conformations and pocket side chains co-adapt to form the bound complex. Here, we introduce Harmony, a harmonic torsional diffusion framework for flexible protein-ligand docking. Harmony parameterizes ligand and side-chain torsional score fields as derivatives of learned harmonic potentials on the circle, which are noised mode by mode in the Fourier basis owing to the properties of variance-exploding diffusion on the torus. This construction makes periodicity explicit and gives the model a frequency-aware inductive bias over rotameric motion. On the PDBBind benchmark, Harmony improves ligand pose accuracy and pocket all-atom reconstruction over recent flexible docking methods. On PoseBusters, it improves the physical validity of generated complexes. Case studies on EBNA1 and KRAS G12D illustrate the method's behavior on a polar and a shallow binding site, respectively. Together, these results indicate that aligning the score parameterization with the geometry of the diffusion process is a simple and effective lever for improving flexible docking.