Markov-State Transition Path Analysis of Electrostatic Channeling.

Electrostatic channeling is a naturally occurring approach to control the flux of charged intermediates in catalytic cascades. Computational techniques have enabled quantitative understanding of such mechanisms, augmenting experimental approaches by modeling molecular interactions in atomic detail. In this work, we report the first utilization of a Markov-state model (MSM) to describe the surface diffusion of a reaction intermediate, glucose 6-phosphate, on an artificially modified cascade where hexokinase and glucose-6-phosphate dehydrogenase are covalently conjugated by a cationic oligopeptide bridge. Conformation space networks are used to represent intermediate transport on enzyme surfaces, along with committor probabilities that assess the desorption probability of the intermediate on each segment of the channeling pathway. For the region between the peptide bridge and downstream active site, the ionic strength dependence of desorption probability by MSM agreed well with that by transition state theory. A kinetic Monte Carlo model integrates parameters from different computational methods to evaluate the contribution of desorption during each step. The approach is validated by calculation of kinetic lag time, which agrees well with experimental results. These results further demonstrate the applicability of molecular simulations and advanced sampling techniques to the design of chemical networks.

Introduction

Multistep catalytic cascades have the potential to increase the catalytic efficiency of chemical/biochemical synthesis via well-controlled reactant and intermediate flux and to enhance the output of catalytic devices such as biofuel cells and biosensors.[1−3] Natural metabolism reveals a variety of cascade reactions confined in single compartments, wherein a wide range of biomolecules are synthesized from few simple precursors, and energy is efficiently extracted by deep oxidation of biofuels.[4,5] A key limitation of these stepwise reactions is the mass transport of reaction intermediates, which in nature is found to be largely facilitated by substrate channeling phenomena, wherein reaction intermediate molecules are directly shuttled to downstream active sites instead of equilibrated to bulk media.[6] Given the complex chemical environment of the cell, channeling increases catalytic efficiency and can prevent intermediates from binding to unproductive sites. The elucidation of these transport processes is therefore of great significance to the understanding of natural channeling mechanisms and application to synthetic cascades.

An effective reaction cascade requires an appropriate spatial organization of catalysts and a favorable pathway between sequential reactive sites. Simple proximal location of such sequential sites has so far proved insufficient to account for the increase of catalytic efficiency, and functional channeling based on molecular interactions is required.[6−8] Widely recognized natural channeling mechanisms include intramolecular tunnels,[9,10] electrostatic channeling,[11−15] chemical swing arms,[16,17] and spatial organization.[18−20]

Electrostatic channeling may be described as diffusion of a charged intermediate molecule within an electrostatic field generated by an oppositely charged pathway. Electrostatic channeling represents a more general mechanism for channeling of ionic substrates, as compared to structure-specific chemical swing arms and sterically bound intramolecular tunneling. So far, the study of natural electrostatic cascades has focused on the bifunctional enzyme thymidylate synthase–dihydrofolate reductase and the tricarboxylic acid (TCA) cycle supercomplex malate dehydrogenase–citrate synthase.[11,12,14,21−24] For both natural cascades, X-ray diffraction and cross-linking mass spectrometry have revealed a positively charged pathway between the two active sites, suggesting an electrostatic channeling process.[11,14] Meanwhile, cascade kinetics have been experimentally studied via lag time analysis and competing bulk reactions, both of which show a significant increase in catalytic efficiency as compared to the nonchanneled enzyme pairs.[21,22] Additionally, the kinetic increase can be largely disabled by increasing the ionic strength (IS) of the surrounding environment[23,25] or by neutralization of the channeling pathway.[12,24] Such mechanisms were validated by our recent report on synthetic cascades with artificially introduced electrostatic channeling, wherein hexokinase (HK) and glucose-6-phosphate dehydrogenase (G6PDH) were covalently conjugated by a cationic oligopeptide bridge.[25]

