Solvent Box Research Articles

Effect of DMSO and Triacetin Solvents on Polyvinylidene Fluoride Polymorphs: Molecular Dynamics Simulations

In this article, molecular dynamics simulations are performed to investigate the interaction between polyvinylidene fluoride (PVDF) polymorphs (β phase has more properties compared with α, such as polarity, higher mechanical strength, and piezoelectric, ferroelectric, and pyroelectric properties), and two solvents, dimethyl sulfoxide (DMSO) and glyceryl triacetate (GTA). Pure solvent boxes are built using two force fields (Compass and Dreiding) to study the stabilization of their density and solubility, which are used with the PVDF surfaces (for each solvent, computations are performed using two surfaces α and β). According to the radial distribution function results, DMSO shows a stronger interaction with two PVDF phases, allowing it to orient chains from the α phase to the β phase, while GTA exhibits a weaker interaction due to specific hydrogen bonds. These findings are confirmed by the experiment, where the fraction F(β) reaches 62% when using DMSO, while it is 36% when using GTA.

Efficient Sampling of Cavity Hydration in Proteins with Nonequilibrium Grand Canonical Monte Carlo and Polarizable Force Fields.

Prediction of the hydration levels of protein cavities and active sites is important to both mechanistic analysis and ligand design. Due to the unique microscopic environment of these buried water molecules, a polarizable model is expected to be crucial for an accurate treatment of protein internal hydration in simulations. Here we adapt a nonequilibrium candidate Monte Carlo approach for conducting grand canonical Monte Carlo simulations with the Drude polarizable force field. The GPU implementation enables the efficient sampling of internal cavity hydration levels in biomolecular systems. We also develop an enhanced sampling approach referred to as B-walking, which satisfies detailed balance and readily combines with grand canonical integration to efficiently calculate quantitative binding free energies of water to protein cavities. Applications of these developments are illustrated in a solvent box and the polar ligand binding site in trypsin. Our simulation results show that including electronic polarization leads to a modest but clear improvement in the description of water position and occupancy compared to the crystal structure. The B-walking approach enhances the range of water sampling in different chemical potential windows and thus improves the accuracy of water binding free energy calculations.

Correlation between chemical denaturation and the unfolding energetics of Acanthamoeba actophorin

The actin filament network is in part remodeled by the action of a family of filament severing proteins that are responsible for modulating the ratio between monomeric and filamentous actin. Recent work on the protein actophorin from the amoeba Acanthamoeba castellani identified a series of site-directed mutations that increase the thermal stability of the protein by 22°C. Here, we expand this observation by showing that the mutant protein is also significantly stable to both equilibrium and kinetic chemical denaturation, and employ computer simulations to account for the increase in thermal or chemical stability through an accounting of atomic-level interactions. Specifically, the potential of mean force (PMF) can be obtained from steered molecular dynamics (SMD) simulations in which a protein is unfolded. However, SMD can be inefficient for large proteins as they require large solvent boxes, and computationally expensive as they require increasingly many SMD trajectories to converge the PMF. Adaptive steered molecular dynamics (ASMD) overcomes the second of these limitations by steering the particle in stages, which allows for convergence of the PMF using fewer trajectories compared with SMD. Use of the telescoping water scheme within ASMD partially overcomes the first of these limitations by reducing the number of waters at each stage to only those needed to solvate the structure within a given stage. In the PMFs obtained from ASMD, the work of unfolding Acto-2 was found to be higher than the Acto-WT by approximately 120 kCal/mol and reflects the increased stability seen in the chemical denaturation experiments. The evolution of the average number of hydrogen bonds and number of salt bridges during the pulling process provides a mechanistic view of the structural changes of the actophorin protein as it is unfolded, and how it is affected by the mutation in concert with the energetics reported through the PMF.

Rheological behavior and molecular dynamics simulation of κ-carrageenan/casein under simulated gastrointestinal electrolyte conditions

