| Introduction

The effects of ligands binding to proteins are numerous and foundational to research efforts in biotechnology and drug discovery, yet efficiently determining the structure, and thereby the function, of such ligand-bound protein complexes has challenged the structural biology community for decades. In the 15th Critical Assessment of Techniques for Structure Prediction (CASP15) [1], template-based approaches to protein-ligand structure determination generally outperformed those based on deep learning (DL). Nonetheless, over the last few years, several new DL methods (mostly diffusion models [2]) for protein-ligand docking and structure prediction have been introduced [3][4][5][6][7], importantly raising the question of which method(s) perform(s) best for a range of diverse prediction targets and how to combine them if they are complementary. We sought to answer this question by designing and evaluating a new DL ensembling method called MULTICOM_ligand for protein-ligand modeling, which we will describe in detail in Section 2.

Our MULTICOM_ligand team (group number 207) participating in the 2024 CASP16 experiment submitted 33 models for 13 incidental ligand pose targets; 1165 models for 233 ligand pose pharma targets; 700 models for 140 ligand affinity pharma targets; and 110 affinity predictions for 110 phase-2 ligand affinity pharma targets, representing a submission for every protein target available in the CASP16 ligand prediction category. To facilitate such a breadth of prediction types, we designed MULTICOM_ligand as a modular software framework for protein-ligand modeling. Originally developed as a DL benchmarking toolkit for protein-ligand docking methods (i.e., PoseBench [8]), we adapted the core predictor modules of this benchmarking pipeline to support the prediction of arbitrary protein-ligand structures with associated confidence and affinity scores from protein sequence and ligand SMILES string inputs.

Moreover, to provide high-accuracy estimates of a protein-ligand complex's binding affinity from only primary sequences, we concurrently developed the new FlowDock generative flow matching model for joint prediction of protein-ligand structure and binding affinity [9]. Notably, our initial development of FlowDock revealed that joint training and prediction of protein-ligand structure and binding affinity yielded top results in various internal affinity prediction benchmarks we used for model prototyping and evaluation. As such, in the CASP16 experiment, we integrated FlowDock into MULTICOM_ligand as an optional add-on of the framework that, as one desires, can use initially provided protein-ligand complex structures as additional model inputs for confidence and affinity estimation. This flexible software design greatly simplified our usage of MULTICOM_ligand for Stage 2 of the CASP16 binding affinity prediction category, in which predictors were given the crystal structure of a protein-ligand complex and asked to estimate the complex's binding affinity using this additional information.

According to the CASP16 experiment's official analysis, in the protein-ligand structure prediction category, MULTICOM_ligand ranked fifth with its predictions' median lDDT-PLI score of 0.58, which denotes a protein-ligand interaction (PLI)focused implementation of the local Distance Difference Test (lDDT) for assessment of biomolecular structure accuracy. Further, in the protein-ligand binding affinity prediction category, MULTICOM_ligand achieved a Kendall's Tau ranking coefficient of 0.32 in Affinity Stage 1, earning it fifth place overall. Notably, MULTICOM_ligand performed better than many CASP16 template-based predictors, demonstrating that deep learning has advanced the state of the art of protein-ligand structure and binding affinity prediction since CASP15.

| Materials and Methods

| Overview of Approach

From primary sequence inputs of a protein and one or more ligands alone, MULTICOM_ligand, visualized in Figure 1, provides users with rank-ordered predicted protein-ligand complex conformations filtered using structural and chemical sanity checks available in the PoseBusters software suite [10] and annotated with estimated per-atom quality scores and binding affinity values produced by our new generative flow matching model FlowDock. This approach is generally summarized in Algorithm 1. The steps of the approach are described in detail in the following subsections.

| Protein-Ligand Inputs

MULTICOM_ligand represents a protein-ligand complex as a pair of single-/multi-chain protein sequence and SMILES string of one or more ligands (S, M). Multiple chains within a protein sequence are delimited using the character ":", whereas multi-ligand SMILES sequences within the same string are separated using the character "." following RDKit's conventions for parsing "fragment" ligands of a single molecule [11]. Certain protein-ligand structure prediction methods employed in MULTICOM_ligand support using predicted protein structures as input to enhance their prediction accuracy. Accordingly, we use ESMFold [12] to provide predicted protein structure inputs to these methods. Note that, during the CASP16 experiment, we instead predicted these protein structures using AlphaFold 3 [13], though the public release of MULTICOM_ligand's source code by default uses the MIT-licensed ESMFold model for these purposes.

