Direct-Space Corrections Enable Fast and Accurate Lorentz-Berthelot Combination Rule Lennard-Jones Lattice Summation

Crystal structure of archaeal IF5A-DHS complex reveals insights into the hypusination mechanism

The translation factor IF5A is highly conserved in Eukarya and Archaea and undergoes a unique post-translational hypusine modification by the deoxyhypusine synthase (DHS) enzyme. DHS transfers the butylamine moiety from spermidine to IF5A using NAD as a cofactor, forming a deoxyhypusine intermediate. IF5A is a key player in protein synthesis, preventing ribosome stalling in proline-rich sequences during translation elongation and facilitating translation elongation and termination. Additionally, human eIF5A participates in various essential cellular processes and contributes to cancer metastasis, with inhibiting hypusination showing anti-proliferative effects. The hypusination pathway of IF5A is therefore an attractive new therapeutic target. We elucidated the 2.0 Å X-ray crystal structure of the archaeal DHS-IF5A complex, revealing hetero-octameric architecture and providing a detailed view of the complex active site including the hypusination loop. This structure, along with biophysical data and molecular dynamics simulations, provides new insights into the catalytic mechanism of the hypusination reaction.

Single-site iron-anchored amyloid hydrogels as catalytic platforms for alcohol detoxification

Constructing effective antidotes to reduce global health impacts induced by alcohol prevalence is a challenging topic. Despite the positive effects observed with intravenous applications of natural enzyme complexes, their insufficient activities and complicated usage often result in the accumulation of toxic acetaldehyde, which raises important clinical concerns, highlighting the pressing need for stable oral strategies. Here we present an effective solution for alcohol detoxification by employing a biomimetic-nanozyme amyloid hydrogel as an orally administered catalytic platform. We exploit amyloid fibrils derived from β-lactoglobulin, a readily accessible milk protein that is rich in coordinable nitrogen atoms, as a nanocarrier to stabilize atomically dispersed iron (ferrous-dominated). By emulating the coordination structure of the horseradish peroxidase enzyme, the single-site iron nanozyme demonstrates the capability to selectively catalyse alcohol oxidation into acetic acid, as opposed to the more toxic acetaldehyde. Administering the gelatinous nanozyme to mice suffering from alcohol intoxication significantly reduced their blood-alcohol levels (decreased by 55.8% 300 min post-alcohol intake) without causing additional acetaldehyde build-up. Our hydrogel further demonstrates a protective effect on the liver, while simultaneously mitigating intestinal damage and dysbiosis associated with chronic alcohol consumption, introducing a promising strategy in effective alcohol detoxification.

Balancing Group I Monatomic Ion-Polar Compound Interactions for Condensed Phase Simulation in the Polarizable Drude Force Field

Molecular dynamics (MD) simulations are a commonly used method for investigating molecular behavior at the atomic level. Achieving reliable MD simulation results necessitates the use of an accurate force field. In the present work, we present a protocol to enhance the quality of group 1 monatomic ions (specifically Li+, Na+, K+, Rb+, and Cs+) with respect to their interactions with common polar model compounds in biomolecules in condensed phases in the context of the Drude polarizable force field. Instead of adjusting preexisting individual parameters for ions, model compounds, and water, we employ atom-pair specific Lennard-Jones (LJ) (known as NBFIX in CHARMM) and through-space Thole dipole screening (NBTHOLE) terms to fine-tune the balance of ion-model compound, ion-water, and model compound-water interactions. This involved establishing a protocol for the optimization of NBFIX and NBTHOLE parameters targeting the difference between molecular mechanical (MM) and quantum mechanical (QM) potential energy scans (PES). It is shown that targeting PES involving complexes that include multiple model compounds and/or ions as trimers and tetramers yields parameters that produce condensed phase properties in agreement with experimental data. Validation of this protocol involved the reproduction of experimental thermodynamic benchmarks, including solvation free energies of ions in methanol and N-methylacetamide, osmotic pressures, ionic conductivities, and diffusion coefficients within the condensed phase. These results show the importance of including more complex ion-model compound complexes beyond dimers in the QM target data to account for many-body effects during parameter fitting. The presented parameters represent a significant refinement of the Drude polarizable force field, which will lead to improved accuracy for modeling ion-biomolecular interactions.

Simulations of Functional Motions of Super Large Biomolecules with a Mixed-Resolution Model

Many large protein machines function through an interplay between large-scale movements and intricate conformational changes. Understanding functional motions of these proteins through simulations becomes challenging for both all-atom and coarse-grained (CG) modeling techniques because neither approach alone can readily capture the full details of these motions. In this study, we develop a multiscale model by employing the popular MARTINI CG model to represent a heterogeneous environment and structurally stable proteins and using the united-atom (UA) model PACE to describe proteins undergoing subtle conformational changes. PACE was previously developed to be compatible with the MARTINI solvent and membrane. Here, we couple the protein descriptions of the two models by directly mixing UA and CG interaction parameters to greatly simplify parameter determination. Through extensive validations with diverse protein systems in solution or membrane, we demonstrate that only additional parameter rescaling is needed to enable the resulting model to recover the stability of native structures of proteins under mixed representation. Moreover, we identify the optimal scaling factors that can be applied to various protein systems, rendering the model potentially transferable. To further demonstrate its applicability for realistic systems, we apply the model to a mechanosensitive ion channel Piezo1 that has peripheral arms for sensing membrane tension and a central pore for ion conductance. The model can reproduce the coupling between Piezo1's large-scale arm movement and subtle pore opening in response to membrane stress while consuming much less computational costs than all-atom models. Therefore, our model shows promise for studying functional motions of large protein machines.

Machine Learning-Based Interfacial Tension Equations for (H2 + CO2)-Water/Brine Systems over a Wide Range of Temperature and Pressure

