QM/CG-MM: Systematic Embedding of Quantum Mechanical Systems in a Coarse-Grained Environment with Accurate Electrostatics.

Quantum Mechanics/Molecular Mechanics (QM/MM) can describe chemical reactions in molecular dynamics (MD) simulations at a much lower cost than ab initio MD. Still, it is prohibitively expensive for many systems of interest because such systems usually require long simulations for sufficient statistical sampling. Additional MM degrees of freedom are often slow and numerous but secondary in interest. Coarse-graining (CG) is well-known to be able to speed up sampling through both reduction in simulation cost and the ability to accelerate the dynamics. Therefore, embedding a QM system in a CG environment can be a promising way of expediting sampling without compromising the information about the QM subsystem. Sinitskiy and Voth first proposed the theory of Quantum Mechanics/Coarse-grained Molecular Mechanics (QM/CG-MM) with a bottom-up CG mapping. Mironenko and Voth subsequently introduced the DFT-QM/CG-MM formalism to couple a Density Functional Theory (DFT) treated QM system and to an apolar environment. Here, we present a more complete theory that addresses MM environments with significant polarity by explicitly accounting for the electrostatic coupling. We demonstrate our QM/CG-MM method with a chloride-methyl chloride SN2 reaction system in acetone, which is sensitive to solvent polarity. The method accurately recapitulates the potential of mean force for the substitution reaction, and the reaction barrier from the best model agrees with the atomistic simulations within sampling error. These models also have generalizability. In two other reactive systems that they have not been trained on, the QM/CG-MM model still achieves the same level of agreement with the atomistic QM/MM models. Finally, we show that in these examples the speed-up in the sampling is proportional to the acceleration of the rotational dynamics of the solvent in the CG system.

Mechanism of phosphate release from actin filaments.

After ATP-actin monomers assemble filaments, the ATP's [Formula: see text]-phosphate is hydrolyzedwithin seconds and dissociates over minutes. We used all-atom molecular dynamics simulations to sample the release of phosphate from filaments and study residues that gate release. Dissociation of phosphate from Mg2+ is rate limiting and associated with an energy barrier of 20 kcal/mol, consistent with experimental rates of phosphate release. Phosphate then diffuses within an internal cavity toward a gate formed by R177, as suggested in prior computational studies and cryo-EM structures. The gate is closed when R177 hydrogen bonds with N111 and is open when R177 forms a salt bridge with D179. Most of the time, interactions of R177 with other residues occlude the phosphate release pathway. Machine learning analysis reveals that the occluding interactions fluctuate rapidly, underscoring the secondary role of backdoor gate opening in Pi release, in contrast with the previous hypothesis that gate opening is the primary event.

Lipid Organization by the Caveolin-1 Complex.

Caveolins are lipid-binding proteins that can organize membrane remodeling and oligomerize into the 8S-complex. The CAV1 8S-complex comprises a disk-like structure, about 15nm in diameter, with a central beta barrel. Further oligomerization of 8S-complexes remodels the membrane into caveolae vessels, with a dependence on cholesterol concentration. However, the molecular mechanisms behind membrane remodeling and cholesterol filtering are still not understood. Performing atomistic Molecular Dynamics simulations in combination with advanced sampling techniques, we describe how the CAV1-8S complex bends the membrane and accumulates cholesterol. Here, our simulations show an enhancing effect by the palmitoylations of CAV1, and we predict that the CAV1-8S complex can extract cholesterol molecules from the lipid bilayer and accommodate them in its beta barrel. Through backmapping to the all-atom level we also conclude that the Martini v2 coarse-grained forcefield overestimates membrane bending, as the atomistic simulations exhibit only very localized bending.

Understanding dynamics in coarse-grained models. IV. Connection of fine-grained and coarse-grained dynamics with the Stokes-Einstein and Stokes-Einstein-Debye relations.

