Constrained Molecular Dynamics Simulations Research Articles

Clarifying the Direct Generation of •OH Radicals in Photocatalytic O2 Reduction: Theoretical Prediction Combined with Experimental Validation.

This work systematically studied thermocatalytic and photocatalytic pathways of formaldehyde degradation and H-assisted O2 reduction over a Pt13/anatase-TiO2(101) composite via DFT calculations together with constrained molecular dynamics (MD) simulations. We show that photocatalytic O2 reduction on Pt/TiO2 can directly generate •OH radicals (*O2 → *OOH → •OH) via two hydrogenation steps with small barriers, and the product selectivity (*H2O2 or •OH) is decided by the relative position between catalyst Fermi level and •OH/*H2O2 redox potential (theoretical determination of 0.07 V referencing to the SHE). Such a novel reaction channel was furthermore validated at the liquid-solid interface via constrained MD simulations and experimental electron paramagnetic resonance detections, and a wide range of H resources, e.g., *HCHO, *HCO, *H (H+ + e-), can always drive the direct •OH generation. The additional portion of e--triggered •OH radicals are prone to diffuse into solution or the TiO2 surface and furthermore cooperate with the conventional h+-driven photooxidations.

Towards realistic structural char models: Generation of stacked graphene-like layers using constrained reactive molecular dynamics simulations

For a molecular understanding of coal combustion, realistic coal and char models are pivotal. We propose a novel top-down method for the construction of char models, consisting of stacked graphene-like layers with defects. We start on the mesoscale with corrugated stacked surfaces as a template. Then the carbon atoms are pre-positioned on these surfaces and surrounded by hydrogen atoms. In a constrained reactive MD simulation the carbon atoms are then rearranged into a chemically reasonable state with dangling bonds saturated by hydrogen atoms. Importantly, the initial shape of the template is maintained to a great extent, leading to a substantial amount of defects in the graphitic planes, especially for high curvatures. In a final step, small molecular fragments and unfavorable arrangements are removed and remaining dangling bonds are saturated with hydrogen atoms. The central advantage over other approaches is the ability to choose the shape and size of the templating planes, with the possibility to use experimental input.The method is validated on single sheets with a protrusion and application examples are shown for stacked graphitic layers with up to 100 Å size and more than 20000 atoms. It represents the first step of a longer term project. In the future, heteroatoms and functional groups will be implemented, and by combining several such stacks, the approach will allow the construction of more realistic char models using the data from HRTEM images as an experimental boundary condition. On this basis important questions like the adsorption and diffusion of volatiles under high temperatures can be studied on a molecular level.

Metal–N4 model single‐atom catalyst with electroneutral quadri‐pyridine macrocyclic ligand for CO2 electroreduction

AbstractMetal–N–C single‐atom catalysts, mostly prepared from pyrolysis of metal‐organic precursors, are widely used in heterogeneous electrocatalysis. Since metal sites with diverse local structures coexist in this type of material and it is challenging to characterize the local structure, a reliable structure–property relationship is difficult to establish. Conjugated macrocyclic complexes adsorbed on carbon support are well‐defined models to mimic the single‐atom catalysts. Metal–N4 site with four electroneutral pyridine‐type ligands embedded in a graphene layer is the most commonly proposed structure of the active site of single‐atom catalysts, but its molecular counterpart has not been reported. In this work, we synthesized the conjugated macrocyclic complexes with a metal center (Co, Fe, or Ni) coordinated with four electroneutral pyridinic ligands as model catalysts for CO2 electroreduction. For comparison, the complexes with anionic quadri‐pyridine macrocyclic ligand were also prepared. The Co complex with the electroneutral ligand expressed a turnover frequency of CO formation more than an order of magnitude higher than that of the Co complex with the anionic ligand. Constrained ab initio molecular dynamics simulations based on the well‐defined structures of the model catalysts indicate that the Co complex with the electroneutral ligand possesses a stronger ability to mediate electron transfer from carbon to CO2.

Clumped isotope reordering kinetics in strontianite and witherite: Experiments and first-principles simulations

The carbonate clumped isotope (∆47) thermometer is a powerful tool for temperature reconstruction based on the thermodynamical 13C—18O bonding preference. Accurate interpretation of temperatures requires preservation of original carbonate ∆47 signatures, which are susceptible to resetting at elevated temperatures during burial. Kinetic models of carbonate 13C—18O bond reordering derived from experimental studies constrain the evaluation of ∆47 preservation and the reconstruction of thermal histories and initial temperatures for reordered carbonates. However, application of these models to different carbonates is limited because kinetic parameters vary with materials and are believed to be affected by factors such as crystal structure, cation, trace element composition, crystal defects, and internal water whose roles in reordering are still not well-understood. Here we present reordering data for strontianite and witherite from heating experiments and evaluate the effects of crystal structure, cation chemistry, cation vacancy and substitution, and internal water on reordering using constrained ab-initio molecular dynamics (cAIMD) simulations. We show that strontianite and witherite reorder significantly faster than other carbonates. Comparing experimental studies yields the following reordering rate sequence: strontianite > witherite > aragonite > calcite > dolomite. The general trend in Arrhenius activation energy for reordering is strontianite < witherite < aragonite < dolomite < calcite (with calcite values based on the average of individual materials). Our cAIMD simulations reveal a similar trend in free energy of activation for oxygen exchange in these carbonates. We suggest that orthorhombic carbonates are more susceptible to reordering than trigonal carbonates because of the shorter distance between carbonate groups and because the weaker M—O interaction in orthorhombic carbonates reduces the structural distortion and the energy cost during transition state formation. For carbonates sharing the same structure, the weak M—O bonding can be compensated by a large unit cell volume, leading to a non-systematic relationship between the free energy of activation and the cation radius. In simulations with lattice-bound water, the free energy of activation is greatly reduced compared with cation vacancy and substitution. The influence of fluid inclusions on measured reordering rates depends on the fraction of carbonates that can exchange with water as suggested by strontianite δ¹⁸O data. Our results serve as a starting point for a more sophisticated but flexible kinetic model containing material-specific factors (e.g., substituted ions, structure water, etc.) for reordering in abiotic and biotic carbonates.

Evaluation of Constrained and Restrained Molecular Dynamics Simulation Methods for Predicting Skin Lipid Permeability.

Recently, molecular dynamics (MD) simulations have been utilized to investigate the barrier properties of human skin stratum corneum (SC) lipid bilayers. Different MD methods and force fields have been utilized, with predicted permeabilities varying by few orders of magnitude. In this work, we compare constrained MD simulations with restrained MD simulations to obtain the potential of the mean force and the diffusion coefficient profile for the case of a water molecule permeating across an SC lipid bilayer. Corresponding permeabilities of the simulated lipid bilayer are calculated via the inhomogeneous solubility diffusion model. Results show that both methods perform similarly, but restrained MD simulations have proven to be the more robust approach for predicting the potential of the mean force profile. Critical to both methods are the sampling of the whole trans-bilayer axis and the following symmetrization process. Re-analysis of the previously reported free energy profiles showed that some of the discrepancies in the reported permeability values is due to misquotation of units, while some are due to the inaccurately obtained potential of the mean force. By using the existing microscopic geometrical models via the intercellular lipid pathway, the permeation through the whole SC is predicted from the MD simulation results, and the predicted barrier properties have been compared to experimental data from the literature with good agreement.

Two-Channel Model for Electron Transfer in a Dye-Catalyst-Dye Supramolecular Complex for Photocatalytic Water Splitting.

To improve the performance of dye‐sensitized photoelectrochemical cell (DS‐PEC) devices for splitting water, the tailoring of the photocatalytic four‐photon water oxidation half‐reaction represents a principle challenge of fundamental significance. In this study, a Ru‐based water oxidation catalyst (WOC) covalently bound to two 2,6‐diethoxy‐1,4,5,8‐diimide‐naphthalene (NDI) dye functionalities provides comparable driving forces and channels for electron transfer. Constrained ab initio molecular dynamics simulations are performed to investigate the photocatalytic cycle of this two‐channel model for photocatalytic water splitting. The introduction of a second light‐harvesting dye in the Ru‐based dye‐WOC‐dye supramolecular complex enables two separate parallel electron‐transfer channels, leading to a five‐step catalytic cycle with three intermediates and two doubly oxidized states. The total spin S=1 is conserved during the catalytic process and the system with opposite spin on the oxidized NDI proceeds from the Ru=O intermediate to the final Ru‐O2 intermediate with a triplet molecular 3O2 ligand that is eventually released into the environment. The in‐depth insight into the proposed photocatalytic cycle of the two‐channel model provides a strategy for the development of novel high‐efficiency supramolecular complexes for DS‐PEC devices with buildup and conservation of spin multiplicity along the reaction coordinate as a design principle.

Unfolding bovine α-lactalbumin with T-jump: Characterizing disordered intermediates via time-resolved x-ray solution scattering and molecular dynamics simulations.

The protein folding process often proceeds through partially folded transient states. Therefore, a structural understanding of these disordered states is crucial for developing mechanistic models of the folding process. Characterization of unfolded states remains challenging due to their disordered nature, and incorporating multiple methods is necessary. Combining the time-resolved x-ray solution scattering (TRXSS) signal with molecular dynamics (MD), we are able to characterize transient partially folded states of bovine α-lactalbumin, a model system widely used for investigation of molten globule states, during its unfolding triggered by a temperature jump. We track the unfolding process between 20 µs and 70 ms and demonstrate that it passes through three distinct kinetic states. The scattering signals associated with these transient species are then analyzed with TRXSS constrained MD simulations to produce protein structures that are compatible with the input signals. Without utilizing any experimentally extracted kinetic information, the constrained MD simulation successfully drove the protein to an intermediate molten globule state; signals for two later disordered states are refined to terminal unfolded states. From our examination of the structural characteristics of these disordered states, we discuss the implications disordered states have on the folding process, especially on the folding pathway. Finally, we discuss the potential applications and limitations of this method.

Echo dephasing and heat capacity from constrained and unconstrained dynamics of triiodothyronine nuclear receptor protein.

The objective of this study is to observe the echo feature curves, vibrational dephasing, and heat capacity of a protein-hormone system taking thyroid hormone receptor-beta (THR-β) as an example. Constrained and unconstrained molecular dynamics simulations are performed by implementing the theory of velocity reassignments to probe the phase coherent state in terms of echo pulses. The constrained vibrations are incorporated by adjusting rigid bonds to all hydrogen atoms with an integrator parameter of 2fs/step in order to reduce the degrees of freedom whereas 1fs/step is used in the free vibrations of the atomic cluster. The nature of temperature auto-correlation functions changes so that echo feature curves also show a distinct nature in the cases of constrained and unconstrained vibrations. There is a large variation in kinetic temperature and internal potential energy in the echo time zone. The temperature rate of change of internal potential energy is the main contributor to the heat capacity of the native state protein-hormone system. The heat capacity of proteins estimated from this technique is in good agreement with the values from experiments. This study shows that triiodothyronine (T3) hormone makes some differences in heat capacity upon binding to the THR-β ligand binding domain (LBD). The physical properties of unliganded THR-β and T3-bound THR-β LBD in the cases of constrained and unconstrained dynamics are observed distinctly under the effect of anharmonicity on the phase coherent state of normal modes and the dephasing time lies in a range of 0.6-0.8ps when the systems are perturbed suddenly.

In-Silico Skin Model: A Multiscale Simulation Study of Drug Transport.

Accurate in-silico models are required to predict the release of drug molecules through skin in order to supplement the in-vivo experiments for faster development/testing of drugs. The upper most layer of the skin, stratum corneum (SC), offers the main resistance for permeation of actives. Most of the SC's molecular level models comprise cholesterol and phospholipids only, which is far from reality. In this study we have implemented a multiscale modeling framework to obtain the release profile of three drugs, namely, caffeine, fentanyl, and naphthol, through skin SC. We report for the first time diffusion of drugs through a realistic skin molecular model comprised of ceramides, cholesterol, and free fatty acid. The diffusion coefficients of drugs in the SC lipid matrix were determined from multiple constrained molecular dynamics simulations. The calculated diffusion coefficients were then used in the macroscopic models to predict the release profiles of drugs through the SC. The obtained release profiles were in good agreement with available experimental data. The partition coefficient exhibits a greater effect on the release profiles. The reported multiscale modeling framework would provide insight into the delivery mechanisms of the drugs through the skin and shall act as a guiding tool in performing targeted experiments to come up with a suitable delivery system.

Conformational Study on Pheromonotropin neuropeptide using NMR and Molecular Dynamics

Pheromonotropic neuropeptide, Pseudaletia pheromonotropin, is an 18 amino acid peptide with the sequence Lys-Leu-Ser-Tyr-Asp-Asp-Lys-Val-Phe-Glu-Asn-Val-Glu-Phe-Thr-Pro-Arg-Leu which is widely distributed in female moths of Bombyx mori. This peptide is structurally related to leucopyrokinin, an insect myotropic neuropeptide, responsible for melanization and reddish coloration hormone (MRCH) activity. A combined NMR and Molecular Dynamics (MD) methods were used to fathom the conformational behaviour in water which was used as a solvent. The structure was investigated using 1D-NMR and 2D-NMR (COSY, TOCSY and ROESY) experiments. The conformation was built by constrained MD simulations using distance and dihedral restraints from NMR data in the GROMACS simulation package. The peptide predominately adopts a β-sheet structure in water.

One-dimensional potential of mean force underestimates activation barrier for transport across flexible lipid membranes

Transport of a fullerene-like nanoparticle across a lipid bilayer is investigated by coarse-grained molecular dynamics (MD) simulations. Potentials of mean force (PMF) acting on the nanoparticle in a flexible bilayer suspended in water and a bilayer restrained to a flat surface are computed by constrained MD simulations. The rate of the nanoparticle transport into the bilayer interior is predicted using one-dimensional Langevin models based on these PMFs. The predictions are compared with the transport rates obtained from a series of direct (unconstrained) MD simulations of the solute transport into the flexible bilayer. It is observed that the PMF acting on the solute in the flexible membrane underestimates the transport rate by more than an order of magnitude while the PMF acting on the solute in the restrained membrane yields an accurate estimate of the activation energy for transport into the flexible membrane. This paradox is explained by a coexistence of metastable membrane configurations for a range of the solute positions inside and near the flexible membrane. This leads to a significant reduction of the contribution of the transition state to the mean force acting on the solute. Restraining the membrane shape ensures that there is only one stable membrane configuration corresponding to each solute position and thus the transition state is adequately represented in the PMF. This mechanism is quite general and thus this phenomenon is expected to occur in a wide range of interfacial systems. A simple model for the free energy landscape of the coupled solute-membrane system is proposed and validated. This model explicitly accounts for effects of the membrane deformations on the solute transport and yields an accurate prediction of the activation energy for the solute transport.

Structural Analysis of Carbon-Added Na–Ga Melts in Na Flux GaN Growth by First-Principles Calculation

We investigated the fundamentals of the effect of C addition on Na flux GaN growth by first-principles calculation. We simulated C-added Na–Ga melts using molecular dynamics (MD) simulations to examine the local melt structure around a N atom. We also calculated C–N bond energy using constrained MD simulations. Results show that a N atom bonded to a C atom and there were no Ga atoms around the N atom because C–N bond energy was larger than Ga–N bond energy. This is the reason for the suppression of heterogeneous nucleation by C addition. It was also found that the C–N bond energy was affected by surrounding Ga atoms and that the C–N atomic distance increased with the Ga coordination number around the N atom.

Collective degrees of freedom involved in absorption and desorption of surfactant molecules in spherical non-ionic micelles

Dynamics of absorption and desorption of a surfactant monomer into and out of a spherical non-ionic micelle is investigated by coarse-grained molecular dynamics (MD) simulations. It is shown that these processes involve a complex interplay between the micellar structure and the monomer configuration. A quantitative model for collective dynamics of these degrees of freedom is developed. This is accomplished by reconstructing a multi-dimensional free energy landscape of the surfactant-micelle system using constrained MD simulations in which the distance between the micellar and monomer centers of mass is held constant. Results of this analysis are verified by direct (unconstrained) MD simulations of surfactant absorption in the micelle. It is demonstrated that the system dynamics is likely to deviate from the minimum energy path on the energy landscape. These deviations create an energy barrier for the monomer absorption and increase an existing barrier for the monomer desorption. A reduced Fokker-Planck equation is proposed to model these effects.

Folding of Small Proteins Using Constrained Molecular Dynamics

The focus of this paper is to examine whether conformational search using constrained molecular dynamics (MD) method is more enhanced and enriched toward "native-like" structures compared to all-atom MD for the protein folding as a model problem. Constrained MD methods provide an alternate MD tool for protein structure prediction and structure refinement. It is computationally expensive to perform all-atom simulations of protein folding because the processes occur on a time scale of microseconds. Compared to the all-atom MD simulation, constrained MD methods have the advantage that stable dynamics can be achieved for larger time steps and the number of degrees of freedom is an order of magnitude smaller, leading to a decrease in computational cost. We have developed a generalized constrained MD method that allows the user to "freeze and thaw" torsional degrees of freedom as fit for the problem studied. We have used this method to perform all-torsion constrained MD in implicit solvent coupled with the replica exchange method to study folding of small proteins with various secondary structural motifs such as, α-helix (polyalanine, WALP16), β-turn (1E0Q), and a mixed motif protein (Trp-cage). We demonstrate that constrained MD replica exchange method exhibits a wider conformational search than all-atom MD with increased enrichment of near-native structures. "Hierarchical" constrained MD simulations, where the partially formed helical regions in the initial stretch of the all-torsion folding simulation trajectory of Trp-cage were frozen, showed a better sampling of near-native structures than all-torsion constrained MD simulations. This is in agreement with the zipping-and-assembly folding model put forth by Dill and co-workers for folding proteins. The use of hierarchical "freeze and thaw" clustering schemes in constrained MD simulation can be used to sample conformations that contribute significantly to folding of proteins.

Heirarchical Constrained Molecular Dynamics Simulations for Proteins

Here we report a constrained molecular dynamics method that allows molecular models ranging from all-torsion to freezing parts of the protein as rigid body connected by flexible hinges, and studying the dynamics of proteins. This method is known as GNEIMO (Generalized Newton-Euler inverse mass operator method). We have derived new algorithms to assign initial velocities in the torsional space that obey the Boltzmann distribution. The effect of assigning initial velocities in the dihedral space (as opposed to transforming the Cartesian velocities to dihedral velocities), on the time taken for equilibrating the molecular system is analyzed.

Solute–solvent friction kernels and solution properties of methyl oxazoline–phenyl oxazoline (MeOx–PhOx) copolymers in binary ethanol–water mixtures

Solvent mixtures often alter the solubility of polymeric substances. Statistical copolymers made from 2-methyl-2-oxazoline (MeOx) and 2-phenyl-2-oxazoline (PhOx) are known for their varying solubilities in pure ethanol, pure water and in binary mixtures of ethanol-water. Constrained Molecular Dynamics (MD) simulations have been carried out with an aim to explain the varying solubilities of the statistical MeOx-PhOx copolymers. The solute-solvent dynamic friction kernels calculated through constrained MD simulations corroborate the solubility pattern in these copolymers. The solvation characteristics have been analyzed in terms of the solute-solvent radial distribution functions (RDFs). The ethanol-soluble MeOx-PhOx copolymers exhibit characteristic solute-composition dependence in the dynamic solute-solvent friction kernels, indicating the strength of the solute-solvent correlations. The aggressive solvation by the ethanol molecules in the binary solvent mixtures has been brought out by the O(solute)-H(ethanol) RDFs which exhibit a characteristic dependence on the ethanol content in the solvent composition. The corresponding O(solute)-H(water) RDFs are devoid of any such composition dependence. For all the MeOx-PhOx copolymers, the O-site solvation is strongly dominated by the water molecules and the N-sites are solvated equally by both ethanol and water molecules.

The Conformation and Assignment of the Proton NMR Spectrum in Water of DX600, a Bioactive Peptide with a Random Coil Conformation

DX600, a small peptide with 26 residues, is a potent, highly selective inhibitor of angiotensin converting enzyme 2 (ACE2). A range of NMR methods including TOCSY and ROESY yield an assignment of its proton spectrum in water and constraints on its conformation. Constrained molecular dynamics simulations of solvated DX600 show that the peptide′s most abundant conformer adopts a predominantly random coil conformation. Constrained by the disulfide bond, its backbone defines an overhand knot with frayed ends.

Hydration Dynamics As Revealed By The Fluorescence Stokes Shift: The Origin Of Slow Hydration Dynamics And Breakdown Of Linear Response

Hydration dynamics in the immediate vicinity of a protein, probed by time-dependent fluorescence Stokes shift experiments, is critical to understand its biological function. Protein Stokes shifts typically exhibit biphasic relaxation following photo-excitation: fast relaxation occurs on a time scale of several picoseconds while slower components indicate additional hydration dynamics on a time scale of tens of picoseconds, or longer. Theoretical studies using both linear response and non-equilibrium molecular dynamics (MD) calculation qualitatively reproduce the observed biphasic behavior of time dependent Stokes shift for Trp-7 (W7) in myoglobin. Comparison with constrained MD simulations with protein frozen at the instant of photo-excitation reveals the molecular mechanism of slow hydration process and establishes the critical role of protein flexibility. Coupled protein-water motion is shown to be necessary for the observation of the slow component of hydration dynamics. Qualitatively similar results are found for a series of additional cases, such as monellin and staph. nuclease. We illustrate why tracking the separate contributions to the Stokes shift without constrained MD studies may not yield an accurate interpretation of protein hydration dynamics. Additionally, we examine the extent to which protein fluctuations obey Gaussian statistics and the linear response approximation to the Stokes shift is valid. Equilibrium fluctuations of the ground-excited energy difference, which control the absorption and fluorescence line shapes, in the ground and excited electronic states are not independent of each other. We illustrate how small differences from Gaussian statistics in one electronic state can be a signature of very significant deviations from linear response theory, such as isomerization, in the other electronic state.

Density functional study of aqueous uranyl(VI) fluoride complexes

Mixed uranyl aquo fluoro complexes [UO 2(H 2O) x F y ] 2− y ( y = 1–4; x + y = 4, 5) have been optimized with BLYP and B3LYP density functionals in vacuo and in a polarizable continuum modeling bulk water, and have been studied at the BLYP level with Car–Parrinello molecular dynamics (MD) simulations. Using constrained MD simulations and thermodynamic integration, the computed free binding energy between aqueous uranyl and fluoride, affording [UO 2(H 2O) 4F] +, is in excellent agreement with experiment. With the same technique, five-coordinate [UO 2F 4(H 2O)] 2− is indicated to be unstable against loss of the water ligand, as the free energy for dissociation is computed to be ca. −7 kcal/mol in aqueous solution.

Density Functional Studies of the Hydrolysis of Aluminum (Chloro)Hydroxide in Water with CPMD and COSMO

Car-Parrinello molecular dynamics (CPMD) and the static density functional method (DFT) with a conductor-like screening model (COSMO) were used to investigate the chemistry of aluminum (chloro)hydroxide in water. With these methods, the stability, reactivity, and acidic nature of the chosen chlorohydrate were able to be determined. Constrained molecular dynamics simulations were used to investigate the binding of chlorine in an aquatic environment. According to the results, aluminum preferred to be 5-fold-coordinated. In addition, the activation energy barriers for the dissociation of chlorine atoms from the original chlorohydrate structure were able to be determined. The actual values for the barriers were 14 +/- 3 and 40 +/- 5 kJ mol (-1). The results also revealed the acidity of the original cationic dimer. DFT with COSMO was used to determine free energy differences for the reactions detected in the molecular dynamic simulations. In conclusion, new results and insight into the aquatic chemistry of the aluminum (chloro)hydroxides are provided.