So far, experimental results provide strong evidence of the occurrence of electrostatic channeling but are also limited because of the lack of in situ insight at a molecular level. In response, various types of models were developed to further study the channeling mechanism and quantify cascade kinetics. Continuum modeling was first applied to study the channeling process, where the spatial distribution of reaction intermediates was represented by a continuous media.[13,26] Using well-defined geometry and transport parameters, the migration of intermediates was shown to strongly control cascade efficiency and product yield, particularly at low IS. However, molecular-level interactions cannot be well represented in the continuum model, and the cascade geometry was coarse-grained for simplicity.

The first molecular simulations of electrostatic channeling were reported by McCammon’s group, using Brownian dynamics to study the probability that simplified but explicit intermediate molecules would reach the vicinity (∼0.7 nm) of a downstream active site.[7,12,15,27] The resulting probabilities were applied to microkinetic models to estimate the overall cascade kinetics.[23] These simulations gave a more detailed description of the channeling process than continuum modeling, but the strong hydrogen bond interactions between charged/polarized intermediate molecules and the cascade surface were not well represented. Short-range electrostatic interactions with cascade surfaces were recently shown by our work to be important in electrostatic channeling, especially at high IS.[25,28] Our study on artificially modified cascades (HK–G6PDH) gave a detailed explanation of the surface diffusion via hydrogen bond interactions and its contribution to cascade kinetics.[25,28] In this work, a surface “hopping” mode was established by molecular dynamics (MD) simulations, and hydrogen-bond-directed surface diffusion was observed on an energy-uniform surface, such as a fully saturated lysine (LYS) and arginine (ARG) oligopeptides. A kinetic Monte Carlo (KMC) model was used to validate the model against experimental results, pointing to the channeling pathway length and desorption during individual hops as key factors to overall channeling efficiency. In addition, an IS study using transition state theory (TST) and umbrella sampling (US) showed that the desorption probability for individual hops on a saturated LYS surface increased from 0.5% at 0 mM to 3.5% at 120 mM IS. During transport over the short (∼2 nm) and nonuniform pathway between the bridge and G6PDH active site, the desorption probability was estimated to be 23% at 120 mM IS.

Markov-state models (MSMs) can be used to integrate large collections of short MD trajectories and predict long time-scale phenomena, based on the assumption that future states depend only on the current state instead of the history.[29−32] Recently, MSM has been widely used to describe protein dynamics and ligand unbinding, where the systems include conformations with multiple degrees of freedom and highly interconnected transition pathways.[33−41] As a result, a conformation space network (CSN) is usually mapped to visualize the dynamic transitions between meta-states of the system. Combined with transition path theory (TPT),[35,39,42−44] the probability of the flux along each explicit transition pathway can be quantitatively calculated to enrich the MSM output, giving a complete energy landscape of the conformation space.

This ability of MSM to describe energy landscapes over multidimensional transition pathways makes it a promising approach to describe intermediate channeling in reaction cascades. Herein, we report an MSM to study the intermediate (glucose 6-phosphate, G6P) transitions on the pathway segment between the upstream HK active site and LYS bridge as well as the region from the bridge to the downstream G6PDH active site. On the basis of parallel MD simulations of intermediate trajectories around the cascade surface, an MSM transition matrix was obtained and a CSN[37] was used to visualize the transition pathways during the channeling process. On the basis of TPT,[35,42−44] committor probabilities were calculated to estimate the desorption probabilities along each channeling segment, which were then integrated in a KMC model and validated against experimentally measured stopped-flow lag times.

Computational Details

Details of MD simulations, the KMC model, and experiments can be found in the Supporting Information and our previous work.[25,28]

Featurization of MD Trajectories

MD trajectories for the bridge–G6PDH segment were taken from our previous work,[28] wherein the G6P molecule was initialized on the surface of G6PDH, between the peptide bridge and binding pocket, such that it had comparable probability to reach either site. Starting at the release point, 500 parallel short simulations (2 ns) were conducted with randomly regenerated atom velocities and counter ions. The resulting 105 MD frames were reanalyzed using an MSM in this work.

