{'fraction_non_alphanumeric': 0.03713814633618942, 'fraction_numerical': 0.013460550855474767, 'mean_word_length': 4.8989158163265305, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Predicting the protein-ligand affinity from molecular dynamics trajectories

Accurate protein-ligand binding affinity prediction is essential in drug design and many other molecular recognition problems. Despite many advances in affinity prediction based on machine learning techniques, they are still limited since the protein-ligand binding is determined by the dynamics of atoms and molecules. To this end, we curated an MD dataset containing 3,218 dynamic protein-ligand complexes and further developed Dynaformer, a graph-based deep learning framework. Dynaformer can fully capture the dynamic binding modes by considering various geometric characteristics of the interaction. Our method shows superior performance on the benchmark dataset over the methods hitherto reported. Moreover, we performed virtual drug-screening on heat shock protein 90 (HSP90) by integrating our model with structure-based docking. We benchmarked our performance against the conventional knowledge-based approach, demonstrating that our method can identify two lead compounds with submicromolar activity (0.207\textmu{}M and 0.517\textmu{}M). We anticipate that large-scale MD dataset and machine learning models will form a new synergy, providing a new route toward accelerated drug discovery.


Introduction
Molecular recognition is of fundamental importance in biology. The interactions between proteins and small molecules regulate biological functions and systems such as enzyme catalysis and signaling pathways [1,2]. Therefore, the research of drug design and discovery has primarily focused on manipulating the molecular binding on proteins. As the binding affinity measures the strength of such interactions, accurate affinity prediction is one of the most critical tasks in computational-aided drug discovery [3].
Recent advancements in computing have enabled protein-ligand binding affinity calculation through physics-and knowledge-based approaches [4]. The molecular docking algorithms are an essential coarse-grained physical model that efficiently predicts and evaluates the binding conformation of the small molecule. However, such models usually incorporate expert knowledge, hand-tuned features, and parameters, leading to biased results [5]. All-atom molecular dynamics (MD) is another popular computational technique to investigate the structural and dynamical properties of the biological system. It is predicated on Newtonian dynamic equations and reveals full atomic detail at a fine spatial-temporal resolution, offering accurate binding affinity prediction (i.e., alchemical free energy calculations) [6,7]. Despite various advantages of MD, it consumes prohibitively high computational resources, rendering MD-based large-scale affinity prediction infeasible [8,9].
Data-driven approaches are another emerging tool in the affinity prediction task [10]. Primarily benefiting from the rapidly growing amount of experimental protein-ligand complex data, data-driven approaches have achieved a high accuracy level by learning the binding modes directly from threedimensional (3D) structures. Nonetheless, current state-of-the-art methods require sophisticated feature engineering while offering limited generalizability. More importantly, various attempts at accuracy improvement, such as topology-, spectral-, sequence-or atom pair-based methods, have mainly centered on feature engineering from the static 3D structure [11][12][13][14][15]. So far, none of which has taken the dynamics of protein-ligand interaction into account. In fact, all proteins are highly dynamic biomolecules, constantly interacting with their surroundings and often adopting various conformational states. As a result, proteins exist in an ensemble of conformational states. Nevertheless, the experimentally resolved structure only represents a temporally averaged snapshot of the most populated state, and the dynamic effects of the structure are inevitably all but lost. Furthermore, the ligand is partially exposed to solvent, therefore, it does not only interact with protein but also with the solvent (i.e., liquid water and buffer ion). With sufficient MD simulation duration, the influence of ligand-solvent interaction and other long-and short-ranged interactions (e.g., Van der Waals interaction, hydrogen bond, and any other polar/apolar interactions) on the ligand conformations are naturally embedded in the MD trajectories. Despite all the benefits, no MD-based protein-ligand dataset is as yet reported.
In this work, we curated a large-scale MD protein-ligand dataset and introduced a dynamic-graph-transformer (Dynaformer) model trained on the MD trajectories. To fully capture the dynamic binding modes, we developed the Dynaformer based on a roto-translation invariant encoding scheme, taking various interaction features into account, i.e., interatomic distances and angles between bonds. Compared to other methods, Dynaformer showed the highest accuracy on the benchmark dataset. Furthermore, we discussed several underlying thermodynamic binding mechanisms that contributed to the improved performance of Dynaformer. We also experimentally demonstrated that, compared to the molecular docking method, our model can identify new lead compounds in a more efficient manner.

Thermodynamics determines the binding affinity
The driving forces that determine the association between protein and ligands are a collective consequence of various complex interactions and energy exchanges among the protein, the ligand, and the aqueous environment [16]. Consequently, as shown in Fig. 1, the state distribution of a biological system is dictated by the thermodynamic laws [17]. Eq. 1 indicates that the actual ligand binding strength is determined by the Gibbs free energy change ∆G, which is a sum of enthalpy change (∆H) and entropy change (T ∆S). Nonetheless, the experimentally resolved state occupies the lowest energy level ∆H 1 , which cannot represent the full picture of the binding mode (Fig. 1). Therefore, the best way to improve a predictive model is to introduce the dynamic effects as inductive bias into the training dataset. It is proved that both ∆H i and ∆S i can be recovered from an MD ensemble of protein-ligand interaction snapshots [18], which is mathematically expressed as: in which T represents the temperature and N indicates the total number of states sampled from MD trajectories. ∆H i and ∆S i represent the enthalpy and entropy at the i-th snapshot (see Methods for detailed derivation). Fig. 1 The state ensemble determines the binding affinity. The Gibbs free energy of protein-ligand complex is a sum of entropy and enthalpy of states, which determines the binding affinity. The state ensemble follows a certain distribution that is governed by thermodynamic laws.

Dynaformer framework
The structure of molecules can be interpreted as graph-structured data, where nodes represent atoms and edges stand for the relationships between atoms, i.e., chemical bonds or non-covalent interactions ( Fig. 2(B, C)). Such graph structure retains the full physical, chemical, and geometric information within the dynamic trajectories. Fig. 2(D) illustrates the basic architecture of Dynaformer, consisting of the multi-head self-attention module and a feed-forward network (see Methods). Dynaformer is directly built upon the Graphormer, which achieves state-of-the-art performances on a wide range of molecule-related tasks [19]. However, Graphormer can only handle isolated molecules with covalent bonds. Whereas in a protein-ligand complex, atoms with non-covalent interactions and their geometry in the binding site play key roles in determining the binding affinity. To capture these intricate features, Dynaformer introduces additional attention bias terms to encode the information regarding the interatomic distances and angles between bonds ( Fig. 2(E)).
To mimic the in vitro aqueous environments, we developed an MD data preprocessing protocol that solvates and neutralize the protein-ligand complex in a simulation box with water molecules and ions ( Fig. 2(A)). In total, 3,218 protein-ligand complexes are selected for MD simulation and the detailed selection criteria are provided in Methods. The ideal duration for the simulation is the time scale at which a biological phenomenon takes place. Nevertheless, the time scale is highly case-specific and could lead to an explosion in computation cost. In this study, we selected 10 nanoseconds (ns) for all systems for the following practical reasons: 1) Most systems reach convergence within 10 ns. 2) Characteristic binding-related events can be observed within 10 ns (discussed in the next section). As shown in Fig. 2(F), each data in the MD database contains 100 uniformly sampled frames from MD simulation. For the molecular simulation packages, we adopt the latest version of NAMD 3.0, which maximizes the performance of moderate-sized (the number of atoms is less than 1 million) MD simulations for modern GPU architectures [20]. To build an efficient interactive graph, we only account for the protein atoms within a certain range, since the strength of interatomic interactions decays rapidly within a fewÅ. Hence, the relative position of the ligand to the protein is important, which can be characterized by root-mean-square deviation (RMSD) from its initial position at the binding site (see Methods). Favorable configurations usually exhibit high stability and the construction of the graph is straightforward; but for unstable configurations, it is often not possible to construct an interactive graph representation as the ligand can leave the original pocket and enters into the solvent. Therefore, to gain comprehensive insights into the configuration stability of the MD dataset, we performed RMSD analysis over all MD simulations (3,218 complexes). Fig.  3(A) demonstrates the mean RSMD (ligand position relative to the original pocket) normalized by the number of atoms of the ligand (# atom) versus the # atom. We characterized the configuration with normalized mean RSMD lower than 0.1 as 'stable' and higher ones as 'unstable'. Fig. 3(A) suggests that 90 % of the configurations are stable, as indicated by their low normalized mean RSMD. Additionally, the normalized mean RSMD quickly rises when the # atom is below 37, suggesting that the configuration stability is highly correlated to the # atom. Such correlation can be explained by the fact that smaller ligands (lower # atom) generally interact less with protein residues than bigger ones and are more susceptible to the impact of the random motion of solvent molecules. Not unexpectedly, Fig. 3(B) shows that the center of the binding affinity distribution of smaller molecules (# atom < 37) is located at pK i = 4.2, which is significantly lower compared to that of the bigger molecules (pK i = 7.8). Fig. 3(C) showcases an example (# atom = 80, PDB ID: 2yge), where the binding is very stable (RSMD = 0.5 ± 0.2Å). While Fig. 3(D) shows a small ligand left the original pocket (# atom < = 24, PDB ID: 4y3j) in which RMSD indicates that the ligand drifts away from the binding site after 6ns. Indeed, the experimental evidence corroborates the MD observation, showing a very low affinity (pK i = 2.17). Meanwhile, Fig. 3(E) shows the Dynaformer predictions of a stable complex (normalized RSMD = 0.08, PDB ID: 4ciw, pK i = 4.82). Although the affinity prediction of individual frames fluctuates from pK i = 6 to 4 due to the dynamic surrounding environment of the ligand, the averaged dynamic affinity prediction (pK i = 4.71) is close to the ground truth. As the affinity predictions using averaged dynamic and experimentally measured structure showed no essential difference on most configurations, we used the experimental measured structure for the model inference on the benchmark dataset. For the performance evaluation and comparison, we used the PDBBind 2016 core set as the benchmark dataset, which is a stand-alone test set for protein-ligand binding affinity prediction (see Methods). We benchmarked Dynaformer against other major methods and the results are shown in Fig. 3(F-H), which include Root Mean Square Error (RMSE), Mean Absolute Error (MAE), and Pearson correlation (Pearson r) of best-performing models in its categories (detailed analysis see Supplementary information). It is evidenced that RSME, MAE, and Pearson r of Dynaformer outperform other baseline methods by a large margin. Additionally, we selected SIGN as our baseline model for prediction comparison in the next section, as it shows the second-best performance in all three metrics.

Molecular dynamics improves affinity prediction
The improved performance is likely due to Dynaformer's ability to learn the dynamic binding modes from MD dataset. Fig.4 displays three cases from the PDBBind 2016 core set, in which Dynaformer shows high accuracy where SIGN deviates considerably. The detailed analysis shows that the dynamic binding process in these cases significantly differs from its static counterpart. These case studies unequivocally suggest that our MD-based dataset contains more information regarding the binding process and Dynaformer is able to learn from it.
The first case (PDB ID: 2v7a) is shown in Fig. 4(A), ITK kinase domain in complex with a ligand, demonstrating the entropy and enthalpy compensation [21]. Here, SIGN underestimated the binding affinity by a large margin of more than 40%. While both static and dynamic data contain enthalpy information, the dynamic data contained extra entropic information (the ligand flexibility), which had been learned by Dynaformer and results in better accuracy. As it can be seen in the zoom-in figure of Fig. 4(A), the head of the ligand (pyrrolo [3,4-c]pyrazole group) was tightly buried in the pocket, whose contacts with protein are highly conserved during the MD simulation. The two nitrogen atoms of the pyrazole group interacted with the carbonyl oxygen of GLU-316 and with the amide nitrogen of MET-318, whereas the nitrogen of the amide group forms a hydrogen bond with the carbonyl oxygen of MET-318. In addition, the pyrrolopyrazole core held hydrophobic contact with ALA-269, VAL-256 and LEU-248. These favorable non-covalent interactions led to a large negative enthalpy change term. Whilst the Root Mean Square Fluctuation (RMSF) plot( Fig. 4(B)), which measures the fluctuation of a group of atoms around their average positions, indicated that the N-methylpiperazine group exhibits less contact with protein and wiggles randomly due to its exposure to the solvent, resulting in a positive entropy change term. Consequently, this ligand showed a strong affinity towards protein (pK i = 8.3).
The second case (PDB ID: 3udh) is demonstrated in Fig. 4(C), where the protease endothiapepsin as the receptor interacted with 1-methyl-L-histidine derived from a novel fragment library [22]. Various conformers were observed during the MD simulation and the ligand drifted away from the original pocket after a few femtoseconds (fs). Such drift-away event was demonstrated by a large deviation (more than 4Å RMSD) of the ligand from its initial position ( Fig. 4(D)). Consequently, such molecules should have a relatively low binding affinity. Indeed, such phenomena were commonly observed in the MD simulation of low-affinity protein-ligand complex (e.g., 4y3j and 3gv9), but SIGN prediction (pK i = 5.9) could deviate as much as 200% from ground truth, indicating that models can only learn such knowledge by training on the MD trajectories. Fig. 4(E) demonstrates the Dynaformer's ability to handle the affinity cliff, by which the ligands with similar conformation possess distinct binding affinity. As shown in Fig. 4(F-H), three thiophenes bind to PTP1B (Protein tyrosine phosphatase 1B) with different potency [23]. 2-phenylthiophene, as the structural core of the three compounds, adopts a nearly identical binding pose to occupy the active site, and tail groups (phenyl ring) entered a new subsite (Supplementary information). The RMSD plot in Fig. 4(I and J) indicates the ligand from 2qbp exhibits high conformational stability at the binding site (RSMD < 1Å), likely due to the Pi-stacking interaction with residue PHE-182 and hydrogen bonds with residue LYS-120 and ARG-221. On the other hand, 2qbq's ligand left the original active site after a very short time (around 1 nanosecond) while the structural core of the ligand of 2qbr still occupies the main active site thanks to the π-stacking interaction with residue TYR-48, albeit with less stability (RSMD > 3Å). These structurally similar but dynamically different ligands provide distinct binding strength (pK i of 2qbp = 8.4, 2qbq = 7.4, 2qbr = 6.3). As shown in Fig. 4(E), Dynaformer provides more accurate predictions (pK i of 2qbp = 7.5, 2qbq = 7.0, 2qbr = 5.9) in comparison with SIGN. More importantly, it differentiates the rank of the binding strength (2qbp > 2qbq > 2qbr), thus providing more reliable guidance toward lead optimization. Notably, the ligand from 2qbp exhibited high conformational stability whereas 2qbr displays high flexibility. such difference was also observed in the interaction diagram. PS, HP, and HB stand for Pi-stacking interaction, hydrophobic interaction, and hydrogen Bond respectively.
Nevertheless, the spatial-temporal patterns of the binding process can be difficult to learn in some cases. Fig. 5(A) shows an example, in which Thrombin as a serine protease is involved in blood coagulation and an inhibitor with a bis(phenyl) methane moiety bind to the active site of thrombin with remarkable potency (PDB ID: 2zda, pK i = 8.4) [24,25]. In this case, both the predictions of SIGN and Dynaformer (pK i = 5.9 and 6.4) did not agree with the experimental result. A plausible explanation might be that the protein with high conformational flexibility violates the assumption that binding does not render significant protein structural change (see Methods), thereby making accurate calculation and prediction of the binding affinity very difficult. This is evidenced in the protein RMSD plot, in which the configuration is not stabilized even after 10 ns (Fig. 5(C)). The RSMF plot also indicates that a few protein side chains show a high degree of freedom (Fig. 5(B)). Consequently, the ligand experienced a drastic change in surroundings, as evidenced by the difference between ligand position at frame 100 and frame 1, and affinity prediction was very difficult in such a scenario. Nevertheless, protein stability could serve as an important criterion for identifying the outliers in the dataset and further improve the accuracy of learning models.

Lead discovery with Dynaformer
As a proof-of-concept, we experimentally validated that Dynaformer is capable of discovering potent lead compounds by re-scoring the docked molecular structures. The original docking score often has a large number of false positives in the top-scored lists [26]. As a result, one needs to perform hundreds of experiments to identify the true lead compounds. Here, we demonstrate that, by re-scoring the docked poses, Dynaformer could more efficiently deliver the lead compounds.
Specifically, we selected Heat Shock Protein 90 (HSP90) as our target protein. Hsp90 is an important chaperone protein involved in pathways of many refractory diseases, including cancer, neurodegenerative diseases and viral infections. A common strategy to inhibit HSP90 functioning is to design small molecules targeting its ATP binding pocket to suppress the ATPase activity in the N-terminal domain. Because studies of inhibitors against HSP90 are relatively abundant, there are 18 protein-ligand structures associated with ATPase activity in the training data with pK i ranging from 3.84 to 8.32. All of these structures could be categorized as the known inhibitor scaffolds summarized in the literature [27,28].
Subsequently, we performed docking on HSP90 against the ChemBridge DIVERSet-EXP Library containing 50,000 small molecules that cover a diverse pharmacophore space. The docked conformers were obtained using Autodock Vina, a widely used open-source docking software [29,30]. The reference protein-ligand structure (PDB ID: 2xdl) was chosen from the PDBBind 2016 core set for identification of the binding pocket. The lead discovery protocol can be summarized as follows: 1) The ordinary docking was performed using the Autodock Vina; 2) molecules were re-scored and re-ranked based on the Dynaformer prediction; 3) The top-ranked binders are then clustered by scaffold and visually inspected for prioritization. To ensure the fairness of prioritization, we created a pool consisting of molecules among the top 10% from both Autodock Vina and Dynaformer. From the pool, 20 molecules are chosen by visual inspection without knowing their origin. Finally, the Surface Plasmon Resonance (SPR) experiments were conducted for each of the 20 molecules to measure the inhibition constant K i . The details and results of both dry lab and wet lab experiments are stated in the Supplementary Information. Fig. 6(A and G) show the structures, experimentally measured K i and the ranking of the top 10 molecules. Numbers in red indicate the Dynaformer rank is better than Autodock Vina and vice versa. The K i comparative statistics are summarized in Fig.6(F). It is evident that Dynaformer can identify the lead molecules more efficiently by assigning them higher ranks. Specifically, the Dynaformer ranked the top 3 molecules at 0.82%, 3.71%, and 0.21% respectively, surpassing Autodock Vina by a large margin. In particular, molecule 1, 3, 4 and 6 show K i of 0.207µM, 1.70µM, 2.50µM and 4.50µM without known precedent in the literature, which can serve as favorable candidates for further lead optimization. The possible binding mode is shown in Fig. 6(B, D and E), which is consistent with findings from rational drug design [31,32]. For example, a previous study shows that 1) the hydrogen bond networks between inhibitors and residues (ASP-93, THR-184, and LYS-58) are crucial for stable binding; 2) the occupancy by the aromatic groups of lipophilic pocket bottom made up of MET-98, PHE-138, TYR-139, VAL-150, and VAL-186 leads to enhanced affinity and selectivity [27,28,33,34]. In addition, we found some molecules resemble previously reported scaffolds. Molecule 2, 7, 10 and 11 are inhibitors containing the resorcinol-like scaffold from natural products, the possible binding mode of 2 is shown in the Fig. 6(C). Molecule 5 and 9 share a similar backbone that belongs to purine scaffold while molecule 8 is similar to an aminotriazine compound reported previously in a fragment-based screening [28,35]. Therefore, the wet lab results showed that Dynaformer is capable of delivering lead compounds with good binding affinity and novel scaffolds. It can be seen that Dynaformer is able to identify the lead compounds more efficiently by assigning them a higher ranks.

