Two-dimensional Free Energy Surface Research Articles

Structures of Th4+ aqueous solutions: insights from AIMD and metadynamics simulations.

Solution chemistry of actinide ions is critical to understanding the solvation behaviors and hydrolysis process. Using tetravalent thorium ion Th4+ as a representative example, we investigate the local structures and dynamic behaviors of hydrated Th4+ ions by ab initio molecular dynamics (AIMD) simulations using the recently developed norm-conserving pseudopotentials and basis sets optimized for actinides (J.-B. Lu et al., J. Chem. Theory Comput. 2021, 17, 3360-3371). AIMD simulations reveal two distinct solvation shells, with the first shell comprising 9 water molecules at approximately rTh-O = 2.50 Å and exhibiting a tricapped trigonal prism geometry. These conclusions are confirmed through metadynamics simulations and further structural analysis. AIMD simulations also show the slight effect of temperature and counterions on the structure of the solution. The structured solvation shells of the highly charged Th4+ ion with the specific geometry, distinct from the structure of liquid water, lead to corresponding structural changes in the hydrogen bond network in water. Additionally, beyond the solvent-shared ion pair (SIP) state observed in the unbiased AIMD simulations, the metadynamics simulations reconstruct a two-dimensional free energy surface that clearly indicates the potential stability of the contact ion pair (CIP) state in the system with Cl- as a counterion. The findings in this work provide insights into the solution chemistry of actinides and serve as a reference for studying other actinide solution systems.

Catalytic Mechanism of Human T-Cell Leukemia Virus Type 1 Protease Investigated by Combined QM/MM Molecular Dynamics Simulations.

Combined quantum mechanical and molecular mechanical (QM/MM) molecular dynamics simulations were performed to investigate the catalytic mechanism of human T-cell leukemia virus type 1 (HTLV-1) protease, a retroviral aspartic protease that is a potential therapeutic target for curing HTLV-1-associated diseases. To elucidate the proteolytic cleavage mechanism, we determined the two-dimensional free energy surfaces of the HTLV-1 protease-catalyzed reactions through various possible pathways. The free energy simulations suggest that the catalytic reactions of the HTLV-1 protease occur in the following sequential steps: (1) a proton is transferred from the lytic water to Asp32', followed by the nucleophilic addition of the resulting hydroxyl to the carbonyl carbon of the scissile bond, forming a tetrahedral oxyanion intermediate, and (2) a proton is transferred from Asp32 to the peptide nitrogen of the scissile bond, leading to the spontaneous breakage of the scissile bond. The rate-limiting step of this catalytic process is the proton transfer from Asp32 to the peptide nitrogen of the scissile bond, with a free energy of activation of 21.1 kcal/mol. This free energy barrier is close to the experimentally determined free energy of activation (16.3 kcal/mol) calculated from the measured catalytic rate constant (kcat). This mechanistic study provides detailed dynamic and structural information that will facilitate the design of mechanism-based inhibitors for the treatment of HTLV-1-associated diseases.

Molecular Dynamics Simulations of Ion Permeation in Human Voltage-Gated Sodium Channels.

The recent determination of cryo-EM structures of voltage-gated sodium (Nav) channels has revealed many details of these proteins. However, knowledge of ionic permeation through the Nav pore remains limited. In this work, we performed atomistic molecular dynamics (MD) simulations to study the structural features of various neuronal Nav channels based on homology modeling of the cryo-EM structure of the human Nav1.4 channel and, in addition, on the recently resolved configuration for Nav1.2. In particular, single Na+ permeation events during standard MD runs suggest that the ion resides in the inner part of the Nav selectivity filter (SF). On-the-fly free energy parametrization (OTFP) temperature-accelerated molecular dynamics (TAMD) was also used to calculate two-dimensional free energy surfaces (FESs) related to single/double Na+ translocation through the SF of the homology-based Nav1.2 model and the cryo-EM Nav1.2 structure, with different realizations of the DEKA filter domain. These additional simulations revealed distinct mechanisms for single and double Na+ permeation through the wild-type SF, which has a charged lysine in the DEKA ring. Moreover, the configurations of the ions in the SF corresponding to the metastable states of the FESs are specific for each SF motif. Overall, the description of these mechanisms gives us new insights into ion conduction in human Nav cryo-EM-based and cryo-EM configurations that could advance understanding of these systems and how they differ from potassium and bacterial Nav channels.

