Feng Zhou,Shiqiu Yin,Yi Xiao,Zaiyun Lin,Weiqiang Fu,Yingsheng J. Zhang
Structure–Kinetic Relationship for Drug Design Revealed by a PLS Model with Retrosynthesis-Based Pre-Trained Molecular Representation and Molecular Dynamics Simulation
Drug design based on kinetic properties is growing in application. Here, we applied retrosynthesis-based pre-trained molecular representation (RPM) in machine learning (ML) to train 501 inhibitors of 55 proteins and successfully predicted the dissociation rate constant (koff) values of 38 inhibitors from an independent dataset for the N-terminal domain of heat shock protein 90α (N-HSP90). Our RPM molecular representation outperforms other pre-trained molecular representations such as GEM, MPG, and general molecular descriptors from RDKit. Furthermore, we optimized the accelerated molecular dynamics to calculate the relative retention time (RT) for the 128 inhibitors of N-HSP90 and obtained the protein–ligand interaction fingerprints (IFPs) on their dissociation pathways and their influencing weights on the koff value. We observed a high correlation among the simulated, predicted, and experimental −log(koff) values. Combining ML, molecular dynamics (MD) simulation, and IFPs derived from accelerated MD helps design a drug for specific kinetic properties and selectivity profiles to the target of interest. To further validate our koff predictive ML model, we tested our model on two new N-HSP90 inhibitors, which have experimental koff values and are not in our ML training dataset. The predicted koff values are consistent with experimental data, and the mechanism of their kinetic properties can be explained by IFPs, which shed light on the nature of their selectivity against N-HSP90 protein. We believe that the ML model described here is transferable to predict koff of other proteins and will enhance the kinetics-based drug design endeavor.
Structure–Kinetic Relationship for Drug Design Revealed by a PLS Model with Retrosynthesis-Based Pre-Trained Molecular Representation and Molecular Dynamics Simulation

Drug design based on kinetic properties is growing in application. Here, we applied retrosynthesis-based pre-trained molecular representation (RPM) in machine learning (ML) to train 501 inhibitors of 55 proteins and successfully predicted the dissociation rate constant (koff) values of 38 inhibitors from an independent dataset for the N-terminal domain of heat shock protein 90α (N-HSP90). Our RPM molecular representation outperforms other pre-trained molecular representations such as GEM, MPG, and general molecular descriptors from RDKit. Furthermore, we optimized the accelerated molecular dynamics to calculate the relative retention time (RT) for the 128 inhibitors of N-HSP90 and obtained the protein–ligand interaction fingerprints (IFPs) on their dissociation pathways and their influencing weights on the koff value. We observed a high correlation among the simulated, predicted, and experimental −log(koff) values. Combining ML, molecular dynamics (MD) simulation, and IFPs derived from accelerated MD helps design a drug for specific kinetic properties and selectivity profiles to the target of interest. To further validate our koff predictive ML model, we tested our model on two new N-HSP90 inhibitors, which have experimental koff values and are not in our ML training dataset. The predicted koff values are consistent with experimental data, and the mechanism of their kinetic properties can be explained by IFPs, which shed light on the nature of their selectivity against N-HSP90 protein. We believe that the ML model described here is transferable to predict koff of other proteins and will enhance the kinetics-based drug design endeavor.