These MD trajectories with full atom coordinates were featurized into a representative time-dependent vector. Reference groups were selected on the G6P molecule and bridge–G6PDH complex to represent the position and orientation of the G6P molecule. Specifically, the phosphorus atom and center of mass (COM) of G6P were selected, along with the COM of 15 residues on the surface of the bridge–G6PDH complex (Figure S1). Additionally, the atom of the bridge–G6PDH complex that was nearest to G6P COM was dynamically chosen to represent the surface distance. As a result, by calculating distances between the two G6P reference groups and 16 groups on the cascade surface, a feature vector of 32 elements, V(t) = [v1, v2, v3, ..., v32], was obtained from the MD trajectories.

For the HK–bridge segment, 500 parallel, 10 ns MD trajectories were generated. In addition to the phosphorus atom and COM of G6P, 20 surface residues and one dynamically chosen nearest atom were used as reference groups. This resulted in a 42-element feature vector, V(t).

Markov-State Model

The Python package MSMBuilder was used to build MSMs.[45] Using the “MiniBatchKMeans” clustering algorithm, the positions and orientations of G6P relative to the cascade surface (HK–bridge and bridge–G6PDH) in each MD frame were grouped into 1000 meta-states, [s1, s2, ..., s1000]. Each meta-state contains many MD frames with similar conformations based on the feature vector. Transitions between these meta-states over a time step of duration Δt = 10 ps were counted, resulting in a 1000 × 1000 count matrix, N(Δt).

As the discrete hopping between each meta-state pair was dynamically reversible, the count matrix, N, may be optionally symmetrized using the relation Nsym = (N + NT)/2, where NT is the transpose of N. Such a symmetric count matrix reflects an equilibrium condition where the counts of forward and reverse transitions between each meta-state pair are equal. Though the symmetrization step was included in this study, results are also presented for the case where the count matrix is not symmetrized (Supporting Information).

Each row of the count matrix may be normalized to obtain the discrete probability distribution of transition from a single meta-state to all other meta-states. The resulting transition probabilities were stored in a transition matrix, M(Δt), as shown in Figure S2. If the system is known to be in a meta-state s, M(Δt) represents the likelihood of making a transition to meta-state, s, after one time step.[29,46] Such transitions can be calculated over n time steps by repeated multiplication by M(Δt)

As n → ∞, the transition matrix, M, converges to a stationary transition matrix, M(∞), representing a stationary distribution of system states that follows a Boltzmann distribution. Each row of M(∞) is identical, and represents the eigenvector, ϵ, associated with the largest eigenvalue of M(Δt).

A CSN of G6P channeling on each segment was rendered using Gephi software with the “Yifan Hu” layout algorithm.[47] Each node of a CSN represents a meta-state, s, with node size proportional to its stationary population, M(∞) = ϵ. Edges connecting each node pair indicate transitions between the meta-states, with thickness representing the average transition probability (M(Δt) + M(Δt))/2. The CSN is thus a visual representation of the entire energy landscape sampled by the ensemble of MD simulations.

Committor Probabilities

Committor probabilities, based on TPT, were employed to quantify the probability of desorption, and therefore the efficiency of the channeling process.[35,42−44] Committor probabilities are calculated by defining groups of meta-states (called basins) and determining the probability that the system, starting from a specific meta-state (node of the CSN), will reach a particular basin first—that is, before any of the other basins. These basins are considered energy minima from which the system will not re-emerge. For any meta-state, the committor probabilities to all basins sum to 1, as the system is ergodic and will ultimately reach one of the basins.