Ethanol exchange between two graphene surfaces in nanoconfined aqueous solution: Rate and mechanism.

We observe, by computer simulations, a remarkable long-distance, rare, but repetitive, exchange of ethanol molecules between two parallel graphene surfaces in nanoconfined, aqueous, ethanol solutions. We compute the rate of exchange as a function of the separation (d) between the two surfaces. We discover that the initiating (or, the launching) step in this exchange is the attainment of an instantaneous orientation of the carbon-oxygen bond vector relative to the graphene surface. This observation led us to construct a two-dimensional free energy surface for this exchange, with respect to two order parameters, namely, (i) the perpendicular distance of ethanol molecule from the graphene surfaces, z, and (ii) the orientation of the O-C bond vector, θ, of the tagged ethanol molecule. For d = 3nm, the rate of exchange is found to be 0.44ns-1 for the force field used. We also vary the force field and determine the sensitivity of the rate. From the free energy landscape, one could determine the minimum energy pathway. We use both, the transition state theory and Kramers' theory, to calculate the rate. The calculated rate agrees well with the simulated value as mentioned above. We find that the rate of exchange phenomenon is sensitive to the interaction strength of graphene and the hydrophobic group of ethanol. The free energy landscape exchange shows dependence on the distance separation of the two hydrophobic surfaces and reveals interesting features.

Strategy to improve Cu-BTC metal-organic frameworks performance in removal of Rhodamine B: MD and WT-MtD simulations assessment

With industry progress, environmental problems have begun to threaten human health. Among them, water pollution is closely related to human life and has attracted researchers’ attention. Hence, coping strategies for these pollutants have become a priority nowadays. Here, we carried out the molecular dynamics (MD) and well-tempered metadynamics simulations to evaluate the interaction of Rhodamine B (Rh B) with a copper (II)-benzene-1,3,5-tricarboxylate metal-organic framework (Cu-BTC/MOF). To design a more efficient dye removal platform, the effect of the -NH2, -OH, and -NO2 functional groups on the efficiency of Cu-BTC/MOF in the adsorption of Rh B is investigated. It is found that the interaction energy of Rh B with -NH2-MOF, -OH-MOF, and -NO2-MOF is about −79.98, −121.87, and −365.55 kJ mol−1, respectively, more than the pristine case. This observation confirms that the functionalization strategy can enhance the Cu-BTC/MOF efficiency. The obtained free energy (FE) values from the metadynamics simulation indicated that for adsorption of Rh B on pristine, -NH2-MOF, -OH-MOF, and -NO2-MOF, the global minimums are located about at −220.47, −234.75, −236.09, and −259.01 kJ mol−1, respectively. The obtained results show that in the two-dimensional FE surfaces, the most stable complex with Rh B belongs to the MOF-NO2 system.

Thermodynamics and Free Energy Landscape of BAR-Domain Dimerization from Molecular Simulations.

Proteins with BAR domains function to bind to and remodel biological membranes, where the dimerization of BAR domains is a key step in this function. These domains can dimerize in solution or after localizing to the membrane surface. Here, we characterize the binding thermodynamics of homodimerization between the LSP1 BAR domain proteins in solution, using molecular dynamics (MD) simulations. By combining the MARTINI coarse-grained protein models with enhanced sampling through metadynamics, we construct a two-dimensional free energy surface quantifying the bound versus unbound ensembles as a function of two distance variables. With this methodology, our simulations can simultaneously characterize the structures and relative stabilities of a range of sampled dimers, portraying a heterogeneous and extraordinarily stable bound ensemble, where the proper crystal structure dimer is the most stable in a 100 mM NaCl solution. Nonspecific dimers that are sampled involve contacts that are consistent with experimental structures of higher-order oligomers formed by the LSP1 BAR domain. Because the BAR dimers and oligomers can assemble on membranes, we characterize the relative alignment of the known membrane binding patches, finding that only the specific dimer is aligned to form strong interactions with the membrane. Hence, we would predict a strong selection of the specific dimer in binding to or assembling when on the membrane. Establishing the pairwise stabilities of homodimer contacts is difficult experimentally when the proteins form stable oligomers, but through the method used here, we can isolate these contacts, providing a foundation to study the same interactions on the membrane.

Free energy surface of initial cap formation in carbon nanotube growth.

