Insight into Molecular Basis and Dynamics of Full-length CRaf Kinase in Cellular Signaling Mechanisms.

Raf kinases play key roles in signal transduction in cells for regulating proliferation, differentiation, and survival. Despite decades of research into functions and dynamics of Raf kinases with respect to other cytosolic proteins, understanding Raf kinases is limited by the lack of their full-length structures at the atomic resolution. Here, we present the first model of the full-length CRaf kinase obtained from Artificial Intelligence/Machine Learning (AI/ML) algorithms with a converging ensemble of structures simulated by large-scale temperature replica exchange simulations. Our model is validated by comparing simulated structures with the latest Cryo-EM structure detailing close contacts among three key domains and regions of the CRaf. Our simulations reveal potentially new epitopes of intra-molecule interactions within the CRaf and reveal a dynamical nature of CRaf kinases, in which the three domains can move back and forth relative to each other for regulatory dynamics. The dynamical conformations are then used in a docking algorithm to shed insight into the paradoxical effect caused by Vemurafenib in comparison with a paradox breaker PLX7904. We propose a model of Raf-heterodimer/KRas-dimer as a signalosome based on the dynamics of the full-length CRaf.

Millisecond molecular dynamics simulations of KRas-dimer formation and interfaces.

Ras dimers have been proposed as building blocks for initiating the extracellular signal-regulated kinase (ERK)/mitogen-activated protein kinase (MAPK) cellular signaling pathway. To better examine the structure of possible dimer interfaces, the dynamics of Ras dimerization, and its potential signaling consequences, we performed molecular dynamics simulations totaling 1 ms of sampling, using an all-atom model of two full-length, farnesylated, guanosine triphosphate (GTP)-bound, wild-type KRas4b proteins diffusing on 29%POPS (1-palmitoyl-2-oleoyl-sn-glycero-3-phospho-L-serine)-mixed POPC (1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine) membranes. Our simulations unveil an ensemble of thermodynamically weak KRas dimers spanning multiple conformations. The most stable conformations, having the largest interface areas, involve helix α2 and a hypervariable region (HVR). Among the dimer conformations, we found that the HVR of each KRas has frequent interactions with various parts of the dimer, thus potentially mediating the dimerization. Some dimer configurations have one KRas G-domain elevated above the lipid bilayer surface by residing on top of the other G-domain, thus likely contributing to the recruitment of cytosolic Raf kinases in the context of a stably formed multi-protein complex. We identified a variant of the α4-α5 KRas-dimer interface that is similar to the interfaces obtained with fluorescence resonance energy transfer (FRET) data of HRas on lipid bilayers. Interestingly, we found two arginine fingers, R68 and R149, that directly interact with the beta-phosphate of the GTP bound in KRas, in a manner similar to what is observed in a crystal structure of GAP-HRas complex, which can facilitate the GTP hydrolysis via the arginine finger of GTPase-activating protein (GAP).

The voltage-sensing domain of a hERG1 mutant is a cation-selective channel.

A cationic leak current known as an "omega current" may arise from mutations of the first charged residue in the S4 of the voltage sensor domains of sodium and potassium voltage-gated channels. The voltage-sensing domains (VSDs) in these mutated channels act as pores allowing nonspecific passage of cations, such as Li+, K+, Cs+, and guanidinium. Interestingly, no omega currents have been previously detected in the nonswapped voltage-gated potassium channels such as the human-ether-a-go-go-related (hERG1), hyperpolarization-activated cyclic nucleotide-gated, and ether-a-go-go channels. In this work, we discovered a novel omega current by mutating the first charged residue of the S4 of the hERG1, K525 to serine. To characterize this omega current, we used various probes, including the hERG1 pore domain blocker, dofetilide, to show that the omega current does not require cation flux via the canonical pore domain. In addition, the omega flux does not cross the conventional selectivity filter. We also show that the mutated channel (K525S hERG1) conducts guanidinium. These data are indicative of the formation of an omega current channel within the VSD. Using molecular dynamics simulations with replica-exchange umbrella sampling simulations of the wild-type hERG1 and the K525S hERG1, we explored the molecular underpinnings governing the cation flow in the VSD of the mutant. We also show that the wild-type hERG1 may form water crevices supported by the biophysical surface accessibility data. Overall, our multidisciplinary study demonstrates that the VSD of hERG1 may act as a cation-selective channel wherein a mutation of the first charged residue in the S4 generates an omega current. Our simulation uncovers the atomistic underpinning of this mechanism.

