Nonconverged Constraints Cause Artificial Temperature Gradients in Lipid Bilayer Simulations

Pathways to a Shiny Future: Building the Foundation for Computational Physical Chemistry and Biophysics in 2050

In the last quarter-century, the field of molecular dynamics (MD) has undergone a remarkable transformation, propelled by substantial enhancements in software, hardware, and underlying methodologies. In this Perspective, we contemplate the future trajectory of MD simulations and their possible look at the year 2050. We spotlight the pivotal role of artificial intelligence (AI) in shaping the future of MD and the broader field of computational physical chemistry. We outline critical strategies and initiatives that are essential for the seamless integration of such technologies. Our discussion delves into topics like multiscale modeling, adept management of ever-increasing data deluge, the establishment of centralized simulation databases, and the autonomous refinement, cross-validation, and self-expansion of these repositories. The successful implementation of these advancements requires scientific transparency, a cautiously optimistic approach to interpreting AI-driven simulations and their analysis, and a mindset that prioritizes knowledge-motivated research alongside AI-enhanced big data exploration. While history reminds us that the trajectory of technological progress can be unpredictable, this Perspective offers guidance on preparedness and proactive measures, aiming to steer future advancements in the most beneficial and successful direction.

Building complex membranes with Martini 3

The Martini model is a popular force field for coarse-grained simulations. Membranes have always been at the center of its development, with the latest version, Martini 3, showing great promise in capturing more and more realistic behavior. In this chapter we provide a step-by-step tutorial on how to construct starting configurations, run initial simulations and perform dedicated analysis for membrane-based systems of increasing complexity, including leaflet asymmetry, curvature gradients and embedding of membrane proteins.

Insertases scramble lipids: Molecular simulations of MTCH2

Scramblases play a pivotal role in facilitating bidirectional lipid transport across cell membranes, thereby influencing lipid metabolism, membrane homeostasis, and cellular signaling. MTCH2, a mitochondrial outer membrane protein insertase, has a membrane-spanning hydrophilic groove resembling those that form the lipid transit pathway in known scramblases. Employing both coarse-grained and atomistic molecular dynamics simulations, we show that MTCH2 significantly reduces the free energy barrier for lipid movement along the groove and therefore can indeed function as a scramblase. Notably, the scrambling rate of MTCH2 in silico is similar to that of voltage-dependent anion channel (VDAC), a recently discovered scramblase of the outer mitochondrial membrane, suggesting a potential complementary physiological role for these mitochondrial proteins. Finally, our findings suggest that other insertases which possess a hydrophilic path across the membrane like MTCH2, can also function as scramblases.

Analyzing lipid distributions and curvature in molecular dynamics simulations of complex membranes

We describe methods to analyze lipid distributions and curvature in membranes with complex lipid mixtures and embedded membrane proteins. We discuss issues involved in these analyses, available tools to calculate curvature preferences of lipids and proteins, and focus on tools developed in our group for visual analysis of lipid-protein interactions and the analysis of membrane curvature.

Assessment of the MARTINI 3 Performance for Short Peptide Self-Assembly

The coarse-grained MARTINI force field, initially developed for membranes, has proven to be an exceptional tool for investigating supramolecular peptide assemblies. Over the years, the force field underwent refinements to enhance accuracy, enabling, for example, the reproduction of protein-ligand interactions and constant pH behavior. However, these protein-focused improvements seem to have compromised its ability to model short peptide self-assembly. In this study, we assess the performance of MARTINI 3 in reproducing peptide self-assembly using the well-established diphenylalanine (FF) as our test case. Unlike its success in version 2.1, FF does not even exhibit aggregation in version 3. By systematically exploring parameters for the aromatic side chains and charged backbone beads, we established a parameter set that effectively reproduces tube formation. Remarkably, these parameter adjustments also replicate the self-assembly of other di- and tripeptides and coassemblies. Furthermore, our analysis uncovers pivotal insights for enhancing the performance of MARTINI in modeling short peptide self-assembly. Specifically, we identify issues stemming from overestimated hydrophilicity arising from charged termini and disruptions in π-stacking interactions due to insufficient planarity in aromatic groups and a discrepancy in intermolecular distances between this and backbone-backbone interactions. This investigation demonstrates that strategic modifications can harness the advancements offered by MARTINI 3 for the realm of short peptide self-assembly.