In this study, we investigated the effects of gastrointestinal electrolytes (including simulated gastric fluid (SGF) and simulated intestinal fluid (SIF)) on the conformational transition and exposed sulfate groups of κ-carrageenan (κ-CGN) and κ-CC (κ-carrageenan/casein) during digestion, which might contribute to the colitis. The conformational transition of κ-CGN and κ-CC were assessed through zeta-potential measurements, rheological analysis, differential scanning calorimetry, Fourier-transform infrared spectroscopy, and scanning electron microscopy. The sulfate group content measurement and molecular dynamics simulation were applied to investigate the exposed quantity of free sulfate groups. The results showed that SGF and SIF promoted the formation of a helical structure of κ-CGN and enhanced the intermolecular interactions of κ-CC. The κ-CGN self-helix formation in the presence of monovalent cations (Na+ and K+), and the helix-helix association of the neighboring κ-CGN chain in the presence of divalent cations (Mg2+ and Ca2+). In κ-CC system, the stability of helix conformation was enhanced by SGF and SIF, a greater number of hydrogen bonds were detected in Na+ and K+ solvent boxes, and Mg2+ and Ca2+ not only improved the number of hydrogen bonds through salt bridges but enhanced the electrostatic interaction of κ-CC, which were in accordance with the lower zeta-potential values, greater thermal hysteresis, and intensified binding energy results. There was an increase in hydrogen bonding sites on –SO3- groups and OH- groups of κ-CGN with the effects of electrolyte cations so that decreased the content of free sulfate groups. In conclusion, SGF and SIF enhanced the binding stability and reduced the contents of free sulfate groups of κ-CC. Our study findings provided value information for the safe processing of carrageenan.

Implementation of Telescoping Boxes in Adaptive Steered Molecular Dynamics.

Long-time dynamical processes, such as those involving protein unfolding and ligand interactions, can be accelerated and realized through steered molecular dynamics (SMD). The challenge has been the extraction of information from such simulations that generalize for complex nonequilibrium processes. The use of Jarzynski's equality opened the possibility of determining the free energy along the steered coordinate, but sampling over the nonequilibrium trajectories is slow to converge. Adaptive steered molecular dynamics (ASMD) and other related techniques have been introduced to overcome this challenge through the use of stages. Here, we take advantage of these stages to address the numerical cost that arises from the required use of very large solvent boxes. We introduce telescoping box schemes within adaptive steered molecular dynamics (ASMD) in which we adjust the solvent box between stages and thereby vary (and optimize) the required number of solvent molecules. We have benchmarked the method on a relatively long α-helical peptide, Ala30, with respect to the potential of mean force and hydrogen bonds. We show that the use of telescoping boxes introduces little numerical error while significantly reducing the computational cost.

Relative Binding Free Energy Calculations for Ligands with Diverse Scaffolds with the Alchemical Transfer Method.

We present an extension of the alchemical transfer method (ATM) for the estimation of relative binding free energies of molecular complexes applicable to conventional, as well as scaffold-hopping, alchemical transformations. Named ATM-RBFE, the method is implemented in the free and open-source OpenMM molecular simulation package and aims to provide a simpler and more generally applicable route to the calculation of relative binding free energies than what is currently available. ATM-RBFE is based on sound statistical mechanics theory and a novel coordinate perturbation scheme designed to swap the positions of a pair of ligands such that one is transferred from the bulk solvent to the receptor binding site while the other moves simultaneously in the opposite direction. The calculation is conducted directly in a single solvent box with a system prepared with conventional setup tools, without splitting of electrostatic and nonelectrostatic transformations, and without pairwise soft-core potentials. ATM-RBFE is validated here against the absolute binding free energies of the SAMPL8 GDCC host-guest benchmark set and against protein-ligand benchmark sets that include complexes of the estrogen receptor ERα and those of the methyltransferase EZH2. In each case the method yields self-consistent and converged relative binding free energy estimates in agreement with absolute binding free energies and reference literature values, as well as experimental measurements.

Observing reorientation dynamics with Time-Resolved fluorescence and molecular dynamics in varying periodic boundary conditions

