Combined Free-Energy Calculation and Machine Learning Methods for Understanding Ligand Unbinding Kinetics.

The determination of drug residence times, which define the time an inhibitor is in complex with its target, is a fundamental part of the drug discovery process. Synthesis and experimental measurements of kinetic rate constants are, however, expensive and time consuming. In this work, we aimed to obtain drug residence times computationally. Furthermore, we propose a novel algorithm to identify molecular design objectives based on ligand unbinding kinetics. We designed an enhanced sampling technique to accurately predict the free-energy profiles of the ligand unbinding process, focusing on the free-energy barrier for unbinding. Our method first identifies unbinding paths determining a corresponding set of internal coordinates (ICs) that form contacts between the protein and the ligand; it then iteratively updates these interactions during a series of biased molecular dynamics (MD) simulations to reveal the ICs that are important for the whole of the unbinding process. Subsequently, we performed finite-temperature string simulations to obtain the free-energy barrier for unbinding using the set of ICs as a complex reaction coordinate. Importantly, we also aimed to enable the further design of drugs focusing on improved residence times. To this end, we developed a supervised machine learning (ML) approach with inputs from unbiased "downhill" trajectories initiated near the transition state (TS) ensemble of the string unbinding path. We demonstrate that our ML method can identify key ligand-protein interactions driving the system through the TS. Some of the most important drugs for cancer treatment are kinase inhibitors. One of these kinase targets is cyclin-dependent kinase 2 (CDK2), an appealing target for anticancer drug development. Here, we tested our method using two different CDK2 inhibitors for the potential further development of these compounds. We compared the free-energy barriers obtained from our calculations with those observed in available experimental data. We highlighted important interactions at the distal ends of the ligands that can be targeted for improved residence times. Our method provides a new tool to determine unbinding rates and to identify key structural features of the inhibitors that can be used as starting points for novel design strategies in drug discovery.

Introduction

A recent paradigm shift in drug design highlighted the importance of long residence time as a key objective in addition to strong binding affinity.[1] The residence time defines the timescale of the ligand bound in the binding pocket.[2,3] It is related to the overall rate of the unbinding process, which could consist of several steps. This therefore requires information about the corresponding high-energy transition states and free-energy barriers, which is challenging to obtain. Even if a drug interacts strongly with its target (high binding affinity), a short residence time can significantly reduce the efficacy of the drug.[4] While many successful drugs have been discovered on the basis of high binding affinity alone, recent studies have shown that for drug efficacy, the kinetics of drug-receptor binding may be, in some targets, more important than affinity.[2] The complexity of the drug-target dissociation may also involve several steps and complex pathways. Accordingly, promising hit candidates with high affinity were discarded for the next step of the drug discovery process due to their low residence time.[5,6]

A major challenge in drug discovery is finding a fast and reliable method to predict the kinetics of ligand–protein interactions.[7] Importantly, for experimental determination of ligand kinetics, ligands first need to be synthesized, which can be expensive and time consuming even for a moderate number of compounds. Different experimental methods have been used to obtain kinetics of ligand-receptor unbinding, such as radioligand binding assays, fluorescence methods, chromatography, isothermal titration calorimetry (ITC), surface plasmon resonance (SPR) spectroscopy, and nuclear magnetic resonance (NMR) spectroscopy.[6,8] Radioligand binding assays and fluorescence binding assays require binding with radiolabelled ligands, where they exploit the physical–chemical characteristics of the ligand between their free and complexed forms with the target. Several successful assays have been used to predict ligand–protein unbinding, for example, fluorescence resonance energy transfer (FRET)[9] or fluorescence correlation spectroscopy (FCS).[10] These methods can suffer from interference (especially fluorescence), lack of accuracy for short residence times, and high cost/hazard in the case of radioligands.[11] SPR is the most widely used assay to measure rate constants associated with (kon and koff) of ligand-receptor unbinding. The receptors are immobilized to a sensor that can distinguish the protein between its ligand-free form and its bound form. This method is label-free; however, the attachment of the protein to the probe may influence the activity of the protein due to conformational changes.[11] To offer a screening approach that alleviates these difficulties, various computational techniques have been proposed as alternatives to estimating the kinetics of unbinding events.[12,13]

Molecular dynamics (MD) is a powerful computational tool to understand at an atomistic level the behavior of biological processes such as protein–ligand interactions.[14] Unbiased MD simulations were successfully used in the initial stage of the drug discovery process, using either multiple independent relatively short simulations[15] or using specialized computer architecture, such as ANTON, where microsecond long simulations are readily accessible. However, due to the limited time scales typically accessible via MD simulations, it is often challenging to obtain sufficient statistical sampling required to calculate kinetic and thermodynamic properties accurately. Drug–protein unbinding processes occur on long time scales, typically ranging from millisecond to hours, depending on the nature and the strength of the interaction between the ligand and target. Some drugs, for example, aclidinium, deoxyconformycin, or tiotropium, have a half-life of hours,[16] requiring prohibitively long time scale simulations and highly demanding computer resources; therefore, enhanced sampling methods are required.[17]