On the bridge–G6PDH segment, the set of meta-states wherein the phosphorus atom of G6P lay less than 0.7 nm from the bridge surface was denoted as the bridge basin. Meta-states with the phosphate atom less than 1.4 nm from the COM of G6PDH’s binding pocket were denoted as the pocket basin, and meta-states with phosphorous atoms more than 0.7 + 1.2 = 1.9 nm above the cascade surface were defined as the desorption basin. For the HK–bridge segment, the cutoff for the binding pocket basin was 2.0 nm. A committor probability matrix, S(∞), was obtained by first zeroing each row of M(Δt) associated with a basin, and setting the diagonal element to one, yielding a sink matrix S(Δt). Then, S(∞) was obtained using S(Δt) for n → ∞.

Using row S(∞) and summing over columns associated with the three basins (bridge, pocket, desorption), a three-component committor probability vector, pc3(i), can be calculated for each meta-state swhere superscripts indicate the probability of transition to a specific basin. Note that for each meta-state, s, the components of pc3(i) sum to 1.

Using pc3(i), an overall three-basin desorption probability, p3bsndes, was obtained. A set of states, s, were chosen which had approximately equal probability of reaching the pocket or the bridge, specifically |pc3poc(k) – pc3brg(k)| < 2%. This can be considered as a set of transition states between the bridge and the pocket. Desorption probability was then calculated according towhere the summation is weighted by the stationary population of each meta-state, given by ϵ.

Alternatively, a two-basin committor probability, pc2(i), was calculated by disabling one basin, typically the source of that segment of the channeling pathway. Using the HK–bridge segment as an example, only the bridge and desorption were defined as basins; two-basin committor probabilities and the resulting desorption probability were calculated for meta-states associated with the pocket by

Results and Discussion

Overview of Cascade Structure and Computational Methodologies

Figure shows the cascade structure (HK–bridge–G6PDH), where HK and G6PDH were covalently linked by a fully saturated peptide bridge. Our previous work showed that the energy and rate terms can be well quantified for individual hops on the energy-uniform surface of the peptide bridge.[25,28] That is, each hop was accompanied by a certain probability of desorption, depending on the energy difference between the hopping and desorption energy barriers (Figure b). However, the pathway between bridge and active sites is a nonuniform 2-dimensional surface composed of flexible hydrophobic, polarized, and charged amino acid residues, necessitating more sophisticated computational approaches.

Figure 1

Computational methodologies to quantify the desorption probability of glucose 6-phosphate on the cascade structure (HK–bridge–G6PDH) composed of HK and G6PDH covalently linked by a LYS oligopeptide bridge. (a) Channeling on a HK–bridge segment by MSM in this work. (b) Discrete hopping/desorption on a peptide bridge studied via TST and US in our previous work.[25,28] (c) Channeling on a bridge–G6PDH segment via TST[28] and MSM.

Specifically, the net charge of G6PDH is slightly positive (+1) and it is even more positive around the binding pocket. On the bridge–G6PDH segment, therefore, the G6P molecule (−2 charge) mostly experiences attractive electrostatic interactions with an oppositely charged surface. As shown in Figure c, the bridge–G6PDH segment is short (2 nm) and a simplified energy barrier can be assumed between the bridge and G6PDH’s binding pocket. Previously, a direct probability measurement was conducted by releasing the G6P molecule around the transition region, which resulted in a comparable probability for G6P to reach the bridge and pocket sites. By utilizing the last MD frame from each parallel simulation, a transition-state-based desorption probability was estimated and applied to a KMC model to quantify the cascade kinetics.[28]

In contrast, the HK monomer possesses a highly negative net charge of −16. Electrostatic channeling in such an environment depends on local surface polarization, with a higher probability of desorption. Moreover, the HK–bridge segment is longer (∼4.1 nm) as shown in Figure a, and there exist multiple energy barriers and minima, making the HK–bridge segment more complex as compared to the bridge–G6PDH segment. As a result, the simplification to a single energy barrier, as previously applied to G6PDH, cannot be well established for G6P transport over HK.

