Full kinetics of CO entry, internal diffusion, and exit in myoglobin from transition-path theory simulations

Abstract

We use Markovian milestoning molecular dynamics (MD) simulations on a tessellation of the collective variable space for CO localization in myoglobin to estimate the kinetics of entry, exit, and internal site-hopping. The tessellation is determined by analysis of the free-energy surface in that space using transition-path theory (TPT), which provides criteria for defining optimal milestones, allowing short, independent, cell-constrained MD simulations to provide properly weighted kinetic data. We coarse grain the resulting kinetic model at two levels: first, using crystallographically relevant internal cavities and their predicted interconnections and solvent portals; and second, as a three-state side-path scheme inspired by similar models developed from geminate recombination experiments. We show semiquantitative agreement with experiment on entry and exit rates and in the identification of the so-called "histidine gate" at position 64 through which ≈90% of flux between solvent and the distal pocket passes. We also show with six-dimensional calculations that the minimum free-energy pathway of escape through the histidine gate is a "knock-on" mechanism in which motion of the ligand and the gate are sequential and interdependent. In total, these results suggest that such TPT simulations are indeed a promising approach to overcome the practical time-scale limitations of MD to allow reliable estimation of transition mechanisms and rates among metastable states.

Figure 1

A sketch of cells and milestones at a portal. P is the portal milestone at HG or Xe1s. S is the spherical boundary comprised of milestones interfacing solvent. (Note that we label each milestone with double-letter indices of cells that share that milestone.)

Figure 2

Overlay of three selected isosurfaces (−12.0, −8.0, and −5.4 kcal/mol) of the 3D FES of the CO COM inside Mb. Discretized MFEP’s computed by direct string method on the full FES are represented by blue balls. Free energies for the cavities are (in kcal/mol): DP, −15.3; Xe4, −14.2; Xe2, −13.1; Xe1, −11.6; ph, −10.9; Xe3, −10.8; NMR, −9.8.

Figure 3

Hitting points (blue) on cell boundaries from all milestoning MD simulations. Discretized MFEP’s computed by direct string method on the full FES are presented by red balls. A full 360° 3D view of the hitting points is available in Supporting Information.

Figure 4

Left panel: Pathways of ligand diffusion in Mb (Magenta). Blue spheres indicate metastable ligand localization sites. Cavity and portal sites are indicated by integer labels which correspond to states of the kinetic networks in the right panel. The sites, PHE33s (9), Xe3s (10), Phs (11), NMRs (12), MIDs (13), Xe2s (15), HG (8), and Xe1s (14) are portal sites, which may be buried in or on the surface of protein, but has channels to solvent. The sites DP (1), Xe4 (2), Xe2 (3), Xe1 (4), Xe3 (5), ph (6), and NMR (7) are cavities. Right panel: the site-specific kinetic network modeling escape of one CO from Mb. Nodes are numbered as “states” with ligand at the site corresponding to the same number in the left-hand figure. The MFPT’s are shown on each edge in ns. Only the MFPT’s on edge 1–8 or 4–14 are between the state (with CO in DP or Xe1) and the dissociated state (with CO in solvent). Dashed arrows indicate that the MFPT is above 2000 ns. The MFPT’s τ_dis‑DP and τ_dis‑Xe1 depend on the concentration of CO in solution (see Table 1 for the values at several different concentrations).

Figure 5

A coarse-grained kinetic network with MFPT’s in ns indicated on the arrows.

Figure 6

6D MFEP for the DP-HG migration projected into 2D subspaces of CO COM and COM of the imidazole ring: (a) x of ring COM vs x of CO COM; (b) y of ring COM vs y of CO COM; (c) z of ring COM vs z of CO COM. (d) Projection of the 6D MFEP into the space of displacement ((x – x₀) + (y − y₀) + (z − z₀))^1/2 of CO COM and of ring COM. (e) Free-energy profile along the MFEP. Circles denote stage 1 and triangles stage 2 (see text). The first image (index 0) corresponds to DP milestone and the last one (index 13) corresponds to the HG milestone.