Large-scale underground hydrogen storage (UHS) plays a vital role in energy transition. H2-brine interfacial tension (IFT) is a crucial parameter in structural trapping in underground geological locations and gas-water two-phase flow in subsurface porous media. On the other hand, cushion gas, such as CO2, is often co-injected with H2 to retain reservoir pressure. Therefore, it is imperative to accurately predict the (H2 + CO2)-water/brine IFT under UHS conditions. While there have been a number of experimental measurements on H2-water/brine and (H2 + CO2)-water/brine IFT, an accurate and efficient (H2 + CO2)-water/brine IFT model under UHS conditions is still lacking. In this work, we use molecular dynamics (MD) simulations to generate an extensive (H2 + CO2)-water/brine IFT databank (840 data points) over a wide range of temperature (from 298 to 373 K), pressure (from 50 to 400 bar), gas composition, and brine salinity (up to 3.15 mol/kg) for typical UHS conditions, which is used to develop an accurate and efficient machine learning (ML)-based IFT equation. Our ML-based IFT equation is validated by comparing to available experimental data and other IFT equations for various systems (H2-brine/water, CO2-brine/water, and (H2 + CO2)-brine/water), rendering generally good performance (with R2 = 0.902 against 601 experimental data points). The developed ML-based IFT equation can be readily applied and implemented in reservoir simulations and other UHS applications.

Mutational analysis in Corynebacterium stationis MFS transporters for improving nucleotide bioproduction

Product secretion from an engineered cell can be advantageous for microbial cell factories. Extensive work on nucleotide manufacturing, one of the most successful microbial fermentation processes, has enabled Corynebacterium stationis to transport nucleotides outside the cell by random mutagenesis; however, the underlying mechanism has not been elucidated, hindering its applications in transporter engineering. Herein, we report the nucleotide-exporting major facilitator superfamily (MFS) transporter from the C. stationis genome and its hyperactive mutation at the G64 residue. Structural estimation and molecular dynamics simulations suggested that the activity of this transporter improved via two mechanisms: (1) enhancing interactions between transmembrane helices through the conserved "RxxQG" motif along with substrate binding and (2) trapping substrate-interacting residue for easier release from the cavity. Our results provide novel insights into how MFS transporters change their conformation from inward- to outward-facing states upon substrate binding to facilitate efflux and can contribute to the development of rational design approaches for efflux improvements in microbial cell factories. KEYPOINTS: • An MFS transporter from C. stationis genome and its mutation at residue G64 were assessed • It enhanced the transporter activity by strengthening transmembrane helix interactions and trapped substrate-interacting residues • Our results contribute to rational design approach development for efflux improvement.

Computing the frequency-dependent NMR relaxation of 1H nuclei in liquid water

We present a computational framework for reliably determining the frequency-dependent intermolecular and intramolecular nuclear magnetic resonance (NMR) dipole-dipole relaxation rates of spin 1/2 nuclei from Molecular Dynamics (MD) simulations. This approach avoids the alterations caused by the well-known finite-size effects of translational diffusion. Moreover, a procedure is derived to control and correct for effects caused by fixed distance-sampling cutoffs and periodic boundary conditions. By construction, this approach is capable of accurately predicting the correct low-frequency scaling behavior of the intermolecular NMR dipole-dipole relaxation rate and thus allows for the reliable calculation of the frequency-dependent relaxation rate over many orders of magnitude. Our approach is based on the utilization of the theory of Hwang and Freed for the intermolecular dipole-dipole correlation function and its corresponding spectral density [L.-P. Hwang and J. H. Freed, J. Chem. Phys. 63, 4017-4025 (1975)] and its combination with data from MD simulations. The deviations from the Hwang and Freed theory caused by periodic boundary conditions and sampling distance cutoffs are quantified by means of random walker Monte Carlo simulations. An expression based on the Hwang and Freed theory is also suggested for correcting those effects. As a proof of principle, our approach is demonstrated by computing the frequency-dependent intermolecular and intramolecular dipolar NMR relaxation rates of 1H nuclei in liquid water at 273 and 298 K based on the simulations of the TIP4P/2005 model. Our calculations are suggesting that the intermolecular contribution to the 1H NMR relaxation rate of the TIP4P/2005 model in the extreme narrowing limit has previously been substantially underestimated.

Computing Accurate True Self-Diffusion Coefficients and Shear Viscosities Using the OrthoBoXY Approach

In a recent paper [Busch, J.; Paschek, D. J. Phys. Chem. B 2023, 127, 7983-7987], we have shown that for molecular dynamics (MD) simulations using orthorhombic periodic boundary conditions with "magic" box length ratios of Lz/Lx = Lz/Ly = 2.7933596497, the self-diffusion coefficients Dx and Dy in x- and y-directions are independent of the system size. They both represent the true self-diffusion coefficient D0 = (Dx + Dy)/2, while the shear viscosity can be calculated from diffusion coefficients in x-, y-, and z-directions, using η = kBT·8.1711245653/[3πLz(Dx + Dy - 2Dz)]. In this contribution, we test this "OrthoBoXY" approach by its application to a variety of different systems: liquid water, dimethyl ether, methanol, triglyme, water/methanol mixtures, water/triglyme mixtures, and imidazolium-based ionic liquids. The chosen systems range from small-sized molecular liquids to complex mixtures and ionic liquids, while spanning a viscosity range of almost 3 orders of magnitude. We assess the efficiency of the method for computing true self-diffusion and viscosity data and provide simple formulas for estimating the required MD simulation lengths and sizes for delivering reliable data with targeted uncertainty levels. Our analysis of the system size dependence of statistical uncertainties for both the viscosity and the self-diffusion coefficient leads us to the conclusion that it is preferable to extend the simulation length instead of increasing the system size. MD simulations consisting of 768 molecules or ion pairs seem to be perfectly adequate.

Biomolecular dynamics in the 21st century

The relevance of motions in biological macromolecules has been clear since the early structural analyses of proteins by X-ray crystallography. Computer simulations have been applied to provide a deeper understanding of the dynamics of biological macromolecules since 1976, and are now a standard tool in many labs working on the structure and function of biomolecules. In this mini-review we highlight some areas of current interest and active development for simulations, in particular all-atom molecular dynamics simulations.

An OrthoBoXY-method for various alternative box geometries

We have shown in a recent contribution [Busch and Paschek, J. Phys. Chem. B, 2023 127, 7983-7987] that for molecular dynamics (MD) simulations of isotropic fluids based on orthorhombic periodic boundary conditions with "magical" box length ratios of Lz/Lx = Lz/Ly = 2.7933596497, the computed self-diffusion coefficients Dx and Dy in x- and y-direction become system size independent. They thus represent the true self-diffusion coefficient D0 = (Dx + Dy)/2, while the shear viscosity can be determined from diffusion coefficients in x-, y-, and z-direction, using the expression η = kBT·8.1711245653/[3πLz(Dx + Dy - 2Dz)]. Here we present a more generalized version of this "OrthoBoXY"-approach, which can be applied to any orthorhombic MD box of any shape. In particular, we would like to test, how the efficiency is affected by using a shape more akin to the cubic form, albeit with different box length ratios Lx/Lz ≠ Ly/Lz and Lx < Ly < Lz. We use NVT and NpT simulations of systems of 1536 TIP4P/2005 water molecules as a benchmark and explore different box geometries to determine the influence of the box shape on the computed statistical uncertainties for D0 and η. Moreover, another "magical" set of box length ratios is discovered with Ly/Lz = 0.57804765578 and Lx/Lz = 0.33413909235, where the self-diffusion coefficient in x-direction becomes system size independent, such that D0 = Dx.