Moreover, the previous TST desorption probability for a bridge–G6PDH segment only utilized the final frame of each parallel MD simulation to determine the channeling results, a small sample of the channeling process. To fully utilize all MD frames and give a completed description of intermediate transport, two separate MSMs are used in this work to integrate full MD trajectories and map the CSN for G6P’s channeling, first from the HK pocket to the LYS bridge (Figure a), and then from the bridge to the G6PDH pocket (Figure c). Transition-path-based committor probabilities[35,42−44] are employed to quantify the desorption probability in each region.

CSN for G6P on a Cascade Surface

Figure shows the CSN of the G6P molecule on the bridge–G6PDH segment in an ionic environment at 120 mM. Three major areas are circled in Figure a, including the “bridge” site, “pocket” site, and transition area. After MD simulations, G6P’s trajectories were collected and the density isosurfaces were mapped as shown in Figure b,c. Figure d–f shows the resulting CSNs, where each node stands for a meta-state corresponding to a group of similar G6P’s positions on a cascade surface. Node size is proportional to meta-state population and edges indicated transitions between nodes. Three basins (Figure d) were first defined to calculate the committor probabilities, including a “bridge” basin in blue, “pocket” basin in green, and “bulk” basin in red. On the basis of the weight of the components of each committor, an RGB color was assigned to each meta-state, as shown in Figure e. Consequently, several major meta-state communities can be identified, such as a blue “bridge” community, green “pocket” community, red “bulk” community, and cyan surface transition communities. Strong connectivity between these CSN communities (Figure e) and the overlap between density isosurfaces (Figure b) indicate that the coverage of MD trajectories was sufficient to map the CSN and estimate the desorption probability during G6P transport from the cationic bridge to G6PDH binding pocket.

Figure 2

Conformational ensemble of the G6P molecule on the bridge–G6PDH segment at 120 mmol/L IS. (a) Top view of the bridge–G6PDH complex structure. (b) Top view of G6P density isosurfaces corresponding to different CSN subsets. (c) Side view of G6P’s isosurfaces. (d) CSN with three pre-defined basins, including “bulk” in red, “pocket” in green, and “bridge” in blue. Node size is proportional to the population and the edges indicate the transitions between nodes. (e) CSN colored by three-basin committor probabilities to each basin. (f) Two-basin CSN using “pocket” and “bulk” basins. (g,h) Two-basin CSN for IS of 40 and 0 mmol/L cases. (i) Triangle color bar where RGB components represent committor probabilities to each basin.

The cyan states, for example, have roughly equal committor probabilities to green and blue basins, as shown by the triangular color bar in Figure i. Meanwhile, these states also have a certain committor probability to the desorption basin, pc3des(i). As stated in Methods, a three-basin desorption probability, p3bsndes, was estimated by taking the weighted desorption probability for the meta-states with comparable committor probabilities (<2%) to bridge and pocket basins (eqs –4). This approach to determining desorption probability represents the basic idea in our previous work,[28] wherein desorption probability was estimated by a direct probability analysis based on TST. This method is effective when no strong energy minima exist outside of the three basins, and the transition meta-states can be assumed to possess high energy levels.

An alternative analysis involving two-basin desorption probability was conducted by disabling the bridge as a basin. As a result, only pocket and bulk were considered as basins when calculating the committor probability. Graphically, the blue color was removed from Figure e and the red and green colors were renormalized to sum to one. This significantly affects the original blue region, as shown by the oval circle in Figure f. In this way, the color of the original bridge basin gives a qualitative visualization of the desorption probability of G6P when starting from the LYS bridge, as compared to the probability of channeling to G6PDH’s active site. For example, the yellow states in Figure f have equal probability of transitioning to either the desorption (red) or pocket (green) basin. As shown in Figure f–h, the color of the previous bridge basin states (oval) turned from green to yellow with increasing IS from 0 to 120 mM, indicating significant desorption under concentrated ionic environments. A two-basin desorption probability, p2bsndes, was calculated according to eqs and 6.