The Single Residue K12 Governs the Exceptional Voltage Sensitivity of Mitochondrial Voltage-Dependent Anion Channel Gating.

The voltage-dependent anion channel (VDAC) is a ß-barrel channel of the mitochondrial outer membrane (MOM) that passively transports ions, metabolites, polypeptides, and single-stranded DNA. VDAC responds to a transmembrane potential by "gating," i.e. transitioning to one of a variety of low-conducting states of unknown structure. The gated state results in nearly complete suppression of multivalent mitochondrial metabolite (such as ATP and ADP) transport, while enhancing calcium transport. Voltage gating is a universal property of ß-barrel channels, but VDAC gating is anomalously sensitive to transmembrane potential. Here, we show that a single residue in the pore interior, K12, is responsible for most of VDAC's voltage sensitivity. Using the analysis of over 40 µs of atomistic molecular dynamics (MD) simulations, we explore correlations between motions of charged residues inside the VDAC pore and geometric deformations of the ß-barrel. Residue K12 is bistable; its motions between two widely separated positions along the pore axis enhance the fluctuations of the ß-barrel and augment the likelihood of gating. Single channel electrophysiology of various K12 mutants reveals a dramatic reduction of the voltage-induced gating transitions. The crystal structure of the K12E mutant at a resolution of 2.6 Å indicates a similar architecture of the K12E mutant to the wild type; however, 60 µs of atomistic MD simulations using the K12E mutant show restricted motion of residue 12, due to enhanced connectivity with neighboring residues, and diminished amplitude of barrel motions. We conclude that ß-barrel fluctuations, governed particularly by residue K12, drive VDAC gating transitions.

Molecular determinants of pro-arrhythmia proclivity of d- and l-sotalol via a multi-scale modeling pipeline.

Drug isomers may differ in their proarrhythmia risk. An interesting example is the drug sotalol, an antiarrhythmic drug comprising d- and l- enantiomers that both block the hERG cardiac potassium channel and confer differing degrees of proarrhythmic risk. We developed a multi-scale in silico pipeline focusing on hERG channel - drug interactions and used it to probe and predict the mechanisms of pro-arrhythmia risks of the two enantiomers of sotalol. Molecular dynamics (MD) simulations predicted comparable hERG channel binding affinities for d- and l-sotalol, which were validated with electrophysiology experiments. MD derived thermodynamic and kinetic parameters were used to build multi-scale functional computational models of cardiac electrophysiology at the cell and tissue scales. Functional models were used to predict inactivated state binding affinities to recapitulate electrocardiogram (ECG) QT interval prolongation observed in clinical data. Our study demonstrates how modeling and simulation can be applied to predict drug effects from the atom to the rhythm for dl-sotalol and also increased proarrhythmia proclivity of d- vs. l-sotalol when accounting for stereospecific beta-adrenergic receptor blocking.

Mapping Ryanodine Binding Sites in the Pore Cavity of Ryanodine Receptors.

Ryanodine (Ryd) irreversibly targets ryanodine receptors (RyRs), a family of intracellular calcium release channels essential for many cellular processes ranging from muscle contraction to learning and memory. Little is known of the atomistic details about how Ryd binds to RyRs. In this study, we used all-atom molecular dynamics simulations with both enhanced and bidirectional sampling to gain direct insights into how Ryd interacts with major residues in RyRs that were experimentally determined to be critical for its binding. We found that the pyrrolic ring of Ryd displays preference for the R4892AGGG-F4921 residues in the cavity of RyR1, which explain the effects of the corresponding mutations in RyR2 in experiments. Particularly, the mutant Q4933A (or Q4863A in RyR2) critical for both the gating and Ryd binding not only has significantly less interaction with Ryd than the wild-type, but also yields more space for Ryd and water molecules in the cavity. These results describe clear binding modes of Ryd in the RyR cavity and offer structural mechanisms explaining functional data collected on RyR blockade.