Nonuniversal impact of cholesterol on membranes mobility, curvature sensing and elasticity

Biological membranes, composed mainly of phospholipids and cholesterol, play a vital role as cellular barriers. They undergo localized reshaping in response to environmental cues and protein interactions, with the energetics of deformations crucial for exerting biological functions. This study investigates the non-universal role of cholesterol on the structure and elasticity of saturated and unsaturated lipid membranes. Our study uncovers a highly cooperative relationship between thermal membrane bending and local cholesterol redistribution, with cholesterol showing a strong preference for the compressed membrane leaflet. Remarkably, in unsaturated membranes, increased cholesterol mobility enhances cooperativity, resulting in membrane softening despite membrane thickening and lipid compression caused by cholesterol. These findings elucidate the intricate interplay between thermodynamic forces and local molecular interactions that govern collective properties of membranes.

Phospholipids are imported into mitochondria by VDAC, a dimeric beta barrel scramblase

Mitochondria are double-membrane-bounded organelles that depend critically on phospholipids supplied by the endoplasmic reticulum. These lipids must cross the outer membrane to support mitochondrial function, but how they do this is unclear. We identify the Voltage Dependent Anion Channel (VDAC), an abundant outer membrane protein, as a scramblase-type lipid transporter that catalyzes lipid entry. On reconstitution into membrane vesicles, dimers of human VDAC1 and VDAC2 catalyze rapid transbilayer translocation of phospholipids by a mechanism that is unrelated to their channel activity. Coarse-grained molecular dynamics simulations of VDAC1 reveal that lipid scrambling occurs at a specific dimer interface where polar residues induce large water defects and bilayer thinning. The rate of phospholipid import into yeast mitochondria is an order of magnitude lower in the absence of VDAC homologs, indicating that VDACs provide the main pathway for lipid entry. Thus, VDAC isoforms, members of a superfamily of beta barrel proteins, moonlight as a class of phospholipid scramblases - distinct from alpha-helical scramblase proteins - that act to import lipids into mitochondria.

Insertases Scramble Lipids: Molecular Simulations of MTCH2

Scramblases play a pivotal role in facilitating bidirectional lipid transport across cell membranes, thereby influencing lipid metabolism, membrane homeostasis, and cellular signaling. MTCH2, a mitochondrial outer membrane protein insertase, has a membrane-spanning hydrophilic groove resembling those that form the lipid transit pathway in known scramblases. Employing both coarse-grained and atomistic molecular dynamics simulations, we show that MTCH2 significantly reduces the free energy barrier for lipid movement along the groove and therefore can indeed function as a scramblase. Notably, the scrambling rate of MTCH2 in silico is similar to that of VDAC, a recently discovered scramblase of the outer mitochondrial membrane, suggesting a potential complementary physiological role for these mitochondrial proteins. Finally, our findings suggest that other insertases which possess a hydrophilic path across the membrane like MTCH2, can also function as scramblases.

CGCompiler: Automated Coarse-Grained Molecule Parametrization via Noise-Resistant Mixed-Variable Optimization

Coarse-grained force fields (CG FFs) such as the Martini model entail a predefined, fixed set of Lennard-Jones parameters (building blocks) to model virtually all possible nonbonded interactions between chemically relevant molecules. Owing to its universality and transferability, the building-block coarse-grained approach has gained tremendous popularity over the past decade. The parametrization of molecules can be highly complex and often involves the selection and fine-tuning of a large number of parameters (e.g., bead types and bond lengths) to optimally match multiple relevant targets simultaneously. The parametrization of a molecule within the building-block CG approach is a mixed-variable optimization problem: the nonbonded interactions are discrete variables, whereas the bonded interactions are continuous variables. Here, we pioneer the utility of mixed-variable particle swarm optimization in automatically parametrizing molecules within the Martini 3 coarse-grained force field by matching both structural (e.g., RDFs) as well as thermodynamic data (phase-transition temperatures). For the sake of demonstration, we parametrize the linker of the lipid sphingomyelin. The important advantage of our approach is that both bonded and nonbonded interactions are simultaneously optimized while conserving the search efficiency of vector guided particle swarm optimization (PSO) methods over other metaheuristic search methods such as genetic algorithms. In addition, we explore noise-mitigation strategies in matching the phase-transition temperatures of lipid membranes, where nucleation and concomitant hysteresis introduce a dominant noise term within the objective function. We propose that noise-resistant mixed-variable PSO methods can both improve and automate parametrization of molecules within building-block CG FFs, such as Martini.