The yellow color in Figure f–h was useful to qualitatively visualize the desorption probability of any transition state. Therefore, we consistently color nodes such that the desorption basin is red, the source of the channeling pathway is blue, and the destination is green. Therefore, the bridge basins are colored green in the CSN of the upstream HK–bridge segment and blue in the downstream bridge–G6PDH CSN.

The HK–bridge segment is longer (∼4.1 nm) as shown in Figure and the length of each parallel simulation was extended to 10 ns. As compared to the CSN of bridge–G6PDH (Figure ), one extra major cluster was found between the pocket and bridge communities, as shown in Figure e. This state community, which appears as cyan, was associated with the LYS-428 and ARG-423 of HK, and represents an electrostatic energy minimum on the HK surface. Such “kinetic traps” justify the MSM approach over the simplified energy barrier previously employed for the bridge–G6PDH segment.[28]

Figure 3

Conformational ensemble of the G6P molecule on the HK–bridge pathway at 0 mmol/L IS. (a) Top view of the HK–bridge complex structure. (b) Top view of G6P density isosurfaces corresponding to different CSN subsets. (c) Side view of G6P’s isosurfaces. (d) CSN with three predefined basins, including “bulk” in red, “pocket” in blue, and “bridge” in green. Node size is proportional to the population and the edges indicate the transitions between nodes. (e) CSN colored by three-basin committor probabilities to each basin. (f) Two-basin CSN using “bridge” and “bulk” basins. (g,h) Two-basin CSN for IS of 40 and 120 mmol/L cases. (i) Triangle color bar where RGB components represent committor probabilities to each basin.

CSNs colored by 2-basin committor probabilities, as shown in Figure f, enable the visualization of IS-dependent desorption probability from the HK to the bridge. As IS increases from 0 to 120 mM, the color of the original pocket basin transitions from light green to orange, indicating an increasing affinity to bulk states (red). For example, the resulting desorption probability by eq was 41% at 0 mM, which agrees well with the CSN colored by 2-basin committor probability shown in Figure f.

IS Dependence and Validation by Direct Probability Measurement

Comparing the CSNs of the HK–bridge and bridge–G6PDH segments, differences include channeling segment length (∼4.1 nm vs ∼2 nm) and additional energy minima on the HK surface (Figures e and 3e). Significantly, the bridge basin is the destination for HK–bridge transport, but the source in the bridge–G6PDH case. As previously discussed regarding Figure , hopping on the LYS bridge can be well represented by TST and US, whereas MSM is employed mainly to study channeling on enzyme surfaces. For this reason, application of a 2-basin desorption probability, p2bsndes, to the bridge–G6PDH segment would introduce significant desorption from the bridge, which was already parameterized by TST and US in the KMC model.

In contrast, a 3-basin desorption probability (p3bsndes) only includes desorption from a single bridge residue closest to G6PDH, which effectively avoids parameter overlap between the bridge and bridge–G6PDH segments. No significant energy minima were observed on the G6PDH surface between the bridge and pocket, indicating that intermediate transport on this area is mostly overcoming transition states at high energy levels. For these reasons, p3bsndes from a 3-basin CSN (Figure e) was chosen to represent transport over the bridge–G6PDH segment.