To accelerate the simulations and sample rare events, different enhanced sampling techniques have been proposed to predict free-energy barriers and uncover the kinetics of biological events.[18,19] These methods include free energy perturbation (FEP),[20,21] metadynamics (MetaD),[22,23] temperature-accelerated MD (TAMD),[24] steered MD (SMD),[25] milestoning,[26] umbrella sampling (US),[27] replica exchange,[28] scaled MD,[29] smoothed potential MD,[30] transition path sampling,[31] τ-random acceleration molecular dynamics simulations (τ-RAMD),[32] and more recently a combination of enhanced MD with machine learning.[33−35] For most of these methods, a key factor is the identification of a collective variable (CV), representing a physical pathway, that allows the calculation of the free energy profile.[36] Hence, correct identification of appropriate CVs becomes a problem, with very few practical ways to build them properly.[37−39] These methods have already been used for ligand unbinding: for example, MetaD was used to predict the ligand–protein unbinding of p38 MAP kinase bound to type II inhibitors,[40] where depending on the set of CVs chosen, different values for koff were obtained, and the closest koff to the experimental data is still one order of magnitude lower. More recently, it was found that using a combination of MetaD and quantum mechanics/molecular mechanics (QM/MM) simulations, a more accurate prediction of the kinetics can be achieved.[41] The residence times of sunitinib and sorafenib in complex with the human endothelial growth factor receptor 2 have been calculated using steered molecular dynamics (SMD).[42] SMD was also used to calculate the unbinding free energy profile for TAK-632 and PLX4720 bound to B-RAF.[43] In both works, the ligands could be distinguished qualitatively to assess shorter or longer residence times; however, the predicted free energy barriers for the unbinding were significantly lower than the experimental data.

To produce accurate free energy profiles using biased simulations with many important degrees of freedom, we need to define an ideal set of CVs that map the full path of the reaction coordinate.[44,45] Usually, the vectors that describe this manifold are selected based on a priori chemical/physical intuition, typically based on the initial binding pose of the ligand. The same set of CVs are then kept constant and used for the full simulation. Considering only CVs from an initial structure implies possibly neglecting essential interactions that occur during the unbinding process, thus significantly affecting the free energy calculation. Additionally, structures resolved by X-ray crystallography or cryo-EM may capture the system in metastable states, which do not always reflect appropriate conformers for ligand binding.

In this work, we introduce a novel enhanced sampling method to obtain accurate free energy barriers for ligand–protein unbinding. Unlike existing methods, we also propose a method that subsequently can identify key molecular features determining the unbinding kinetics. We suggest an iterative way of assigning our CVs during the unbinding trajectory and then use these CVs as the driving force to pull the ligand out from the pocket and to perform the sampling for accurate free energy calculations. Similarly to, e.g., τ-RAMD (which, however, does not provide a free energy profile), there is no need to a priori select CVs; these naturally arise from unbinding trajectories that build a reliable path of unbinding taking the flexibility and dynamics of the system into consideration.

The CVs extracted from our trajectories sufficiently describe a full pathway for the unbinding process. Subsequently, we optimize this path in the space of the identified CVs to obtain a minimum free energy profile using the finite-temperature string method.[46] While different unbinding trajectories may lead to slightly different variations due to multiple local minima along the paths, we typically expect that the main transition state ensembles would be captured by all of these paths similarly after the convergence to the minimum free energy pathway. This is the main underlying assumption behind the finite-temperature string method, which was proven to work very well even for complex systems.[47,48] Our results accordingly show little variations in the unbinding free energy barriers using different starting pathways for free energy calculations.

In addition to determining unbinding rates, we also aim to identify key molecular descriptors that provide guidance for further design of drugs based on improved residence times. We propose a systematic approach to identify key low-dimensional sets of internal coordinates using machine learning (ML) approaches. Machine learning methods have been widely successful in multidimensional data-driven problems, which are also applied to biomolecular simulations to determine key CVs.[49−51] Here, we develop a novel approach making use of our obtained string unbinding pathway and, within that, the knowledge of the transition state (TS) ensemble. We explored two different ML methods in this study: neural networks (NN),[52,54] which provide efficient training on complex high-dimensional data, and gradient boosting decision trees (GBDT),[53] which allow straightforward evaluation of feature importances (FI). We generate unbiased “downhill” trajectories initiated at our TS and used these to train a ML model that predicts the fate of binding or unbinding.