Applying an excess entropy scaling formalism to the coarse-grained (CG) dynamics of liquids, we discovered that missing rotational motions during the CG process are responsible for artificially accelerated CG dynamics. In the context of the dynamic representability between the fine-grained (FG) and CG dynamics, this work introduces the well-known Stokes-Einstein and Stokes-Einstein-Debye relations to unravel the rotational dynamics underlying FG trajectories, thereby allowing for an indirect evaluation of the effective rotations based only on the translational information at the reduced CG resolution. Since the representability issue in CG modeling limits a direct evaluation of the shear stress appearing in the Stokes-Einstein and Stokes-Einstein-Debye relations, we introduce a translational relaxation time as a proxy to employ these relations, and we demonstrate that these relations hold for the ambient conditions studied in our series of work. Additional theoretical links to our previous work are also established. First, we demonstrate that the effective hard sphere radius determined by the classical perturbation theory can approximate the complex hydrodynamic radius value reasonably well. Furthermore, we present a simple derivation of an excess entropy scaling relationship for viscosity by estimating the elliptical integral of molecules. In turn, since the translational and rotational motions at the FG level are correlated to each other, we conclude that the "entropy-free" CG diffusion only depends on the shape of the reference molecule. Our results and analyses impart an alternative way of recovering the FG diffusion from the CG description by coupling the translational and rotational motions at the hydrodynamic level.

Cracked actin filaments as mechanosensitive receptors.

Actin filament networks are exposed to mechanical stimuli, but the effect of strain on actin filament structure has not been well established in molecular detail. This is a critical gap in understanding because the activity of a variety of actin-binding proteins has recently been determined to be altered by actin filament strain. We therefore used all-atom molecular dynamics simulations to apply tensile strains to actin filaments and find that changes in actin subunit organization are minimal in mechanically strained, but intact, actin filaments. However, a conformational change disrupts the critical D-loop to W-loop connection between longitudinal neighboring subunits, which leads to a metastable cracked conformation of the actin filament whereby one protofilament is broken prior to filament severing. We propose that the metastable crack presents a force-activated binding site for actin regulatory factors that specifically associate with strained actin filaments. Through protein-protein docking simulations, we find that 43 evolutionarily diverse members of the dual zinc-finger-containing LIM-domain family, which localize to mechanically strained actin filaments, recognize two binding sites exposed at the cracked interface. Furthermore, through its interactions with the crack, LIM domains increase the length of time damaged filaments remain stable. Our findings propose a new molecular model for mechanosensitive binding to actin filaments.

The SARS-CoV-2 nucleoprotein associates with anionic lipid membranes.

Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) is a lipid-enveloped virus that acquires its lipid bilayer from the host cell it infects. SARS-CoV-2 can spread from cell to cell or from patient to patient by undergoing assembly and budding to form new virions. The assembly and budding of SARS-CoV-2 is mediated by several structural proteins known as envelope (E), membrane (M), nucleoprotein (N), and spike (S), which can form virus-like particles (VLPs) when co-expressed in mammalian cells. Assembly and budding of SARS-CoV-2 from the host ER-Golgi intermediate compartment is a critical step in the virus acquiring its lipid bilayer. To date, little information is available on how SARS-CoV-2 assembles and forms new viral particles from host membranes. In this study, we used several lipid binding assays and found the N protein can strongly associate with anionic lipids including phosphoinositides and phosphatidylserine. Moreover, we show lipid binding occurs in the N protein C-terminal domain, which is supported by extensive in silico analysis. We demonstrate anionic lipid binding occurs for both the free and the N oligomeric forms, suggesting N can associate with membranes in the nucleocapsid form. Based on these results, we present a lipid-dependent model based on in vitro, cellular, and in silico data for the recruitment of N to assembly sites in the lifecycle of SARS-CoV-2.

Making the cut: Multiscale simulation of membrane remodeling.

Biological membranes are dynamic heterogeneous materials, and their shape and organization are tightly coupled to the properties of the proteins in and around them. However, the length scales of lipid and protein dynamics are far below the size of membrane-bound organelles, much less an entire cell. Therefore, multiscale modeling approaches are often necessary to build a comprehensive picture of the interplay of these factors, and have provided critical insights into our understanding of membrane dynamics. Here, we review computational methods for studying membrane remodeling, as well as passive and active examples of protein-driven membrane remodeling. As the field advances towards the modeling of key aspects of organelles and whole cells - an increasingly accessible regime of study - we summarize here recent successes and offer comments on future trends.