LipIDens: simulation assisted interpretation of lipid densities in cryo-EM structures of membrane proteins

Cryo-electron microscopy (cryo-EM) enables the determination of membrane protein structures in native-like environments. Characterising how membrane proteins interact with the surrounding membrane lipid environment is assisted by resolution of lipid-like densities visible in cryo-EM maps. Nevertheless, establishing the molecular identity of putative lipid and/or detergent densities remains challenging. Here we present LipIDens, a pipeline for molecular dynamics (MD) simulation-assisted interpretation of lipid and lipid-like densities in cryo-EM structures. The pipeline integrates the implementation and analysis of multi-scale MD simulations for identification, ranking and refinement of lipid binding poses which superpose onto cryo-EM map densities. Thus, LipIDens enables direct integration of experimental and computational structural approaches to facilitate the interpretation of lipid-like cryo-EM densities and to reveal the molecular identities of protein-lipid interactions within a bilayer environment. We demonstrate this by application of our open-source LipIDens code to ten diverse membrane protein structures which exhibit lipid-like densities.

A high-throughput structural dynamics approach for identification of potential agonists of FFAR4 for type 2 diabetes mellitus therapy

Diabetes mellitus is a metabolic disorder that persists as a global threat to the world. A G-protein coupled receptor (GPCR), free fatty acid receptor 4 (FFAR4), has emerged as a potential target for type 2 diabetes mellitus (T2DM) and obesity-related disorders. The current study has investigated the FFAR4, deploying 3-dimensional structure modeling, molecular docking, machine learning, and high-throughput virtual screening methods to unravel the receptor's crucial and non-crucial binding site residues. We screened four lakh compounds and shortlisted them based on binding energy, stereochemical considerations, non-bonded interactions, and pharmacokinetic profiling. Out of the screened compounds, four compounds were selected for ligand-bound simulations. The molecular dynamic simulations were carried out for 1µs for native FFAR4 and 500 ns each for complexes of FFAR4 with compound 1, compound 2, compound 3, and compound 4. Our findings showed that in addition to reported binding site residues ARG99, ARG183, and VAL98 in known agonists like TUG-891, the amino acids ARG22, ARG24, THR23, TRP305, and GLU43 were also critical binding site residues. These amino acids impart stability to the FFAR4 complexes and contribute to the stronger binding affinity of the compounds. The study also indicated that aromatic residues like PHE211 are crucial for recognizing the active site's pi-pi and C-C double bonds. Since FFAR4 is a membrane protein, the simulation studies give an insight into the mechanisms of the crucial protein-lipid and lipid-water interactions. The analysis of the molecular dynamics trajectories showed all four compounds as potential hit molecules that can be developed further into potential agonists for T2DM therapy. Amongst the four compounds, compound 4 showed relatively better binding affinity, stronger non-bonded interactions, and a stable complex.Communicated by Ramaswamy H. Sarma.

Investigation of nano- and microdomains formed by ceramide 1 phosphate in lipid bilayers