To test this approach on a simple analytical model system, we generated trajectory data using a collection of one-dimensional (1D) model potentials, including one selected double-well potential. Our results demonstrate that our novel ML analysis can identify the key features correlated with this selected double-well potential to define the end states and thus can be used for key feature selection successfully. To demonstrate the applicability and accuracy of this approach on challenging complex biomolecular systems, we obtained free energy barriers for two ligands bound to CDK2 with PDB IDs of 3sw4 (18K) and 4fkw (62K) (Figure ).[55] Cyclin-dependent kinase 2 (CDK2) is a crucial regulator in eukaryotic cell growth: deregulation of CDK2 has been associated with unscheduled cell proliferation resulting in cancer progression and aggressiveness.[56,57] Selective inhibition of this protein makes it an appealing target in treating multiple tumors of specific genotypes.[58] Several molecules are currently under clinical evaluation as CDK2 inhibitors for cancer treatment, such as AT759,[59] AG-024322,[60] dinaciclib,[61] roniciclib,[62] milciclib.[63] Furthermore, CDK2 is an ideal benchmark system with its relatively small size and well-documented kinetic data for the binding of a range of different molecules.[55]

Figure 1

Illustration of the simulation system. (a) CDK2 bound to two different ligands: (b) thiazolyl-pyrimidine derivative (18K) and (c) carboxylate oxindole derivative (62K), originated from PDB structures 3sw4 and 4fkw, respectively. Structural details of ATP pockets are shown for the two systems (bottom), with the ligands in the bound (green sticks), unbound (red sticks), and transition states (gray sticks). Dashed lines depict key interactions.

Methods

All MD simulations were carried out in NAMD 2.12,[64] using the AMBER ff14SB force field for the protein,[65] and using the general Amber force field (GAFF) for the ligands.[66] The MD simulation setup is detailed in SI Section 1.

Unbinding Simulations

Our unbinding method is illustrated algorithmically in Figure . An explorational unbiased MD simulation of at least 20 ns was performed to identify the initial interactions between the protein and the ligand in the bound state. These initial simulations allow us to define the first set of CVs describing all distances between the heavy atoms of the ligand and the heavy atoms of the protein smaller than din = 3.5 Å, our interaction cutoff. The identified interactions will generate a single one-dimensional CV as the sum of these M distances, d, and will be used for iteratively biasing the simulations to observe an unbinding trajectory.

Figure 2

Flowchart illustrating the steps for the unbinding protocol.

At every iteration, we will define our bias as a harmonic restraint: , where D = D0 + (Mdtar). Here, we aim to reach the target value D for our 1D CV starting from the initial value at the beginning of the nth iteration D0. dtar is the incremental factor, set to 1 Å, representing the average increase we aim to achieve per distance for the next iteration. The targeted D value will be reached progressively within the next 10 ns long MD simulation for every iteration. The force constant, k, was set to 20 kcal/(mol Å2).

At the end of each iteration, the biased trajectory was analyzed, and novel interactions were identified, within din of the ligand, that are present for more than half of the total simulation time (i.e., 5 ns). These novel interactions are then added to the list of interactions that define the main CV for the next iteration. Additionally, we also re-evaluate existing interactions. If a distance during the last 5 ns of the trajectory exceeds dout = 6 Å or its variance exceeds dvar = 1 Å, then the distance is removed from the main CV in the next iteration. This exclusion factor will ensure that once a protein–ligand atom pair distance has exceeded dout, and therefore, there is no significant interaction between these atoms, this interaction is no longer biased. Similarly, loosely interacting atom pairs have higher distance fluctuations, and thus the corresponding weak interaction does not need to be included in the bias.

To reduce the number of interactions between the ligand and the protein and to remove redundancies, we combine atoms that are part of an equivalent group where a rotational degree of freedom can interconvert the atoms from one to the other (for example, benzene ring or carboxylic groups, Figure S1). Here, we considered the center of mass of that functional group and not the individual atoms.

The iterative process will end when no more distances are present in the main CV from the last iteration n; thus, there are no more stable interactions between the ligand and the protein, suggesting that the ligand is outside the binding pocket. Figure S2.I–VI represent the distances included in the unbinding trajectories.

Free Energy Calculations