Interaction of MRI Contrast Agent [Gd(DOTA)]- with Lipid Membranes: A Molecular Dynamics Study.

Contrast agents are important imaging probes in clinical MRI, allowing the identification of anatomic changes that otherwise would not be possible. Intensive research on the development of new contrast agents is being made to image specific pathological markers or sense local biochemical changes. The most widely used MRI contrast agents are based on gadolinium(III) complexes. Due to their very high charge density, they have low permeability through tight biological barriers such as the blood-brain barrier, hampering their application in the diagnosis of neurological disorders. In this study, we explore the interaction between the widely used contrast agent [Gd(DOTA)]- (Dotarem) and POPC lipid bilayers by means of molecular dynamics simulations. This metal complex is a standard reference where several chemical modifications have been introduced to improve key properties such as bioavailability and targeting. The simulations unveil detailed insights into the agent's interaction with the lipid bilayer, offering perspectives beyond experimental methods. Various properties, including the impact on global and local bilayer properties, were analyzed. As expected, the results indicate a low partition coefficient (KP) and high permeation barrier for this reference compound. Nevertheless, favorable interactions are established with the membrane leading to moderately long residence times. While coordination of one inner-sphere water molecule is maintained for the membrane-associated chelate, the physical-chemical attributes of [Gd(DOTA)]- as a MRI contrast agent are affected. Namely, increases in the rotational correlation times and in the residence time of the inner-sphere water are observed, with the former expected to significantly increase the water proton relaxivity. This work establishes a reference framework for the use of simulations to guide the rational design of new contrast agents with improved relaxivity and bioavailability and for the development of liposome-based formulations for use as imaging probes or theranostic agents.

Molecular Dynamics Simulation of Complex Reactivity with the Rapid Approach for Proton Transport and Other Reactions (RAPTOR) Software Package.

Simulating chemically reactive phenomena such as proton transport on nanosecond to microsecond and beyond time scales is a challenging task. Ab initio methods are unable to currently access these time scales routinely, and traditional molecular dynamics methods feature fixed bonding arrangements that cannot account for changes in the system's bonding topology. The Multiscale Reactive Molecular Dynamics (MS-RMD) method, as implemented in the Rapid Approach for Proton Transport and Other Reactions (RAPTOR) software package for the LAMMPS molecular dynamics code, offers a method to routinely sample longer time scale reactive simulation data with statistical precision. RAPTOR may also be interfaced with enhanced sampling methods to drive simulations toward the analysis of reactive rare events, and a number of collective variables (CVs) have been developed to facilitate this. Key advances to this methodology, including GPU acceleration efforts and novel CVs to model water wire formation are reviewed, along with recent applications of the method which demonstrate its versatility and robustness.

Prediction of the essential intermolecular contacts for side-binding of VASP on F-actin.

Vasodilator-stimulated phosphoprotein (VASP) family proteins play a crucial role in mediating the actin network architecture in the cytoskeleton. The Ena/VASP homology 2 (EVH2) domain in each of the four identical arms of the tetrameric VASP consists of a loading poly-Pro region, a G-actin-binding domain (GAB), and an F-actin-binding domain (FAB). Together, the poly-Pro, GAB, and FAB domains allow VASP to bind to sides of actin filaments in a bundle, and recruit profilin-G-actin to processively elongate the filaments. The atomic resolution structure of the ternary complex, consisting of the loading poly-Pro region and GAB domain of VASP with profilin-actin, has been solved over a decade ago; however, a detailed structure of the FAB-F-actin complex has not been resolved to date. Experimental insights, based on homology of the FAB domain with the C region of WASP, have been used to hypothesize that the FAB domain binds to the cleft between subdomains 1 and 3 of F-actin. Here, in order to develop our understanding of the VASP-actin complex, we first augment known structural information about the GAB domain binding to actin with the missing FAB domain-actin structure, which we predict using homology modeling and docking simulations. In earlier work, we used mutagenesis and kinetic modeling to study the role of domain-level binding-unbinding kinetics of Ena/VASP on actin filaments in a bundle, specifically on the side of actin filaments. We further look at the nature of the side-binding of the FAB domain of VASP at the atomistic level using our predicted structure, and tabulate effective mutation sites on the FAB domain that would disrupt the VASP-actin complex. We test the binding affinity of Ena with mutated FAB domain using total internal reflection fluorescence microscopy experiments. The binding affinity of VASP is affected significantly for the mutant, providing additional support for our predicted structure.