Protein–drug interactions can be described by both thermodynamic properties (i.e., equilibrium constant Ki, KD, and half-maximal inhibitory concentration IC50) and kinetic properties (i.e., association kon and dissociation rate constant, koff). Ligands with the same affinity may have different association and dissociation rates, so kinetic parameters often provide more useful information for drug design. Compounds having a short retention time in active sites may lead to low drug occupancy at targets and poor efficacy. In contrast, for compounds with a longer retention time or slower dissociation rate, their doses can be reduced to achieve better selectivity and fewer side effects. Studies show that kinetic properties are strongly correlated with the pharmacological activity of drugs in cells and in vivo. Heitman et al. studied the relationship among intracellular efficacy, affinity, and residence time of 10 A2A adenosine receptor agonists, finding that there is a stronger correlation between intracellular efficacy and residence time (R2 = 0.90) than equilibrium inhibition constant Ki (R2 = 0.13) of these agonists. One strategy to increase selectivity of drugs is to optimize the kinetic properties to obtain drugs with longer retention time. Applying kinetics-based drug design techniques, Khanna et al. optimized the kinetic property of EZH2 inhibitors and achieved a 10 times increase of retention time only by changing −OCH3 to −SCH3. Zhou et al. designed and synthesized 44 novel derivatives and discovered an ACE inhibitor as a clinical candidate with much lower koff than donepezil. Among them, compound 12 demonstrated a much better efficacy and a lower effective dose than that of donepezil. Using ensemble molecular docking and relative binding free-energy calculation, Miller et al. found a novel compound TDI-11861, which shows higher Sac (ADCY10) binding affinity and longer residence time than its precursor (3181 versus 25 s). Combining induced fit docking and MM/GBSA calculation cubic by cubic, Bai and Jiang constructed a sophisticated energy landscape of ligand dissociation process and found a potential lead compound, which gives better binding kinetic and thermodynamic properties with either TcAChE or mAChE. Because the wet lab experiments to determine the kinetic properties of protein–drug interactions are expensive and time-consuming, new in silico methods such as molecular dynamics (MD) simulation and ML are imperative. Free MD without any bias is the most straightforward way to calculate kinetic properties like kon and koff. For example, using the supercomputer ANTON, Shaw et al. first reported the binding process of a molecule directly from outside of the pocket by unbiased MD. With 400 μs MD, Pantsar et al. found the important role of water and protein conformational change on the kinetic properties of two p38α MARK inhibitors with short and long residence times but nearly identical activities. Buch et al. ran 495 short free MD and used MSM to construct the complete binding process of the trypsin/benzamidine complex, in which 187 trajectories were consistent with the binding mode found in the crystal. In general, the dissociation time of molecules from their binding pockets is often on the time scale of seconds to hours. Most computing tasks cannot afford such expensive MD simulations. Since 2000, many enhanced sampling methods have been developed to calculate koff. Capelli et al. used the infrequent metadynamics (InMetaD) to calculate the koff of the M2 receptor/iperoxo complex, and the result is only an order of magnitude away from the experimental value. InMetaD has also been used to calculate the koff of trypsin-benzamidine, kinase, biotin–streptavidin, and FKBP. Umbrella sampling (US) is useful for studying the dissociation pathway and calculating the free energy of different metastable states. For example, US was used to study the dissociation process and binding free energy of various protein complexes such as benzamidine-trypsin, acetylcholinesterase, cathepsin K, type I dehydrogenase, HSP90, and factor Xa. For more complex systems in which the dissociation pathway is straight, the steered MD (SMD) was used together with the US to find the optimal dissociation path and to calculate the potential of mean force (PMF). For example, US+SMD was used to calculate the binding free energy between maltose-binding protein and maltose. You et al. have applied US+MM/GBSA (molecular mechanics with generalized Born and surface area solvation) to study the mechanism of dissociation of two small molecules from ATP allosteric channels of p38 MAP kinase. PMF can also be used to estimate the retention time (1/koff) by calculating the mean free passage times (MFPTs) from the initial unbinding barrier. SMD was also used to calculate the relative residence times of p38α kinase and FAK. Potterton et al. predicted the relative residence time of 17 small molecules on the A2A adenosine receptor from changes in water–ligand interaction energy. Mollica et al. used scaled MD (sMD) to calculate the relative residence time of the same series of ligands on proteins such as HSP90, glucose regulation (Grp78), and the A2A GPCR. Zhou et al. combined sMD and InMetaD to calculate the residence time and dissociation mechanism of ASEM. In 2018, using random accelerated MD (RAMD), Kokh et al. performed large-scale koff calculation to calculate the relative residence time of 70 diverse drug-like ligands of N-HSP90 protein and discussed the effect of different substituents on the residence time of small molecules. The number of inhibitors of N-HSP90 was extended to 94 in a following paper. They found that there is a strong correlation between the calculated retention time and the experimental value. By combining RAMD and experimental results, Berger et al. found that PF-562271 is more selective for focal adhesion kinase (FAK) kinase than proline-rich tyrosine kinase 2. RAMD was also used to calculate the kinetic behaviors of a series of inhibitors of the T4 lysozyme, and the results showed that the ligands with a longer residence time have more intermediate metastable states along the dissociation pathways. Recently, RAMD was successfully applied to studying the relative residence times, dissociation mechanisms, and the allosteric effects for the two important membrane-embedded drug targets: β2-adrenergic receptor and muscarinic acetylcholine receptor M2. They found that the dissociation mechanisms observed in the relatively cheap RAMD simulations are consistent with the much more computationally expensive free MD and MetaD simulations. They also uncovered the relationship between the residence time and the allosteric modulation and associated changes in the ligand dissociation pathways. Although there are many successful cases of MD in the calculation of retention time and koff, due to the intrinsic drawback of the force field and insufficient sampling, the calculated retention time and koff often deviate from the ground true experimental values especially when the ligands have diverse scaffolds. Additionally, most enhanced sampling methods can only produce relative retention times. In recent years, through the development of artificial intelligence (AI), more accurate predictions of koff were seen by combining molecular simulation and ML. In 2016, a multi-target machine learning (MTML) with interaction fingerprints was used to predict the koff of HIV-1 protease. Based on electrostatic interaction terms and conformational dynamics, a classification model was used to train and classify 39 inhibitors of HIV-1 protease into four classes, predicting effects with 74% accuracy. This paper argues that electrostatic interactions contribute more to the koff value than van der Waals interactions. Similarly, Zhang et al. used a partial least squares (PLS) regression model to learn interaction features from position-restrained MD to predict the koff of p38 MAPK class II inhibitors and obtained good results. A 3D QSAR (quantitative structure–activity relationship) model was also applied to predict koff, and Qu et al. used VolSurf descriptors as features to predict the koff of HIV-1 protease inhibitors using PLS regression models. The Wade group predicted koff with high accuracy using the comparative binding energy (COMBINE) method for N-HSP90 and HIV-1 protease inhibitors. The same group also used interaction fingerprints (IFPs) from a large number of RAMD trajectories to predict the koff of N-HSP90. They found that information along the dissociation pathway is more important than the bound state. Similarly, Huang et al. predicted the koff of 37 inhibitors of HIV-1 using IFPs derived from SMD. Based on graph methods like contact principal component and pairwise mean Euclidean distance analysis, Bray et al. developed an ML model for clustering the unbinding trajectories and path reaction coordinate detection. A recent literature reported the prediction of koff for a large array of proteins and ligands using the random forest method with RF-score scoring functions and protein secondary structures as features. Although the correlation coefficient rP reached 0.78, the prediction power of this model became very poor for proteins that had not been seen in the training set and the model was not tested on another dataset. Compared to the traditional molecular simulation methods, the ML method to predict koff still faces many challenges such as proper feature representation and model selection as well as methods to increase its interpretability to guide drug design. To enhance ML predictability and increase model interpretability, we trained a new molecular representation named RPM to predict the kinetic koff and compared its performance to other popular pre-trained molecular representations, i.e., GEM (geometry-enhanced molecular representation), MPG (molecular pre-training graph-based deep learning framework), and a commonly used molecular predictor from RDKit. In our previous work, we trained a molecular pre-training model (G2GT) based on retrosynthesis reaction data. By cutting off the connections between encoder and decoder stages of the model, we used the output of the encoder as a representation (RPM) of the molecule to predict the koff property. The model was successfully tested in an independent dataset and showed a high correlation with experimental data. The mechanisms of the kinetic properties can be explained using IFPs. As a case study, we applied our ML/MD protocol to two N-HSP90 inhibitors, which were not in our model training dataset, and the predicted koff values showed a consistent trend with experimental values. By combining ML and MD simulation, coupled with the molecular representation features and IFP analysis, we unveil the features that have the most impact on the kinetic properties of the protein target and the interpretable details on how the molecules achieve their long retention time.