Once the ligand is outside of the binding pocket, to determine the minimum free energy path for the unbinding trajectory, we use the finite-temperature string method.[46,68] The initial path and the full set of distances (CVs) are taken from the obtained unbinding trajectory. We extract these CV values for each interatomic distance along the initial unbinding path to construct the minimum free energy unbinding pathway iteratively, building a string of 100 windows in the coordinate space. For each window and each CV, we apply a position restrain equidistantly along the initial fitted string, using a force constant of 20 kcal/(mol Å2). We perform biased simulations using these restraints for a total time of 5 ns per window. From the obtained set of trajectories, a high-order (8) polynomial fitting is applied using the average values for each collective coordinate to build the subsequent set of refined CV positions. The procedure is carried out iteratively until the convergence of free energy profiles and the pathway. This is verified by ensuring that the maximal change of each CV between subsequent iterations is below 7% (or 0.3 Å) from the previous iteration. By adding multiple overlapping biasing potentials along the dissociation pathways, which are parameterized via the identified CVs, string simulations can sufficiently sample the high-dimensional path describing the full unbinding trajectory in detail. Finally, to obtain the corresponding potential of mean force (PMF), we unbias the simulations using the binless implementation[46] of the weighted histogram analysis method (WHAM).[69]

We note that our method does not aim to calculate binding free energies or kon rates. These would require simulations of a completely dissociated ligand and protein system, for which the string method is not an efficient algorithm. To this aim, routinely used efficient and accurate FEP calculations can be combined with our method to determine binding free energies and koff rates, respectively, from which the kon rates can be derived.[14,70]

Machine Learning Transition State Analysis (MLTSA)

We developed a machine learning transition state analysis (MLTSA) method to identify novel descriptors that determine the fate of a trajectory from the TS, which is applicable to unbinding simulations but also suitable for other applications as a low-dimensional feature selection method for highly complex processes where a TS region is identified. In our case, the novel molecular interactions between the drug molecule and the protein for unbinding provide key signatures that determine the unbinding kinetics.

To test the validity of the MLTSA, we created an analytical model and compared the ability of two ML approaches to detect correlated features: a multilayer perceptron (MLP) architecture NN model and gradient boosting decision trees (GBDT), a common ML approach in feature selection.

The analytical model was based on using multidimensional trajectories generated via a set of one-dimensional (1D) free energy potentials (SI, Section 5). Two types of potentials were used, both a set of single-well (SW) and double-well (DW) potentials. We used all but one of the DW potentials as “noise” and one of the DW potentials to define the outcome of the process, as the decisive coordinate to classify trajectories as “IN” or “OUT”. We generated trajectories using Langevin dynamics along 25 1D potentials. We used these trajectories to define 180 input features analogously to our observable CVs by computing linear combinations of the original coordinates (SI, Section 5). In our example, 11 of these 180 contained the selected DW potential with some nonzero coefficient (Tables S1.I and S1.II). We used these set of CVs to train the ML methods to predict the trajectory outcomes. Importantly, we aimed to identify the CVs that had the largest coefficients for our key selected DW potential.

We trained the MLP to analyze the model data sets of the downhill trajectories and predict their possible outcome from early on data, i.e., at 30–60 steps of the downhill trajectory for the analytical model. The training was performed using the Scikit-learn library.[71] We trained a simple model with an MLP classifier architecture, using three main layers (input, hidden, and output) with as many input nodes as input features depending on the system of study (for the analytical model 180 were used, for CDK2 see Table S2.II), fully connected to a hidden layer with 100 hidden neurons and ending in an output layer with one output node each for IN or OUT classifications. The model was optimized using the Adam solver[72] and using the ReLu[73] function as an activation function for the hidden layer. The training was done with a learning rate of 0.001, iterating over data until convergence or upon reaching the maximum number of iterations (500 epochs). Convergence is determined by the tolerance and the number of epochs with no change in loss. When there are 10 consecutive epochs with less than 0.0001 improvement on the given loss, the training stops, and convergence is reached. The same parameters were used for both the analytical model and CDK2 data.

We also tested the GBDT model using the Scikit-learn library as a comparison to the MLP approach. This method provides feature importances (FI) that enable the ranking and identification of relevant features. We trained 500 decision stumps as weak learners for GBDT minimizing a logistic loss function, with a learning rate of 0.1. The criterion for the quality of the splits was the Friedman Mean Squared Error (MSE), with a minimum of 2 samples to split an internal node, and a minimum of 1 sample to be at a leaf node. The maximum depth of the individual regression estimators was 3, without a limit on the maximum number of features to consider as the best split, without maximum on leaf nodes and using a validation fraction of 0.1. The same parameters were used for both the analytical model system and the CDK2 simulations.