Molecular Tuning of Reactivity of Zeolite Protons in HZSM-5.

In acidic HZSM-5 zeolite, the reactivity of a methanol molecule interacting with the zeolite proton is amenable to modification via coadsorbing a stochiometric amount of an electron density donor E to form the [(E)(CH3OH)(HZ)] complex. The rate of the methanol in this complex undergoing dehydration to dimethyl ether was determined for a series of E with proton affinity (PA) ranging from 659 kJ mol-1 for C6F6 to 825 kJ mol-1 for C4H8O and was found to follow the expression: Ln(Rate) - Ln(RateN2) = ß(PA - PAN2)γ, where E = N2 is the reference and ß and γ are constants. This trend is probably due to the increased stability of the solvated proton in the [(E)(CH3OH)(HZ)] complex with increasing PA. Importantly, this is also observed in steady-state flow reactions when stoichiometric quantities of E are preadsorbed on the zeolite. As demonstrated with E being D2O, the effect on methanol reactivity diminishes when E is present in excess of the [(E)(CH3OH)(HZ)] complex. It is proposed that the methanol dehydration reaction involves [(E)(CH3OH)(CH3OH)(HZ)] as the transition state, which is supported by the isotopologue distribution of the initial dimethyl ether formed when a flow of CH3OH was passed over ZSM-5 containing one CD3OH per zeolite proton. The implication of this on the mechanism of catalytic methanol dehydration on HZSM-5 is discussed.

Kinetic network modeling with molecular simulation inputs: A proton-coupled phosphate symporter.

Phosphate, an essential metabolite involved in numerous cellular functions, is taken up by proton-coupled phosphate transporters of plants and fungi within the major facilitator family. Similar phosphate transporters have been identified across a diverse range of biological entities, including various protozoan parasites linked to human diseases, breast cancer cells with increased phosphate requirements, and osteoclast-like cells engaged in bone resorption. Prior studies have proposed an overview of the functional cycle of a proton-driven phosphate transporter (PiPT), yet a comprehensive understanding of the proposed reaction pathways necessitates a closer examination of each elementary reaction step within an overall kinetic framework. In this work, we leverage kinetic network modeling in conjunction with a "bottom-up" molecular dynamics approach to show how such an approach can characterize the proton-phosphate co-transport behavior of PiPT under different pH and phosphate concentration conditions. In turn, this allows us to reveal the prevailing reaction pathway within a high-affinity phosphate transporter under different experimental conditions and to uncover the molecular origin of the optimal pH condition of this transporter.

Cooperative Membrane Binding of HIV-1 Matrix Proteins.

The HIV-1 assembly process begins with a newly synthesized Gag polyprotein being targeted to the inner leaflet of the plasma membrane of the infected cells to form immature viral particles. Gag-membrane interactions are mediated through the myristoylated (Myr) N-terminal matrix (MA) domain of Gag, which eventually multimerize on the membrane to form trimers and higher order oligomers. The study of the structure and dynamics of peripheral membrane proteins like MA has been challenging for both experimental and computational studies due to the complex transient dynamics of protein-membrane interactions. Although the roles of anionic phospholipids (PIP2, PS) and the Myr group in the membrane targeting and stable membrane binding of MA are now well-established, the cooperative interactions between the MA monomers and MA-membrane remain elusive in the context of viral assembly and release. Our present study focuses on the membrane binding dynamics of a higher order oligomeric structure of MA protein (a dimer of trimers), which has not been explored before. Employing time-lagged independent component analysis (tICA) to our microsecond-long trajectories, we investigate conformational changes of the matrix protein induced by membrane binding. Interestingly, the Myr switch of an MA monomer correlates with the conformational switch of adjacent monomers in the same trimer. Together, our findings suggest complex protein dynamics during the formation of the immature HIV-1 lattice; while MA trimerization facilitates Myr insertion, MA trimer-trimer interactions in the immature lattice can hinder the same.