| Structure Prediction Methods

Based on the results of our previous benchmark of DL-based protein-ligand docking methods [8], MULTICOM_ligand employed four representative DL methods to predict the structure for a protein-ligand sequence input: DiffDock-L [5], DynamicBind [4], RoseTTAFold-All-Atom [6], and NeuralPLexer [7]. We then grouped these methods into one of two groups, DL docking methods (i.e., DiffDock-L and DynamicBind) and DL co-folding methods (i.e., RoseTTAFold-All-Atom and NeuralPLexer), where the former group uses a predicted protein structure to perform DL-based molecular docking and the latter group predicts full protein-ligand complex conformations from primary sequence inputs. Last, we further subdivided these DL docking and DL cofolding groups into protein-fixed/protein-flexible categories (i.e., DiffDock-L/DynamicBind) and predictive/generative categories (i.e., RoseTTAFold-All-Atom/NeuralPLexer), respectively.

| Ranking Heuristics

One of the primary hypotheses driving this work is that geometrically similar ligand poses predicted by different DL methods should largely coincide with an accurate protein-ligand binding ALGORITHM 1 | MULTICOM_ligand for protein-ligand structure and affinity prediction.

Notation: (X: intermediate protein or protein-ligand structure; X: final protein-ligand structure; B: binding affinity, Ĉ : confidence score)

1: Input: Protein sequence and ligand SMILES string (S, M) 2: Predict X init ← ESMFold(S) 3: Sample X dd ← DiffDock -L S, M, X init 4: Sample X db ← DynamicBind S, M, X init 5: Sample X np ← NeuralPLexer S, M, X init 6: Predict X rfaa ← RoseTTAFold -All -Atom(S, M) 7: Rank X con ← StructureConsensus X dd,db,np,rfaa 8: Bust X bust ← PoseBustersFilters(X con ) 9: if Is Multi-Ligand then 10: Clash Bust X bust ← ClashFilters X bust 11: end if 12: Finalize X , Ĉ, B ← FlowDockAssess S, M, X bust 13: Output: Sampled top-5 heavy-atom structures X with confidence scores Ĉ and binding affinities B pocket and pose prediction overall. That is, when all DL methods have predicted nearly the same binding pocket and ligand pose for a given ligand molecule, they have, in essence, reached a "structural consensus" on the location and orientation of the crystal ligand pose. Based on this consensus (n.b., which may be misled if the majority of methods predict a similar incorrect binding pocket), we formulate an unsupervised ranking metric that calculates the pairwise root mean square deviation (RMSD) of all ligand poses predicted by each DL method and rank-orders the poses according to their average pairwise RMSD to each other. This provides a simple, computationally efficient heuristic (similar to that of [14] for protein complex structure ranking) for selecting our "best guess" of the location and orientation of a ligand pose given a pool of predictions produced by various DL prediction methods.

| Ligand Pose Filters

An important component of MULTICOM_ligand's design is that it not only curates a list of rank-ordered protein-ligand complex conformations produced by some of the latest DL prediction methods but also re-ranks (i.e., down-weights) its top-5 predicted conformations if any prediction fails to pass each of the standardized structural and chemical validity tests available in the PoseBusters software suite [10]. This provides an additional layer of filtering to ensure that MULTICOM_ligand's top predictions are ordered according to a secondary heuristic that posits that accurate ligand poses must not only be identified through a consensus of different prediction methods but must also not contain any violations of known ligand biochemistry such as non-planar ring conformations or steric clashes with protein heavy atoms.

| Selected Poses

During MULTICOM_ligand's initial stage of development, we discovered the need to add another layer of ranking: the possibility of encountering multi-ligand prediction targets for which accurate poses can be identified but may contain undesirable (and unrealistic) inter-ligand steric clashes between ligand heavy atoms. Notably, this phenomenon frequently occurs with DL methods such as DiffDock-L and DynamicBind which were originally trained on only single-ligand protein complexes, necessitating a stopgap measure to prevent such (clashing) poses from being selected as MULTICOM_ligand's top-ranked pose.

Consequently, for multi-ligand prediction targets, MULTICOM_ ligand automatically assigns predictions made by the DL method NeuralPLexer (n.b., which was trained on multi-ligand protein complexes with inter-ligand steric clash penalties) a higher rank than any other method's predictions, to discourage (potentially) clashing poses produced by the other (single-ligand) DL methods from being selected as MULTICOM_ligand's top pose prediction.

| Confidence and Affinity Prediction

