Energy Barriers , Molecular Motions , and Residence Time in Ligand Dissociation : A Computational Study on Type II Inhibitors Binding to CDK 8 / CycC

It is important but unclear how protein-ligand system motions affect free energy profiles, create energy barriers, and lead to slow residence time. We computed residence time (RT) and potential of mean force (PMF) of the dissociations of five structure-kinetic relationship (SKRs) series inhibitors of CDK8/CycC using metadynamics and milestoning theory. By using a novel way to represent the reaction coordinate based on principal component analysis (PCA), we found one-to-one mapping of hydrogen bond (H-bond) breakage and protein domain motions with the energy barriers in the PMFs. This work provides a novel and well-defined approach to study the dissociation kinetics of large and flexible systems. peer-reviewed) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprint (which was not . http://dx.doi.org/10.1101/169607 doi: bioRxiv preprint first posted online Jul. 28, 2017;

2][3][4] koff and residence time (RT) are particularly important as determinants of the efficacy of a drug candidate. 5It is desirable to optimize these properties for drug discovery but, binding kinetics and its determinants in the ligand-receptor structures are not yet fully understood.
It is also important to predict these quantities computationally when experimental measurements are not available.Therefore, computational methods, which are able to provide an atomistic description of temporal and spatial dissociation pathway details of ligand-receptor, can be employed to unravel the secret behind the RT and kinetics.
Cyclin-dependent kinase 8 (CDK8) is a kinase protein from the CDK family.7] Abnormal activities of CDK8 and its partner CycC are related to a variety of human cancers. 8CDK8 has an allosteric binding site adjacent to the ATP binding controlled by a DMG motif (Asp-Met-Gly) that characterizes the DFG-in/DFG-out conformations as in other protein kinases. 91] A structure-kinetic relationship (SKR) series compounds were developed and targeted at the allosteric binding site of CDK8. 12Up to date, a few computational works were performed to study CDK8.
Molecular dynamics (MD) simulation was used to study the structural stability of CDK8/CycC complex. 13The RT of SKR candidate compounds series 12 has been studied by metadynamics stimulations 14 using reweighting technique, but the timescale is fast on nanosecond level. 15lestoning theory is useful method for studying kinetics. 16Using rates and transition time from unbiased samplings between milestones dividing the transitions, milestoning computes the lifetime of each milestone and the kernel matrix which includes the transition probabilities between milestones.The potential of mean force (PMF) and the time the system takes to travel between milestones can be computed from the kernel matrix and the lifetime.Milestoning has witnessed a few applications and succeeded in reproducing the free energy surface and transition time. 17How to describe the reaction coordinate using collective variables (CVs) faithfully is a challenging problem for applying methods including milestoning, metadynamics, umbrella sampling 18 and other methods.It is critical to include degrees of freedom involved in the reaction coordinate as more as possible.With this concern in mind, it is nature to solve this issue by using principal component analysis (PCA), a mathematical method that extracts the major motions from a collection of data.Therefore, it is desirable to take the advantage of PCA to suggest a more straightforward but robust way to define reaction coordinate for calculation using milestoning theory.
In this work, we used metadynamics to provide dissociation pathways and computed PMF and RT for five ligands from the SKR compounds series 12  In our metadynamics trajectories, ligands completely dissociate at about 5 to 30 ns.Instead of reweighting the kinetics from the metadynamics trajectories, we used milestoning theory to more accurately evaluate RTs.RTs were computed as the time it takes to bring the ligand from the bound state milestone to the solvated state milestone (where the gray horizontal lines intersect with the PMFs in Figure 1).The computed RTs are listed in Table 1.The relative ranking distinguishes ligand1 and 2 that have much slower RTs than the rest ligands.
Among the ligands, ligand 1 has a RT on the millisecond level which reproduces the experimental value the best.The RTs of other ligands are on the micro-or submicro-level.
Interestingly, the computed RT of ligand 2, whose experimental RT is slightly slower, is hundred times faster than ligand 1.We noticed that ligand 2 forms half less H-bonds with CDK8 and requires minor protein motions for it to dissociate compared to ligand 1.
Therefore, the faster kinetics may be an intrinsic property underlying in the ligand 2 dissociation pathway found by metadynamics.The same reversed rank of ligand 1 and 2 was reported previously. 15e PMF plots along milestones describing ligand dissociations was computed by using milestoning theory and are shown in Figure 1.In the figure, the binding affinities and the energy barriers are labeled correspondingly.The computed binding affinities agree with experimental values reasonably well with deviations smaller than 3 kcal/mol.The PMF of ligand 1 suggests a binding affinity of roughly -12 kcal/mol, which is comparable to the experimental value -10 kcal/mol (Table SI 1).Interestingly, the computed energy barrier is higher than the experimental value but the computed RT is faster.The energy barriers shown in PMF correspond to observable conformational changes (Figure 2) and trespassing the energy barriers requires significant amount of time (Table 2).β1-2 on the N-lobe of CDK8 moves up for one or two Å to open the ATP binding site, and the R-group of the ligand 1 is unlocked (L1-A).Then, the R-group of the ligand moves up towards the N-lobe of CDK8 (L1-B).The huge energy barrier C comes when the R-group of the ligand shrinks into the ATP by passing β1-2 binding site (L1-C).This energy barrier takes 15 ns to cross, and is slowest in all observed energy barriers.The energy barrier D comes from the motion when the R-group of ligand 1 passes LYS52.Between D and E, the R-group of ligand 1 finds a local comfortable position between β1-3 and αC helix.The bulky part of the ligand passes through the narrow tunnel between αC helix and activation loop and results in the energy barrier E. Finally, the ligand gets solvated on the surface of CDK8 (L1-F and L1-G).
The motions and interaction indicated by the PMF can be further quantified by measuring distances, angles and H-bonds along the dissociation milestones.We computed the distances between centers of mass of heavy atoms of β1-2 and β8, αC helix and β8, C-lobe of CDK8 and CycC, and angle between centers of mass of heavy atoms of N-Lobe, hinge region and C-lobe of CDK8 for ligand 1 system, and projected these profiles onto the corresponding path and milestones (Figure SI 2).We also computed the H-bond between ligand 1 and CDK8, and projected the existence of H-bond onto the milestones (Figure SI 3).These plots correlate with the energy barriers and system motions perfectly.At energy barrier A, the opening of the CDK8 binding site is indicated by the β1-2/β8 and αC helix/β8 distances (Figure SI  The ligand 2 binding affinity predicted by milestoning is roughly -10 kcal/mol.Although the computed binding affinity of ligand 2 is similar to experimental value (Table SI 1), the computed RT is orders faster.Unlike Ligand 1, the PMF of ligand 2 only suggests three major relevant barriers (Figures 1 and 2); the R-group of ligand 2 passes through α-C helix (L2-A); ligand 2 completely moves out of the allosteric binding site (L2-B); ligand 2 gets solvated and desolvated (L2-C, D, E, F).Moreover, ligand 2 starts to encounter a huge energy barrier (about 3 kcal/mol) at milestone 5, while ligand 1 starts to encounter the first major energy barrier at milestone 20.This is because the stretching motion and preparation for dissociation of ligand 1 can be clearly distinguished by PC1/PC2 space path, but for ligand 2 the motions in the first 20 ns metadynamics simulations are dominated by protein and ligand vibrations near the same bound state so that the track of these regions in PC1/PC2 space are convoluted and we cannot separate the motions using our approach (Figure SI 4).The PMF suggests that three major events are involved in the dissociation of ligand 2, and this is less complicated than ligand 1.The PMFs of the other ligands indicate roughly -6 kcal/mol binding affinities for ligand 5, 10 and 11, which are 2 to 3 kcal/mol shallower than the experimental values (Table SI 1).Similar to ligand 2, the tracks in PC1/PC2 space of these three ligands are also convoluted in the first 10-20 ns.Therefore, the dynamics included in these periods of simulations are projected onto only a few milestones (1-5), and the PMFs of these ligands also start to encounter the major energy barrier within 10 milestones.The first distinguishable energy barriers for these ligands are also from the motion when the ligand passes through αC helix.The other major energy barriers in the PMFs of these ligands are from desolvation of the ligands.
The average lifetime of the milestones for the systems are about 1 to 7 ps.In our calculation, the time resolution is 0.1 ps, and this means the timestep resolution is able to capture the transition between the milestones.The lifetime of the milestones is affected by two factors.
First, although our milestones in the PC1/PC2 space is evenly distributed along the path, the areas divided into each milestone are not exactly identical.The first 20 milestones for ligand 1 system are located close to each other, so that the life time of these milestones can be as low as 0.5 ps (Figure SI 1 and Figure SI 5).The other factor is the nature of the population in that milestone.The lifetime of the milestones corresponding to the sloped regions near the energy barriers of the PMF is generally shorter than those corresponding to the flat regions on the PMF (Figure SI 5).The computed RT of the systems are several orders faster than experimental values and this is similar to other works.A millisecondlevel substrate translocation study using milestoning showed the energy barrier was 4-5 kcal/mol and transition time was 1.6 μs. 19Another milestoning study of CO entry and exit in myoglobin also suggested minimal energy barrier path with energy barrier and transition time on the same scale. 20This discrepancy implies that MD simulations tend to overestimate the rate of slow motions.
Although we used distance of centers of mass, a naïve way to define the CVs, in metadynamics, we computed PMFs and RTs using a novel approach.This approach extracts the most remarkable motions from the metadynamics trajectory by choosing the first two PC modes, provides a 2D projection of the complex motions of the diffusion of ligand and ligand/protein internal motions.The path combines related motions from multiple degrees of freedom, not limiting to a few CVs like distance or angles, although the ligand-protein relative motions dominate in these degrees of freedom.The advantage of this approach is that it reduces the dimension of the spatial motions, but keeps the most important dimensions of the motions.This approach can be further improved.First, the 2D projection of the smoothed track must be able to suggest a clear path without convoluted loops in the important regions of the dissociation.For example, as shown in Figure SI 6, the conformations involved with dissociation of ligand 11 is projected onto the PC1/PC2 space as a folded curve indicated by the black arrow, and this region cannot be effectively distinguished by the 2D projection.For the same reason, the early period of the metadynamics trajectories for ligand 2, 5, 10 and 11 are more or less projected onto a small convoluted region on the PC1/PC2 space (Figure SI 4), and this partly results in the overestimated dissociation kinetics and underestimated PMF for these ligands.An optimization of the metadynamcis path using methods like string method 21 may further improve the performance by removing the convoluted regions in the PC1/PC2 space.Also, only the first two PC mode may not be able to describe all major motions in the system.As shown in Figure SI 7, PC1 and PC2 captures 85% (sum of eigenvalue of PC1 and PC2) of the motions in ligand 1 complex, and roughly 75% for the other ligand complexes.This implies that this approach can be expanded to 3-dimensional cases in which the first three PC modes can be used to generate a 3D PC space to include slightly more major motions of the system.A 3D curve in the space can be used to define milestones that resemble a series of Domino along a thread in the same way.Although the path will be 3D, but the milestones will still be one dimensional as in the 2D space.
We used metadynamics and milestoning theory to study the dissociations of five SKR series compound from CDK8/CycC complex and suggested a novel way based on PCA to define the reaction coordinate of the dissociation.By using this novel reaction coordinate definition, we yielded RTs on milli-or micro-second level, and explained the origins of energy barriers indicated on PMF by using H-bonds and protein motions that were projected to the same milestone coordinate.peer-reviewed) is the author/funder.All rights reserved.No reuse allowed without permission.