Water Exchange from the Buried Binding Sites of Cytochrome P450 Enzymes 1A2, 2D6, and 3A4 Correlates with Conformational Fluctuations

Human cytochrome P450 enzymes (CYPs) are critical for the metabolism of small-molecule pharmaceuticals (drugs). As such, the prediction of drug metabolism by and drug inhibition of CYP activity is an important component of the drug discovery and design process. Relative to the availability of a wide range of experimental atomic-resolution CYP structures, the development of structure-based CYP activity models has been limited. To better characterize the role of CYP conformational fluctuations in CYP activity, we perform multiple microsecond-scale all-atom explicit-solvent molecular dynamics (MD) simulations on three CYP isoforms, 1A2, 2D6, and 3A4, which together account for the majority of CYP-mediated drug metabolism. The MD simulations employ a variety of positional restraints, ranging from keeping all CYP atoms close to their experimentally determined coordinates to allowing full flexibility. We find that, with full flexibility, large fluctuations in the CYP binding sites correlate with efficient water exchange from these buried binding sites. This is especially true for 1A2, which, when restrained to its crystallographic conformation, is unable to exchange water between the binding site and bulk solvent. These findings imply that, in addition to crystal structures, a representative ensemble of conformational states ought to be included when developing structure-based CYP activity models.

Comprehensive mapping of mutations in TDP-43 and α-Synuclein that affect stability and binding

Abnormal aggregation and amyloid inclusions of TAR DNA-binding protein 43 (TDP-43) and α-Synuclein (α-Syn) are frequently co-observed in amyotrophic lateral sclerosis, Parkinson's disease, and Alzheimer's disease. Several reports showed TDP-43 C-terminal domain (CTD) and α-Syn interact with each other and the aggregates of these two proteins colocalized together in different cellular and animal models. Molecular dynamics simulation was conducted to elucidate the stability of the TDP-43 and Syn complex structure. The interfacial mutations in protein complexes changes the stability and binding affinity of the protein that may cause diseases. Here, we have utilized the computational saturation mutagenesis approach including structure-based stability and binding energy calculations to compute the systemic effects of missense mutations of TDP-43 CTD and α-Syn on protein stability and binding affinity. Most of the interfacial mutations of CTD and α-Syn were found to destabilize the protein and reduced the protein binding affinity. The results thus shed light on the functional consequences of missense mutations observed in TDP-43 associated proteinopathies and may provide the mechanisms of co-morbidities involving these two proteins.Communicated by Ramaswamy H. Sarma.

Phosphorylation motif dictates GPCR C-terminal domain conformation and arrestin interaction

Arrestin-dependent G protein-coupled receptor (GPCR) signaling pathway is regulated by the phosphorylation state of GPCR's C-terminal domain, but the molecular bases of arrestin:receptor interaction are to be further illuminated. Here we investigated the impact of phosphorylation on the conformational features of the C-terminal region from three rhodopsin-like GPCRs, the vasopressin V2 receptor (V2R), the growth hormone secretagogue or ghrelin receptor type 1a (GHSR), and the β2-adernergic receptor (β2AR). Using phosphomimetic variants, we identified pre-formed secondary structure elements, or short linear motifs (SLiMs), that undergo specific conformational transitions upon phosphorylation. Of importance, such conformational transitions appear to favor arrestin-2 binding. Hence, our results suggest a model in which the phosphorylation-dependent structuration of the GPCR C-terminal regions would modulate arrestin binding and therefore signaling outcomes in arrestin-dependent pathways.

OrthoBoXY: A Simple Way to Compute True Self-Diffusion Coefficients from MD Simulations with Periodic Boundary Conditions without Prior Knowledge of the Viscosity

Recently, an analytical expression for the system size dependence and direction-dependence of self-diffusion coefficients for neat liquids due to hydrodynamic interactions has been derived for molecular dynamics (MD) simulations using orthorhombic unit cells. Based on this description, we show that for systems with a "magic" box length ratio of Lz/Lx = Lz/Ly = 2.7933596497 the computed self-diffusion coefficients Dx and Dy in the x- and y-direction become system-size independent and represent the true self-diffusion coefficient D0 = (Dx + Dy)/2. Moreover, by using this particular box geometry, the viscosity can be determined with a reasonable degree of accuracy from the difference of components of the diffusion coefficients in x-, y-, and z-directions using the simple expression η = kBT × 8.1711245653/[3πLz(Dx + Dy - 2Dz)], where kB denotes Boltzmann's constant and T represents the temperature. MD simulations of TIP4P/2005 water for various system sizes using both orthorhombic and cubic box geometries are used to test the approach.

Calculations of Absolute Solvation Free Energies with Transformato─Application to the FreeSolv Database Using the CGenFF Force Field

We recently introduced transformato, an open-source Python package for the automated setup of large-scale calculations of relative solvation and binding free energy differences. Here, we extend the capabilities of transformato to the calculation of absolute solvation free energy differences. After careful validation against the literature results and reference calculations with the PERT module of CHARMM, we used transformato to compute absolute solvation free energies for most molecules in the FreeSolv database (621 out of 642). The force field parameters were obtained with the program cgenff (v2.5.1), which derives missing parameters from the CHARMM general force field (CGenFF v4.6). A long-range correction for the Lennard-Jones interactions was added to all computed solvation free energies. The mean absolute error compared to the experimental data is 1.12 kcal/mol. Our results allow a detailed comparison between the AMBER and CHARMM general force fields and provide a more in-depth understanding of the capabilities and limitations of the CGenFF small molecule parameters.

Identifying and overcoming the sampling challenges in relative binding free energy calculations of a model protein:protein complex