Computational membrane biophysics: From ion channel interactions with drugs to cellular function.

The rapid development of experimental and computational techniques has changed fundamentally our understanding of cellular-membrane transport. The advent of powerful computers and refined force-fields for proteins, ions, and lipids has expanded the applicability of Molecular Dynamics (MD) simulations. A myriad of cellular responses is modulated through the binding of endogenous and exogenous ligands (e.g. neurotransmitters and drugs, respectively) to ion channels. Deciphering the thermodynamics and kinetics of the ligand binding processes to these membrane proteins is at the heart of modern drug development. The ever-increasing computational power has already provided insightful data on the thermodynamics and kinetics of drug-target interactions, free energies of solvation, and partitioning into lipid bilayers for drugs. This review aims to provide a brief summary about modeling approaches to map out crucial binding pathways with intermediate conformations and free-energy surfaces for drug-ion channel binding mechanisms that are responsible for multiple effects on cellular functions. We will discuss post-processing analysis of simulation-generated data, which are then transformed to kinetic models to better understand the molecular underpinning of the experimental observables under the influence of drugs or mutations in ion channels. This review highlights crucial mathematical frameworks and perspectives on bridging different well-established computational techniques to connect the dynamics and timescales from all-atom MD and free energy simulations of ion channels to the physiology of action potentials in cellular models. This article is part of a Special Issue entitled: Biophysics in Canada, edited by Lewis Kay, John Baenziger, Albert Berghuis and Peter Tieleman.

Stereoselective recognition of morphine enantiomers by µ-opioid receptor.

Stereospecific recognition of chiral molecules plays a crucial role in biological systems. The µ-opioid receptor (MOR) exhibits binding affinity towards (-)-morphine, a well-established gold standard in pain management, while it shows minimal binding affinity for the (+)-morphine enantiomer, resulting in a lack of analgesic activity. Understanding how MOR stereoselectively recognizes morphine enantiomers has remained a puzzle in neuroscience and pharmacology for over half-a-century due to the lack of direct observation techniques. To unravel this mystery, we constructed the binding and unbinding processes of morphine enantiomers with MOR via molecular dynamics simulations to investigate the thermodynamics and kinetics governing MOR's stereoselective recognition of morphine enantiomers. Our findings reveal that the binding of (-)-morphine stabilizes MOR in its activated state, exhibiting a deep energy well and a prolonged residence time. In contrast, (+)-morphine fails to sustain the activation state of MOR. Furthermore, the results suggest that specific residues, namely D1142.50 and D1473.32, are deprotonated in the active state of MOR bound to (-)-morphine. This work highlights that the selectivity in molecular recognition goes beyond binding affinities, extending into the realm of residence time.

Taming Conformational Heterogeneity on Ion Racetrack to Unveil Principles that Drive Membrane Permeation of Cyclosporines.

Our study reveals the underlying principles governing the passive membrane permeability in three large N-methylated macrocyclic peptides (N-MeMPs): cyclosporine A (CycA), Alisporivir (ALI), and cyclosporine H (CycH). We determine a series of conformers required for robust passive membrane diffusion and those relevant to other functions, such as binding to protein targets or intermediates, in the presence of solvent additives. We investigate the conformational interconversions and establish correlations with the membrane permeability. Nuclear magnetic resonance (NMR) and cyclic ion-mobility spectrometry-mass spectrometry (cIMS-MS) are employed to characterize conformational heterogeneity and identify cis-amides relevant for good membrane permeability. In addition, ion mobility selected cIMS-MS and infrared (IR) multiple-photon dissociation (IRMPD) spectroscopy experiments are conducted to evaluate the energy barriers between conformations. We observe that CycA and ALI, both cyclosporines with favorable membrane permeabilities, display multiple stable and well-defined conformers. In contrast, CycH, an epimer of CycA with limited permeability, exhibits fewer and fewer stable conformers. We demonstrate the essential role of the conformational shift from the aqueous cis MeVal11-MeBmt1 state (A1) to the closed conformation featuring cis MeLeu9-MeLeu10 (C1) in facilitating membrane permeation. Additionally, we highlight that the transition from A1 to the all-trans open conformation (O1) is specifically triggered by the presence of CaCl2. We also capture a set of conformers with cis Sar3-MeLeu4, MeLeu9-MeLeu10, denoted as I. Conformationally selected cIMS-MS and IRMPD data of [CycA+Ca]2+ show immediate repopulation of the original population distribution, suggesting that CaCl2 smooths out the energy barriers. Finally, our work presents an improved sampling molecular dynamics approach based on a refined force field that not only consistently and accurately captures established conformers of cyclosporines but also exhibits strong predictive capabilities for novel conformers.