Molecular Dynamics
We obtained the initial structures of CDK8/CycC-ligand 1, 2, 5, and 11 complexes from PDB database (PDB ID: 4F6W, 4F7L, 4F6U, 4F7N 1 ), and manually modified PDB structure of CDK8/CycC-ligand 11 (4F7N) to get the initial structure of CDK8/CycCligand 10.We built the missing activation loop of the structures by using Swiss Model 2 based on p38 DFG-out crystals structure (PDB ID: 1W82 3 ) as a template.The other missing residues of CDK8 residue 1 to 360 and CycC residue 361 to 619 were also added by using Swiss-Model.We determined the protonation states of histidine residues in the CDK8/CycC-ligand complexes by using MCCE 4 , and all histidine residues were in the single protonated state.Then we added missing hydrogen atoms to the structures by using tleap module from AmberTool14. 5BER FF14SB force field and GAFF 6 were used for CDK8/CycC complex and the ligands.We performed minimizations on the hydrogen atoms, sidechain atoms, and the whole structures to relax the complex conformations.Then we solvated the five complexes with TIP3P and a water buffer size of 12 Å and added 6 Cl -ions to neutralize the formal charges of the system.After we minimized the water and the solvated systems, we equilibrated the water molecules with the solute fixed at 298K for 20 ns in NPT (isothermal-isobaric) ensemble.Then the systems were heated from 200 K to 298K gradually.Production runs of the five systems were performed at 298K for 500 ns in NPT ensemble and saved every 2 ps with a 2-fs time step.All MD simulations were performed using pmemd.cudafrom AMBER 14 package. 5The trajectories were resaved every 20 ps for post-analysis.