Relative alchemical binding free energy calculations are routinely used in drug discovery projects to optimize the affinity of small molecules for their drug targets. Alchemical methods can also be used to estimate the impact of amino acid mutations on protein:protein binding affinities, but these calculations can involve sampling challenges due to the complex networks of protein and water interactions frequently present in protein:protein interfaces. We investigate these challenges by extending a GPU-accelerated opensource relative free energy calculation package (Perses) to predict the impact of amino acid mutations on protein:protein binding. Using the well-characterized model system barnase:barstar, we describe analyses for identifying and characterizing sampling problems in protein:protein relative free energy calculations. We find that mutations with sampling problems often involve charge-changes, and inadequate sampling can be attributed to slow degrees of freedom that are mutation-specific. We also explore the accuracy and efficiency of current state-of-the-art approaches-alchemical replica exchange and alchemical replica exchange with solute tempering-for overcoming relevant sampling problems. By employing sufficiently long simulations, we achieve accurate predictions (RMSE 1.61, 95% CI: [1.12, 2.11] kcal/mol), with 86% of estimates within 1 kcal/mol of the experimentally-determined relative binding free energies and 100% of predictions correctly classifying the sign of the changes in binding free energies. Ultimately, we provide a model workflow for applying protein mutation free energy calculations to protein:protein complexes, and importantly, catalog the sampling challenges associated with these types of alchemical transformations. Our free open-source package (Perses) is based on OpenMM and available at https://github.com/choderalab/perses .

Structure-based virtual screening to identify potential lipase inhibitors to reduce lipid storage in Wolman disorder

Wolman disorder (WD) was first described in Iranian-Jewish (IJ) children, and it is caused by a deficiency of the lysosomal acid lipase (LAL). Newborns with WD are healthy and active at birth but soon develop severe malnutrition symptoms and often die before 1 year. In particular, spleens, livers, bone marrows, intestines, adrenal glands, and lymph nodes accumulate harmful amounts of lipids. G87V mutation in LIPA is responsible for Wolman disorder. Some reports suggest that δ-tocopherol can reduce lipid accumulation in cholesterol storage disorders. Hence, we used δ-tocopherol for the virtual screening process in this study. Initially, the lead compounds were docked with native and G87V mutant LIPA. Subsequently, the ADME and toxicity parameters for screened compounds were determined to ensure the safety profiles. Finally, the molecular dynamics simulations result indicated that dl-alpha-Tocopherol-13C3, a molecule obtained from the PubChem database, is identified as a potential and stable lead molecule that could be effective against the G87V mutant form of LIPA.

Discovering potential inhibitors against SARS-CoV-2 by targeting Nsp13 Helicase

The rise in the incidence of COVID-19 as a result of SARS-CoV-2 infection has threatened public health globally. Till now, there have been no proper prophylactics available to fight COVID-19, necessitating the advancement and evolution of effective curative against SARS-CoV-2. This study aimed at the nonstructural protein 13 (nsp13) helicase as a promising target for drug development against COVID-19. A unique collection of nucleoside analogs was screened against the SARS-CoV-2 helicase protein, for which a molecular docking experiment was executed to depict the selected ligand's binding affinity with the SARS-CoV-2 helicase proteins. Simultaneously, molecular dynamic simulations were performed to examine the protein's binding site's conformational stability, flexibility, and interaction with the ligands. Key nucleoside ligands were selected for pharmacokinetic analysis based on their docking scores. Selected ligands (cordycepin and pritelivir) showed excellent pharmacokinetics and were well stabilized at the proteins' binding site throughout the MD simulation. We have also performed binding free energy analysis or the binding characteristics of ligands with Nsp13 by using MM-PBSA and MM-GBSA. Free energy calculation by MM-PBSA and MM-GBSA analysis suggests that pritelivir may work as viable therapeutics for efficient drug advancement against SARS-CoV-2 Nsp13 helicase, potentially arresting the SARS-CoV-2 replication.Communicated by Ramaswamy H. Sarma.

Comparison of the United- and All-Atom Representations of (Halo)alkanes Based on Two Condensed-Phase Force Fields Optimized against the Same Experimental Data Set

The level of accuracy that can be achieved by a force field is influenced by choices made in the interaction-function representation and in the relevant simulation parameters. These choices, referred to here as functional-form variants (FFVs), include for example the model resolution, the charge-derivation procedure, the van der Waals combination rules, the cutoff distance, and the treatment of the long-range interactions. Ideally, assessing the effect of a given FFV on the intrinsic accuracy of the force-field representation requires that only the specific FFV is changed and that this change is performed at an optimal level of parametrization, a requirement that may prove extremely challenging to achieve in practice. Here, we present a first attempt at such a comparison for one specific FFV, namely the choice of a united-atom (UA) versus an all-atom (AA) resolution in a force field for saturated acyclic (halo)alkanes. Two force-field versions (UA vs AA) are optimized in an automated way using the CombiFF approach against 961 experimental values for the pure-liquid densities ρ_liq and vaporization enthalpies ΔH_vap of 591 compounds. For the AA force field, the torsional and third-neighbor Lennard-Jones parameters are also refined based on quantum-mechanical rotational-energy profiles. The comparison between the UA and AA resolutions is also extended to properties that have not been included as parameterization targets, namely the surface-tension coefficient γ, the isothermal compressibility κ_T, the isobaric thermal-expansion coefficient α_P, the isobaric heat capacity c_P, the static relative dielectric permittivity ϵ, the self-diffusion coefficient D, the shear viscosity η, the hydration free energy ΔG_wat, and the free energy of solvation ΔG_che in cyclohexane. For the target properties ρ_liq and ΔH_vap, the UA and AA resolutions reach very similar levels of accuracy after optimization. For the nine other properties, the AA representation leads to more accurate results in terms of η; comparably accurate results in terms of γ, κ_T, α_P, ϵ, D, and ΔG_che; and less accurate results in terms of c_P and ΔG_wat. This work also represents a first step toward the calibration of a GROMOS-compatible force field at the AA resolution.

Orientation and Conformation of Hydrophobin at the Oil-Water Interface: Insights from Molecular Dynamics Simulations

