Exploring ligand binding pathways on proteins using hypersound-accelerated molecular dynamics

Abstract

Capturing the dynamic processes of biomolecular systems in atomistic detail remains difficult despite recent experimental advances. Although molecular dynamics (MD) techniques enable atomic-level observations, simulations of "slow" biomolecular processes (with timescales longer than submilliseconds) are challenging because of current computer speed limitations. Therefore, we developed a method to accelerate MD simulations by high-frequency ultrasound perturbation. The binding events between the protein CDK2 and its small-molecule inhibitors were nearly undetectable in 100-ns conventional MD, but the method successfully accelerated their slow binding rates by up to 10-20 times. Hypersound-accelerated MD simulations revealed a variety of microscopic kinetic features of the inhibitors on the protein surface, such as the existence of different binding pathways to the active site. Moreover, the simulations allowed the estimation of the corresponding kinetic parameters and exploring other druggable pockets. This method can thus provide deeper insight into the microscopic interactions controlling biomolecular processes.

Fig. 1. Schematic illustration of the modeling of hypersound shock waves in MD simulations and hypersound-perturbed water dynamics at 298 K.

A (Top) Generation of hypersound shock waves in six directions (+X, +Y, +Z, −X, −Y, and −Z), propagating from the X₀, Y₀, Z₀, X₁, Y₁, and Z₁ surfaces, respectively, of the simulation box. (Bottom) Each shock wave consisted of 5 cycles and involved 80 (=16 × 5 cycles) velocity pulses (indicated by vertical bars) applied every N MD time steps (see “Methods” for details). B–D Spatial variation of B mass density, C pressure in the +X direction (p_x), and D X component of kinetic energy (k_x), measured at different simulation times. E, F Time dependence of E mass density and F the pressure, measured at different X positions; the corresponding positions are shown in (B) and (C). Shock waves were generated in the X = 0–1 nm region [X₀ surface in (A)] of the simulation box.

Fig. 2. Microscopic binding pathways of CDK2 inhibitors.

A, B Representative binding pathways of A CS3 and B CS242 ligands to the ATP-binding pocket of CDK2. (Top) Projections of binding conformations observed in the whole set of MD trajectories (colored dots) and of a representative binding pathway (black line) onto the first and second principal components (PC1 and PC2) calculated from principal component analysis (PCA). Ten (CS3) and 7 (CS242) representative binding poses (magenta sticks) on CDK2 (gray surfaces) are shown alongside the crystallographic pose (green sticks), the closest conformation to which was assigned as Pose 1. (Bottom) Potential energy (black) and free energy (red) trajectories corresponding to the pathway shown in the PCA map. The highest-energy transition state is indicated by a black (potential energy) or red (free energy) arrow. The upper panel shows an enlarged view of these trajectories close to the highest-energy transition state. Note that transition states occur A immediately before/after the ligand enters the CDK2 pocket and B during conformational rearrangements taking place after pocket entry. C Schematic illustration of microscopic and macroscopic kinetic models. The conventional kinetic model assumes a single binding pathway with a single transition state. However, at the single-molecule level, the ligand binds to the protein through multiple pathways with different highest-energy transition state conformations.

Fig. 3. Ligand-dependent binding site selectivity.

A–D Specific binding sites of A CS3, B CS242, C 2AN, and D 9YZ ligands on the CDK2 surface. The chemical structure of each ligand is shown at the top of the figure. After eliminating residues frequently accessed by all ligands, the backbone (ribbons) and side chains (thin sticks) of CDK2 were colored gradually by the interaction frequency of the individual ligands (e.g., residues with higher and lower frequencies are indicated by red and blue colors, respectively). Enlarged views of the ATP pocket, allosteric site 1, and allosteric site 2 (surfaces) are shown, with the crystallographic poses of CS3/CS242, 2AN, and 9YZ inhibitors (white sticks), respectively. The percentages shown in the models indicate the probabilities of capturing the ligand-binding event in the hypersound-perturbed simulations (Supplementary Table 4).