To summarize, a 2-basin desorption probability, p2bsndes, was used for the HK–bridge segment and 3-basin desorption probability, p3bsndes, was used for the bridge–G6PDH segment. Figure summarizes the IS dependence of desorption probability for G6P transport on the HK–bridge and bridge–G6PDH segments and for individual hops on the LYS bridge. Generally, all types of desorption increase with increasing IS values because of increased shielding of electrostatic interactions. Specifically, individual hops on the LYS bridge reflect the lowest desorption probability, because of the strong hydrogen bond interactions between G6P’s phosphate and LYS’s ε-ammonium groups. As discussed in our previous work,[28] such strong interactions were able to maintain a desorption probability as low as 3.5% at high IS (120 mM). As for the bridge–G6PDH segment, the bridge is strongly charged and G6PDH has a relatively neutral net charge (+1.1). As a result, the long-range electrostatic impact of the LYS bridge was able to fully cover the transition area and even reach the binding pocket of G6PDH, providing a good electrostatic confinement of surface diffusion. As a result, desorption probability on the bridge–G6PDH segment at 0 mM can be as low as 5% as shown in Figure . At 120 mM, that probability increased to 25%. Through the entire IS range, the desorption probabilities as calculated by the 3-basin CSN, p3bsndes, agree with previous results from TST analysis within error.[28]

Figure 4

IS dependence of desorption probabilities for different segments of cascade surface, including the HK–bridge segment, bridge–G6PDH segment, and single hops on the bridge. Desorption on the bridge–G6PDH segment was quantified both by MSM and TST, each using the same MD trajectories.[28] Comparison with the results of a nonsymmetrized MSM can be found in Figure S3.

On the HK–bridge segment, desorption probability increased to 40–62%. One reason is the net charge of −16 on the HK subunit, which introduces long-range repulsion that frustrates electrostatic interactions with the LYS bridge. The Debye length, a measure of the range of electrostatic interactions, decreases from 9.8 nm at 1 mM IS, to 2.2 nm at 20 mM and 0.9 nm at 120 mM. Given the 4.1 nm length of the HK–bridge segment, long-range electrostatic confinement was likely negligible at higher IS, and channeling relies on local electrostatic interactions, with higher probability of desorption. Only when the intermediate approaches the bridge, (light blue cluster in Figure e), did desorption probability significantly decrease.

It should be noted that the channeling efficiency at low IS may be underestimated, because meta-states of more than 1.9 nm from the cascade surface were defined as the desorption basin, even though Debye lengths were greater than this distance (e.g., 9.8 nm at 1 mM). Therefore, the bridge–HK complex still exerts certain electrostatic forces on intermediate molecules in “desorption” states.

KMC Model on Cascade Kinetics

To further quantify the cascade kinetics, a previously developed KMC model was used to integrate probabilities generated by the MSM and other computational models.[28]Figure a shows the schematic diagram of KMC model. Specifically, the two sequential active sites, E1 and E2, were connected by discrete hopping sites numbered from 1 to 4. Explicit rate constants were assigned to individual events, including hopping, desorption, adsorption, and reaction. A universal bulk reservoir was employed to modulate the concentration of reaction intermediates. All event sites were allowed to exchange molecules to the bulk environment. Hopping on the bridge is reversible, but is irreversible between active sites (E1, E2) and the bridge.

Figure 5

Cascade kinetics via the KMC model. (a) Schematic diagram of the KMC model. (b) Comparison of experimental lag time and KMC lag time, including intermediate desorption on various segments. “Free” represents the nonconjugated cascade (complete desorption on the KMC bridge); “bridge” includes desorption during individual hops on the peptide bridge: “bridge–G6PDH” adds desorption on the bridge–G6PDH segment, using either the MSM or TST as shown in Figure . “HK–bridge–G6PDH” reflects the cascade kinetics involving desorption on all segments.

MSM-derived desorption probability (Figure ) was used to parameterize the KMC model. As discussed in our previous work,[28] the intermediate surface transition rate (>1 μs–1) was much faster than the turnover frequency of active sites (0.6–7 s–1). Therefore, instead of their absolute values, the ratio of the hopping and desorption rate constants were considered. The MSM desorption probabilities, p2bsndes and p3bsndes, were used to parameterize the ratio khop/kdes according to:

Rate constants derived from MSM and other methodologies are summarized in Table S1. Parameterization details for other KMC rate constants can be found in our previous work.[28]