This work presents a combined study of time-resolved fluorescence spectroscopy and all-atom molecular dynamics simulation to investigate periodic boundary conditions’ and water models’ influence on the orientation dynamics and translational and rotational diffusion of peptides in solution. We have characterized the effects of solvent box size and water model choice on the dynamics of two peptide systems, NATA and WK5. Computationally, translational, and rotational diffusion and internal fluctuations are investigated through all-atom molecular dynamics simulation with two water models and different box sizes. These results are compared with time-resolved fluorescence anisotropy decay (FAD) measurements. The associated time constant and orientation dynamics from FAD measurement along the 1Lb axis provided baseline data to validate molecular dynamics simulation. The modeling results show that diffusion rates vary roughly in inverse proportion to water model viscosity, as one would expect. Corrections for finite box size are significant for translational diffusion and insignificant for rotational diffusion. This study also finds that internal dynamics described by autocorrelation functions and kinetic network models are relatively insensitive to both box size and water model properties. Our observation suggests that different peptide properties respond differently to a change in simulation conditions. Communicated by Ramaswamy H. Sarma

CityApps: A bioinformatics tool for predicting the key residues of enzymes weakly interacting with monovalent metal ions

Monovalent metal ions play an essential role in some enzyme reactions. Due to the loose interactions, identifying the critical residues interacting with monovalent metal ions based on a molecular dynamic (MD) analysis is still a challenge. In this study, we first defined the enrichment factor (Eaa) and the dwell time (td) to describe the interactions with metal ions in each enzyme residue during MD analysis. A "sandwich-like" solvent box for MD analysis was then constructed to avoid the interference of the initial position of metal ions. Finally, a bioinformatics tool, CityApps, was developed to facilitate the automatization of MD analysis and the calculation of the interaction parameter. Using the tool, the key residues (D191, Q272, and N329) responsible for K+-mediated activation of a methyl parathion hydrolase (MPH) from Pseudomonas sp. WBC-3 were identified based on Eaa, td, and RMSF differences. Ala mutation showed that the mutant D191A/Q272A/N329A was difficult to be activated by 1 M K+. In addition, Q272A showed a 176.19 % increase in specific activity compared to the wild MPH type. Our results indicate that CityApps developed here could be an effective tool for identifying enzyme residues interacting with monovalent metal ions and those residues that are critical for catalytic reaction.

Stability of Lipid Coatings on Nanoparticle-Decorated Surfaces.

Lipid membranes supported on solid surfaces and nanoparticles find multiple applications in industrial and biomedical technologies. Here, we explore in silico the mechanisms of the interactions of lipid membranes with nanostructured surfaces with deposited nanoparticles and explain the characteristic particle size dependence of the uniformity and stability of lipid coatings observed in vitro. Simulations are performed to demonstrate the specifics of 1,2-dimyristoyl-sn-glycero-3-phosphocholine (DMPC) lipid membrane adhesion to hydrophilic and hydrophobic nanoparticles ranging in size from 1.5 to 40 nm using an original coarse-grained molecular dynamics model with implicit solvent and large simulation boxes (scales up to 280 × 154 × 69 nm3). We find that one of the major factors that affects the uniformity and stability of lipid coatings is the disjoining pressure in the water hydration layer formed between the lipid membrane and hydrophilic solid surface. This effect is accounted for by introducing a special long-range lipid-solid interaction potential that mimics the effects of the disjoining pressure in thin water layers. Our simulations reveal the physical mechanisms of interactions of lipid bilayers with solid surfaces that are responsible for the experimentally observed nonmonotonic particle size dependence of the uniformity and stability of lipid coatings: particles smaller than the hydration layer thickness (<2-3 nm) or larger than ∼20 nm are partially or fully enfolded by a lipid bilayer, whereas particles of the intermediate size (5-20 nm) cause membrane perforation and pore formation. In contrast, hydrophobic nanoparticles, which repel the hydration layer, tend to be encapsulated within the hydrophobic interior of the membrane and coated by the lipid monolayer. The proposed model can be further extended and applied to a wide class of systems comprising nanoparticles and nanostructured substrates interacting with lipid and surfactant bilayers and monolayers.

The Gibbs free energy of cavity formation in a diverse set of solvents.