Metadynamics
We obtained the equilibrated conformations at 298K, and conformations after 100, 400, 500 ns MD simulations as initial structures, and performed twelve metadynamics simulations starting from these initial structures for each of the CDK8/CycC-ligand complexes.The distance between centers of mass of the heavy atoms on the ligands and the heavy atoms of residue 26 to 38 96 to 106 356 to 359 on CDK8 were chosen as the collective variables (CVs), and the forces from the Gaussian functions were only applied on the atoms in the ligands.A Gaussian of 0.02 kcal/mol was deposited on to the CV coordinate every 0.1 ps.A 1-fs time step were used in the simulations, and the metadynamics run were performed until the ligand dissociation event occurs (roughly 25-30 ns).The metadynamics simulations were performed using NAMD. 7The metadynamics simulations were saved every 5 ps as initial conformations for milestoning calculations.The sum of percentages of eigenvalues of PC1 and PC2 are labeled.
using milestoning theory with a novel way for defining the milestones.The details of methods are included in SI.Using PC1 and PC2 of the metadynamics trajectories, we unambiguously defined the dissociation pathways and set up milestoning calculations using this definition (Figure SI 1).The computed PMF plots indicate clear energy barriers that correspond to and originate from the hydrogen bond (H-bond) formations and important interactions between ligand and CDK8, and also protein motions featured by relative motions of N-lobe, C-lobe of CDK8 and CycC.The computed RTs are on the time scale of milliseconds and microseconds and correlate with the compute PMF reasonably well.
2 A and B).At energy barrier C, ligand 1 starts to move through the ATP and allosteric binding sites, where the N-lobe of CDK8 needs to move away to open up the binding site again (Figure SI 2 A and B), and H-bonds between ligand 1 and Val27, Arg29, Asp103, Ala155, Asn156 and Arg356 break right at this energy barrier (Figure SI 3).As ligand 1 moves further toward allosteric dissociation, H-bonds between ligand 1 and Glu66 and Asp173, which lock ligand 1 in the binding site, breaks right at energy barrier D and E (Figure SI 3).Meantime, the β1-2/β8 distance starts to restore to native state value, as ligand 1 has already leave this region (Figure SI 2 A).On the contrary, the αC helix/β8 and C-lobe/CycC distances do not restore until ligand 1 completely dissociates at energy barrier F (Figure SI 2 B and C).This is also the moment when all bound state and intermediate state H-bonds break (Figure SI 3).These observations indicate that the ligandprotein interactions and protein domain motions contribute specifically to the energy barriers shown in the PMF.