G2GT is a transformer model that uses a large amount of retrosynthetic reaction data as input to train. The model combines the encoder and decoder architecture of Transformer with a molecular graph. The encoder structure of G2GT extracts the structural information, the interaction relationship between molecules and output vectors, which the decoder uses to predict possible synthons. In this stage, we obtain the intrinsic relationship between input and output structures such as functional groups and intermolecular relationships. By cutting off the connections between encoder and decoder stages, we use the latent vector output as a representation of the molecule (RPM). The RPM was used with the PLS regression method to predict the koff of 501 ligands of 55 proteins. This new molecular representation takes consideration of molecular reactivity, function of chemical groups, and intermolecular relationships, thus showing superior performance in property prediction such as koff (complete process of koff prediction in Figure 1).

A common method to measure a pre-trained model’s performance is to test and evaluate it on downstream tasks. In order to show that our model can learn a rich and generalized representation, especially for downstream koff tasks, we compare our RPM with two mainstream pre-trained models: MPG and GEM, which are currently recognized as two state-of-the-art pre-trained models for molecular property prediction. MPG is a pre-trained molecular representation based on molecular graphs. It combines the classical Neural Message Passing (MPNN) for the Quantum Chemistry framework with the powerful transformer block to learn molecular representations. GEM is also a GNN-based molecular pre-training representation. In order to solve the problem of traditional GNN that three-dimensional molecular information will partially lose in graph representation, GEM introduces molecular graphs and bond angle graphs in which bond-angle information is calculated and shared to jointly model molecular representation. Our RPM and the two aforementioned models are all self-supervised pre-trained molecular representations. We also added the traditional RDKit molecular descriptors and then use a unified predictor (i.e., PLS) to predict the koff value. The training data is from the dataset of 501 inhibitors (dataset A) collected by Amangeldiuly et al. We adopt the leave-one-out strategy to train and predict these samples (Table 1). The results show that our RPM molecular representation gives the best prediction for the koff task, which refers to a higher Pearson correlation coefficient (rP = 0.76), compared with GEM (rP = 0.47) and MPG (rP = 0.65), and the lowest MSE (0.74).