Conclusion
In this work, we presented a graph transformer model to learn the underlying physico-chemical principles of protein-ligand interaction from the MD trajectory. For the first time, we constructed a large-scale dynamic protein-ligand dataset, which encompasses various dynamic effects. We also designed a graph model that encodes geometric-chemical characteristics of the binding process.
As a result, our model outperforms the current state-of-the-art on the benchmark dataset significantly. Here, we manually curated an MD protein-ligand dataset with 3,218 trajectories, but it represents only a small fraction of the available database. Besides, not all protein types are included using the current study.i.e., the membrane protein requires a bi-lipid layer rather than an aqueous environment and many metalloprotein simulations are beyond the capacity of any MD simulation software due to the highly complicated nature of protein-metal interaction. In the future, one might develop a more efficient way to automate the large-scale MD simulation pipeline. Moreover, the MD itself suffers from some drawbacks, such as the inherently approximate nature of the force field calculation and insufficient sampling of MD trajectories. Further progress on the force field accuracy of MD might help improve the quality of the MD dataset.
By re-scoring the docked poses with our model, we showcased the effectiveness of machine learning-based lead compound discovery. We validated that our model can considerably reduce the experimental cost compared to the traditional knowledge-based method. Nevertheless, as the current virtual drug-screening relies on molecular docking, and the ligand database, it is limited both in scoring power and search space. One possible method for further development is to replace structure-based docking with the generative model to establish a de novo virtual screening pipeline, as the generative model could generative model potentially provide new molecules which do not exist in the ligand library. Finally, by fusing molecular dynamics and machine learning, we create not only a symbiotic relationship between the two scientific fields but also provide better tools for describing biological host-guest problems.