Figure 1 .
Figure 1.PMF plots of the five ligand dissociations along milestones.Important energy barriers are label by A-G.The desolvated states used for computing binding affinities and RTs are indicated by the gray horizontal line.Binding affinities are labeled next to the black double arrows.

Figure 2 .
Figure 2. Conformational changes of CDK8/CycC-ligand 1 and 2 complexes along milestones.Ligands are shown in licorice and CDK8/CycC is shown in new cartoon.Indices of milestones are labeled on left top corner.First, second and third (if any) conformations are shown in gray, cyan and yellow.Conformational changes are indicated by green arrows.

Figure SI 1 .
Figure SI 1.The projection of metadynamics trajectory of ligand 1 onto PC1/PC2 space and definition of milestones.A: the tracks of original trajectory (dark and light blue dots) and smoothed trajectory (red), and the manually defined path (green dots).B: the optimized path (purple) and milestones (black lines).C: the path and milestone (black lines) and all projection of short replicas of unbiased MD simulations used for milestoning theory (the color areas).D: the area of projections in C within 10 from the path (black line) used for computing PMF and RT.

Figure SI 5 .Figure SI 6 .Figure SI 7 .
Figure SI 5.The PMF and lifetime of milestones.PMF is shown as black dotted line and the lifetime is shown as blue dotted line.

Table 1 .
The experimental and computed residence time (RT) of the five ligands.

Table 2 .
The time required to cross energy barriers for ligand 1 and ligand 2. The time required to cross one energy barrier is in column 3, and the time required to travel from bound state to pass one energy barrier is in column 4.
The experimental data of CDK8/CycC-ligand complexes taken from ref1.peer-reviewed) is the author/funder.All rights reserved.No reuse allowed without permission.