Initial cap formation is an important process of carbon nanotube (CNT) growth where a hexagonal carbon network is lifted off from the catalyst surface. In this study, free energy surface (FES) of initial cap formation in the CNT growth is investigated by metadynamics simulation. A two-dimensional collective variable (CV) space is newly developed to examine the complicated formation process of the cap structure, which consists of the formation of a hexagonal carbon network and lift-off of the network from the catalyst surface. States before and after the lift-off of the carbon network are clearly distinguished in the two-dimensional FES. Therefore, free energy difference before and after the lift-off can be directly derived from the two-dimensional FES. It was revealed that the cap structure is stable at a high temperature due to the entropy effect, while the carbon network covering the catalyst surface is energetically stable. The new insight in this study is achieved owing to metadynamics simulation in conjunction with a newly developed two-dimensional CV space since it is impossible to explore FES for such complicated processes in the framework of conventional molecular dynamics simulation.

Efficient computation of free energy surfaces of chemical reactions using ab initio molecular dynamics with hybrid functionals and plane waves.

Ab initio molecular dynamics (AIMD) simulations employing density functional theory (DFT) and plane waves are routinely carried out using density functionals at the level of generalized gradient approximation (GGA). AIMD simulations employing hybrid density functionals are of great interest as it offers a more accurate description of structural and dynamic properties than the GGA functionals. However, the computational cost for carrying out calculations using hybrid functionals and plane wave basis set is at least two orders of magnitude higher than that using GGA functionals. Recently, we proposed a strategy that combined the adaptively compressed exchange operator formulation and the multiple time step integration scheme to reduce the computational cost by an order of magnitude [J. Chem. Phys. 151, 151102 (2019)]. In this work, we demonstrate the application of this method to study chemical reactions, in particular, formamide hydrolysis in an alkaline aqueous medium. By actuating our implementation with the well-sliced metadynamics scheme, we can compute the two-dimensional free energy surface of this reaction at the level of hybrid-DFT. This work also investigates the accuracy of the PBE0 (hybrid) and the PBE (GGA) functionals in predicting the free energetics of this chemical reaction.

Combining LE4PD Normal Modes and Markov State Modeling to Elucidate the Fluctuation Dynamics of Ubiquitin

Local fluctuations are important for protein binding and molecular recognition because they provide conformational states that can be trapped through a conformational selection mechanism of binding. Thus, an accurate characterization of local fluctuations can be important for modeling the mechanisms leading to the biological activity of a protein. In this presentation we quantify the fluctuation dynamics of the regulatory protein ubiquitin and propose a novel theoretical approach to modelling its fluctuations. A coarse-grained, diffusive, mode-dependent description of the protein's fluctuation dynamics is accomplished using the Langevin Equation for Protein Dynamics (LE4PD); this equation decomposes the dynamics of a protein, modeled through a molecular dynamics (MD) simulation, into dynamical pathways that explore mode-dependent free energy surfaces. Using as input to the theory a one-microsecond, atomistic MD simulation, we calculate the time scales of the slow, high-amplitude fluctuations by modeling the kinetics of barrier crossing in the two-dimensional free-energy surfaces using a Markov state model (MSM) approach. We find that the LE4PD predicts slow fluctuations in three biologically significant binding regions in ubiquitin: the C-terminal tail, the Lys11 loop, and the 50 s loop. The combination of the LE4PD modes with the MSM analysis can provide useful information on the role of fluctuations in the process of molecular recognition regulating the biological activity of ubiquitin, with the LE4PD modes describing the location and magnitude of high-amplitude conformational changes and the MSM giving the associated timescales.

Kinetics analysis of ubiquitin local fluctuations with Markov state modeling of the LE4PD normal modes.

Local fluctuations are important for protein binding and molecular recognition because they provide conformational states that can be trapped through a selection mechanism of binding. Thus, an accurate characterization of local fluctuations may be important for modeling the kinetic mechanism that leads to the biological activity of a protein. In this paper, we study the fluctuation dynamics of the regulatory protein ubiquitin and propose a novel theoretical approach to model its fluctuations. A coarse-grained, diffusive, mode-dependent description of fluctuations is accomplished using the Langevin Equation for Protein Dynamics (LE4PD). This equation decomposes the dynamics of a protein, simulated by molecular dynamics, into dynamical pathways that explore mode-dependent free energy surfaces. We calculate the time scales of the slow, high-amplitude fluctuations by modeling the kinetics of barrier crossing in the two-dimensional free energy surfaces using Markov state modeling. We find that the LE4PD predicts slow fluctuations in three important binding regions in ubiquitin: the C-terminal tail, the Lys11 loop, and the 50 s loop. These results suggest that the LE4PD can provide useful information on the role of fluctuations in the process of molecular recognition regulating the biological activity of ubiquitin.