Biological membranes are renowned for their intricate complexity, with the formation of membrane domains being pivotal to the successful execution of numerous cellular processes. However, due to their nanoscale characteristics, these domains are often understudied, as the experimental techniques required for quantitative investigation present significant challenges. In this study we employ spot-variation z-scan fluorescence correlation spectroscopy (svzFCS) tailored for artificial lipid vesicles of varying composition and combine this approach with high-resolution imaging. This method has been harnessed to examine the lipid-segregation behavior of distinct types of ceramide-1-phosphate (C1P), a crucial class of signaling molecules, within these membranes. Moreover, we provide a quantitative portrayal of the lipid membranes studied and the domains induced by C1P at both nano and microscales. Given the lack of definitive conclusions from the experimental data obtained, it was supplemented with comprehensive in silico studies-including the analysis of diffusion coefficient via molecular dynamics and domain populations via Monte Carlo simulations. This approach enhanced our insight into the dynamic behavior of these molecules within model lipid membranes, confirming that nano- and microdomains can co-exist in lipid vesicles.

Martini 3 Coarse-Grained Force Field for Cholesterol

Cholesterol plays a crucial role in biomembranes by regulating various properties, such as fluidity, rigidity, permeability, and organization of lipid bilayers. The latest version of the Martini model, Martini 3, offers significant improvements in interaction balance, molecular packing, and inclusion of new bead types and sizes. However, the release of the new model resulted in the need to reparameterize many core molecules, including cholesterol. Here, we describe the development and validation of a Martini 3 cholesterol model, addressing issues related to its bonded setup, shape, volume, and hydrophobicity. The proposed model mitigates some limitations of its Martini 2 predecessor while maintaining or improving the overall behavior.

An implementation of the Martini coarse-grained force field in OpenMM

We describe a complete implementation of Martini 2 and Martini 3 in the OpenMM molecular dynamics software package. Martini is a widely used coarse-grained force field with applications in biomolecular simulation, materials, and broader areas of chemistry. It is implemented as a force field but makes extensive use of facilities unique to the GROMACS software, including virtual sites and bonded terms that are not commonly used in standard atomistic force fields. OpenMM is a flexible molecular dynamics package widely used for methods development and is competitive in speed on GPUs with other commonly used packages. OpenMM has facilities to easily implement new force field terms, external forces and fields, and other nonstandard features, which we use to implement all force field terms used in Martini 2 and Martini 3. This allows Martini simulations, starting with GROMACS topology files that are processed by custom scripts, with all the added flexibility of OpenMM. We provide a GitHub repository with test cases, compare accuracy and performance between GROMACS and OpenMM, and discuss the limitations of our implementation in terms of direct comparison with GROMACS. We describe a use case that implements the Modeling Employing Limited Data method to apply experimental constraints in a Martini simulation to efficiently determine the structure of a protein complex. We also discuss issues and a potential solution with the Martini 2 topology for cholesterol.

Crosstalk between regulatory elements in disordered TRPV4 N-terminus modulates lipid-dependent channel activity

Intrinsically disordered regions (IDRs) are essential for membrane receptor regulation but often remain unresolved in structural studies. TRPV4, a member of the TRP vanilloid channel family involved in thermo- and osmosensation, has a large N-terminal IDR of approximately 150 amino acids. With an integrated structural biology approach, we analyze the structural ensemble of the TRPV4 IDR and the network of antagonistic regulatory elements it encodes. These modulate channel activity in a hierarchical lipid-dependent manner through transient long-range interactions. A highly conserved autoinhibitory patch acts as a master regulator by competing with PIP2 binding to attenuate channel activity. Molecular dynamics simulations show that loss of the interaction between the PIP2-binding site and the membrane reduces the force exerted by the IDR on the structured core of TRPV4. This work demonstrates that IDR structural dynamics are coupled to TRPV4 activity and highlights the importance of IDRs for TRP channel function and regulation.

Partitioning into phosphatidylcholine-cholesterol membranes: liposome measurements, coarse-grained simulations, and implications for bioaccumulation

Membrane-water partitioning is an important physical property for the assessment of bioaccumulation and environmental impact. Here, we advance simulation methodology for predicting the partitioning of small molecules into lipid membranes and compare the computational predictions to experimental measurements in liposomes. As a step towards high-throughput screening, we present an automated mapping and parametrization procedure to produce coarse-grained models compatible with the Martini 3 force field. The methodology is general and can also be used for other applications where coarse-grained simulations are appropriate. This article addresses the effect on membrane-water partitioning of adding cholesterol to POPC (1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine) membranes. Nine contrasting neutral, zwitterionic and charged solutes are tested. Agreement between experiment and simulation is generally good, with the most challenging cases being permanently charged solutes. For all solutes, partitioning is found to be insensitive to membrane cholesterol concentration up to 25% mole fraction. Hence, for assessment of bioaccumulation into a range of membranes (such as those found in fish), partitioning data measured in pure lipid membranes are still informative.

