Adaptive Multilevel Splitting Method for Molecular Dynamics Calculation of Benzamidine-Trypsin Dissociation Time

Abstract

Adaptive multilevel splitting (AMS) is a rare event sampling method that requires minimal parameter tuning and allows unbiased sampling of transition pathways of a given rare event. Previous simulation studies have verified the efficiency and accuracy of AMS in the calculation of transition times for simple systems in both Monte Carlo and molecular dynamics (MD) simulations. Now, AMS is applied for the first time to an MD simulation of protein-ligand dissociation, representing a leap in complexity from the previous test cases. Of interest is the dissociation rate, which is typically too low to be accessible to conventional MD. The present study joins other recent efforts to develop advanced sampling techniques in MD to calculate dissociation rates, which are gaining importance in the pharmaceutical field as indicators of drug efficacy. The system investigated here, benzamidine bound to trypsin, is an example common to many of these efforts. The AMS estimate of the dissociation rate was found to be (2.6 ± 2.4) × 10(2) s(-1), which compares well with the experimental value.

Figure 1

Schematic of basic AMS algorithm for a dissociation process of a ligand in an initially bound state from a binding site (blue-green). In this case, N = 3 replicas are used and the reaction coordinate is defined as the radius about the initial state. In (a), an initial trajectory segment is generated for each replica. The suprema of the segments are compared, and the replica with the segment of lowest supremum (red) is killed, as shown in (b). Subsequently in (c), a surviving replica is randomly picked (green in this case); its trajectory segment up to the supremum of the killed replica is cloned into the killed replica and simulation is restarted until the trajectory returns to the initial state. Once again, the replica which has the least progress along the reaction coordinate (blue) is identified and killed, as shown in (d). The process is repeated until all replicas have surpassed z_max (not shown).

Figure 2

To avoid being stuck in small trajectory loop segments near the starting point, the first step of the original algorithm, shown in (a), is altered to incorporate a starting point z_min, shown in (b), a small distance away from z₀. Replica trajectory segments now begin at z_min, but still terminate at z₀.

Figure 3

Trypsin (light blue) and benzamidine in a randomly picked equilibration frame at reaction coordinate z_eq (red), AMS initial (yellow) bound states, and AMS final unbound state (purple). The reaction coordinate z is defined as the center-of-mass distance between non-hydrogen benzamidine atoms and the Cα atoms of 16 residues near the binding site (blue spheres). AMS initial and final state structures were extracted from trajectory frames during AMS simulation.

Figure 4

(a) Normalized sample distributions of z values within initial (green) and other (red) metastable portions of the trajectory during equilibrium simulation. The initial conditions of the AMS simulation are defined by z₀ and z_min levels, indicated by dashed lines. The inset shows the time evolution of the z value. (b) Simulation trajectory frames of benzamidine corresponding to initial (green) and other (red) metastable state portions, each 20 ns long, within the binding site on trypsin (grey).

Figure 5

Benzamidine is gradually pulled away from trypsin (details in Section S1 of the Supporting Information). The force profile height reflects the amount of resistance against the pulling force, which drops to near zero when benzamidine is far enough to escape the influence of trypsin. While being a crude measurement of the potential of mean force, this calculation adequately serves as a quick and simple means of locating the unbound state along the reaction coordinate. Note the correspondence of the force peak around z = 2 Å to the potential of mean force barrier of the initial metastable state.

Figure 6

Loop times are defined as the time taken for the system, starting from z_min, to reach z₀ and return to z_min without reaching z_max first. The distribution of loop times shown above was obtained from the portion of the equilibrium trajectory corresponding to the initial metastable state. The average loop time is extracted from the distribution, as described in Section S3 of the Supporting Information, as a necessary step in calculating the dissociation rate.

Figure 7

AMS trajectories are histogrammed by branching point z coordinates. Displayed with a logarithmic scale in the y axis, the histogram clearly shows a large concentration of loops about the initial metastable state.