SCC-DFTB-PIMD Method To Evaluate a Multidimensional Quantum Free-Energy Surface for a Proton-Transfer Reaction.

The self-consistent charge density functional tight binding method was combined with the path-integral molecular dynamics method for the first time to evaluate the two-dimensional free-energy surface including nuclear quantum effects of a proton-transfer reaction in a 2,4-dichlorophenol-trimethylamine complex. A statistically converged two-dimensional quantum free-energy surface was evaluated by the multidimensional blue moon ensemble method. The accuracy was guaranteed by optimizing the repulsive potential between the sp3-hybridized nitrogen and hydrogen atoms in a SCC-DFTB3 parameter set for the system to reproduce high-level quantum chemical calculations. The present study illustrates the usefulness of this new approach to investigate nuclear quantum effects in various realistic proton-transfer reactions.

Computational simulation of mechanism and isotope effects on acetal heterolysis as a model for glycoside hydrolysis

Abstract DFT calculations for the equilibrium isotope effect for deuterium substitution at the anomeric centre Cα in 2-(p-nitrophenoxy)tetrahydropyran with continuum solvation show significant variation in the range of relative permittivity 2 ≤ ε ≤ 10. One-dimensional scans of potential energy (with implicit solvation by water) or of free energy (from QM/MM potentials of mean force with explicit aqueous solvation with a hybrid AM1/OPLS method) for heterolysis of the bond between Cα and the nucleofuge do not show a transition state. A two-dimensional free-energy surface that considers also the distance between Cα and a nucleophilic water indicates a pre-association DN*ANint ‡ mechanism with a transition state involving nucleophilic attack upon an ion-pair intermediate, and this is supported by good agreement between the mean values of the calculated and experimental α-D KIEs. However, the magnitudes of the standard deviations about the mean values for the making and breaking C–O bonds suggest that the transition state is rather plastic, with Cα–Onu≈2 ± 0.4 Å and Cα–Olg≈3 ± 0.5 Å. Not only is nucleophilic solvent assistance necessary, but there is also evidence for electrophilic assistance through specific hydrogen bonding to the nucleofuge.

Molecular dynamics simulations of the interaction of phospholipid bilayers with polycaprolactone

ABSTRACTThe molecular interaction between common polymer chains and the cell membrane is unknown. Molecular dynamics simulations offer an emerging tool to characterise the nature of the interaction between common degradable polymer chains used in biomedical applications, such as polycaprolactone, and model cell membranes. Herein we characterise with all-atomistic and coarse-grained molecular dynamics simulations the interaction between single polycaprolactone chains of varying chain lengths with a phospholipid membrane. We find that the length of the polymer chain greatly affects the nature of interaction with the membrane, as well as the membrane properties. Furthermore, we next utilise advanced sampling techniques in molecular dynamics to characterise the two-dimensional free energy surface for the interaction of varying polymer chain lengths (short, intermediate, and long) with model cell membranes. We find that the free energy minimum shifts from the membrane-water interface to the hydrophobic core of the phospholipid membrane as a function of chain length. Finally, we perform coarse-grained molecular dynamics simulations of slightly larger membranes with polymers of the same length and characterise the results as compared with all-atomistic molecular dynamics simulations. These results can be used to design polymer chain lengths and chemistries to optimise their interaction with cell membranes at the molecular level.

Effect of Intercalated Water on Potassium Ion Transport through Kv1.2 Channels Studied via On-the-Fly Free-Energy Parametrization.

We introduce a two-dimensional version of the method called on-the-fly free energy parametrization (OTFP) to reconstruct free-energy surfaces using Molecular Dynamics simulations, which we name OTFP-2D. We first test the new method by reconstructing the well-known dihedral angles free energy surface of solvated alanine dipeptide. Then, we use it to investigate the process of K+ ions translocation inside the Kv1.2 channel. By comparing a series of two-dimensional free energy surfaces for ion movement calculated with different conditions on the intercalated water molecules, we first recapitulate the widely accepted knock-on mechanism for ion translocation and then confirm that permeation occurs with water molecules alternated among the ions, in accordance with the latest experimental findings. From a methodological standpoint, our new OTFP-2D algorithm demonstrates the excellent sampling acceleration of temperature-accelerated molecular dynamics and the ability to efficiently compute 2D free-energy surfaces. It will therefore be useful in large variety complex biomacromolecular simulations.