Hydrophobins, a new class of potential protein emulsifiers, have been extensively employed in the food, pharmaceutical, and chemical industries. However, the knowledge of the underlying molecular mechanism of protein adsorption at the oil-water interface remains elusive. In this study, all-atom molecular dynamics simulations were performed to probe the adsorption orientation and conformation change of class II hydrophobin HFBI at the cyclohexane-water interface. It was proposed that a hydrophobic dipole of the protein could be used to quantitatively predict the orientation of the adsorbed HFBI. Simulation results revealed that HFBI adsorbed at the interface with the patch-up orientation toward the oil phase, regardless of its initial orientations. HFBI's secondary structure was maintained to be intact in the course of simulations despite relatively significant variations in the tertiary structure observed, which could well preserve the bioactivity of HFBI. From the energy analysis, the driving force for interface adsorption was primarily determined by van der Waals interactions between HFBI and cyclohexane. Further analysis indicated that the adsorption orientation and conformation of HFBI at the oil-water interface were typically regulated by the hydrophobic patch and some key residues. This study provides some insights into the orientation, conformation, and adsorption mechanism of proteins at the oil-water interface and theoretical guidelines for the design and development of novel biological emulsifiers involved in the food, pharmaceutical, and chemical industries.

Computational Insights of Unfolding of N-Terminal Domain of TDP-43 Reveal the Conformational Heterogeneity in the Unfolding Pathway

TDP-43 proteinopathies is a disease hallmark that characterizes amyotrophic lateral sclerosis (ALS) and frontotemporal lobar degeneration (FTLD). The N-terminal domain of TDP-43 (NTD) is important to both TDP-43 physiology and TDP-43 proteinopathy. However, its folding and dimerization process is still poorly characterized. In the present study, we have investigated the folding/unfolding of NTD employing all-atom molecular dynamics (MD) simulations in 8 M dimethylsulfoxide (DMSO) at high temperatures. The MD results showed that the unfolding of the NTD at high temperature evolves through the formation of a number of conformational states differing in their stability and free energy. The presence of structurally heterogeneous population of intermediate ensembles was further characterized by the different extents of solvent exposure of Trp80 during unfolding. We suggest that these non-natives unfolded intermediate ensembles may facilitate NTD oligomerization and subsequently TDP-43 oligomerization, which might lead to the formation of irreversible pathological aggregates, characteristics of disease pathogenesis.

Molecular simulations of proteins: From simplified physical interactions to complex biological phenomena

Molecular dynamics simulation is the most popular computational technique for investigating the structural and dynamical behaviour of proteins, in search of the molecular basis of their function. Far from being a completely settled field of research, simulations are still evolving to best capture the essential features of the atomic interactions that govern a protein's inner motions. Modern force fields are becoming increasingly accurate in providing a physical description adequate to this purpose, and allow us to model complex biological systems under fairly realistic conditions. Furthermore, the use of accelerated sampling techniques is improving our access to the observation of progressively larger molecular structures, longer time scales, and more hidden functional events. In this review, the basic principles of molecular dynamics simulations and a number of key applications in the area of protein science are summarized, and some of the most important results are discussed. Examples include the study of the structure, dynamics and binding properties of 'difficult' targets, such as intrinsically disordered proteins and membrane receptors, and the investigation of challenging phenomena like hydration-driven processes and protein aggregation. The findings described provide an overall picture of the current state of this research field, and indicate new perspectives on the road ahead to the upcoming future of molecular simulations.

Accurate Simulations of Lipid Monolayers Require a Water Model with Correct Surface Tension

Lipid monolayers provide our lungs and eyes their functionality and serve as proxy systems in biomembrane research. Therefore, lipid monolayers have been studied intensively including using molecular dynamics simulations, which are able to probe their lateral structure and interactions with, e.g., pharmaceuticals or nanoparticles. However, such simulations have struggled in describing the forces at the air-water interface. Particularly, the surface tension of water and long-range van der Waals interactions have been considered critical, but their importance in monolayer simulations has been evaluated only separately. Here, we combine the recent C36/LJ-PME lipid force field that includes long-range van der Waals forces with water models that reproduce experimental surface tensions to elucidate the importance of these contributions in monolayer simulations. Our results suggest that a water model with correct surface tension is necessary to reproduce experimental surface pressure-area isotherms and monolayer phase behavior. The latter includes the liquid expanded and liquid condensed phases, their coexistence, and the opening of pores at the correct area per lipid upon expansion. Despite these improvements of the C36/LJ-PME with certain water models, the standard cutoff-based CHARMM36 lipid model with the 4-point OPC water model still provides the best agreement with experiments. Our results emphasize the importance of using high-quality water models in applications and parameter development in molecular dynamics simulations of biomolecules.

All-Atom Modeling of Complex Cellular Membranes

Cell membranes are composed of a variety of lipids and proteins where they interact with each other to fulfill their roles. The first step in modeling these interactions in molecular simulations is to have reliable mimetics of the membrane's lipid environment. This Feature Article presents our recent efforts to model complex cellular membranes using all-atom force fields. A short review of the CHARMM36 (C36) lipid force field and its recent update to incorporate the long-range dispersion is presented. Key examples of model membranes mimicking various species and organelles are given. These include single-celled organisms such as bacteria (E. coli., chlamydia, and P. aeruginosa) and yeast (plasma membrane, endoplasmic reticulum, and trans-Golgi network) and more advanced ones such as plants (soybean and Arabidopsis thaliana) and mammals (ocular lens, stratum corneum, and peripheral nerve myelin). Leaflet asymmetry in composition has also been applied to some of these models. With the increased lipid diversity in the C36 lipid FF, these complex models can better reflect the structural, mechanical, and dynamic properties of realistic membranes and open an opportunity to study biological processes involving other molecules.

Molecular understanding of charge effect on desalination performance in lamellar MoS2 membranes

The effect of atomic charge information on the desalination performance of lamellar MoS2 membranes was investigated at the molecular level.

Hydrophobicity of Self-Assembled Monolayers of Alkanes: Fluorination, Density, Roughness, and Lennard-Jones Cutoffs

The interplay of fluorination and structure of alkane self-assembled monolayers and how these affect hydrophobicity are explored via molecular dynamics simulations, contact angle goniometry, and surface-enhanced infrared absorption spectroscopy. Wetting coefficients are found to grow linearly in the monolayer density for both alkane and perfluoroalkane monolayers. The larger contact angles of monolayers of perfluorinated alkanes are shown to be primarily caused by their larger molecular volume, which leads to a larger nearest-neighbor grafting distance and smaller tilt angle. Increasing the Lennard-Jones force cutoff in simulations is found to increase hydrophilicity. Specifically, wetting coefficients scale like the inverse square of the cutoff, and when extrapolated to the infinite cutoff limit, they yield contact angles that compare favorably to experimental values. Nanoscale roughness is also found to reliably increase monolayer hydrophobicity, mostly via the reduction of the entropic part of the work of adhesion. Analysis of depletion lengths shows that droplets on nanorough surfaces partially penetrate the surface, intermediate between Wenzel and Cassie-Baxter states.