Thermodynamics of binding
In the following, we show that the free energy ∆G for protein-ligand binding, which determines the binding affinity, can be more accurately described using molecular dynamics trajectories than the 3D static counterpart. As it was shown in Ref. [16,36], ∆G can be expressed as the sum of two components, the vacuum binding free energy and the solvation free energy difference: where R and T stand for the universal gas constant and temperature respectively; K i is the thermodynamic equilibrium constant of the binding process; ∆G vacuum is the binding free energy between protein and ligand in the vacuum. this energy can be further split into two parts: where E bind pl is the ensemble-averaged protein-ligand interaction energy and the term −T ∆S is the entropic contribution. Based on the interaction entropy scheme [36], ∆G vacuum can be expressed as: in which β is equal to 1/RT and E bind pl stands for the interaction energy of protein and ligand (e.g., hydrogen and hydrophobic interaction, Van der Waals interaction, Pi-Pi interactions, long-range electrostatic interactions, etc.); E pw , E lw are the interaction energies of protein-water, ligand-water; E p , E l , E w are internal energies of the protein, ligand, and waters, respectively. The enthalpy terms ∆H i equals to E pl i , E pl i is the interaction energy between protein and ligand at state i. Meanwhile, the entropy term ∆S i equals e ∆E pl i /kT , where ∆E pl i represents the deviation of state i and the average interaction energy (k is the Boltzmann constant). Equation. 4 proves that the binding free energy can be evaluated using the ensemble average. Such ensemble calculation can be performed by averaging over the MD simulations: where N is the number of MD frames. Finally, ∆G solv is the free energy difference between the protein-ligand complex and the sum of stand-alone protein and ligand systems. In principle, computing ∆G solv would require three independent MD simulations of protein, ligand and complex. In order to simplify the calculation, one typical assumption is that no significant conformational changes occur upon binding, so that ∆G solv can be approximated using one single frame. Therefore, compared to the 3D static structure, which is essentially a snapshot of the most populated states of the complex. An ensemble of uncorrelated snapshots collected from an equilibrated molecular dynamics simulation can take both enthalpy and entropy into account, presenting a more accurate input feature for predicting protein-ligand binding affinity.