A final component of MULTICOM_ligand's design is its ability to annotate its top-5 predicted protein-ligand structure conformations with estimated per-atom confidence scores and perligand binding affinity values. This is made possible by our new FlowDock generative model, a version of NeuralPLexer finetuned with geometric flow matching for joint protein-ligand structure and binding affinity prediction. Notably, the original NeuralPLexer model was trained as a denoising diffusion probabilistic model [2,[15][16][17] that predicts protein(-multi)-ligand complex structures and their confidence scores from primary sequence inputs, whereas FlowDock generalizes NeuralPLexer's diffusion generation framework with the emerging generative modeling framework of conditional flow matching [18][19][20] to enable generative (multi-ligand) structure predictions starting from biophysics-informed and empirical prior distributions [21][22][23].

At a high level, flow matching (n.b., as a generalization of denoising diffusion) has a DL model learn to solve an ordinary differential equation (ODE) that transforms data points derived from an easy-to-sample prior distribution X 0 (e.g., a Gaussian distribution) to another empirical distribution X 1 (e.g., the distribution of crystal structures in the RCSB Protein Data Bank (PDB) [24]). A DL model learns a solution to such an ODE by repeatedly "denoising" an interpolative noising schedule whereby, for a random time step t ∈ [0, 1] sampled during training, an input data point x 1 ∈ X 1 (e.g., a 3D biomolecular crystal structure) is "noised" according to time step t typically using simple linear interpolation such as x t = (1 -t) ⋅ x 0 + t ⋅ x 1 , and the model is then tasked with predicting the original version of this data point x 1 . Once trained, for example, in the context of structure prediction, such a DL model can be run iteratively to sample multiple 3D biomolecular structures for a primary input sequence starting from time step t = 0 representing a fully random point cloud x 0 ∈ X 0 . Note that additional (e.g., Gaussian) noise is typically injected into these training and sampling processes to ensure the model produces more than a trivial mapping between point masses [21].

Importantly, the primary novelty of flow matching is that one's prior distribution can be arbitrarily chosen, in contrast to denoising diffusion, for which typically only a Gaussian prior distribution can be used. This makes the modeling of 3D biomolecules, in particular, much more flexible in that, with flow matching, one can specify a prior distribution informed by known biophysical properties such as a harmonic prior [21] or one derived from the outputs of another DL structure prediction model such as ESMFold [12]. As such, for a given protein sequence and ligand SMILES string, FlowDock takes precise advantage of this modeling flexibility by sampling an initial protein structure using ESMFold and an initial molecule-like ligand conformation from a harmonic prior distribution at the start of its structure prediction sampling processes, which considerably reduces its training and prediction dynamics for arbitrary protein-ligand complexes.

In addition to introducing bespoke prior distributions for structure sampling, FlowDock repurposes NeuralPLexer's frozen (i.e., non-trainable) confidence estimation module as an additional (trainable) binding affinity prediction module, which was then fine-tuned for binding affinity estimation using the wellknown PDBBind 2020 dataset [25,26]. Overall, FlowDock's model design provides a simple add-on module within

MULTICOM_ligand to report (when requested) confidence scores for (predicted) protein and ligand heavy atom coordinates and binding affinity values for each ligand based on their (predicted) heavy atom coordinates (n.b., Pearson's correlation between the two: -0.127).

| Results and Discussions

The blind structure prediction benchmarking results of MULTICOM_ligand in the CASP16 experiment, as illustrated in Figure 2, demonstrate that our DL ensembling approach to protein-ligand structure modeling (n.b., ranked fifth among 34 predictor groups) reliably produces structurally accurate ligand-bound poses (~2.5 average (Model: 1) RMSD) of the diverse, pharmaceutically relevant protein complexes available in this experiment. Furthermore, Figure 3a,b (for affinity Stages 1 and 2, respectively) illustrate that MULTICOM_ligand's predicted protein-ligand binding affinities are modestly correlated (Pearson's R values of 0.30 and 0.31, respectively) with their ground-truth values (n.b., ranking fifth among 28 predictor groups), highlighting the real-world utility of our approach to estimating binding affinities for virtual screening in drug discovery [27].

As Figures 2 and 3a,b showcase, many of MULTICOM_ligand's pose and binding affinity predictions are highly accurate, yielding several predicted poses with an RMSD less than 1 and estimated affinities nearly identical to their corresponding true values. Nonetheless, in several other cases, MULTICOM_ligand's pose predictions yielded RMSDs above 4, indicating that in these cases our approach failed to identify the correct protein-ligand binding pockets for DL-based docking. Moreover, the gaps between MULTICOM_ligand's predicted affinities and their true counterparts were occasionally large, suggesting that in these cases MULTICOM_ligand was unsuccessful in differentiating weak from strong binding.