Improving Estimation of the Koopman Operator with Kolmogorov-Smirnov Indicator Functions.

It has become common to perform kinetic analysis using approximate Koopman operators that transform high-dimensional timeseries of observables into ranked dynamical modes. The key to the practical success of the approach is the identification of a set of observables that form a good basis on which to expand the slow relaxation modes. Good observables are, however, difficult to identify a priori and suboptimal choices can lead to significant underestimations of characteristic time scales. Leveraging the representation of slow dynamics in terms of Hidden Markov Models (HMM), we propose a simple and computationally efficient clustering procedure to infer surrogate observables that form a good basis for slow modes. We apply the approach to an analytically solvable model system as well as on three protein systems of different complexities. We consistently demonstrate that the inferred indicator functions can significantly improve the estimation of the leading eigenvalues of Koopman operators and correctly identify key states and transition time scales of stochastic systems, even when good observables are not known a priori.

Correction: Identifying key determinants and dynamics of SARS-CoV-2/ACE2 tight interaction.

[This corrects the article DOI: 10.1371/journal.pone.0257905.].

Identifying key determinants and dynamics of SARS-CoV-2/ACE2 tight interaction.

SARS-CoV-2 virus, the causative agent of Covid-19, has fired up a global pandemic. The virus interacts with the human receptor angiotensin-converting enzyme 2 (ACE2) for an invasion via receptor binding domain (RBD) on its spike protein. To provide a deeper understanding of this interaction, we performed microsecond simulations of the RBD-ACE2 complex for SARS-CoV-2 and compared it with the closely related SARS-CoV discovered in 2003. We show residues in the RBD of SARS-CoV-2 that were mutated from SARS-CoV, collectively help make the RBD anchor much stronger to the N-terminal part of ACE2 than the corresponding residues on RBD of SARS-CoV. This would result in a reduced dissociation rate of SARS-CoV-2 from human receptor protein compared to SARS-CoV. The phenomenon was consistently observed in simulations beyond 500 ns and was reproducible across different force fields. Altogether, our study adds more insight into the critical dynamics of the key residues at the virus spike and human receptor binding interface and potentially aids the development of diagnostics and therapeutics to combat the pandemic efficiently.

How Anionic Lipids Affect Spatiotemporal Properties of KRAS4B on Model Membranes.

RAS proteins are small membrane-anchored GTPases that regulate key cellular signaling networks. It has been recently shown that different anionic lipid types can affect the spatiotemporal properties of RAS through dimerization/clustering and signaling fidelity. To understand the effects of anionic lipids on key spatiotemporal properties of RAS, we dissected 1 ms of data from all-atom molecular dynamics simulations for KRAS4B on two model anionic lipid membranes that have 30% of POPS mixed with neutral POPC and 8% of PIP2 mixed with POPC. We unveiled the orientation space of KRAS4B, whose kinetics were slower and more distinguishable on the membrane containing PIP2 than the membrane containing POPS. Particularly, the PIP2-mixed membrane can differentiate a third kinetic orientation state from the other two known orientation states. We observed that each orientation state may yield different binding modes with an RAF kinase, which is required for activating the MAPK/ERK signaling pathway. However, an overall occluded probability, for which RAF kinases cannot bind KRAS4B, remains unchanged on the two different membranes. We identified rare fast diffusion modes of KRAS4B that appear coupled with orientations exposed to cytosolic RAF. Particularly, on the membrane having PIP2, we found nonlinear correlations between the orientation states and the conformations of the cationic farnesylated hypervariable region, which acts as an anchor in the membrane. Using diffusion coefficients estimated from the all-atom simulations, we quantified the effect of PIP2 and POPS on the KRAS4B dimerization via Green's function reaction dynamics simulations, in which the averaged dimerization rate is 12.5% slower on PIP2-mixed membranes.