Datasets preparation
PDBBind dataset. PDBBind is a public protein-ligand binding benchmark dataset curated from Protein Data Bank (PDB) [37,38]. PDBBind is a comprehensive crystal structure collection of target proteins, ligands and corresponding experimentally measured binding affinity data. PDBBind dataset consists of three parts: general set, refined set and core set. Structures from PDBBind version 2019 are used to conduct our study. Explicitly, The general set contains totally 21,382 complexes. The refined set is a subset of the general set containing 4,852 complexes with better quality. The core set is a stand-alone test set consisting of 285 high-quality complexes for testing the accuracy of binding affinity prediction. The general set and refined set are updated annually, while the core set is not. Therefore, we evaluate our model and all baseline models on the latest version, PDBBind core set 2016.
MD dataset. We choose complexes with relatively high quality in PDB-Bind refined set as the initial state of MD. We used web-based CHARMM-GUI to build complex systems and prepare their inputs for MD simulation [39,40]. The MD is performed using NAMD, which is a parallel molecular dynamics program designed for high-performance simulation of large biomolecular systems. The MD simulation is based on the CHARMM36 force field. Each complex is solvated in a truncated periodic TIP3P water box, and the minimum distances from the surfaces of the box to the complex atoms are set to 10 A. Counter ions are added to neutralize systems, and the initial configuration is decided using a short Monte Carlo simulation. CHARMM ligand topology and parameter file is generated using ParamChem service. The missing residue is modelled using GalxyFill. We used A semi-isotropic Langevin-piston method (piston period of 50 fs) with piston decay of 25 fs and Langevin temperature coupling with a friction coefficient of 1 ps −1 to control the temperature and pressure, respectively. The simulation time step was set to 2 fs in conjunction with the SHAKE algorithm to constrain the covalent bonds involving hydrogen atoms. Simulation temperature is maintained at 303.15 K. The simulation time step was set to 2 fs in conjunction with the SHAKE algorithm to constrain the covalent bonds involving hydrogen atoms. Simulation temperature is maintained at 303.15 K. For a stable protein-ligand complex, 2 ns is sufficient for computing free energy calculation. On the other hand, ligand-drift-away events (or other events) become more obvious if MD duration is longer. Therefore, 0.5 ns NVT (constant volume and temperature) is set for equilibrating stage before 10 ns NPT (constant pressure and temperature) dynamics is performed. In practice, such procedure can fail on several occasions: 1) Complex small molecules that lead to failed force field parametrization (e.g., Glycans etc). 2) The protein structure contains too many missing residues. 3) specific type of metallic ion, due to its complex electronic structure (e.g., Fe, Mn, Co etc.). After removing failed cases, totally 3,218 MD simulation are performed, each containing 100 frames.