Optimal Bond Constraint Topology for Molecular Dynamics Simulations of Cholesterol

We recently observed artificial temperature gradients in molecular dynamics (MD) simulations of phase-separating ternary lipid mixtures using the Martini 2 force field. We traced this artifact to insufficiently converged bond length constraints with typical time steps and default settings for the linear constraint solver (LINCS). Here, we systematically optimize the constraint scaffold of cholesterol. With massive virtual sites in an equimomental arrangement, we accelerate bond constraint convergence while preserving the original cholesterol force field and dynamics. The optimized model does not induce nonphysical temperature gradients even at relaxed LINCS settings and is at least as fast as the original model at the strict LINCS settings required for proper thermal sampling. We provide a python script to diagnose possible problems with constraint convergence for other molecules and force fields. Equimomental constraint topology optimization can also be used to boost constraint convergence in atomistic MD simulations of molecular systems.

How Does the Study MD of pH-Dependent Exposure of Nanoparticles Affect Cellular Uptake of Anticancer Drugs?

The lack of knowledge about the uptake of NPs by biological cells poses a significant problem for drug delivery. For this reason, designing an appropriate model is the main challenge for modelers. To address this problem, molecular modeling studies that can describe the mechanism of cellular uptake of drug-loaded nanoparticles have been conducted in recent decades. In this context, we developed three different models for the amphipathic nature of drug-loaded nanoparticles (MTX-SS-γ-PGA), whose cellular uptake mechanism was predicted by molecular dynamics studies. Many factors affect nanoparticle uptake, including nanoparticle physicochemical properties, protein-particle interactions, and subsequent agglomeration, diffusion, and sedimentation. Therefore, the scientific community needs to understand how these factors can be controlled and the NP uptake of nanoparticles. Based on these considerations, in this study, we investigated for the first time the effects of the selected physicochemical properties of the anticancer drug methotrexate (MTX) grafted with hydrophilic-γ-polyglutamic acid (MTX-SS-γ-PGA) on its cellular uptake at different pH values. To answer this question, we developed three theoretical models describing drug-loaded nanoparticles (MTX-SS-γ-PGA) at three different pH values, such as (1) pH 7.0 (the so-called neutral pH model), (2) pH 6.4 (the so-called tumor pH model), and (3) pH 2.0 (the so-called stomach pH model). Exceptionally, the electron density profile shows that the tumor model interacts more strongly with the head groups of the lipid bilayer than the other models due to charge fluctuations. Hydrogen bonding and RDF analyses provide information about the solution of the NPs with water and their interaction with the lipid bilayer. Finally, dipole moment and HOMO-LUMO analysis showed the free energy of the solution in the water phase and chemical reactivity, which are particularly useful for determining the cellular uptake of the NPs. The proposed study provides fundamental insights into molecular dynamics (MD) that will allow researchers to determine the influence of pH, structure, charge, and energetics of NPs on the cellular uptake of anticancer drugs. We believe that our current study will be useful in developing a new model for drug delivery to cancer cells with a much more efficient and less time-consuming model.

Reorganization of the outer layer of a model of the plasma membrane induced by a neuroprotective aminosterol