Molecular Mechanism of Conductance Enhancement in Narrow Cation-Selective Membrane Channels.

Membrane proteins known as ryanodine receptors (RyRs) display large conductance of ∼1 nS and nearly ideal charge selectivity. Both properties are inversely correlated in other large-conductance but nonselective biological nanopores (i.e., α-hemolysin) used as industrial biosensors. Although recent cryo-electron microscopy structures of RyR2 show similarities to K+- and Na+-selective channels, it remains unclear whether similar ion conduction mechanisms occur in RyR2. Here, we combine microseconds of all-atom molecular dynamics (MD) simulations with mutagenesis and electrophysiology experiments to investigate large K+ conductance and charge selectivity (cation vs anion) in an open-state structure of RyR2. Our results show that a water-mediated knock-on mechanism enhances the cation permeation. The polar Q4863 ring may function as a confinement zone amplifying charge selectivity, while the cytoplasmic vestibule can contribute to the efficiency of the cation attraction. We also provide direct evidence that the rings of acidic residues at the channel vestibules are critical for both conductance and charge discrimination in RyRs.

Improved QM/MM Linear-Interaction Energy Model for Substrate Recognition in Zinc-Containing Metalloenzymes.

One of the essential challenges in the description of receptor-drug interactions in the presence of various polyvalent cations (such as zinc, magnesium, or iron) is the accurate assessment of the electronic effects due to cofactor binding. The effects can range from partial electronic polarization of the proximal atoms in a receptor and bound substrate to long-range effects related to partial charge transfer and electronic delocalization effects between the cofactor and the drug. Here, we examine the role of the explicit account for electronic effects for a panel of small-molecule inhibitors binding to the zinc-aminopeptidase PfA-M1, an essential target for antimalarial drug development. Our study on PfA-M1:inhibitor interactions at the QM level reveals that the partial charge and proton transfer due to bound zinc ion are important mechanisms in the inhibitors' recognition and catalysis. The combination of classical MD simulations with a posteriori QM/MM corrections with novel DFTB parameters for the zinc cation and the linear-interaction energy (LIE) approach offers by far the most accurate estimates for the PfA-M1:inhibitor binding affinities, opening the door for future inhibitor design.

Estimation of Potentials of Mean Force from Nonequilibrium Pulling Simulations Using Both Minh-Adib Estimator and Weighted Histogram Analysis Method.

Nonequilibrium pulling simulations have been a useful approach for investigating a variety of physical and biological problems. The major target in the simulations is to reconstruct reliable potentials of mean force (PMFs) or unperturbed free-energy profiles for quantitatively addressing both equilibrium mechanistic properties and contributions from nonequilibrium processes. While several current nonequilibrium methods were shown to be accurate in computing free-energy profiles in systems with relatively simple dynamics, they have proved to be unsuitable in more complicated systems. To extend the applicability of nonequilibrium sampling, we demonstrate a novel method that combines Minh-Adib's bidirectional estimator with nonlinear WHAM equations to reconstruct and assess PMFs from relatively fast pulling trajectories. We test the method in a one-dimensional model system and in a system of an antibiotic gramicidin-A (gA) channel, which is considered a significant challenge for nonequilibrium sampling. We identify key parameters for efficiently performing pulling simulations to improve and ensure the convergence and accuracy of estimated PMFs. We show that a few pulling trajectories of a relatively fast pulling speed v = 10 Å/ns can return a fair estimate of the PMF of a single potassium ion in gA.

Is the G-quadruplex an effective nanoconductor for ions?