The concept of the formation of a solute-sized cavity in a solvent is widely used in the theories of solvation processes; however, most of the studies of cavity formation using atomistic simulations were limited to water and hydrocarbon models. We calculated the Gibbs free energy of cavity formation ΔcavG for a structurally diverse set of 23 common organic solvents. For the calculation, molecular dynamics simulations of solvent boxes were conducted, and the Widom particle insertion method was applied. The results obtained with two different force fields for the same solvent were in good agreement with each other in most cases. The obtained cavity size dependences of ΔcavG allowed ranking the solvents by the free energy cost of creation of a cavity with a certain size. Surprisingly, this cost was somewhat higher in glycerol, formamide, and ethylene glycol than in water. In general, higher values of ΔcavG are observed for the solvents with a branched network of intermolecular hydrogen bonds and strongly polar aprotic solvents. The numerical results can be used to improve the accuracy of the calculation of the cavity term in non-aqueous continuum solvation models.

Straightforward Incorporation of Multiple Ligand Types into Molecular Dynamics Simulations for Efficient Binding Site Detection and Characterization.

Binding site identification and characterization is an important initial step in structure-based drug design. To account for the effects of protein flexibility and solvation, several cosolvent molecular dynamics (MD) simulation methods that incorporate small organic molecules into the protein's solvent box to probe for binding sites have been developed. However, most of these methods are highly inefficient, as they allow for the use of only one probe type at a time, which means that multiple sets of simulations have to be performed to map different types of binding sites. The high probe concentrations used in some of these methods also necessitate the use of artificial repulsive forces to prevent the probes from aggregating. Here, we present multiple-ligand-mapping MD (mLMMD), a method that incorporates multiple types of probes for simultaneous and efficient mapping of different types of binding sites without the need for introduction of artificial forces that may cause unintended mapping artifacts. We validate the method on a diverse set of 10 proteins and show that the mLMMD probes are able to reliably identify hydrophobic, hydrogen-bonding, charged, and cryptic binding sites in all of the test cases. Our results also highlight the potential utility of mLMMD for virtual screening and rational drug design.

Effect of solvent on protein structure and dynamics

Understanding how much solvents influence the structures and dynamics of proteins is important to understand functional mechanisms of solvated proteins. We propose a solvated potential model that approximates the potential energy of a solvated protein by projecting the solvent information into the protein structure. Using the model, we derive three properties of the solvent. First, the influence of the solvent on protein structure and dynamics, mostly by the bulk solvent, decays drastically (near-exponentially) as the distances of the solvent from the protein increase. Using this decay pattern, we suggest the economical size of solvent boxes in molecular dynamics simulations. Second, the hydration shell regulates the protein dynamics by effecting extra interactions within the protein structure. Lastly, the lowest frequency modes are determined mostly by protein structures.

Molecular dynamics simulation on the dissolution process of Kaempferol cluster

Abstract The solid-liquid extraction is an effective method to extract bioactive components from crude medicinal herbs. The dissolution of bioactive molecules is the key to efficient solid-liquid extraction. In this work, a five-nanosecond molecular dynamics simulation was carried out to study the dissolution mechanism of a cluster, including 15 kaempferol molecules in a cubic solvent box of 6 nm. The theoretical and experimental infrared and Raman spectra were taken to confirm the low-energy kaempferol conformers for building the cluster. The effects of ethanol concentration and temperature on the kaempferol dispersion efficiency are investigated. The simulation results show that different concentration of ethanol solution has different dissolution behavior. Because of the intermolecular interaction between kaempferol and solvent, the kaempferol cluster can be decomposed gradually. In order to quantitatively describe the diffusion process of kaempferol cluster, the radial distribution function of the barycenter of the cluster around kaempferol molecules, the number of H-bond between kaempferol molecules, the solvent accessible area of the cluster, the dissolution free energy and the diffusion constant are analyzed. The results show that the optimal dispersion condition is at 80–90% of the ethanol concentration and 323.15 K of the system temperature, which agree with the experiment of extraction. This research provides new insight into the dissolution process of kaempferol on the nanosecond scale. Also, the result of molecular dynamics (MD) can effectively predict the extraction conditions of the actual experiment.

Investigating the solid-liquid extraction process of puerarin by molecular dynamics simulations