Figure b shows IS-dependent lag time plots for the HK–bridge–G6PDH cascade with calculated results compared to the experiment. Detailed discussion of these lag time experiments as well as the TST results and calculation for the bridge alone can be found in our previous work.[25,28] Briefly, the final experimental solutions (pH = 7.0) contained 1.4 mM citrate, 0.5 mM ATP, 5.6 mM MgCl2, and 0.8 mM NADP+, adding up to an approximately 20 mM background IS. In the stop-flow lag time experiments, substrate-saturated kinetics at HK and observable lag times were obtained by using concentrations of 278 mM for glucose and 8 nM for the HK–bridge–G6PDH complex. Product evolution at the downstream enzyme (G6PDH) was determined spectrophotometrically by monitoring the production of NADPH at a wavelength of 340 nm.

The turnover frequency of G6PDH (∼6.2 s–1) was much faster than that of upstream HK (∼0.7 s–1), maximizing the capability for channeling as compared to desorption. In a reaction-limited process, such a fast downstream reaction rate results in an empty channel as described by queuing theory.[48]

The contribution of a bridge comprising five LYS residues was studied by disallowing desorption between the bridge and active sites. The slope of the of KMC lag time (bridge) agrees well with the experiment, but the absolute values were consistently lower. By including desorption in the bridge–G6PDH region, as calculated by TST (bridge–G6PDH–TST), the calculated lag time shifts upward, closer to the experimental data, whereas the slope becomes more comparable to the experiment. Replacing the TST calculation with the MSM result (bridge–G6PDH–MSM) did not significantly change the lag time plot, consistent with the desorption probability plots in Figure . Finally, MSM-based desorption from the HK–bridge segment is considered to yield a KMC analysis of the complete pathway (HK–bridge–G6PDH). The calculated lag time shifts further upward to the same level as the experimental lag time, with best agreement at high IS value, and underpredicted lag time at low IS.

As discussed above, the 1.9 nm definition of the desorption basin likely explains discrepancies at low IS. This boundary is controlled largely by limitations in the size of the MD box (∼2–4 nm gap between cascade surface and box boundary). In the present MD simulations with full atoms, it is impractical to specify a distance between the cascade surface and box boundary greater than 10 nm. This limitation may be overcome by using a multiscale approach; for example, MD simulations could be used to build the MSM for the proximal surface diffusion and Brownian dynamics or implicit solvent MD simulations used to consider long-range electrostatics. Such a combination of techniques in a single CSN is a topic for future study.

Design Rules for Future Cascades

On the basis of the present analysis, cascade designs incorporating electrostatic channeling must still seek to maximize long-range electrostatic interactions as well as site–site proximity. Mutation of linker sites could be an option to realize the former goal and segment length should be well balanced to enhance adsorption but not greatly exceed the Debye length in the working environment. Where possible, any strong local energy minima on channeling surfaces should be eliminated, in order to create a continuous, favorable transition pathway and avoid significant desorption.

Conclusions

A Markov-state approach was applied to the simulation of surface diffusion via electrostatic channeling. The multidimensional pathways over the surface of a two-enzyme cascade were mapped from millions of physical states to a smaller distribution of meta-states in a human-readable CSN. MSM-derived intermediate transport efficiency from the LYS bridge to the binding pocket of G6PDH was found to agree well with a previous, simplified TST analysis. The longer pathway between HK and the LYS bridge introduced additional energy minima as well as overall electrostatic repulsion, resulting in significantly increased desorption.

A KMC model was used to integrate desorption probabilities as calculated by varying computational methods, and allowed comparison to experimental results. Desorption probabilities from different segments of the cascade (bridge, bridge–G6PDH, and HK–bridge) were introduced independently, elucidating the contribution of each segment to the overall cascade efficiency. Comparison with the experiment was poorest at low IS, where MD simulations are limited to length scales lower than the Debye length. A combination of MD and another, more coarse-grained molecular simulation could be a possible way to address this issue.