Validating the CHARMM36m protein force field with LJ-PME reveals altered hydrogen bonding dynamics under elevated pressures

The pressure-temperature phase diagram is important to our understanding of the physics of biomolecules. Compared to studies on temperature effects, studies of the pressure dependence of protein dynamic are rather limited. Molecular dynamics (MD) simulations with fine-tuned force fields (FFs) offer a powerful tool to explore the influence of thermodynamic conditions on proteins. Here we evaluate the transferability of the CHARMM36m (C36m) protein force field at varied pressures compared with NMR data using ubiquitin as a model protein. The pressure dependences of J couplings for hydrogen bonds and order parameters for internal motion are in good agreement with experiment. We demonstrate that the C36m FF combined with the Lennard-Jones particle-mesh Ewald (LJ-PME) method is suitable for simulations in a wide range of temperature and pressure. As the ubiquitin remains stable up to 2500 bar, we identify the mobility and stability of different hydrogen bonds in response to pressure. Based on those results, C36m is expected to be applied to more proteins in the future to further investigate protein dynamics under elevated pressures.

Evaluation of nine condensed-phase force fields of the GROMOS, CHARMM, OPLS, AMBER, and OpenFF families against experimental cross-solvation free energies

Experimental solvation free energies are nowadays commonly included as target properties in the validation of condensed-phase force fields, sometimes even in their calibration. In a previous article [Kashefolgheta et al., J. Chem. Theory. Comput., 2020, 16, 7556-7580], we showed how the involved comparison between experimental and simulation results could be made more systematic by considering a full matrix of cross-solvation free energies . For a set of N molecules that are all in the liquid state under ambient conditions, such a matrix encompasses N×N entries for considering each of the N molecules either as solute (A) or as solvent (B). In the quoted study, a cross-solvation matrix of 25 × 25 experimental value was introduced, considering 25 small molecules representative for alkanes, chloroalkanes, ethers, ketones, esters, alcohols, amines, and amides. This experimental data was used to compare the relative accuracies of four popular condensed-phase force fields, namely GROMOS-2016H66, OPLS-AA, AMBER-GAFF, and CHARMM-CGenFF. In the present work, the comparison is extended to five additional force fields, namely GROMOS-54A7, GROMOS-ATB, OPLS-LBCC, AMBER-GAFF2, and OpenFF. Considering these nine force fields, the correlation coefficients between experimental values and simulation results range from 0.76 to 0.88, the root-mean-square errors (RMSEs) from 2.9 to 4.8 kJ mol-1, and average errors (AVEEs) from -1.5 to +1.0 kJ mol-1. In terms of RMSEs, GROMOS-2016H66 and OPLS-AA present the best accuracy (2.9 kJ mol-1), followed by OPLS-LBCC, AMBER-GAFF2, AMBER-GAFF, and OpenFF (3.3 to 3.6 kJ mol-1), and then by GROMOS-54A7, CHARM-CGenFF, and GROMOS-ATB (4.0 to 4.8 kJ mol-1). These differences are statistically significant but not very pronounced, and are distributed rather heterogeneously over the set of compounds within the different force fields.

Considerations of Recent All-Atom Lipid Force Field Development

Molecular simulations of biological molecules require an accurate description of molecular interactions through a force field (FF). The focus of this Perspective is on all-atom lipid FFs. Recent additions to the CHARMM36 lipid FF continue to expand a researcher's ability to probe membrane structure and function with a wide variety of biologically important lipids. Currently, there is an effort to reduce the assumptions in all-atom lipid FFs. The inclusion of long-range dispersion interaction through particle-mesh Ewald is allowing for more accurate descriptions of lipid bilayer and monolayer properties without additional computational cost. Soon, simulations with lipid FFs will no longer depend on short-range cutoffs and will accurately represent long-range dispersion. This requires efficient FF parametrization with an automated approach due to FF complexity. In addition, polarizable FFs for lipids will be important for the next generation of simulations that accurately represent how molecule interactions respond to a varied environment.

An exact a posteriori correction for hydrogen bond population correlation functions and other reversible geminate recombinations obtained from simulations with periodic boundary conditions. Liquid water as a test case

The kinetics of breaking and re-formation of hydrogen bonds (HBs) in liquid water is a prototype of reversible geminate recombination. HB population correlation functions (HBPCFs) are a means to study the HB kinetics. The long-time limiting behavior of HBPCFs is controlled by translatoric diffusion and shows a t^-3/2 time-dependence, which can be described by analytical expressions based on the HB acceptor density and the donor-acceptor inter-diffusion coefficient. If the trajectories are not properly "unwrapped," the presence of periodic boundary conditions (PBCs) can perturb this long-time limiting behavior. Keeping the trajectories "wrapped," however, allows for a more efficient calculation of HBPCFs. We discuss the consequences of PBCs in combination with "wrapped" trajectories following from the approximations according to Luzar-Chandler and according to Starr, each deviating in a different fashion from the true long-time limiting behavior, but enveloping the unperturbed function. A simple expression is given for estimating the maximum time up to which the computed HBPCFs reliably describe the long-time limiting behavior. In addition, an exact a posteriori correction for systems with PBCs for "wrapped" trajectories is derived, which can be easily computed and which is able to fully recover the true t^-3/2 long-time behavior. For comparison, HBPCFs are computed from MD simulations of TIP4P/2005 model water for varying system sizes and temperatures of 273 and 298 K using this newly introduced correction. Implications for the computations of HB lifetimes and the effect of the system-size are discussed.

Engineered MATE multidrug transporters reveal two functionally distinct ion-coupling pathways in NorM from Vibrio cholerae

Multidrug and toxic compound extrusion (MATE) transport proteins confer multidrug resistance on pathogenic microorganisms and affect pharmacokinetics in mammals. Our understanding of how MATE transporters work, has mostly relied on protein structures and MD simulations. However, the energetics of drug transport has not been studied in detail. Many MATE transporters utilise the electrochemical H+ or Na+ gradient to drive substrate efflux, but NorM-VC from Vibrio cholerae can utilise both forms of metabolic energy. To dissect the localisation and organisation of H+ and Na+ translocation pathways in NorM-VC we engineered chimaeric proteins in which the N-lobe of H+-coupled NorM-PS from Pseudomonas stutzeri is fused to the C-lobe of NorM-VC, and vice versa. Our findings in drug binding and transport experiments with chimaeric, mutant and wildtype transporters highlight the versatile nature of energy coupling in NorM-VC, which enables adaptation to fluctuating salinity levels in the natural habitat of V. cholerae.