The solid-liquid extraction is an efficient method to obtain bioactive components from natural medicinal herbs. The diffusion of the bioactive molecule is the key to efficient solid-liquid extraction. The current research aims to investigate the dissolution mechanism and the dispersal behavior of puerarin in the solid-liquid extraction process by molecular dynamics simulations. A cluster including 20 puerarin molecules immersed into the solvent box with 80 Å thicknesses and 5 ns of MD simulations are performed to all systems. The effects of ethanol concentration and extraction temperature on the molecular dispersion efficiency of puerarin are investigated. The simulation results show various dissolution behaviors in different concentration of ethanol solution. Due to the H-bonding interaction between puerarin and solvents, the puerarin molecules can move gradually out of the cluster. To quantitatively describe the diffusion process of puerarin cluster, the RDFs for oxygen-oxygen (O-O) pair of puerarin cluster, the diffusion coefficient of the cluster, and the number of hydrogen bonds and lifetime of hydrogen bonds between puerarin molecules are calculated. The results show that the best condition of dispersion is at 60% of the ethanol concentration and 323.15 K of the system temperature, which can be fitted with the experiment of extraction. This research provides clear images of the dissolution process of drug molecular dispersions at the molecular level. Also, the MD can effectively predict the extraction conditions of the actual experiment.

Comment on 'Valid molecular dynamics simulations of human hemoglobin require a surprisingly large box size'.

A recent molecular dynamics investigation into the stability of hemoglobin concluded that the unliganded protein is only stable in the T state when a solvent box is used in the simulations that is ten times larger than what is usually employed (El Hage et al., 2018). Here, we express three main concerns about that study. In addition, we find that with an order of magnitude more statistics, the reported box size dependence is not reproducible. Overall, no significant effects on the kinetics or thermodynamics of conformational transitions were observed.

General Purpose Water Model Can Improve Atomistic Simulations of Intrinsically Disordered Proteins.

Unconstrained atomistic simulations of intrinsically disordered proteins and peptides (IDP) remain a challenge: widely used, "general purpose" water models tend to favor overly compact structures relative to experiment. Here we have performed a total of 93 μs of unrestrained MD simulations to explore, in the context of IDPs, a recently developed "general-purpose" 4-point rigid water model OPC, which describes liquid state of water close to experiment. We demonstrate that OPC, together with a popular AMBER force field ff99SB, offers a noticeable improvement over TIP3P in producing more realistic structural ensembles of three common IDPs benchmarks: 55-residue apo N-terminal zinc-binding domain of HIV-1 integrase ("protein IN"), amyloid β-peptide (Aβ42) (residues 1-42), and 26-reside H4 histone tail. As a negative control, computed folding profile of a regular globular miniprotein (CLN025) in OPC water is in appreciably better agreement with experiment than that obtained in TIP3P, which tends to overstabilize the compact native state relative to the extended conformations. We employed Aβ42 peptide to investigate the possible influence of the solvent box size on simulation outcomes. We advocate a cautious approach for simulations of IDPs: we suggest that the solvent box size should be at least four times the radius of gyration of the random coil corresponding to the IDP. The computed free energy landscape of protein IN in OPC resembles a shallow "tub" - conformations with substantially different degrees of compactness that are within 2 kB T of each other. Conformations with very different secondary structure content coexist within 1 kB T of the global free energy minimum. States with higher free energy tend to have less secondary structure. Computed low helical content of the protein has virtually no correlation with its degree of compactness, which calls into question the possibility of using the helicity as a metric for assessing performance of water models for IDPs, when the helicity is low. Predicted radius of gyration ( R g) of H4 histone tail in OPC water falls in-between that of a typical globular protein and a fully denatured protein of the same size; the predicted R g is consistent with two independent predictions. In contrast, H4 tail in TIP3P water is as compact as the corresponding globular protein. The computed free energy landscape of H4 tail in OPC is relatively flat over a significant range of compactness, which, we argue, is consistent with its biological function as facilitator of internucleosome interactions.

Natural payload delivery of the doxorubicin anticancer drug from boron nitride oxide nanosheets