The flowchart of the MLTSA method is illustrated in Figure S3. For the analytical model, we run 180 trajectories for the ML training and a separate validation set with 50 additional unseen trajectories. Following the flowchart, after labeling them as “IN” or “OUT” using the decisive coordinate, we created a data set for the ML algorithms containing 180 features per frame (Figure S4). We trained the ML models at different time frames (Figure S5) to observe the evolution of the accuracy throughout the simulations. The accuracy and number of epochs used in training are given in Table S3. This allows us to find a time range in the simulations where the classification problem is neither hard nor too trivial. Using this range, we trained the MLP model to analyze the importance of the features with our novel method. In a similar fashion to feature permutation,[74,75] or other model inspection techniques,[76−79] the MLTSA uses the Global Mean (GM) approach,[77] which swaps the value of each feature, one at a time with the mean value of the feature across all data used for training. This altered data set is used for prediction again expecting to get the same accuracy as the training on noncorrelated features and an accuracy drop on the correlated features, which depends on the level of correlation. For the comparison with GBDT and its FI, we trained the model at the same time and fetched the FI from the model to compare it with the accuracy drop analysis (Figure S6).

For the application of the MLTSA on CDK2, first, we identified the approximate TS location by selecting the last simulation frames from the highest energy five windows near the TS point of the obtained PMF. From each of these five starting coordinates, we then ran 50 independent unbiased MD simulations, each 5 ns long. We classified and labeled these short “downhill” trajectories by considering a combination of two key distances (Table S2.I) to identify which simulations finish either in a ligand-bound position (IN) or in a ligand unbound position (OUT). We then selected the starting structure (i.e., our TS) that provides the closest to a 1:1 ratio of IN and OUT events among these trajectories, and we ran 200 additional 5 ns long unbiased MD simulations with this starting point. We considered all interatomic distances (heavy atoms only) between the ligand and the protein within 6 Å at the TS starting position and determined the values of these distances along downhill trajectories (Table S2.II). These constitute a data set of distances for each simulation trajectory, and we aimed to select the most important features from these with our MLTSA method.

The number of epochs and convergence of the loss function for each model can be found in Tables S4.I–S4.II and Figure S7. Thus, using the frames coming from the multiple short, unbiased MD simulation trajectories starting from our TS, we provided a data set of distances extracted along the trajectory, as well as the future outcome of the IN or OUT events as the desired answer/classification. We performed the ML training at several different time ranges of the trajectories (Figure S8), to observe the predicted accuracy at different time ranges along the simulations. From all of the available trajectories for each system, we reserve a part for further validation to avoid the overfitting of our model. The rest is used for training, with all frames from the trajectories concatenated and randomly mixed, then split in different fractions as training (0.7) and test (0.3) sets. The trained model is additionally verified to have a similar prediction accuracy on the unseen trajectories.

Using our trained model, we assess which features are the most important for the model to predict whether the simulation is classified as bound (IN) or unbound (OUT). To do so, we apply our own feature reduction approach (FR), in which every single distance (i.e., feature) is excluded one-by-one from the analysis, and we calculate the drop in accuracy compared with the full set of distances present. Different from the standard approach,[79] where the real value of each excluded feature is replaced with a zero, here we replace the value for each excluded feature with the global mean of that selected feature across the simulations, thus canceling the variance of the aforementioned feature.

Results and Discussion

MLTSA Analytical Test and Validation

ML training on the model potential-derived trajectories was performed with both MLP and GBDT ML methods. We performed the MLP training at different time frames and trajectory lengths, from the 0th time step to the 500th step in intervals ranging from 10 to 150 frames at a time to assess the accuracy through time (Figure S5). Using a suitable time range consisting of the 30th–60th simulation steps from each trajectory, the trained ML methods found the classification problem accurately solvable but not too trivial. We replicated the complete process 100 times by generating 180 new independent simulations for each replica and performing the ML training. The MLP achieved an average test accuracy of over 94% and an average validation accuracy of over 93%, whereas the GBDT achieved over 99% on the test set and 91% on the validation set.

To identify the selected DW potential and its highest correlated features from the data set, we calculated the accuracy drop (MLTSA as in Figure S3) using the trained MLP and compared this approach to the FI using GBDT. Training accuracies for both ML models at 1 and 5 DW potentials can be found in Table S3. Results of both feature analysis methods are found in Figure for the 1 DW data set and in Figure S6 for the 5 DW potential data set.

Figure 3

Comparison between GBDT (top) and MLTSA with NN (bottom) feature analysis methods for the 1 DW data set. Correlated features are marked from blue (0%) to red (100%) depending on the mixing coefficient, α (× symbols, color scale on the right, five highest mixing coefficients are also displayed for the data points). Uncorrelated features (small black symbols) are at 0 FI for GBDT and show no loss of accuracy for MLTSA with MLPs. Correlated features all show a significant accuracy drop for the MLP, while only the top correlated features have high FI using GBDT.