Kirkwood-Buff-Derived Force Field for Peptides and Proteins: Philosophy and Development of KBFF20

A new classical nonpolarizable force field, KBFF20, for the simulation of peptides and proteins is presented. The force field relies heavily on the use of Kirkwood-Buff theory to provide a comparison of simulated and experimental Kirkwood-Buff integrals for solutes containing the functional groups common in proteins, thus ensuring intermolecular interactions that provide a good balance between the peptide-peptide, peptide-solvent, and solvent-solvent distributions observed in solution mixtures. In this way, it differs significantly from other biomolecular force fields. Further development and testing of the intermolecular potentials are presented here. Subsequently, rotational potentials for the ϕ/ψ and χ dihedral degrees of freedom are obtained by analysis of the Protein Data Bank, followed by small modifications to provide a reasonable balance between simulated and observed α and β percentages for small peptides. This, the first of two articles, describes in detail the philosophy and development behind KBFF20.

Interactions Between Nucleosomes: From Atomistic Simulation to Polymer Model

The organization of genomes in space and time dimension plays an important role in gene expression and regulation. Chromatin folding occurs in a dynamic, structured way that is subject to biophysical rules and biological processes. Nucleosomes are the basic unit of chromatin in living cells, and here we report on the effective interactions between two nucleosomes in physiological conditions using explicit-solvent all-atom simulations. Free energy landscapes derived from umbrella sampling simulations agree well with recent experimental and simulation results. Our simulations reveal the atomistic details of the interactions between nucleosomes in solution and can be used for constructing the coarse-grained model for chromatin in a bottom-up manner.

CHARMM36 Lipid Force Field with Explicit Treatment of Long-Range Dispersion: Parametrization and Validation for Phosphatidylethanolamine, Phosphatidylglycerol, and Ether Lipids