We use a stepwise pulling protocol in molecular dynamics simulations to identify how a G-quadruplex selects and conducts Na(+), K(+), and NH4(+) ions. By estimating the minimum free-energy changes of the ions along the central channel via Jarzynski's equality, we find that the G-quadruplex selectively binds the ionic species in the following order: K(+) > Na(+) > NH4(+). This order implies that K(+) optimally fits the channel. However, the features of the free-energy profiles indicate that the channel conducts Na(+) best. These findings are in fair agreement with experiments on G-quadruplexes and reveal a profoundly different behavior from the prototype potassium-ion channel KcsA, which selects and conducts the same ionic species. We further show that the channel can also conduct a single file of water molecules and deform to leak water molecules. We propose a range for the conductance of the G-quadruplex.

Parallel-pulling protocol for free-energy evaluation.

Jarzynski's equality (JE) allows us to compute free-energy differences from distributions of work. In molecular dynamics simulations, the traditional way of constructing work distributions is to perform as many pulling simulations as possible. But reliable work distributions are not always produced in a finite number of simulations. The computational cost of using JE is not less than other commonly used methods such as thermodynamic integration and umbrella sampling methods. Here we first show a different proof of JE based on the idea of stepwise pulling procedures that is efficient in computing free energies by using JE. The key point in our proof is that the processes of turning on or off a harmonic potential to perform work are described by double Heaviside functions of time. We then show that the distributions of work performed by the potential can be easily generated from the distributions of a reaction coordinate along a pathway. Based on the proof, we propose sequential and parallel stepwise pulling protocols for generating work distributions that require suitable relaxation time at each pulling step. The criterion for reliable work distributions is that there must be sufficient mutual overlaps between the adjacent distributions of the reaction coordinate along the pathway. We arrive at an alternative formula (besides JE) to compute free-energy differences from the averaged values of the reaction coordinate. The combination of JE and the alternative formula provides a viable way to determine the accuracy of computed free-energy differences. For the stretching of a deca-alanine molecule, our approach requires 21 parallel simulations and relaxation time as small as 0.4 ns for each simulation to estimate free-energy differences with an uncertainty of about 13%.

Demonstration of Jarzynski's equality in open quantum systems using a stepwise pulling protocol.

We present a generalization of Jarzynski's equality, applicable to quantum systems, that is related to discretized mechanical work and free-energy changes. The theory is based on a stepwise pulling protocol. We find that work distribution functions can be constructed from fluctuations of a reaction coordinate along a reaction pathway in the stepwise pulling protocol. We also propose two sets of equations to determine the two possible optimal pathways that provide the most significant contributions to free-energy changes. We find that the transitions along these most optimal pathways, satisfying both sets of equations, follow the principle of detailed balance. We then test the theory by explicitly computing the free-energy changes for a one-dimensional quantum harmonic oscillator. This approach suggests a feasible way of measuring the fluctuations to experimentally test Jarzynski's equality in many-body systems, such as Bose-Einstein condensates.

Molecular mechanism of flip-flop in triple-layer oleic-acid membrane: correlation between oleic acid and water.

We perform all-atom molecular dynamics simulations to study a pure oleic acid (OA) membrane in water that results in a triple-layer structure. We compute the pressure profiles to examine the hydrophobic and hydrophilic regions, and to estimate the surface tension (≈34.5 mN/m), which is similar to those of lipid membranes. We observe that the membrane of OAs having a large diffusion coefficient (0.4 × 10(-7) cm(2)/s) along the normal to the membrane is an ideal model to study oleic acid flip-flop. In the model, the membrane contains a middle layer serving as an intermediate for water and OAs to easily migrate (flip-flop) from one to other leaflets. Water molecules surrounding OA head-groups help to reduce the barriers at the hydrophobic interface to trigger flip-flop events. Within 500 ns, we observe 175 flip-flop events of OAs and 305 events of water traversing the membrane. The ratio of water passing rate (k(H(2)O) = 0.673 ns(-1)) to OA flip-flop rate (k(OA) = 0.446 ns(-1)) is 3/2. The ratio of the totally correlated water-OA events to the totally uncorrelated water-OA events, n(cor)/n(uncor), is also 3/2. The probability of the totally and partially correlated events is 69%. The results indicate that the trans-membrane movement of water and OAs is cooperative and correlated, and agrees with experimentally measured absorption rates. They support the idea that OA flip-flop is more favorable than transport by means of functional proteins. This study might provide further insight into how primitive cell membranes work, and how the interplay and correlation between water and fatty acids may occur.