Model architecture
Dynaformer is a Transformer-based model that learns the protein-ligand binding features from molecular dynamics data. On top of Graphormer [19], we propose a roto-translation invariant encoding scheme to learn the interaction energy E bind pl which associates with the binding affinity. On a higher-level, the Dynaformer predicts the binding affinity given the collective interaction patterns (from molecule dynamics) of a ligand at the active binding site. Given the protein-ligand interaction graph as G = (V, E) with V = {v 0 , v 1 , v 2 , . . . , v n } and E = {e i,j | i, j ∈ V}, where n represents the number of atoms and v i , e ij represents the input atom or edge features, the design of Dynaformer is described as follows. Note that v 0 is a virtual node for aggregating information from all other nodes to predict the binding affinity. Enthalpic motivation of Dynaformer. As mentioned before, The interaction energy of protein and ligand E bind pl is directly associated with the binding affinity through Equation 4. E bind pl mainly consist of the non-covalent bonds between them, which can be fully represented by distances between atoms and angles between bonds. To encode the E bind pl , we exploit Gaussian Basis Function (GBF), which can be defined as Equation 7, where µ and σ are learnable scalars and x are embedded scalars representing the encoded distances d χ(i,j) or angles s ψ(i,j) . Two functions d χ(i,j) and s φ(i,j) have the form of ·(i, j) : V × V → R 1 that encode the distances or angles and the corresponding atomic features. Hence, the inter-atomic energy (i.e., the enthalpic contribution) can be encoded by Equation 8,9: where W χ and W ψ are the parameter matrices for linear transformation. Dynaformer layer. In each encoder layer, Dynaformer has two components, a multi-head self-attention module (MHA) and a feed-forward network (FFN). The encoded distances and angles are fused into MHA module which is used to update atom representation and the energy contribution of each atom. In each head of MHA, the attention is calculated as Equation 10,11, where Q = HW Q , K = HW K , V = HW V are linear transformed atom representations denoted by H = [h 0 , · · · , h n ]. The attention module is used for updating atomic representation by a weighted summation from all atoms. Similar to d χ(i,j) and a ψ(i,j) , s φ(i,j) is an additional attention bias, which are learnable scalars indexed by function φ(i, j) = [v i |SPD ij |v j ] W ψ where SPD ij is the shortest path distance between node i and j on the interaction graph. The other term c ij = 1 N N n=1 e n w T n encodes the chemical bonds on the shortest path between node i and j, where w n denotes the weight that projects e n to a scalar.
Predictions through the virtual node. A virtual node is a meta node that collects information for affinity prediction from all other nodes. In every layer of Dynaformer, other nodes (atoms) gradually aggregate information to the virtual node through Equation 10, and at the last layer, the representation of the virtual node is exploited to predict the binding affinity by a linear layer. Further, to extend the receptor field of the binding pocket, three fingerprint vectors are calculated and fused into the virtual node [41][42][43]. By design, such a scheme allows the Dynaformer to automatically learn the effective structural factors that strongly associate with binding free energy through dynamic trajectories, and the frames in the trajectories also serve as a rational data augmentation approach in the model training.
Ablation study. Details can be found in the Supplemental Information.