Long-range Lennard-Jones (LJ) interactions have been incorporated into the CHARMM36 (C36) lipid force field (FF) using the LJ particle-mesh Ewald (LJ-PME) method in order to remove the inconsistency of bilayer and monolayer properties arising from the exclusion of long-range dispersion [Yu, Y.; Semi-automated Optimization of the CHARMM36 Lipid Force Field to Include Explicit Treatment of Long-Range Dispersion. J. Chem. Theory Comput. 2021, 10.1021/acs.jctc.0c01326. (preceding article in this issue)]. The new FF is denoted C36/LJ-PME. While the first optimization was based on three phosphatidylcholines (PCs), this work extends the validation and parametrization to more lipids including PC, phosphatidylethanolamine (PE), phosphatidylglycerol (PG), and ether lipids. The agreement with experimental structure data is excellent for PC, PE, and ether lipids. C36/LJ-PME also compares favorably with scattering data of PG bilayers but less so with NMR deuterium order parameters of 1,2-dimyristoyl-sn-glycero-3-phospho-(1'-rac-glycerol) (DMPG) at 303.15 K, indicating a need for future optimization regarding PG-specific parameters. Frequency dependence of NMR T1 spin-lattice relaxation times is well-described by C36/LJ-PME, and the overall agreement with experiment is comparable to C36. Lipid diffusion is slower than C36 due to the added long-range dispersion causing a higher viscosity, although it is still too fast compared to experiment after correction for periodic boundary conditions. When using a 10 Å real-space cutoff, the simulation speed of C36/LJ-PME is roughly equal to C36. While more lipids will be incorporated into the FF in the future, C36/LJ-PME can be readily used for common lipids and extends the capability of the CHARMM FF by supporting monolayers and eliminating the cutoff dependence.

Semi-automated Optimization of the CHARMM36 Lipid Force Field to Include Explicit Treatment of Long-Range Dispersion

The development of the CHARMM lipid force field (FF) can be traced back to the early 1990s with its current version denoted CHARMM36 (C36). The parametrization of C36 utilized high-level quantum mechanical data and free energy calculations of model compounds before parameters were manually adjusted to yield agreement with experimental properties of lipid bilayers. While such manual fine-tuning of FF parameters is based on intuition and trial-and-error, automated methods can identify beneficial modifications of the parameters via their sensitivities and thereby guide the optimization process. This work introduces a semi-automated approach to reparametrize the CHARMM lipid FF with consistent inclusion of long-range dispersion through the Lennard-Jones particle-mesh Ewald (LJ-PME) approach. The optimization method is based on thermodynamic reweighting with regularization with respect to the C36 set. Two independent optimizations with different topology restrictions are presented. Targets of the optimizations are primarily liquid crystalline phase properties of lipid bilayers and the compression isotherm of monolayers. Pair correlation functions between water and lipid functional groups in aqueous solution are also included to address headgroup hydration. While the physics of the reweighting strategy itself is well-understood, applying it to heterogeneous, complex anisotropic systems poses additional challenges. These were overcome through careful selection of target properties and reweighting settings allowing for the successful incorporation of the explicit treatment of long-range dispersion, and we denote the newly optimized lipid force field as C36/LJ-PME. The current implementation of the optimization protocol will facilitate the future development of the CHARMM and related lipid force fields.

Investigating a Combined Stochastic Nucleation and Molecular Dynamics-Based Equilibration Approach for Constructing Large-Scale Polycrystalline Films

The morphology of small-molecule organic semiconducting materials can vary from single crystals via polycrystalline films with varying grain sizes to amorphous structures, depending on the process conditions. This structural variety affects the electronic properties and, thus, the performance of organic electronic devices. A nucleation-equilibration approach is investigated, whose focus is on the construction of morphologies with controlled variations in the average grain size. Its computational requirements are low because nucleation is purely based on geometrical considerations, thus allowing the construction of model systems of experimentally relevant sizes. Its application is demonstrated for C₆₀ and pentacene by generating single-component films that vary from amorphous to crystalline structures. It is further generalized to two-component films and applied to C₆₀: pentacene blends as well as dilute n-doped C₆₀ structures. When combined with electronic structure calculations in the future, the nucleation-equilibration approach can offer insights into the impact of polycrystallinity on electronic and charge-transport properties in the absence of any knowledge about the growth mechanism and for a broad set of systems.

Determination of van der Waals Parameters Using a Double Exponential Potential for Nonbonded Divalent Metal Cations in TIP3P Solvent

A double exponential (DE) functional form for Lennard-Jones (LJ) interactions, proposed in our previous study, has many advantages over LJ potentials including a natural softcore characteristic for the convenience of the pathway-based free-energy calculations, fast convergence, and flexibility in use. In this work, we put the first step on the application of the DE functional form by identifying a DE potential, coined DE-TIP3P, for molecular simulations using the TIP3P water model. The developed DE-TIP3 potential was better than LJ potential in reproducing the experimental water properties. Afterward, we developed the nonbonded models of 15 divalent metal ions, which frequently appear and play vital roles in biological systems, to be consistent with the DE-TIP3P potential and TIP3P water model. Our nonbonded models were as good as the complicated nonbonded dummy cationic models by Jiang et al. and the nonbonded 12-6-4 LJ models by Li and Merz in reproducing the experimental properties of those ions. Moreover, our nonbonded models achieved a better performance than the compromise (CM) LJ models and 12-6-4 LJ models, developed by Li and Merz, in reproducing the properties of MgCl₂ in aqueous solution.

A compression strategy for particle mesh Ewald theory

Particle Mesh Ewald (PME) has become a standard method for treating long-range electrostatics in molecular simulations. Although the method has inferior asymptotic computational complexity to its linear scaling competitors, it remains enormously popular due to its high efficiency, which stems from the use of fast Fourier transforms (FFTs). This use of FFTs provides great challenges for scaling the method up to massively parallel systems, in large part because of the need to transfer large amounts of data. In this work, we demonstrate that this data transfer volume can be greatly reduced as a natural consequence of the structure of the PME equations. We also suggest an alternative algorithm that supplants the FFT with a linear algebra approach, which further decreases communication costs at the expense of increased asymptotic computational complexity. This linear algebra based approach is demonstrated to have great potential for latency hiding by interleaving communication and computation steps of the short- and long-range electrostatic terms.

Enhanced Selective Hydrogen Permeation through Graphdiyne Membrane: A Theoretical Study

Graphdiyne (GDY), with uniform pores and atomic thickness, is attracting widespread attention for application in H2 separation in recent years. However, the challenge lies in the rational design of GDYs for fast and selective H2 permeation. By MD and DFT calculations, several flexible GDYs were constructed to investigate the permeation properties of four pure gas (H2, N2, CO2, and CH4) and three equimolar binary mixtures (H2/N2, H2/CO2, and H2/CH4) in this study. When the pore size is smaller than 2.1 Å, the GDYs acted as an exceptional filter for H2 with an approximately infinite H2 selectivity. Beyond the size-sieving effect, in the separation process of binary mixtures, the blocking effect arising from the strong gas-membrane interaction was proven to greatly impede H2 permeation. After understanding the mechanism, the H2 permeance of the mixtures of H2/CO2 was further increased to 2.84 × 105 GPU by reducing the blocking effect with the addition of a tiny amount of surface charges, without sacrificing the selectivity. This theoretical study provides an additional atomic understanding of H2 permeation crossing GDYs, indicating that the GDY membrane could be a potential candidate for H2 purification.

Update of the CHARMM36 United Atom Chain Model for Hydrocarbons and Phospholipids

Accurate lipid force field (FF) parameters used in molecular dynamics (MD) simulations are crucial for understanding the properties of lipid-containing systems and biological processes related to lipids. The last update of the CHARMM36 united atom chain model (C36UA) was in 2013 [Lee, S. J. Phys. Chem. B 2014, 118, 547 556]; it utilized CHARMM36 (C36) lipid FF parameters for headgroups and OPLS-UA Lennard-Jones (LJ) parameters for tails. Simulations with the FF were able to reproduce many experimental observables of lipid bilayers accurately, but to be more applicable for a wide range of lipids, additional FF parameter optimization was needed. In this work, we present an update of the model, named C36UAr. The parameterization included the LJ parameters for hydrocarbons and related dihedrals. Bulk liquid properties (density, heat of vaporization, isothermal compressibility, and diffusion constant) of model compounds were used as targets for the LJ parameter fitting, and dihedrals were fit to either quantum mechanical (QM) or potential of mean force (PMF) calculations using C36. Thermodynamic reweighting was used to further improve the parameters. Bilayer simulations of various lipid headgroups (phosphatidylcholine, phosphatidylethanolamine, and phosphatidylglycerol) and tails (saturated, monounsaturated, and polyunsaturated) were performed to validate the model, and significant improvements were seen in bilayer properties, including surface area, membrane thicknesses, NMR deuterium order parameters, and density profiles. C36UAr was also compared to the hydrogen mass repartitioning (HMR) method. The high accuracy and competitive efficiency shown in this study make C36UAr one of the best choices for studies of membrane structure and membrane-associated proteins.

A New Lipid Force Field (FUJI)

To explore inhomogeneous and anisotropic systems such as lipid bilayers, the Lennard-Jones particle mesh Ewald (LJ-PME) method has been applied without a conventional isotropic dispersion correction. As the popular AMBER and CHARMM lipid force fields were developed using a cutoff scheme, their lipid bilayers unacceptably shrink when using the LJ-PME method. In this study, a new all-atom lipid force field (FUJI) was developed on the basis of the AMBER force-field scheme including the Lipid14 van der Waals parameters. Point charges were calculated using the restrained electrostatic potentials of many lipid conformers. Further, torsion energy profiles were calculated using high-level ab initio molecular orbitals (LCCSD(T)/Aug-cc-pVTZ//LMP2/Aug-cc-pVTZ), following which the molecular mechanical dihedral parameters were derived through a fast Fourier transform. By incorporation of these parameters into a new lipid force field without fitting experimental data, the desired lipid characteristics such as the area per lipid and lateral diffusion coefficients were obtained through GROMACS molecular dynamics simulations using the LJ-PME method and virtual hydrogen sites. The calculated area per lipid and lateral diffusion coefficients showed satisfactory agreement with experimental data. Furthermore, the electron-density profiles along the membrane normal were calculated for pure lipid bilayers, and the resulting membrane thicknesses agreed well with the experimental values. As the new lipid force field is compatible with FUJI for protein and small molecules, the new FUJI force field will offer accurate modeling for complex systems consisting of various membrane proteins and lipids.