Molecular dynamics simulations of HIV-1 matrix-membrane interactions at different stages of viral maturation.

Although the structural rearrangement of the membrane-bound matrix (MA) protein trimers upon HIV-1 maturation has been reported, the consequences of MA maturation on the MA-lipid interactions are not well understood. Long-timescale molecular dynamics simulations of the MA multimeric assemblies of immature and mature virus particles with our realistic asymmetric membrane model have explored MA-lipid interactions and lateral organization of lipids around MA complexes. The number of stable MA-phosphatidylserine and MA-phosphatidylinositol 4,5-bisphosphate (PIP2) interactions at the trimeric interface of the mature MA complex is observed to be greater compared to that of the immature MA complex. Our simulations identified an alternative PIP2-binding site in the immature MA complex where the multivalent headgroup of a PIP2 lipid with a greater negative charge binds to multiple basic amino acid residues such as ARG3 residues of both the MA monomers at the trimeric interface and highly basic region (HBR) residues (LYS29, LYS31) of one of the MA monomers. Our enhanced sampling simulations have explored the conformational space of phospholipids at different binding sites of the trimer-trimer interface of MA complexes that are not accessible by conventional unbiased molecular dynamics. Unlike the immature MA complex, the 2' acyl tail of two PIP2 lipids at the trimeric interface of the mature MA complex is observed to sample stable binding pockets of MA consisting of helix-4 residues. Together, our results provide molecular-level insights into the interactions of MA trimeric complexes with membrane and different lipid conformations at the specific binding sites of MA protein before and after viral maturation.

On the Key Influence of Amino Acid Ionic Liquid Anions on CO2 Capture.

Amino acid ionic liquids (AAILs) are promising green materials for CO2 capture and conversion due to their large chemical structural tunability. However, the structural understanding of the AAILs underlying the CO2 reaction dynamics remains uncertain. Herein, we examine the steric effects of AAIL anions with various chemical structures on CO2 capture behavior. Based on ab initio free-energy sampling, we assess reaction mechanisms for carbamate formation via a two-step reaction pathway with a zwitterion intermediate undergoing dynamic proton transfer. Our results show that free-energy barriers for carbamate formation can be significantly reduced as the degree of steric hindrance of the anions decreases. Further analyses reveal that reduced steric hindrance of anions causes markedly stronger intermolecular interactions between zwitterion and anions, leading to an increased kinetically favorable intermolecular proton transfer for carbamate production. We also describe the correlation strength between intramolecular interactions within the zwitterion and intermolecular interactions between the zwitterion and anions. We conclude that the favored structural flexibility due to the less steric hindrance of the zwitterion leads to enhanced intermolecular interactions, facilitating proton transfer to nearby AAIL anions for carbamate formation. Our study provides invaluable insight into the influence of various degrees of steric hindrance of the AAIL anions governing CO2 chemisorption. These findings may aid in the design of optimal AAIL solvents for the CO2 capture process.

HIV-1 capsid shape, orientation, and entropic elasticity regulate translocation into the nuclear pore complex.

Nuclear import and uncoating of the viral capsid are critical steps in the HIV-1 life cycle that serve to transport and release genomic material into the nucleus. Viral core import involves translocating the HIV-1 capsid at the nuclear pore complex (NPC). Notably, the central channel of the NPC appears to often accommodate and allow passage of intact HIV-1 capsid, though mechanistic details of the process remain to be fully understood. Here, we investigate the molecular interactions that operate in concert between the HIV-1 capsid and the NPC that regulate capsid translocation through the central channel. To this end, we develop a "bottom-up" coarse-grained (CG) model of the human NPC from recently released cryo-electron tomography structure and then construct composite membrane-embedded CG NPC models. We find that successful translocation from the cytoplasmic side to the NPC central channel is contingent on the compatibility of the capsid morphology and channel dimension and the proper orientation of the capsid approach to the channel from the cytoplasmic side. The translocation dynamics is driven by maximizing the contacts between phenylalanine-glycine nucleoporins at the central channel and the capsid. For the docked intact capsids, structural analysis reveals correlated striated patterns of lattice disorder likely related to the intrinsic capsid elasticity. Uncondensed genomic material inside the docked capsid augments the overall lattice disorder of the capsid. Our results suggest that the intrinsic "elasticity" can also aid the capsid to adapt to the stress and remain structurally intact during translocation.