Transition paths in single-molecule force spectroscopy.

In a typical single-molecule force spectroscopy experiment, the ends of the molecule of interest are connected by long polymer linkers to a pair of mesoscopic beads trapped in the focus of two laser beams. At constant force load, the total extension, i.e., the end-to-end distance of the molecule plus linkers, is measured as a function of time. In the simplest systems, the measured extension fluctuates about two values characteristic of folded and unfolded states, with occasional transitions between them. We have recently shown that molecular (un)folding rates can be recovered from such trajectories, with a small linker correction, as long as the characteristic time of the bead fluctuations is shorter than the residence time in the unfolded (folded) state. Here, we show that accurate measurements of the molecular transition path times require an even faster apparatus response. Transition paths, the trajectory segments in which the molecule (un)folds, are properly resolved only if the beads fluctuate more rapidly than the end-to-end distance of the molecule. Therefore, over a wide regime, the measured rates may be meaningful but not the transition path times. Analytic expressions for the measured mean transition path times are obtained for systems diffusing anisotropically on a two-dimensional free energy surface. The transition path times depend on the properties both of the molecule and of the pulling device.

Molecular Dynamics of Oxazole Yellow Dye in its Ground and First Excited Electronic States in Solution and when Intercalated in dsDNA.

Oxazole yellow (YOPRO), a cyanine dye consisting of benzoxazole and quinoline rings connected by a linker, is almost nonfluorescent in water, but its fluorescence is greatly enhanced after intercalation in double-stranded DNA, forming the basis of DNA concentration assays. To explore this difference, new potential energy surfaces for the two linker dihedral angles in the ground S0 and first excited S1 electronic states are developed. Umbrella sampling molecular dynamics is used to obtain the free energy of rotation around the two dihedral angles of the linker. The two-dimensional free energy surface of the S1 state, spanning the Franck-Condon transition point from the S0 electronic state minimum (dihedral 1 around 180°, dihedral 2 around 0°) to the S1 state minimum (∼90, ∼0), is obtained in water and when intercalated. In water, YOPRO's S1 free energy surface is completely downhill from the Franck-Condon point, whereas when intercalated, there is a barrier on the path. Thus, when intercalated in DNA, S1 YOPRO is more constrained than in water, supporting the hypothesis that intercalation does inhibit ring rotational motion around the linker and therefore strongly reduces the nonradiative relaxation, resulting in higher fluorescence intensity.

Conformational Dynamics of the Lipopolysaccharide from Escherichia coli O91 Revealed by Nuclear Magnetic Resonance Spectroscopy and Molecular Simulations.

The outer leaflet of the outer membrane in Gram-negative bacteria contains lipopolysaccharides (LPS) as a major component, and the outer membrane provides a physical barrier and protection against hostile environments. The enterohemorrhagic Escherichia coli of serogroup O91 has an O-antigen polysaccharide (PS) with five sugar residues in the repeating unit (RU), and the herein studied O-antigen PS contains ∼10 RUs. 1H-13C HSQC-NOESY experiments on a 1-13C-labeled PS were employed to deduce 1H-1H cross-relaxation rates and transglycosidic 3JCH related to the ψ torsional angles were obtained by 1H-1H NOESY experiments. Dynamical parameters were calculated from the molecular dynamics (MD) simulations of the PS in solution and compared to those from 13C nuclear magnetic resonance (NMR) relaxation studies. Importantly, the MD simulations can reproduce the dynamical behavior of internal correlation times along the PS chain. Two-dimensional free energy surfaces of glycosidic torsion angles delineate the conformational space available to the O-antigen. Although similar with respect to populated states in solution, the O-antigen in LPS bilayers has more extended chains as a result of spatial limitations due to close packing. Calcium ions are highly abundant in the phosphate-containing core region mediating LPS-LPS association that is crucial for maintaining bilayer integrity, and the negatively charged O-antigen promotes a high concentration of counterbalancing potassium ions. The ensemble of structures present for the PS in solution is captured by the NMR experiments, and the similarities between the O-antigen on its own and as a constituent of the full LPS in a bilayer environment make it possible to realistically describe the LPS conformation and dynamics from the MD simulations.

