The microscopic observation of biomolecular processes such as protein folding, protein interactions, and enzyme reactions, most of which occur on timescales ranging from microseconds to seconds 1, is of great interest to the molecular biology community. Although molecular dynamics (MD) simulations enable atomic-level observations, they are limited to several microseconds on standard high-performance computers and are thus normally applicable only to relatively fast processes 2. Recently, the kinetics of slower protein interaction processes were explored through long MD simulations spanning timescales of tens of microseconds to milliseconds 3-6, which were achieved through the development of special-purpose supercomputers for high-speed simulations (e.g., ANTON 7) and/or algorithms to aggregate many short simulations (e.g., Markov state models (MSMs) 8). Unfortunately, MD-specific supercomputers such as ANTON are accessible only to a limited number of researchers owing to their high costs. While MSMs have a lower requirement for simulation power, this method is very sensitive to the choice of hyperparameters 9, which makes MSM approaches less than straightforward to use.
To overcome these problems, we have developed a new MD simulation method that utilizes high-frequency ultrasound (hereafter denoted as hypersound) shock waves to accelerate the dynamics. This method falls into the category of nonequilibrium MD simulations under external field perturbation 10, 11. Its key advantage is that it allows naturally “slow” processes such as those mentioned above to be frequently and directly captured in a series of single MD trajectories performed on standard high-performance computers. In the experimental field, ultrasound irradiation procedures have been applied to accelerate various kinds of chemical reactions 12, 13 and to synthesize nanoparticles 14. This acceleration is considered to be induced by acoustic cavitation (i.e., the repeated growth and collapse of cavitation bubbles formed by the ultrasound waves, which generate local high-temperature/pressure regions in solution) 15. Inspired by these results, in this study, we first analyzed the hypersound-dependent behavior of a liquid water model to test the effect of shock waves with a protein-size wavelength. Next, to assess the effect of hypersound acceleration on biomolecular processes, we performed short (100–200 ns) simulations to capture the slow binding of small-molecule inhibitor compounds (CS3 and CS242) to cyclin-dependent kinase 2 (CDK2) 16, as a representative system in which the binding event would be nearly undetectable in standard MD. The simulations showed a significant acceleration of the binding process under hypersound irradiation compared to standard MD simulations. The accelerated simulations revealed the existence of various conformationally and energetically diverse binding pathways, suggesting that the assumption of a single pathway/transition state made in conventional kinetic models may be inaccurate. The present method allowed not only the estimation of kinetic parameters of slow-binding inhibitors, but also the full exploration of druggable sites. This approach would thus be helpful for efficiently understanding the microscopic mechanism of slow biomolecular processes.
To simulate hypersound shock waves with protein-size wavelengths, their frequency was set to 625 GHz (corresponding to a period of 1.6 ps ) (Fig. 1A), which is more than 100 times higher than that of currently used ultrasound waves 17. The hypersound-perturbed MD simulation of liquid water at 298 K showed the generation and propagation of a high-density region (Fig. 1B and Supplementary Movie 1). Then, we analyzed the MD trajectory focusing on wave propagation along the X direction as a representative example (Fig. 1C–F). As the X coordinate of the first wave reached 4 nm at a simulation time of 1.7 ps after passing through X = 2 nm at 0.7 ps (Fig. 1E), the propagation speed of the shock wave could be estimated to be 2,000 m/s, which is similar to the speed of sound in water (~ 1500 m s-1) 18. The wavelength of the simulated shock wave was estimated to be 3.2 nm (2000 m s-1 × 1.6 ps), which corresponds to the hydrodynamic radii of globular proteins consisting of 300–400 residues 19. This confirms the successful generation of hypersound waves, which is appropriate to perturb biomolecular processes. The pressure (px) and kinetic energy (kx) of the liquid water model also exhibited periodic fluctuations with the same phase as the density, reaching ~ 2000 atmospheres and 0.4–0.5 kcal/mol (corresponding to 400–500 K) at the center (X = 4 nm) of the simulation box (Figs. 1C, D, and F). In contrast, the macroscopic properties of liquid water were not affected by the hypersound irradiation, except for the diffusion constant, which slightly increased to 6.30 ± 0.10 × 10–5 cm2 s-1, equivalent to the corresponding parameter of bulk water at 305 K (Supplementary Table 1). These results demonstrate that hypersound irradiation of a liquid solvent generates local higher pressure and temperature regions appropriate for promoting chemical processes 15 without altering the macroscopic properties of the liquid.
We next assessed the influence of shock wave perturbation on the association between the CDK2 protein and its slow-binding ATP-competitive inhibitors CS3 and CS242. The probability of observing the ligand-binding event in conventional MD simulations of 100 ns was estimated to be only 0.7% (2/283, corresponding to 2 out of 283 MD runs resulting in binding) for CS3 and 0.5% (2/369) for CS242 (Supplementary Table 2). On the other hand, higher probabilities were attained in 100-ns long hypersound-perturbed MD simulations. Hypersound irradiation with a higher amplitude or frequency resulted in an increase in the probability of observing ligand binding without a significant decrease in computation speed (Supplementary Table 3). When using a parameter set of (N=50 and vmax=400m/s), the CS3 and CS242 binding probabilities were 12.4% (22/177) and 4.8% (11/227), respectively (Supplementary Table 2), showing that the perturbation successfully increased the association rate by 17.7 times (12.4/0.7%) for CS3 and 9.6 times (4.8/0.5%) for CS242. MD trajectories that exhibited stable ligand binding were extended to 200 ns to observe the behavior of the bound ligand, and based on their percentages, the association rate constants (kon) under hypersound irradiation, using the same parameters as above, were estimated to be 3.68 × 106 for CS3 and 1.92 × 106 M-1 s-1 for CS242 (Table 1). This analysis proved the effectiveness of the hypersound-perturbed simulations in enhancing the sampling of infrequent binding events; this approach can thus be applied to extract further atomic-level information on these processes, as follows.
Accelerated MD simulations revealed that multiple transitions between different conformations took place within each individual binding pathway (see Fig. 2A and Supplementary Movie 2 for CS3 and Fig. 2B and Supplementary Movie 3 for CS242) This emerges from the inspection of the 67 (CS3) and 14 (CS242) binding pathways observed in the hypersound-perturbed MD simulations, a few representative cases of which are shown in Supplementary Figs. 1 and 2. It should be noted that these pathways contain those observed in conventional MD simulations (Supplementary Fig. 3). The potential energy trajectories (also displayed in the figures) reveal the occurrence of multiple energy barriers along each binding pathway, and show that the position and height of the highest-energy transition state depend on the binding pathway (Fig. 2C). The trajectories indicate that the ligand tends to adopt energetically unstable configurations upon (i) entry into the CDK2 pocket (Fig. 2A, and Supplementary Figs. 1A and 2A) or (ii) conformational rearrangement in the pocket interior (Fig. 2B, and Supplementary Figs. 1B and 2B). These effects have not been previously captured by ensemble-averaged kinetic experiments 16, 20 or existing generalized-ensemble MD simulations (Supplementary Fig. 3) 21, which predict a plausible pathway by efficiently exploring the conformational space. Ligand unbinding was also observed in some of these trajectories, most of which also exhibited different binding and unbinding pathways (Supplementary Figs. 1C and 2C). This suggests that the conventional kinetic model based on identical binding/unbinding pathways is not always valid at the single-molecule level. The trajectories of individual ligand molecules captured by the hypersound perturbation approach revealed the complex microscopic nature of the CDK2-inhibitor binding kinetics, highlighting the effectiveness of this approach in exposing effects not accessible by other experimental and computational techniques.
By averaging the energy barriers observed in the 9 (CS3) and 6 (CS242) trajectories that exhibited stable ligand binding under hypersound irradiation with the parameters predominantly used in the simulations (Supplementary Table 2), the activation energies for CS3 and CS242 binding to CDK2 were estimated to be 3.9 ± 1.8 and 6.7 ± 2.4 kcal/mol, respectively (p=0.02, one-sided Student’s t-test, Table 1), consistent with the relative height of the energy barriers estimated from the free energy landscape (see the Supplementary section 1.7), suggesting that the slower CDK2 association rate of CS242 than CS3 can be attributed to a higher energy barrier. The calculated Arrhenius parameters describing the kon dependence on the temperature (see Supplementary section 1.7) indicate that hypersound irradiation increased both the frequency factor (i.e., from 2.2 × 108 to 8.1 × 108 M-1 s-1 for CS3 and from 2.3 × 109 to 3.6 × 109 M-1 s-1 for CS242) and the effective temperature (from 298 to 362 K (CS3) and 445 K (CS242), see Table 1). The increase in the frequency factor could be attributed to enhanced ligand diffusion (Table 1) 17, which would increase the collision frequency of the ligand with the protein pocket, while enhancing the thermal motions of the solvent molecules did not result in increased ligand diffusion and an acceleration of the ligand binding process (the Supplementary section 1.8). In addition, the generation of local high-energy/pressure regions in the solvent leads to an increase in the effective temperature 15; however, the “macroscopic” temperature of the system remained unchanged at 298 K (Supplementary Table 1). As shown in Supplementary Fig. 4, the native interactions in the CDK2 structure were stably maintained during a series of hypersound-perturbed simulations with different frequencies and amplitudes of shock waves, confirming that the local high-energy regions do not induce thermal denaturation of CDK2. These results suggest that hypersound perturbation accelerates the protein-ligand association process by enhancing the cooperative local motions of the solvent molecules without affecting the native structure of the biomolecules, highlighting the general applicability of this approach for the acceleration of molecular processes in solution.
Table 1 Kinetic parameters of the CDK2-ligand binding process determined by conventional and hypersound-perturbed MD simulations. See also the Supplementary section 1.7. The values are presented as the mean ± SD.
|
|
Association rate constant
kon (M-1 s-1)
|
Activation energy
E (kcal mol-1)
|
Diffusion constant
D (×10-5 cm2 s-1)
|
Steric factor
log(P)
|
Frequency factor
log(A (M-1 s-1))
|
Effective temperature
T (K)
|
|
CS3
(conventional MD)
|
3.35 × 105 a
|
3.9 ± 1.8
|
0.17 ± 0.05
|
‒0.76
|
8.35 ± 0.13
|
298
|
|
CS3
(hypersound MD) b
|
3.68 × 106
|
3.9 ± 1.8
|
0.61± 0.16
|
‒0.76
|
8.91 ± 0.12
|
362
|
|
CS242
(conventional MD)
|
3.21 × 104 a
|
6.7 ± 2.4
|
0.19 ± 0.07
|
0.20
|
9.36 ± 0.17
|
298
|
|
CS242
(hypersound MD) b
|
1.92 × 106
|
6.7 ± 2.4
|
0.30 ± 0.10
|
0.20
|
9.56 ± 0.15
|
445
|
a Experimentally determined kon values retrieved from the Community Structure-Activity Resource (CSAR) database 16.
b Hypersound shock waves with N = 50 steps, vmax = 400 m/s, and Tint = 2,400 N were used to determine the kon, D, A, and T parameters, which vary depending on the hypersound parameters (Supplementary Table 3).
The identification of previously undiscovered druggable pockets on the protein surface plays a key role in expanding the therapeutic target range of the protein 22. The above hypersound-perturbed binding simulations allowed us to explore the properties of other sites beyond the ATP-binding pocket on the CDK2 surface. We analyzed the binding simulation data of the ATP-competitive inhibitors CS3 and CS242 and two allosteric inhibitors, 2AN and 9YZ, whose binding sites are distinct from the ATP-binding pocket 23 24. Hypersound irradiation accelerated the binding of all ligands to both the ATP-binding site and the two allosteric sites 1 and 2 (Fig. 3 and Supplementary Table 4). Allosteric site 2 was frequently accessed by all ligands (Fig. 3A-D), suggesting that this site is remarkably nonspecific because of its shallow pocket shape 24. In contrast, the CS3/CS242 (Fig. 3A and B) and 2AN (Fig. 3C) ligands showed a more specific preference for the ATP-binding pocket and allosteric site 1, respectively, compared to their nonspecific association with allosteric site 2, supporting the suggestion that these ligands prefer to associate with individual binding sites, as observed in their cocrystal structures 16, 23. The analysis of specific and nonspecific sites based on binding simulations of multiple ligands may thus be useful for the prediction of ligand-dependent binding site selectivity and for the exploration of novel druggable cryptic sites that can allosterically regulate enzymatic activity 22.
This study shows that hypersound-stimulated MD simulations have the potential to accelerate protein-ligand binding kinetics through a solvent-mediated mechanism without collapse of the protein structure, thus enabling atomic-scale observation of ligand-binding processes within time scales accessible by standard MD (~ 100 ns). In contrast to other advanced MD methods that accelerate biomolecular processes 25, 26, this method does not require prior knowledge of the protein-ligand complex structure. In this way, the simulations can provide significant insights into fundamental biological mechanisms (such as the discovery of microscopic ligand binding pathways involving various bound conformations) and facilitate drug discovery (as illustrated by the present exploration of druggable binding sites on the protein surface). The present acceleration method code is publicly available (see the Code Availability section), can be implemented on standard high-performance computers, and is suitable for parallel computing because performing multiple independent simulations in parallel enables the sampling of a higher number of binding events and thus produces an improved statistical description of the process under study. Furthermore, the hypersound irradiation modeled in the simulations is not a fictitious computational procedure, but a real physical process, even though hypersound waves of molecular (several nanometers) wavelength have not yet been realized 17, which currently hampers the experimental assessment of its impact. Further applications of the present technique to model other biomolecular (e.g., protein conformational changes and protein-protein interactions) and non-biomolecular (e.g., phase transitions of materials) processes are required to assess its general effectiveness in modeling slow dynamic events.