The highest correlated features (colored depending on the correlation level, a color bar in Figure right panel) were correctly identified by both MLP and GBDT models. For GBDT, only the top three features show a high FI value (labels added to data points in Figure ), whereas the rest of the correlated features ranging from α∼34% up to ∼60% do not show a significant FI value. In addition, despite three features (#48, #89, and #136) having 40.34, 34.80, and 35.48% mixing coefficients, respectively, GBDT did not capture their correlation, showing values very close to 0. For the MLP, the top three distances are similarly captured as in the FI with the highest accuracy drops. Importantly, all correlated features have a nonzero accuracy drop, showing that they are correctly identified.

Using the data set with increased complexity consisting of 5 DW potentials and 15 correlated features, we observed a similar performance of the two ML methods (Figure S6). GBDT correctly captured and ranked the top three features (#8, #25, and #35). However, most other important features scored a FI value very close to 0. Out of 15 correlated features, GBDT did not identify 12 of them with high FI, whereas the MLP captured all of them. However, the MLP accuracy drop did not rank the top four features in the correct order, scoring the third most correlated feature with the biggest accuracy drop.

Considering both analytical models, we found that whereas GBDT has a higher specificity to rank the top correlated features in the correct order, MLP has a higher sensitivity and captures all correlated features but cannot necessarily identify the highest-ranked ones quantitatively using the accuracy drop as the measure. Therefore, a combination of the two ML methods can further help identify the most important features. In more complex systems, this performance might not be directly generalizable, however, due to the simple linear correlation of the CVs of this model.

CDK Kinase Unbinding Free Energy Calculations

For each system, we performed three independent simulation replicas starting from the respective equilibrated system. For each replica, we performed the initial unbiased MD simulation, followed by our unbinding trajectory determination procedure, and subsequently calculated the minimum free-energy path and the corresponding free energy profile using the finite-temperature string method (Figure ).

Figure 4

PMF of the unbinding path for 18K (a) and 62K (b). The free energy profile is obtained from a representative replica; the standard error, shown as a shaded area, is obtained by dividing the full data set into four subgroups.

Figure shows a representative result of the unbinding process for selected interactions. The first distance (blue line) is identified from the initial unbiased bound simulation as being shorter than 3.5 Å. Later during the biased unbinding process at 30 ns, a new interaction is found (orange line) and at 90, 120, and 130 ns, more distances are included in the main CV (green, red, purple, and brown). Additionally, interactions are progressively being removed as they are breaking (above 6 Å). Details of the selected CVs during the unbinding iterations are in the panels of Figure S2.I–VI for every replica.

Figure 5

(a) Unbinding trajectory of ligand 62K represented as selected snapshots along the trajectory at 0, 70, 141, and 219 ns from left to right, respectively. Representative distances used for the bias are shown as colored dashed lines (for the full set of distances, please refer to Figure S2.I–VI). Representative distances included in the CV along the unbinding trajectory are shown in (b) and the corresponding distance values plotted in (c). The lower dashed line at 3.5 Å is the cutoff value below which an interaction is included in the main CV; the upper cutoff at 6 Å is the value above which the distance is excluded from the CV.

Overall, while the identified CVs in different replicas vary, a few common key CVs are present in all unbinding trajectories within all replicas (Figure ). Even if the actual unbinding pathways have differences among the replicas, as seen by looking at the distances found along the paths (Figure S2.I–VI), they are all expected to pass through the same TS ensemble and show generally the same mechanism (see animated gifs for the final minimum free energy paths, SI, Section 11 and SI, Section 8—60K/4FKU system). This can also be confirmed from the consistent free energy profiles (Figure S9).

Figure 6

CVs obtained from the unbinding of 18K (a) and 62K (b); representative distances shown in dashed lines (yellow: interaction from the initial structure, cyan: interaction found during the unbinding trajectory); red sticks represent the coordinate of the ligand when it is outside the pocket. These distances appear in each of the three replicas for each system.

Additionally, we also performed the unbinding calculations for a third ligand, 60K, that is, analogous to 62K (Figure S10). Interestingly, we identified that all three replicate string pathways originating from three distinct unbinding simulations present a rotation of the hydrazineyl N=C bond, leading to a cis(Z)-trans(E) isomerization of the ligand near the TS (Figures S11 and S12). This is due to, on one hand, initial strong forces in string simulations that could be avoided in the future, and, on the other hand, to force field inaccuracies with a too low energy of the transform and too low barrier for the related dihedral angle rotation as determined by density functional theory (DFT) calculations (Figure S13). When compared with 62K, which does not exhibit this behavior in any of the three replicas, we can observe a lower energy for the 60K trans state that enables it to avoid the TS bottleneck. Correspondingly, all three distinct replicas result in a consistently too low unbinding free energy barrier when compared with the experiment (Figure S9). Animated trajectories along the string simulations for all replicas are provided and can be accessed through SI Section 11—additional resources.

The energy barrier extracted from the PMF of our simulations agrees closely with the experimental koff rates and is very well reproducible within the same system (Table and Figure S9). The shape of free energy profiles is also consistent among the replicas; however, the exact shape depends on the CVs identified in that replica (Figure S9 and Table S5). Generally, a higher number of CVs results in a broader TS region (e.g., Figure S9, ligand 62K). In addition, results for the third ligand, 60K, are also presented, demonstrating a consistent underestimation of the free energy barrier due to the discontinuity of the dihedral angle along the minimum free energy paths.[80]

Table 1

Ligand Binding Kinetic and Thermodynamic Values of 3sw4 and 4fkw Systems from Dunbar et al.[55] and Calculated Results Obtained from Our Simulationsa

PDB	ligand	KD (M)	kon [M–1 s–1]	koff [s–1]	ΔGexp‡ (kcal/mol)	ΔGcalc‡ (kcal/mol)
3sw4	18K	9.61 × 10–7	1.00 × 105	0.0823	18.93(±0.17)	16.29(±0.21)
4fkw	62K	4.73 × 10–8	6.49 × 104	0.00261	20.97(±0.05)	20.27(±1.06)

ΔGcalc‡ was calculated using the Eyring–Polanyi equation: k = kBT/h exp(−ΔG‡/kBT) at 298 K.[81]

Importantly, comparing the same ligand within the three different replicas in all systems provide very similar free energy barriers, expressed with a low standard error. Our energy barriers consistently reproduce high-energy barriers also seen experimentally thanks to the introduction of numerous key CVs that are not only taken from the initial ligand-bound conformation but, instead, introduced along the unbinding paths (Figure ).

We observe only one main barrier corresponding to the breaking of drugs with the His84 H-bonding contact (Figure ), suggesting that the different replicas do indeed share the same TS ensemble, despite the slightly varying pathways and identified CVs along the path. This H-bond was reported as a key interaction in many ligands in complex with CDK2/CDK5.[82,83] This interaction was included in the initial unbiased simulation in bound systems at the beginning of the unbinding procedure. However, during the unbinding trajectories, once this important H-bond between His84 and the ligand is broken, new interactions are formed, for varying time scales. For 18K, in all of the three replicas, H-bonds are formed with the exocyclic amino group of the ligand (N5) and the backbone oxygen of Glu81 and subsequently with the backbone oxygen of His84. 62K presents a sulfonamide terminal group, which, during the trajectory, interacts with Val163 and His84 of CDK2.

To analyze which distances are the most important in the TS region, we implemented our MLTSA method. Starting with two data sets of 139 (62K) and 148 (18K) independent downhill trajectories for each system, and an initial set of CVs of over 170 (Table S2.II), we obtained key distances for each system that are major determinants for the prediction of whether a molecule ends up in the bound or unbound states (Figure ). By training with trajectory data from up to 0.3 ns of each downhill simulation, the model can predict with high accuracy the IN or OUT outcome of the trajectories, more specifically: 80.11% for 18K and 93.83% for 62K. To confirm the effectiveness of the ML training, we compared the ML prediction accuracy using optimal thresholds of our main string CVs (Figure ) to determine the outcome at 5 ns of downhill simulations (Figure S14.I–II). Importantly, the ML model predicts the outcome more accurately at early times (before ∼0.3 ns) than using the best possible prediction via the string reaction coordinate: with above around 80–94% accuracy versus ∼55–61%, respectively, for the ML and the main CV (Figure S14.I,II).

Figure 7

Representation of the PMF of ligand 62K along the string coordinate and the path of multiple downhill trajectories started at the TS (in green) for further analysis. From the TS coordinate as a starting point, a set of simulations leading to both an IN position (blue) and an OUT position (red) are represented as lines. The green dots illustrate the free energy profile data points obtained from the WHAM calculation using the string window as string coordinate. The green line represents the fitting obtained from the green dots. The yellow shade represents the simulation time portion used for analysis during our machine learning-based approach.

Using the trained model, we then performed a feature reduction analysis to identify which CV features affect the overall prediction ability of the ML model the most. For both molecules, we were able to select the most important structural features (Figure ) that lead to the significant reduction of the prediction accuracy when such features were eliminated (these were kept as a constant value and fed to the ML, see Figure S3 for details), while other features did not affect the overall accuracy of the predictions.

Figure 8

Identification of the essential distances (feature reduction) from the largest accuracy drop using the last 50% (yellow), 25% (red), and 10% (blue) of the frames up to the first 0.3 ns of the simulations for (a) 18K and (c) 62K. The different shades in the background group the different features according to the atom of the ligand involved. Features presenting a significant decrease in accuracy are labeled (see Table S2.II) and portrayed as a three-dimensional (3D) representation on the right side of each plot: (b) 18K and (d) 62K.

We also compared the validity of the feature reduction approach with GBDT to identify FIs from the GBDT model. The results obtained show broad similarity with our main MLTSA approach (Figure S15.I,II), and they both outperform the baseline approach without ML. This suggests that alternative ML models may also be used successfully and further validate our results.

The MLTSA is significantly less computationally intensive than either the unbinding simulations or the string calculations. The short downhill trajectories can be trivially parallelized, which constitute the main cost of the MLTSA analysis. The ML training and accuracy drop calculations have a negligible cost compared with these; therefore, MLTSA could be a quick and effective approach to understand key CVs at the TS.

Conclusions

Optimizing ligand unbinding kinetics is a very challenging problem for small-molecule drug discovery and design that can lead to the development of drugs with superior efficacy. To tackle this, we have developed a new method, which allows us to calculate the free energy barrier for the ligand unbinding process, therefore providing quantitative information about the residence time of a specific ligand. Our method involves an exploration step, where a ligand unbinding path is determined together with key collective variables that describe this path. Subsequently, we performed accurate free energy calculations using the complete set of identified interactions as CVs along the unbinding path via the finite-temperature string method. This provides us with free energy barriers and an ensemble of structures at the transition state of the ligand unbinding process. The novelty of the method lies in the combination of automated iterative addition and removal of the collective variables determining an unbinding trajectory, which allows us to discover novel interactions not available a priori, based on the interactions from the bound structure. We found that while unbinding trajectories show different paths between different replicas for the same system, our method nevertheless identifies the key interactions important during the unbinding process and provides consistent free energy barriers. The combination of generating an initial path and identifying the important CVs for the unbinding process with the string method for accurate free energy calculations using high-dimensional reaction coordinates provide an efficient way to obtain quantitative kinetics of ligand unbinding.

We tested this method using a well-studied cancer drug target, CDK2, using two drug molecules with measured kinetic profiles. We obtained energy barriers in agreement with experiments using our method, which demonstrates the fundamental importance of determining a well-selected, high-dimensional set of CVs for the correct description of the process and kinetics results.

We explored analytical 180-dimensional systems using one or multiple DW potentials. We performed the ML analysis both with GBDT and MLP methods. Our results demonstrate for simple linear mixing models that they both can capture correctly the most important correlated features. The MLP is a faster approach and is more sensitive to correlated features; however, sometimes it could not rank the top features in their correct order. On the contrary, the GBDT feature importances could miss lowly correlated features in a data set but can more accurately rank the top features. The average training time using a single core was around ∼3.5 min/model to converge, whereas the GBDT training took about ∼5 min/model. Thus, we suggest that a joint approach with both models, which may complement one another, could be used to identify relevant CVs. Nonetheless, future studies with nonlinear correlated time series can further help to explore the performances of these and other ML methods. Importantly, analogous analysis can be performed for various complex processes, including ones with multiple states as possible outcomes.

To aid the kinetics-based design of novel compounds, we also developed a novel method, MLTSA, to identify the most important features involved at the TS of the unbinding. We generated multiple trajectories initiated at the TS, which either terminated in the bound state or in the unbound state. We then trained a multilayer perceptron ML algorithm to predict the outcome of the trajectories using a set of CVs and data drawn from the initial segment of the trajectories only. By doing so, we were able to demonstrate that the ML was able to predict the trajectory outcomes with much higher accuracy than using the original set of CVs used for the free energy calculations. A feature importance analysis was further employed to then identify the key CVs and the corresponding structural features that determined the fate of the trajectories, which therefore are the most important descriptors of the TS.

In addition to binding rates, we also aimed to identify specific molecular features and interactions with the target protein that allows us to design kinetic properties of the ligand. Using our ML methods, we identified multiple interactions between the protein and specific parts of the ligands that were of major importance for trajectories to cross the TS. Important protein–ligand interactions at the TS-bound poses for CDK2 correspond to functional groups of the distal ends of the ligands. Besides His84, a known key residue for interaction with multiple CDK2/4 inhibitors, here we also identified additional common interactions within CDK2 across the ligands, for example, between Lys89 and the sulfonamide groups or between Asp145 and the carboxylic group and the ester group for 62K, respectively. Importantly, to perform this analysis, we require the approximate knowledge of the TS structures as well as the MLTSA approach generating a set of downhill unbiased trajectories from these starting points. Our algorithms enable us to uncover novel design objectives for a kinetics-based lead optimization process.