Nonlinear reconstruction of single-molecule free-energy surfaces from univariate time series.

The stable conformations and dynamical fluctuations of polymers and macromolecules are governed by the underlying single-molecule free energy surface. By integrating ideas from dynamical systems theory with nonlinear manifold learning, we have recovered single-molecule free energy surfaces from univariate time series in a single coarse-grained system observable. Using Takens' Delay Embedding Theorem, we expand the univariate time series into a high dimensional space in which the dynamics are equivalent to those of the molecular motions in real space. We then apply the diffusion map nonlinear manifold learning algorithm to extract a low-dimensional representation of the free energy surface that is diffeomorphic to that computed from a complete knowledge of all system degrees of freedom. We validate our approach in molecular dynamics simulations of a C(24)H(50) n-alkane chain to demonstrate that the two-dimensional free energy surface extracted from the atomistic simulation trajectory is - subject to spatial and temporal symmetries - geometrically and topologically equivalent to that recovered from a knowledge of only the head-to-tail distance of the chain. Our approach lays the foundations to extract empirical single-molecule free energy surfaces directly from experimental measurements.

Nonequilibrium Relaxation of Conformational Dynamics Facilitates Catalytic Reaction in an Elastic Network Model of T7 DNA Polymerase.

Nucleotide-induced conformational closing of the finger domain of DNA polymerase is crucial for its catalytic action during DNA replication. Such large-amplitude molecular motion is often not fully accessible to either direct experimental monitoring or molecular dynamics simulations. However, a coarse-grained model can offer an informative alternative, especially for probing the relationship between conformational dynamics and catalysis. Here we investigate the dynamics of T7 DNA polymerase catalysis using a Langevin-type elastic network model incorporating detailed structural information on the open conformation without the substrate bound. Such a single-parameter model remarkably captures the induced conformational dynamics of DNA polymerase upon dNTP binding, and reveals its close coupling to the advancement toward transition state along the coordinate of the target reaction, which contributes to significant lowering of the activation energy barrier. Furthermore, analysis of stochastic catalytic rates suggests that when the activation energy barrier has already been significantly lowered and nonequilibrium relaxation toward the closed form dominates the catalytic rate, one must appeal to a picture of two-dimensional free energy surface in order to account for the full spectrum of catalytic modes. Our semiquantitative study illustrates the general role of conformational dynamics in achieving transition-state stabilization, and suggests that such an elastic network model, albeit simplified, possesses the potential to furnish significant mechanistic insights into the functioning of a variety of enzymatic systems.

Calculating Watson-Crick to Hoogsteen Transition Kinetics in DNA with Langevin Dynamics and Fokker-Planck Diffusion in Reduced Configuration Space

Previous studies have shown that an important aspect of designing drugs for RNA drug targets, in particular, the transactivation response element (TAR) from HIV type 1 (HIV-1), is the determination of the thermodynamics and kinetics for WC to HG (Watson-Crick to Hoogsteen) conformation transitions. Previous work by Eunae Kim and Ioan Andricioaei successfully implemented umbrella sampling to produce a two-dimensional free energy surface for WC to HG transitions, where the two independent variables are the “flip-over” angle and “flip-out” angle of the adenine moiety undergoing the transition. The current study builds upon the work of Kim and Andricioaei by exploring the kinetics that arise from the time evolution of the probability density function of configurations, subjected to this potential energy surface. A position-dependent diffusion matrix was first calculated from the umbrella sampling data, in order to accurately describe diffusive motions in this reduced space. Next, a 2D polynomial of order 20 was fit to the energy surface, in order to obtain an analytical expression for the 2D surface. The diffusion matrix and analytical function describing the 2D surface were then used to carry out Langevin dynamics. Numerous Langevin trajectories between states were then analyzed in multiple ways, including calculating first passage time distributions between WC and HG configurations, and determination of regions in configuration space of maximum flux. In order to make the system numerically tractable for a Fokker-Planck approach to diffusion, the energy surface was also fit to a simpler 2D analytical function with a 20 term Fourier series solely in the “flip-over” dimension and a separate harmonic term approximating the “flip-out” dimension. Both approaches were utilized to determine favorable pathways for WC/HG transitions, and their respective time scales.