To demonstrate the generalizability of our pre-trained molecular representations, we tested 38 ligands on an independent dataset collected by Liu et al. (dataset B). The 38 ligands are all from the same target protein, N-HSP90. HSP90 is a chaperone protein that assists in protein folding, blocking its N-terminal domain by small-molecule inhibitors leading to degradation of the client proteins and diminished tumor growth. Brough et al. reported the new HSP90 inhibitor (VER-52296/NVP-AUY922), which inhibits tumor growth by ∼50%. Specifically, a significant reduction of tumor growth in HCT116 colon cancer xenografts was observed with VER-52296/NVP-AUY922. In comparison with other proteins, N-HSP90 consists of more kinetic data and there are lots of literatures on koff simulation and prediction for N-HSP90.Figure 2 shows that our prediction value has a similar high correlation and lower MSE than the leave-one-out validation. This suggests that our pre-trained molecular representation has superior quality for cross-dataset representations and exhibits enormous potential for other molecular property prediction task.

The koff prediction by pre-trained molecular representations and the PLS model is based only on the molecular features and lack of the protein pocket information. For example, the same molecule may have different binding affinities and koff to different targets. To increase the selectivity of the target, one needs to consider both the features of molecules and the pocket information and protein–ligand interactions. RAMD is an efficient method to rapidly obtain the relative retention time, the dissociation pathways, and the IFPs. We performed 90 independent RAMD simulations starting from six different initial structures for each complex to avoid stochastic error. The initial structures of the complex are the cocrystal structures either from the Protein Data Bank (PDB) or from docking poses by Amangeldiuly et al. Short, unbiased MD simulations were independently performed six times, 10 ns each run, for each ligand to obtain the equilibrium structures of the bound state, followed by 90 RAMDs to calculate the retention time. The detailed MD method is in the Materials and Methods. The calculated −log(koff) was plotted against experimental values (Figure 3). The Pearson correlation coefficient between the calculated −log(koff) from RAMD and the experimental values is not high (rP = 0.534) for the 100 inhibitors from dataset A. The outliers, which were also observed by Kokh et al., can be attributed to the structures of the ligands, the binding modes, and the force field. However, the correlation between the calculated and the measured log(koff) is significantly improved when grouping the 100 inhibitors by scaffolds. For example, the calculated Pearson correlation coefficient for scaffold I, which contains the hydroxy-indazole structure, is 0.76. Kokh et al. also found that the retention time for the 10 compounds of the amino-quinazoline and amino-pyrrolopyrimidine was systematically underestimated. Consistent with the experiment, we also found that ligands bound to the helical conformation display slow dissociation rates in comparison with those bound to the loop conformation (Figure 3c). The predicted and simulated −log(koff) values have similar populations as the experiment values and the largest population in the range from 2.0 to 3.0 for helical conformation. We also performed RAMD simulation on the inhibitors of N-HSP90 protein from dataset B. The initial structures of the complexes are the cocrystal structures from PDB or from docking poses by Liu et al. The inhibitors that do not share a common scaffold with others were excluded, and the remaining 28 new inhibitors were used for RAMD simulation. Both the predicted −log(koff) values and simulated ones have a high correlation with the experimental values (rP = 0.73) (figure in the Supporting Information). As mentioned above, the ML model is good at predicting the koff but the factors behind the ML model that contribute to the prediction can only be understood by analyzing the IFPs. Therefore, we used the IFPs extracted from RAMD as features to study their relationship to the kinetic properties of N-HSP90 inhibitors. The IFPs that are generated from the 90 RAMD trajectories using the python script from Kokh et al. consist of hydrogen bond donor (HD), hydrogen bond acceptor (HA), hydrophobic (HY), salt bridges (IN), π–π stacking (AR), and halogen−π interaction (HL). From the PLS coefficient (Figure 4), one can analyze the most important IFPs, which have a strong correlation with the kinetic properties of ligands. These IFPs involve residues like ASN51, SER52, ASP54, ASP93, GLY97, ASP102, and LEU107, most of which are polar and surrounding the ATP binding pocket, TYR139 and THR184 inside the hydrophobic sub-pocket, and PHE138 at the entrance of the hydrophobic sub-pocket. The results are in line with the findings by Kokh et al. A distal residue ILE110 at the exit of the hydrophobic sub-pocket does not interact with ligands in the bound state but shows the importance during the dissociation process. Figure 5 outlines these important residues and their positions on N-HSP90. Figure 6 shows the schematic visualization of the RAMD dissociation trajectories of the four compounds with the slowest (Figure 6a,b) and fastest (Figure 6c,d) dissociation rates among the 28 simulated inhibitors from dataset B. There are more interactions of 5j6m_ligand than 6ei5_ligand along their dissociation pathways. Aforementioned important residues, ASN51, ASP93, GLY97, PHE138, and THR184, are all existent and retained during the dissociation process for 5j6m_ligand. However, for 6ei5_ligand, these interactions are either absent or transient during the dissociation process. Furthermore, the dissociation pathways show that for 6ei5_ligand, the transition flow between each intermediate is small while the transition to the unbound state has a larger flow (thicker gray arrow in Figure 6), which suggests a faster dissociation. On the other hand, for 5j6m_ligand, there is a large transition flow between clusters 1 to 6 before it transits into cluster 7, which is the last intermediate before fully unbinding. That means it costs more time for 5j6m_ligand to transit between different intermediates and the transition barrier to the unbound state is higher.