Mechanism of Phosphate Release from Actin Filaments.

After ATP-actin monomers assemble filaments, the γ-phosphate is hydrolyzed from ATP within seconds and dissociates from the filament over several minutes. We used all-atom well-tempered metadynamics molecular dynamics simulations to sample the release of phosphate from filaments along with unbiased molecular dynamics simulations to study residues that gate release. Dissociation of phosphate from Mg2+ is rate limiting and associated with an energy barrier of 20 kcal/mol, consistent with experimental rates of phosphate release. Phosphate then diffuses in an internal cavity toward a gate formed by R177 suggested in prior computational studies and cryo-EM structures. The gate is closed when R177 hydrogen bonds with N111 and is open when R177 forms a salt bridge with D179. Most of the time interactions of R177 with other residues occludes the phosphate release pathway. Machine learning analysis reveals that the occluding interactions fluctuate rapidly. These occluded states have not been documented in cryo-EM reconstructions.

Conformational transitions of the HIV-1 Gag polyprotein upon multimerization and gRNA binding.

During the HIV-1 assembly process, the Gag polyprotein multimerizes at the producer cell plasma membrane, resulting in the formation of spherical immature virus particles. Gag-genomic RNA (gRNA) interactions play a crucial role in the multimerization process, which is yet to be fully understood. We performed large-scale all-atom molecular dynamics simulations of membrane-bound full-length Gag dimer, hexamer, and 18-mer. The inter-domain dynamic correlation of Gag, quantified by the heterogeneous elastic network model applied to the simulated trajectories, is observed to be altered by implicit gRNA binding, as well as the multimerization state of the Gag. The lateral dynamics of our simulated membrane-bound Gag proteins, with and without gRNA binding, agree with prior experimental data and help to validate our simulation models and methods. The gRNA binding is observed to affect mainly the SP1 domain of the 18-mer and the matrix-capsid linker domain of the hexamer. In the absence of gRNA binding, the independent dynamical motion of the nucleocapsid domain results in a collapsed state of the dimeric Gag. Unlike stable SP1 helices in the six-helix bundle, without IP6 binding, the SP1 domain undergoes a spontaneous helix-to-coil transition in the dimeric Gag. Together, our findings reveal conformational switches of Gag at different stages of the multimerization process and predict that the gRNA binding reinforces an efficient binding surface of Gag for multimerization, and also regulates the dynamic organization of the local membrane region itself.

Can a coarse-grained water model capture the key physical features of the hydrophobic effect?

Coarse-grained (CG) molecular dynamics can be a powerful method for probing complex processes. However, most CG force fields use pairwise nonbonded interaction potentials sets, which can limit their ability to capture complex multi-body phenomena such as the hydrophobic effect. As the hydrophobic effect primarily manifests itself due to the nonpolar solute affecting the nearby hydrogen bonding network in water, capturing such effects using a simple one CG site or "bead" water model is a challenge. In this work, we systematically test the ability of CG one site water models for capturing critical features of the solvent environment around a hydrophobe as well as the potential of mean force (PMF) of neopentane association. We study two bottom-up models: a simple pairwise (SP) force-matched water model constructed using the multiscale coarse-graining method and the Bottom-Up Many-Body Projected Water (BUMPer) model, which has implicit three-body correlations. We also test the top-down monatomic (mW) and the Machine Learned mW (ML-mW) water models. The mW models perform well in capturing structural correlations but not the energetics of the PMF. BUMPer outperforms SP in capturing structural correlations and also gives an accurate PMF in contrast to the two mW models. Our study highlights the importance of including three-body interactions in CG water models, either explicitly or implicitly, while in general highlighting the applicability of bottom-up CG water models for studying hydrophobic effects in a quantitative fashion. This assertion comes with a caveat, however, regarding the accuracy of the enthalpy-entropy decomposition of the PMF of hydrophobe association.