Trodusquemine is an amphipathic aminosterol that has recently shown therapeutic benefit in neurodegenerative diseases altering the binding of misfolded proteins to the cell membrane. To unravel the underlying mechanism, we studied the interactions between Trodusquemine (TRO) and lipid monolayers simulating the outer layer of the plasma membrane. We selected two different compositions of dioleoylphosphatidylcholine (DOPC), sphingomyelin (SM), cholesterol (Chol) and monosialotetrahexosylganglioside (GM1) lipid mixture mimicking either a lipid-raft containing membrane (L_d+S_o phases) or a single-phase disordered membrane (L_d phase). Surface pressure-area isotherms and surface compressional modulus-area combined with Brewster Angle Microscopy (BAM) provided the thermodynamic and morphological information on the lipid monolayer in the presence of increasing amounts of TRO in the monolayer. Experiments revealed that TRO forms stable spreading monolayers at the buffer-air interface where it undergoes multiple reversible phase transitions to bi- and tri-layers at the interface. When TRO was spread at the interface with the lipid mixtures, we found that it distributes in the lipid monolayer for both the selected lipid compositions, but a maximum TRO uptake in the rafts-containing monolayer was observed for a Lipid/TRO molar ratio equal to 3:2. Statistical analysis of BAM images revealed that TRO induces a decrease in the size of the condensed domains, an increase in their number and in the thickness mismatch between the L_d and S_o phase. Experiments and MD simulations converge to indicate that TRO adsorbs preferentially at the border of the S_o domains. Removal of GM1 from the lipid L_d+S_o mixture resulted in an even greater TRO-mediated reduction of the size of the S_o domains suggesting that the presence of GM1 hinders the localization of TRO at the S_o domains boundaries. Taken together these observations suggest that Trodusquemine influences the organization of lipid rafts within the neuronal membrane in a dose-dependent manner whereas it evenly distributes in disordered expanded phases of the membrane model.

Nanoscale membrane curvature sorts lipid phases and alters lipid diffusion

The precise spatiotemporal control of nanoscale membrane shape and composition is the result of a complex interplay of individual and collective molecular behaviors. Here, we employed single-molecule localization microscopy and computational simulations to observe single-lipid diffusion and sorting in model membranes with varying compositions, phases, temperatures, and curvatures. Supported lipid bilayers were created over 50-nm-radius nanoparticles to mimic the size of naturally occurring membrane buds, such as endocytic pits and the formation of viral envelopes. The curved membranes recruited liquid-disordered lipid phases while altering the diffusion and sorting of tracer lipids. Disorder-preferring fluorescent lipids sorted to and experienced faster diffusion on the nanoscale curvature only when embedded in a membrane capable of sustaining lipid phase separation at low temperatures. The curvature-induced sorting and faster diffusion even occurred when the sample temperature was above the miscibility temperature of the planar membrane, implying that the nanoscale curvature could induce phase separation in otherwise homogeneous membranes. Further confirmation and understanding of these results are provided by continuum and coarse-grained molecular dynamics simulations with explicit and spontaneous curvature-phase coupling, respectively. The curvature-induced membrane compositional heterogeneity and altered dynamics were achieved only with a coupling of the curvature with a lipid phase separation. These cross-validating results demonstrate the complex interplay of lipid phases, molecular diffusion, and nanoscale membrane curvature that are critical for membrane functionality.

Small ionic radii limit time step in Martini 3 molecular dynamics simulations

Among other improvements, the Martini~3 coarse-grained force field provides a more accurate description of the solvation of protein pockets and channels through the consistent use of various bead types and sizes. Here, we show that the representation of Na+ and Cl- ions as ``tiny' (TQ5) beads limits the accessible time step to 25~fs. By contrast, with Martini~2, time steps of 30-40~fs were possible for lipid bilayer systems without proteins. This limitation is relevant for systems that require long equilibration times. We derive a quantitative kinetic model of time-integration instabilities in molecular dynamics as a function of time step, ion concentration and mass, system size, and simulation time. We demonstrate that ion-water interactions are the main source of instability at physiological conditions, followed closely by ion-ion interactions. We show that increasing the ionic masses makes it possible to use time steps up to 40~fs with minimal impact on static equilibrium properties and on dynamical quantities such as lipid and ion diffusion coefficients. Increasing the size of the bead representing the ions (and thus changing their hydration) also permits longer time steps. For a soluble protein, we find that increasing the mass of tiny beads also on the protein permits simulations with 30-fs time step. The use of larger time steps in Martini~3 results in a more efficient exploration of configuration space. The kinetic model of MD simulation crashes can be used to determine the maximum allowed time step upfront for an efficient use of resources and whenever sampling efficiency is critical.