To conduct a case study, we searched the PDBbind database (http://www.pdbbind.org.cn/) and found two N-HSP90 inhibitors (PDB ID: 6fcj(65) and 6hhr(66)), which have experimental koff data and have not been included in our training datasets. These two inhibitors will be ideal to test our ML/MD protocol performance. 6fcj has a pyrazole scaffold and 6hhr has a triazolone one. The experimental −log(koff) values for 6fcj and 6hhr are 1.47 and 0.95, respectively. This suggests that 6fcj has a retention time four times longer than that of 6hhr. Our model’s predicted −log(koff) values for 6fcj and 6hhr are 2.12 and 1.55, respectively, which also show a four times difference in retention time. These results suggest that our model shows promise in predicting −log(koff) values. Furthermore, we plan to expand our model to include other scaffolds and explore different scenarios. Additionally, we intend to verify our model’s predictions through prospective experimental validation in the future. To decipher the mechanism behind the improved kinetic properties, we analyzed the trajectories and plot the IFPs along the dissociation path (Figure 7). The detailed RMSD with respect to the bound state of the 90 trajectories and the dissociation pathway in 2D graph representation are in the Supporting Informations. 6hhr_ligand and 6fcj_ligand are both inhibitors for the loop conformation of N-HSP90. There are two structural differences between the two molecules. First, 6hhr_ligand has a triazolone and 6fcj_ligand a pyrazole scaffold. Second, 6hhr_ligand has fluorine and sulfur as the substitution group and 6fcj_ligand has chlorine as a substitution group. Comparing the IFPs between the two ligands, we found that 6fcj has extra IFPs like HD_ASN51, HA_SER52, HD_LYS58, AR_LYS58, and HY_LEU107, while 6hhr has extra interactions with HA_GLY97, and HY_THR109. As shown in Figure 4, HD_ASN51, HA_SER52, and HY_LEU107 are all important residues. The fact that HA_GLY97 is not as important as HD_GLY97 means that a longer retention time is expected when GLY97 is a hydrogen bond donor not acceptor. Furthermore, both LEU107 and THR109 are at the entrance of the pocket, but LEU107 has more impact on the kinetic properties than THR109, likely due to forming hydrophobic interactions with chlorine.

The RPM-based PLS ML model was successfully applied to predict the koff of 138 inhibitors of N-HSP90 protein. The relative retention time is calculated using RAMD simulation. Both methods can produce the −log(koff), which has a high correlation to the corresponding experimental data. The molecular features from the pre-trained model and the IFPs from RAMD together give more information for the drug design. Using our ML koff prediction approach, we studied two new N-HSP90 inhibitors [6fcj_ligand and 6hhr_ligand]. The predicted −log(koff)s show consistent trends with experiments, and the mechanism behind is explained by IFPs. The resulting molecular mechanisms were elucidated in both molecular features and IFPs along the dissociation pathways. Our proposed protocol offers a feasible way for the kinetics-driven drug design approach and increases the success of finding molecules with proper kinetic profiles to the interested targets.

We present a novel deep learning model using numerous authentic reaction data for molecular representation learning. Based on our previous work for retrosynthesis prediction, we cut off the connections between the encoder and decoder stages and output the vectors as the representation of the input molecule. Because our model utilizes the data of the paired reactant product of chemical reactions, rather than the traditional single-molecule properties, we believe our model extracts more information on molecular functional groups and reactivities, especially certain interactions between molecules that current existing models do not take advantage of during training. Therefore, as a representation that takes consideration of molecular reactivity and intermolecular relationships, our model shows superior performance in property prediction such as koff (Figure 1).

As an extension of linear regression, multiple linear regression (MLR), and principal component regression (PCR), PLS regression is used to investigate multicollinearity between two groups of variables. For big-p, little-n problems, using PLS to build a model has advantages over MLR and other methods. Moreover, PLS has been successfully applied in koff prediction for HSP90, HIV-1, p38-MARK, and adenosine receptor with a small kinetics dataset. In this paper, we use PLS regression and RPM molecular representation to predict koff.

The relative retention times for the inhibitors were calculated from RAMD, developed by Lüdemann et al. Short, unbiased MD simulations were independently performed six times, 10 ns each run, for each ligand to obtain the equilibrium ensemble of the bound state. A random force was then applied to the ligand to accelerate the dissociation rate of the ligands. The force was set to be 20 kcal/mol/Å, and the calculated retention time can be controlled to be at a nanosecond scale. The retention time was calculated as the average value for the 90 independent RAMDs. In all MD simulations, the amber99sb-star-ildn force field and TIP3P model were used for protein and water. The ligand force field was parametrized with Gaff2 using ACPYE with AM1-BCC charges. The overall temperature of the system was kept constant, pairing independently for each group at 300 K with a Nosé–Hoover thermostat. The pressure was coupled to a Parrinello–Rahman barostat at 1 atm separately in every dimension. The temperature and pressure time constants of the coupling were 1 and 5 ps, respectively. The integration of the equations of motion was performed by using a leapfrog algorithm with a time step of 2 fs. Periodic boundary conditions were implemented in all systems. A cutoff of 1.2 nm was implemented for the Lennard-Jones and the direct space part of the Ewald sum for Coulombic interactions. The Fourier space part of the Ewald splitting was computed by using the particle mesh Ewald method, with a grid length of 0.12 nm on the side and a cubic spline interpolation. The MD simulation package used is Gromacs 2019.