We studied the behavior of doxorubicin (DOX; an anticancer drug) molecules loaded on a boron nitride oxide nanosheet (BNO-NS) using the density functional theory (DFT), time-dependent density functional theory (TDDFT), and molecular dynamic (MD) simulation methods. We found that DOX molecules in π-π or covalent interaction with BNO-NS preserve their optical properties in water. Moreover, the BNO-NS vector allowed stabilizing the DOX molecules on a cellular membrane contrary to isolated DOX that randomly moved in the solvent box without any interaction with the cell membrane. From these results, we conclude that hydrophilic BNO-NS represents a good candidate for DOX molecule transport and stabilization near a cell membrane. In this drug delivery system, the choice of BNO-NS as nanovector is important because it allows delivering an elevated therapeutic dose directly on the cancer cell target without hindrance of the DOX payload.

In Silico Study of Phospholipids as An Oral Insulin Delivery System

Phospholipids have been applied as a material for an oral insulin delivery system that either in the form of solid lipid nanoparticles, liposomes, or microemulsions. However, the administration of this delivery system for diabetic treatment requires twenty times higher of dosage than the injection one. In the present study, we performed molecular dynamics simulation to evaluate how phospholipids can protect insulin at the molecular level and why its efficiency for the diabetic treatment is low. Phosphatydilcholine group, such as DPPC, OPPC, LPC and combine of DPPC-LPC, was used as the phospholipid material to model the delivery system in our simulations because in vitro and in vivo studies of these phospholipids as the delivery system have been well-conducted. In these simulations, the waters-phospholipid-insulin system was prepared by randomly mixing phospholipids molecules inside the solvent box containing an insulin molecule. NPT ensemble was used in all simulations with the setting temperature at 310 K and 1 atm for the pressure. The simulation results showed that DPPC and LPC delivery systems were able to maintain both secondary and tertiary structures of the encapsulated insulin within the time range of simulation. However, all phospholipid delivery system models failed to fully encapsulate insulin surface within 150 ns simulation. This study may explain the reason for the low success rate of using phospholipid as the delivery system of oral insulin in the experiment.

Valid molecular dynamics simulations of human hemoglobin require a surprisingly large box size.

Recent molecular dynamics (MD) simulations of human hemoglobin (Hb) give results in disagreement with experiment. Although it is known that the unliganded (T[Formula: see text]) and liganded (R[Formula: see text]) tetramers are stable in solution, the published MD simulations of T[Formula: see text] undergo a rapid quaternary transition to an R-like structure. We show that T[Formula: see text] is stable only when the periodic solvent box contains ten times more water molecules than the standard size for such simulations. The results suggest that such a large box is required for the hydrophobic effect, which stabilizes the T[Formula: see text] tetramer, to be manifested. Even in the largest box, T[Formula: see text] is not stable unless His146 is protonated, providing an atomistic validation of the Perutz model. The possibility that extra large boxes are required to obtain meaningful results will have to be considered in evaluating existing and future simulations of a wide range of systems.

QM/MM Protocol for Direct Molecular Dynamics of Chemical Reactions in Solution: The Water-Accelerated Diels-Alder Reaction.

We describe a solvent-perturbed transition state (SPTS) sampling scheme for simulating chemical reaction dynamics in condensed phase. The method, adapted from Truhlar and Gao's ensemble-averaged variational transition state theory, includes the effect of instantaneous solvent configuration on the potential energy surface of the reacting system (RS) and allows initial conditions for the RS to be sampled quasiclassically by TS normal mode sampling. We use a QM/MM model with direct dynamics, in which QM forces of the RS are computed at each trajectory point. The SPTS scheme is applied to the acceleration of the Diels-Alder reaction of cyclopentadiene (CP) + methyl vinyl ketone (MVK) in water. We explored the effect of the number of SPTS and of solvent box size on the distribution of bond lengths in the TS. Statistical sampling of the sampling was achieved when distribution of forming bond lengths converged. We describe the region enclosing the partial bond lengths as the transition zone. Transition zones in the gas phase, SMD implicit solvent, QM/MM, and QM/MM+QM (3 water molecules treated by QM) vary according to the ability of the medium to stabilize zwitterionic structures. Mean time gaps between formation of C-C bonds vary from 11 fs for gas phase to 25 fs for QM/MM+QM. Mean H-bond lengths to O(carbonyl) in QM/MM+QM are 0.14 Å smaller at the TS than in MVK reactant, and the mean O(carbonyl)-H(water)-O(water) angle of H-bonds at the TS is 10° larger than in MVK reactant.