Like all other CASP16 ligand predictor groups, MULTICOM_ligand's affinity predictions given crystal protein-ligand structures as additional inputs in affinity Stage 2 were not statistically significant in their differences to those of Stage 1 (Pearson's R of 0.31 vs. 0.30), highlighting that, to make accurate binding affinity predictions, FlowDock's representations derived from primary sequence inputs were generally more useful to the model rather than additional structural context proved to be. In the following subsections, we examine a subset of MULTICOM_ ligand's CASP16 pose predictions to study its relative strengths and weaknesses revealed by the experiment.

| L1004

As one of the first pharma targets released to predictors for CASP16, ligand target L1004 represents a globular protein with a well-defined binding pocket for molecular docking. As such, the challenge presented by this target is largely in modeling the most accurate pose of this novel ligand within the pocket rather than locating the pocket itself. MULTICOM_ligand's topranked prediction for this target (Figure 4a) yielded a precise ligand RMSD of 0.510 and an lDDT-PLI of 0.963 (ranking 2nd overall). Interestingly, MULTICOM_ligand's rank-3 prediction achieved an even lower ligand RMSD of 0.483 (ranking 1st overall), suggesting that our structural consensus ranking heuristic mislabeled our most accurate pose for this target yet still ranked it among the ensemble's top-5 predictions. The top prediction's AlphaFold 3 protein structure for L1004 had a protein backbone RMSD (BB-RMSD) of 0.224, highlighting that our DL ensemble methods each had access to a highly structurally accurate (hololike) protein structure for ligand docking or pose prediction for this target, which contributed to their success in this case.

| L1009

Due to the hierarchical naming structure of CASP16's pharma ligand targets, target L1009 contains the same binding pocket as target L1004 yet asks predictors to provide poses for a new and conformationally distinct ligand. MULTICOM_ligand's topranked prediction for this target (Figure 4b) achieved a ligand RMSD of 0.603 and an lDDT-PLI of 0.950 (ranking 2nd overall), comparable to its predictions for L1004 with high overall accuracy. Again of interest, MULTICOM_ligand's rank-2 prediction yielded even better results with a ligand RMSD of 0.525 (ranking 1st overall), further emphasizing the importance in future work of identifying efficient ways of augmenting our structural consensus ranking heuristic (e.g., with FlowDock's predicted confidence scores).

| T1214

CASP16 incidental ligand target T1214 represents a beta barrel membrane protein structure interacting with a single PQQ ligand molecule. Figure 4c shows MULTICOM_ligand's (failed) top-ranked prediction of this target, which achieved a modest ligand heavy atom RMSD of 4.244 (n.b., 3.994 with the crystal protein structure) and lDDT-PLI of 0.357 (ranking 34th overall). As our initial AlphaFold 3 prediction of this target's beta barrel protein structure yielded a reasonable BB-RMSD of 1.687 (n.b., compared to the BB-RMSD of 0.822 achieved by the top-ranking group for this target), one possible explanation for the difficulties MULTICOM_ligand faced for this target is that membrane proteins constitute approximately only 5% of the PDB's composition [24]. Consequently, we posit that deep learning-based docking methods trained on common subsets of the PDB such as PDBBind [25] are likely to underperform for such targets, since their predictions are primarily optimized for docking with more common types of (e.g., helical) proteins. This suggests that MULTICOM_ligand's performance may be improved as new deep learning methods (in particular co-folding methods) trained on more balanced mixtures of biomolecular data are introduced.

| Conclusions

In CASP16, we introduced MULTICOM_ligand, a deep learning-based ensembling method for protein-ligand structure prediction combined with flow matching for joint structure and binding affinity prediction. Its blind assessment results in the CASP16 experiment demonstrate its efficacy and utility for real-world drug discovery efforts. Future work could include investigating whether FlowDock's predicted confidence scores could enhance the ranking performance of MULTICOM_ligand's structural consensus heuristic and whether the latest DL co-folding methods such as AlphaFold 3 [13], Chai-1 [28], and NeuralPLexer 3 [29] may benefit from a DL ensembling approach like our new MULTICOM_ligand method or if their predictions may be augmented with additional rank-ordering and binding affinity estimations provided by lightweight generative models such as our new FlowDock model [9].