Evaluation Metrics
Binding affinity. In this work, the binding affinity of a protein-ligand pair is the dissociation constant (K d ) or more narrowly defined, inhibition constant (K i ). Both K d and K i are equilibrium constants that describe how potent a ligand is binding to a protein, as shown below. In this paper, we use the negative logarithm of K d or K i , i.e. − log 10 K d or − log 10 K i (pK d or pK i ) to represent the binding affinity.
[P L] K d /Ki [P ] + [L] Pearson r. Also known as the Pearson correlation coefficient, Pearson r measures the linear relationship between two sequences of data {x i } and {y i }. It is defined as: where N is the number of samples,x,ȳ are the mean values of x and y, respectively. r varies between [−1, 1] with 0 indicating no correlation and −1 or 1 indicating exact negative or positive linear relationship. This metric is wildly used in the virtual screen where rank order is crucial. Note that Pearson r is invariant to the scale nor location of the two variables, which means it does not suggest the error between prediction and groundtruth. MAE, MSE and RMSE. As a complementary part of Pearson r, we use MAE, MSE and RMSE to determine the average error of binding affinity prediction. They are defined as: where N is the number of samples, {x i } and {y i } are the predictions and the groundtruth binding affinity. RMSD and RMSF. The RMSD measures the average deviation of a matching set of atoms between two structures. It can be expressed as: where δ i is the euclidean distance between the atom i and the reference position. Whereas the RMSF is the fluctuation around an average, per atom or residue, over a sequence of structures (e.g., from a trajectory).
where x i is the coordinates of the particle i, andx i is the ensemble average position of i. Where the RMSD quantifies how much a structure diverges from a reference over time, the RSMF can measure the mobility of parts of the system. While RMSD is frequently calculated using the initial state, the RMSF should be calculated to an average structure of the simulation. An area of the structure with high RMSF values frequently diverges from the average, indicating high mobility. When RMSF analysis is carried out on proteins, it is typically restricted to backbone or alpha-carbon atoms; these are more characteristic of conformational changes than the more flexible side-chains.