- Graduated (1985) from St.Petersburg (Leningrad) State University, Physical Faculty , Dept. Quantum Field Theory & Statistical Physics.
- Ph.D.(1988) St.Petersburg State University, Physical Faculty, Molecular Biophysics Department
- 1988-1992: researcer at Scientific Research Institute of Physics, St.Petersburg State University, Molecular Biophysics Department; Molecular simulations group
- 1993-2001: research associate Division of Physical Chemistry, Stockholm University
- 2001-2007: Associate Professor (docent) Div. of Physical Chemistry, Stockholm University
- from 2007: Professor Div. of Physical Chemistry, Departnment of Materials and Environmental Chemistry, Stockholm University
I am teaching the following courses:
- Physical Chemistry, 15 hp: a ground level course for 2nd year students
- Chemical modeling, 7.5 hp: An advanced course for 3-rd year bachelor or master students
Students of physical chemistry and related specialities (material chemistry, quantum chemistry, biophysics, etc) are welcomed to apply for graduate research projects (examensarbete) in the field of molecular simulations. The content of the work may be discussed in each individual case.
Molecular Computer Simulations
All material things around us (and we ourselves) consist of molecules. Molecules consist of atoms, atoms consist of electrons and atomic nuclei. Behavior of all these small parts of materia is determined by the fundamental laws of quantum mechanics. Molecular motion at "normal" conditions can be described by the laws of Newtonian mechanics. By knowing these laws, it is not difficult to predict what would happens with two atoms taken separately from the rest of the world. This can be done even with the help of pen and paper, as it was done before the 20th century, for example in calculations of the planetary motion. When the number of atoms in the studied system increases, calculations on a paper become no longer feasible. Here computers come to help. With modern computers, it is possible to compute the force acting on each atom in a system of many thousand atoms, and using the second law of Newton predict the motion of these atoms and thus behavior of the whole system. In cases when precision of the classical mechanics is not satisfactory, laws of quantum mechanics can be used though these laws require more extensive calculations. For larger molecular systems, some more simplified description than the atomistic level can be used. Important is that if computer simulations use only fundamental laws of quantum, classical and statistical mechanics, they are able to produce very accurate and detailed information about the studied molecular system, information which is often impossible or unacceptable expensive to extract from an experiment. Also, computer simulations are indispensable in visualization of molecular motion and in interpretation of experimental data. It is not surprisingly that the number of works in the area of computer simulations grows much faster that the average, and now computer simulations are routinely used even in the applied research, such as drug design or development of technological processes. Because of constantly growing range of applications, further development of simulation methods, algorithms and software is necessary to address to appearing demands. This is what research of the Molecular Computer Simulations group are about.
Development of molecular modeling methodology
Force field development for:
Lipid bilayers: Slipids Force Field
Inorganic surfaces in aqqueous media: TiO2, ZnO
Adsorption of gas molecules in zeolites
Inverse Monte Carlo
Advanced Sampling Methods and free energy calculations:
MDynamix: molecular dynamics simulation software
MagiC : software package for systematic multiscale coarse-graining
Ab-initio simulations of metal oxide surfaces
Modeling of biomembranes / lipid bilayers
DNA and DNA folding in chromatine
Interactions of biomolecules with inorganic surfaces
Modelling of nanotoxicity
Gas adsorbtion in zeolites
Polymorphic forms of drugs
I urval från Stockholms universitets publikationsdatabas
Coarse-Grained Simulation of Rodlike Higher-Order Quadruplex Structures at Different Salt Concentrations
2017. Matúš Rebič (et al.). ACS Omega 2 (2), 386-396Artikel
We present a coarse-grained (CG) model of a rodlike higher-order quadruplex with explicit monovalent salts, which was developed from radial distribution functions of an underlying reference atomistic molecular dynamics simulation using inverse Monte Carlo technique. This work improves our previous CG model and extends its applicability beyond the minimal salt conditions, allowing its use at variable ionic strengths. The strategies necessary for the model development are clearly explained and discussed. The effects of the number of stacked quadruplexes and varied salt concentration on the elasticity of the rodlike higher-order quadruplex structures are analyzed. The CG model reproduces the deformations of the terminal parts in agreement with experimental observations without introducing any special parameters for terminal beads and reveals slight differences in the rise and twist of the G-quartet arrangement along the studied biopolymer. The conclusions of our study can be generalized for other G-quartet-based structures.
Extension of the Slipids Force Field to Polyunsaturated Lipids
2016. Inna Ermilova, Alexander P. Lyubartsev. Journal of Physical Chemistry B 120 (50), 12826-12842Artikel
The all-atomic force field Slipids (Stockholm Lipids) for lipid bilayers simulations has been extended to polyunsaturated lipids. Following the strategy adopted in the development of previous versions of the Slipids force field, the parametrization was essentially based on high-level ab initio calculations. Atomic charges and torsion angles related to polyunsaturated lipid tails were parametrized using structures of dienes molecules. The new parameters of the force field were validated in simulations of bilayers composed of seven polyunsaturated lipids. An overall good agreement was found with available experimental data on the areas per lipids, volumetric properties of bilayers, deuterium order parameters, and scattering form factors. Furthermore, simulations of bilayers consisting of highly polyunsaturated lipids and cholesterol molecules have been carried out. The majority of cholesterol molecules were found in a position parallel to bilayer normal with the hydroxyl group directed to the membrane surface, while a small fraction of cholesterol was found in the bilayer center parallel to the membrane plane. Furthermore, a tendency of cholesterol molecules to form chain-like clusters in polyunsaturated bilayers was qualitatively observed.
Equilibrium properties of 3-arm star-shaped polyions
2017. I. A. Silanteva (et al.). Nanosystems: Physics, Chemistry, Mathematics 8 (1), 108-120Artikel
The entropic sampling Monte Carlo method within Wang-Landau algorithm is applied to investigate properties of a lattice model of strongly charged flexible 3-arm star-shaped polyelectrolyte. The density of states is calculated, from which the canonical properties of the system in a wide temperature range are obtained by simple integration. The effects of the arm length and the short-range monomer-monomer potential on the thermal and structural properties of star polyions are studied. We calculate such characteristics as mean square radius of gyration and its components, the radius vector of the center of mass, components of the tensor of inertia and parameters characterizing the shape of the polyion. In this work, we focus on how these characteristics are influenced by the change of the reduced temperature which, within the considered model, is a parameter combining the effect of real temperature, linear charge density and solvent dielectric permittivity. The coil-globule transition is observed in most of the considered cases, and for the polyions with the longest length of arms ( 24), the transition from a liquid globule to a solid-like state is observed. Comparison of polyelectrolyte models with neutral ones is given.
Base sequence specificity of counterion binding to DNA
2016. Alessio Atzori (et al.). Canadian journal of chemistry (Print) 94 (12), 1181-1188Artikel
Nucleic acids are highly charged biopolymers whose secondary structure is strongly dependent on electrostatic interactions. Solvent molecules and ions are also believed to play an important role in mediating and directing both sequence recognition and interactions with other molecules, such as proteins and a variety of ligands. Therefore, to fully understand the biological functions of DNA, it is necessary to understand the interactions with the surrounding counterions. It is well known that monovalent counterions can bind to the minor groove of DNA with consecutive sequences of four, or more, adenine and thymine (A-tracts) with relatively long residence times. However, much less is known about their binding to the backbone and to the major groove. In this work, we used molecular dynamics simulations to both investigate the interactions between the backbone and major groove of DNA and one of its physiological counterions (Na+) and evaluate the relationship between these interactions and the nucleotide sequence. Three dodecamers, namely CGAAAATTTTCG, CGCTCTAGAGCG, and CGCGAATTCGCG, were simulated using the Toukan-Rahman flexible SPC water model and Smith and Dang parameters for Na+, revealing a significant sequence dependence on the ion binding to both backbone and major groove. In the absence of experimental data on the atomistic details of the studied interactions, the reliability of the results was evaluated performing the simulations with additional sets of potential parameters for ions and solvent, namely the A. qvist or the Joung and Cheatham ion parameters and the TIP3P water model. This allowed us to evaluate the results by verifying which features are preserved independently from the parameters adopted.
Force Field Development for Lipid Membrane Simulations
2016. Alexander P. Lyubartsev, Alexander L. Rabinovich. Biochimica et Biophysica Acta - Biomembranes 1858 (10), 2483-2497Artikel
With the rapid development of computer power and wide availability of modelling software computer simulations of realistic models of lipid membranes, including their interactions with various molecular species, polypeptides and membrane proteins have become feasible for many research groups. The crucial issue of the reliability of such simulations is the quality of the force field, and many efforts, especially in the latest several years, have been devoted to parametrization and optimization of the force fields for biomembrane modelling. In this review, we give account of the recent development in this area, covering different classes of force fields, principles of the force field parametrization, comparison of the force fields, and their experimental validation. This article is part of a Special Issue entitled: Biosimulations edited by Ilpo Vattulainen and Tomasz Rog.
Multiscale coarse-grained modelling of chromatin components
2016. Nikolay Korolev, Lars Nordenskiöld, Alexander P. Lyubartsev. Advances in Colloid and Interface Science 232, 36-48Artikel
To model large biomolecular systems, such as cell and organelles an atomistic description is not currently achievable and is not generally practical. Therefore, simplified coarse-grained (CG) modelling becomes a necessity. One of the most important cellular components is chromatin, a large DNA-protein complex where DNA is highly compacted. Recent progress in coarse graining modelling of the major chromatin components, double helical DNA and the nucleosome core particle (NCP) is presented. First, general principles and approaches allowing rigorous bottom-to-top generation of interaction potentials in the CG models are presented. Then, recent CG models of DNA are reviewed and their adequacy is benchmarked against experimental data on the salt dependence of DNA flexibility (persistence length). Furthermore, a few recent CG models of the NCP are described and their application for studying salt-dependent NCP-NCP interaction is discussed. An example of a multiscale approach to CG modelling of chromatin is presented where interactions and self-assembly of thousands of NCPs in solution are observed.
Reactive wetting properties of TiO2 nanoparticles predicted by ab initio molecular dynamics simulations
2016. Erik G. Brandt, Lorenzo Agosta, Alexander P. Lyubartsev. Nanoscale 8 (27), 13385-13398Artikel
Small-sized wet TiO2 nanoparticles have been investigated by ab initio molecular dynamics simulations. Chemical and physical adsorption of water on the TiO2-water interface was studied as a function of water content, ranging from dry nanoparticles to wet nanoparticles with monolayer coverage of water. The surface reactivity was shown to be a concave function of water content and driven by surface defects. The local coordination number at the defect was identified as the key factor to decide whether water adsorption proceeds through dissociation or physisorption on the surface. A consistent picture of TiO2 nanoparticle wetting at the microscopic level emerges, which corroborates existing experimental data and gives further insight into the molecular mechanisms behind nanoparticle wetting. These calculations will facilitate the engineering of metal oxide nanoparticles with a controlled catalytic water activity.
Transferable force-field for modelling of CO2, N-2, O-2 and Ar in all silica and Na+ exchanged zeolites
2016. Bojan Vujić, Alexander P. Lyubartsev. Modelling and Simulation in Materials Science and Engineering 24 (4)Artikel
In this work we propose a new force field for modelling of adsorption of CO2, N-2, O-2 and Ar in all silica and Na+ exchanged Si-Al zeolites. The force field has a standard molecular-mechanical functional form with electrostatic and Lennard-Jones interactions satisfying Lorentz-Berthelot mixing rules and thus has a potential for further extension in terms of new molecular types. The parameters for the zeolite framework atom types are optimized by an iterative procedure minimizing the difference with experimental adsorption data for a number of different zeolite structures and Si:Al ratios. The new force field shows a good agreement with available experimental data including those not used in the optimization procedure, and which also shows a reasonable transferability within different zeolite topologies. We suggest a potential usage in screening of different zeolite structures for carbon capture and storage process, and more generally, for separation of other gases.
Binding energy calculations for hevein-carbohydrate interactions using expanded ensemble molecular dynamics simulations
2015. Chaitanya A. K. Koppisetty (et al.). Journal of Computer-Aided Molecular Design 29 (1), 13-21Artikel
Accurate estimation of protein-carbohydrate binding energies using computational methods is a challenging task. Here we report the use of expanded ensemble molecular dynamics (EEMD) simulation with double decoupling for estimation of binding energies of hevein, a plant lectin with its monosaccharide and disaccharide ligands GlcNAc and (GlcNAc)(2), respectively. In addition to the binding energies, enthalpy and entropy components of the binding energy are also calculated. The estimated binding energies for the hevein-carbohydrate interactions are within the range of +/- 0.5 kcal of the previously reported experimental binding data. For comparison, binding energies were also estimated using thermodynamic integration, molecular dynamics end point calculations (MM/GBSA) and the expanded ensemble methodology is seen to be more accurate. To our knowledge, the method of EEMD simulations has not been previously reported for estimating biomolecular binding energies.
Molecular Dynamics Simulations of Adsorption of Amino Acid Side Chain Analogues and a Titanium Binding Peptide on the TiO2 (100) Surface
2015. Erik G. Brandt, Alexander P. Lyubartsev. The Journal of Physical Chemistry C 119 (32), 18126-18139Artikel
Adsorption profiles and adsorption free energies were determined for the side chain analogues of the 20 naturally occurring amino acids and a titanium binding peptide on the TiO2 (100) surface. Microsecond simulations with umbrella sampling and metadynamics were used to sample the free energy barriers associated with desolvation of strongly bound water molecules at the TiO2 surface. Polar and aromatic side chain analogues that hydrogen bond either to surface waters or directly to the metal oxide surface were found to be the strongest binders. Further, adsorption simulations of a 6 residue titanium binding peptide identified two binding modes on TiO2 (100). The peptide structure with lowest free energy was shown to be stabilized by a salt bridge between the end termini. A comparison between the free energies of the side chain analogues of the peptide sequence and the peptide itself shows that the free energy contributions are not additive. The simulations emphasize that tightly bound surface waters play a key role for peptide and protein structures when bound to inorganic surfaces in biological environments.
Multiscale modelling of nucleosome core particle aggregation
2015. Alexander P. Lyubartsev (et al.). Journal of Physics 27 (6)Artikel
The nucleosome core particle (NCP) is the basic building block of chromatin. Under the influence of multivalent cations, isolated mononucleosomes exhibit a rich phase behaviour forming various columnar phases with characteristic NCP-NCP stacking. NCP stacking is also a regular element of chromatin structure in vivo. Understanding the mechanism of nucleosome stacking and the conditions leading to self-assembly of NCPs is still incomplete. Due to the complexity of the system and the need to describe electrostatics properly by including the explicit mobile ions, novel modelling approaches based on coarse-grained (CG) methods at the multiscale level becomes a necessity. In this work we present a multiscale CG computer simulation approach to modelling interactions and self-assembly of solutions of NCPs induced by the presence of multivalent cations. Starting from continuum simulations including explicit three-valent cobalt(III) hexammine (CoHex(3+)) counterions and 20 NCPs, based on a previously developed advanced CG NCP model with one bead per amino acid and five beads per two DNA base pair unit (Fan et al 2013 PLoS One 8 e54228), we use the inverse Monte Carlo method to calculate effective interaction potentials for a 'super-CG' NCP model consisting of seven beads for each NCP. These interaction potentials are used in large-scale simulations of up to 5000 NCPs, modelling self-assembly induced by CoHex(3+). The systems of 'super-CG' NCPs form a single large cluster of stacked NCPs without long-range order in agreement with experimental data for NCPs precipitated by the three-valent polyamine, spermidine(3+).
Systematic Optimization of a Force Field for Classical Simulations of TiO2-Water Interfaces
2015. Erik G. Brandt, Alexander P. Lyubartsev. The Journal of Physical Chemistry C 119 (32), 18110-18125Artikel
Atomistic force field parameters were developed for the TiO2-water interface by systematic optimization with respect to experimentally determined crystal structures (lattice parameters) and surface thermodynamics (water adsorption enthalpy). Optimized force field parameters were determined for the two cases where TiO2 was modeled with or without covalent bonding. The nonbonded TiO2 model can be used to simulate different TiO2 phases, while the bonded TiO2 model is particularly useful for simulations of nanosized TiO2 and biomatter, including protein-surface and nanoparticle-biomembrane simulations. The procedure is easily generalized to parametrize interactions between other inorganic surfaces and biomolecules.
Systematic hierarchical coarse-graining with the inverse Monte Carlo method
2015. Alexander P. Lyubartsev (et al.). Journal of Chemical Physics 143 (24)Artikel
We outline our coarse-graining strategy for linking micro-and mesoscales of soft matter and biological systems. The method is based on effective pairwise interaction potentials obtained in detailed ab initio or classical atomistic Molecular Dynamics (MD) simulations, which can be used in simulations at less accurate level after scaling up the size. The effective potentials are obtained by applying the inverse Monte Carlo (IMC) method [A. P. Lyubartsev and A. Laaksonen, Phys. Rev. E 52(4), 3730-3737 (1995)] on a chosen subset of degrees of freedom described in terms of radial distribution functions. An in-house software package MagiC is developed to obtain the effective potentials for arbitrary molecular systems. In this work we compute effective potentials to model DNA-protein interactions (bacterial LiaR regulator bound to a 26 base pairs DNA fragment) at physiological salt concentration at a coarse-grained (CG) level. Normally the IMC CG pair-potentials are used directly as look-up tables but here we have fitted them to five Gaussians and a repulsive wall. Results show stable association between DNA and the model protein as well as similar position fluctuation profile.
A Coarse-Grained DNA Model Parameterized from Atomistic Simulations by Inverse Monte Carlo
2014. Nikolay Korolev (et al.). Polymers 6 (6), 1655-1675Artikel
Computer modeling of very large biomolecular systems, such as long DNA polyelectrolytes or protein-DNA complex-like chromatin cannot reach all-atom resolution in a foreseeable future and this necessitates the development of coarse-grained (CG) approximations. DNA is both highly charged and mechanically rigid semi-flexible polymer and adequate DNA modeling requires a correct description of both its structural stiffness and salt-dependent electrostatic forces. Here, we present a novel CG model of DNA that approximates the DNA polymer as a chain of 5-bead units. Each unit represents two DNA base pairs with one central bead for bases and pentose moieties and four others for phosphate groups. Charges, intra-and inter-molecular force field potentials for the CG DNA model were calculated using the inverse Monte Carlo method from all atom molecular dynamic (MD) simulations of 22 bp DNA oligonucleotides. The CG model was tested by performing dielectric continuum Langevin MD simulations of a 200 bp double helix DNA in solutions of monovalent salt with explicit ions. Excellent agreement with experimental data was obtained for the dependence of the DNA persistent length on salt concentration in the range 0.1-100 mM. The new CG DNA model is suitable for modeling various biomolecular systems with adequate description of electrostatic and mechanical properties.
A new AMBER-compatible force field parameter set for alkanes
2014. Alexei M. Nikitin, Yury V. Milchevskiy, Alexander P. Lyubartsev. Journal of Molecular Modeling 20 (3), 2143Artikel
We present a new force field parameter set for simulating alkanes. Its functional form and parameters are chosen to make it directly compatible with the AMBER94/99/12 family of force fields implemented in the available software. The proposed parameterization enables universal description of both the conformational and thermodynamic properties of linear, branched, and cyclic alkanes. Such unification is achieved by using two essential principles: (1) reduction of the Lennard-Jones radius for all sp3 carbons to 1.75 angstrom; (2) separate optimization of Lennard-Jones well depths for carbons with different degree of substitution. The new parameter set may prove to be optimal for description of alkyl residues in a broad range of biomolecules, from amino acids to lipids with their extended linear tails.
Bond orientation properties in lipid molecules of membranes
2014. Alexander L. Rabinovich, Alexander P. Lyubartsev. 25th IUPAP Conference on Computational Physics (CCP2013)Konferens
Atomistic molecular dynamics simulations have been carried out for 16 different fully hydrated phosphatidylcholine lipid bilayers, having 16 or 18 carbon atoms in fully saturated sn-1 chain and from 18 to 22 carbon atoms in sn-2 chain with different degree of unsaturation, with the purpose to investigate the effect of unsaturation on physical properties of lipid bilayers. Special attention has been paid to profiles of C-C and C-H bond order parameters of lipid molecules and the orientational fluctuations of these bond vectors. It was shown that the study of anisotropy degree of bond orientations probability distributions allows distinguishing extended regions with different types of angular fluctuations of bonds in a membrane formed by lipid molecules with unsaturated chains.
Molecular Dynamics Simulations Demonstrate the Regulation of DNA-DNA Attraction by H4 Histone Tail Acetylations and Mutations
2014. Nikolay Korolev (et al.). Biopolymers 101 (10), 1051-1064Artikel
The positively charged N-terminal histone tails play a crucial role in chromatin compaction and are important modulators of DNA transcription, recombination, and repair. The detailed mechanism of the interaction of histone tails with DNA remains elusive. To model the unspecific interaction of histone tails with DNA, all-atom molecular dynamics (MD) simulations were carried out for systems of four DNA 22-mers in the presence of 20 or 16 short fragments of the H4 histone tail (variations of the 16-23 a. a. KRHRKVLR sequence, as well as the unmodified fragment a. a. 13-20, GGAKRHRK). This setup with high DNA concentration, explicit presence of DNA-DNA contacts, presence of unstructured cationic peptides (histone tails) and K+ mimics the conditions of eukaryotic chromatin. A detailed account of the DNA interactions with the histone tail fragments, K+ and water is presented. Furthermore, DNA structure and dynamics and its interplay with the histone tail fragments binding are analysed. The charged side chains of the lysines and arginines play major roles in the tailmediated DNA-DNA attraction by forming bridges and by coordinating to the phosphate groups and to the electronegative sites in the minor groove. Binding of all species to DNA is dynamic. The structure of the unmodified fully-charged H4 16-23 a. a. fragment KRHRKVLR is dominated by a stretched conformation. The H4 tail a. a. fragment GGAKRHRK as well as the H4 Lys16 acetylated fragment are highly flexible. The present work allows capturing typical features of the histone tail-counterionDNA structure, interaction and dynamics.
Molecular Dynamics Studies of Liposomes as Carriers for Photosensitizing Drugs
2014. Joakim P. M. Jämbeck (et al.). Journal of Chemical Theory and Computation 10 (1), 5-13Artikel
Liposomes are proposed as drug delivery systems and can in principle be designed so as to cohere with specific tissue types or local environments. However, little detail is known about the exact mechanisms for drug delivery and the distributions of drug molecules inside the lipid carrier. In the current work, a coarse-grained (CG) liposome model is developed, consisting of over 2500 lipids, with varying degrees of drug loading. For the drug molecule, we chose hypericin, a natural compound proposed for use in photodynamic therapy, for which a CG model was derived and benchmarked against corresponding atomistic membrane bilayer model simulations. Liposomes with 21-84 hypericin molecules were generated and subjected to 10 microsecond simulations. Distribution of the hypericins, their orientations within the lipid bilayer, and the potential of mean force for transferring a hypericin molecule from the interior aqueous droplet through the liposome bilayer are reported herein.
Systematic Implicit Solvent Coarse Graining of Dimyristoylphosphatidylcholine Lipids
2014. Alexander Mirzoev, Alexander P. Lyubartsev. Journal of Computational Chemistry 35 (16), 1208-1218Artikel
We have used systematic structure-based coarse graining to derive effective site-site potentials for a 10-site coarse-grained dimyristoylphosphatidylcholine (DMPC) lipid model and investigated their state point dependence. The potentials provide for the coarse-grained model the same site-site radial distribution functions, bond and angle distributions as those computed in atomistic simulations carried out at four different lipid-water molar ratios. It was shown that there is a non-negligible dependence of the effective potentials on the concentration at which they were generated, which is also manifested in the properties of the lipid bilayers simulated using these potentials. Thus, effective potentials computed at low lipid concentration favor to more condensed and ordered structure of the bilayer with lower average area per lipid, while potentials obtained at higher lipid concentrations provide more fluid-like structure. The best agreement with the reference data and experiment was achieved using the set of potentials derived from atomistic simulations at 1:30 lipid:water molar ratio providing fully saturated hydration of DMPC lipids. Despite theoretical limitations of pairwise coarse-grained potentials expressed in their state point dependence, all the resulting potentials provide a stable bilayer structure with correct partitioning of different lipid groups across the bilayer as well as acceptable values of the average lipid area, compressibility and orientational ordering. In addition to bilayer simulations, the model has proven its robustness in modeling of self-aggregation of lipids from randomly dispersed solution to ordered bilayer structures, bicelles, and vesicles.
Update to the General Amber Force Field for Small Solutes with an Emphasis on Free Energies of Hydration
2014. Joakim P. M. Jämbeck, Alexander P. Lyubartsev. Journal of Physical Chemistry B 118 (14), 3793-3804Artikel
An approach to a straightforward reparametrization of the general AMBER force field (GAFF) for small organic solutes and druglike compounds is presented. The parametrization is based on specific pair interactions between the solvent and the solute, namely, the interactions between the constituting atoms of the solute and the oxygen of water were tuned in order to reproduce experimental hydration free energies for small model compounds. The key of the parametrization was to abandon the Lorentz-Berthelot combination rules for the van der Waals interactions. These parameters were then used for larger solutes in order to validate the newly derived pair interactions. In total close to 600 hydration free energies are computed, ranging from simple alkanes to multifunctional drug compounds, and compared to experimental data. The results show that the proposed parameters work well in describing the interactions between the solute and the solvent and that the agreement in absolute numbers is good. This modified version of GAFF is a good candidate for computing and predicting hydration free energies on a large scale, which has been a long-sought goal of computational chemists and can be used in rational drug design.
2D to 3D crossover of the magnetic properties in ordered arrays of iron oxide nanocrystals
2013. Bertrand Faure (et al.). Nanoscale 5 (3), 953-960Artikel
The magnetic 2D to 3D crossover behavior of well-ordered arrays of monodomain gamma-Fe2O3 spherical nanoparticles with different thicknesses has been investigated by magnetometry and Monte Carlo (MC) simulations. Using the structural information of the arrays obtained from grazing incidence small-angle X-ray scattering and scanning electron microscopy together with the experimentally determined values for the saturation magnetization and magnetic anisotropy of the nanoparticles, we show that MC simulations can reproduce the thickness-dependent magnetic behavior. The magnetic dipolar particle interactions induce a ferromagnetic coupling that increases in strength with decreasing thickness of the array. The 2D to 3D transition in the magnetic properties is mainly driven by a change in the orientation of the magnetic vortex states with increasing thickness, becoming more isotropic as the thickness of the array increases. Magnetic anisotropy prevents long-range ferromagnetic order from being established at low temperature and the nanoparticle magnetic moments instead freeze along directions defined by the distribution of easy magnetization directions.
An Advanced Coarse-Grained Nucleosome Core Particle Model for Computer Simulations of Nucleosome-Nucleosome Interactions under Varying Ionic Conditions
2013. Yanping Fan (et al.). PLoS ONE 8 (2), e54228Artikel
In the eukaryotic cell nucleus, DNA exists as chromatin, a compact but dynamic complex with histone proteins. The first level of DNA organization is the linear array of nucleosome core particles (NCPs). The NCP is a well-defined complex of 147 bp DNA with an octamer of histones. Interactions between NCPs are of paramount importance for higher levels of chromatin compaction. The polyelectrolyte nature of the NCP implies that nucleosome-nucleosome interactions must exhibit a great influence from both the ionic environment as well as the positively charged and highly flexible N-terminal histone tails, protruding out from the NCP. The large size of the system precludes a modelling analysis of chromatin at an all-atom level and calls for coarse-grained approximations. Here, a model of the NCP that include the globular histone core and the flexible histone tails described by one particle per each amino acid and taking into account their net charge is proposed. DNA wrapped around the histone core was approximated at the level of two base pairs represented by one bead (bases and sugar) plus four beads of charged phosphate groups. Computer simulations, using a Langevin thermostat, in a dielectric continuum with explicit monovalent (K+), divalent (Mg2+) or trivalent (Co(NH3)(6)(3+)) cations were performed for systems with one or ten NCPs. Increase of the counterion charge results in a switch from repulsive NCP-NCP interaction in the presence of K+, to partial aggregation with Mg2+ and to strong mutual attraction of all 10 NCPs in the presence of CoHex(3+). The new model reproduced experimental results and the structure of the NCP-NCP contacts is in agreement with available data. Cation screening, ion-ion correlations and tail bridging contribute to the NCP-NCP attraction and the new NCP model accounts for these interactions.
Another piece of the membrane puzzle
2013. Joakim P. M. Jämbeck, Alexander P. Lyubartsev. Journal of Chemical Theory and Computation 9 (1), 774-784Artikel
To be able to model complex biological membranes in a more realistic manner, the force field Slipids (Stockholm lipids) has been extended to include parameters for sphingomyelin (SM), phosphatidylglycerol (PG), phosphatidylserine (PS) lipids, and cholesterol. Since the parametrization scheme was faithful to the scheme used in previous editions of Slipids, all parameters are consistent and fully compatible. The results of careful validation of a number of key structural properties for one and two component lipid bilayers are in excellent agreement with experiments. Potentials of mean force for transferring water across binary mixtures of lipids and cholesterol were also computed in order to compare water permeability rates to experiments. In agreement with experimental and simulation studies, it was found that the permeability and partitioning of water is affected by cholesterol in lipid bilayers made of saturated lipids to the largest extent. With the extensions of Slipids presented here, it is now possible to study complex systems containing many different lipids and proteins in a fully atomistic resolution in the isothermic-isobaric (NPT) ensemble, which is the proper ensemble for membrane simulations.
Computer Simulation of Lipid Membranes
2013. A. L. Rabinovich, Alexander P. Lyubartsev. POLYMER SCIENCE SERIES C 55 (1), 162-180Artikel
Rapid development of computer power during the last decade has made molecular simulations of lipid bilayers feasible for many research groups, which, together with the growing general interest in investigations of these very important biological systems has lead to tremendous increase of the number of research on the computational modeling of lipid bilayers. In this review, we give account of the recent progress in computer simulations of lipid bilayers covering mainly the period of the last 7 years, and covering only several selected subjects: methodological (development of the force fields for lipid bilayer simulations, use of coarse-grained models) and scientific (studies of the role of lipid unsaturation, and the effect of cholesterol and other inclusions on properties of the bilayer).
Exploring the Free Energy Landscape of Solutes Embedded in Lipid Bilayers
2013. Joakim P. M. Jämbeck, Alexander P. Lyubartsev. Journal of Physical Chemistry Letters 4 (11), 1781-1787Artikel
Free energy calculations are vital for our understanding of biological processes on an atomistic scale and can offer insight to various mechanisms. However, in some cases, degrees of freedom (DOFs) orthogonal to the reaction coordinate have high energy barriers and/or long equilibration times, which prohibit proper sampling. Here we identify these orthogonal DOFs when studying the transfer of a solute from water to a model membrane. Important DOFs are identified in bulk liquids of different dielectric nature with metadynamics simulations and are used as reaction coordinates for the translocation process, resulting in two- and three-dimensional space of reaction coordinates. The results are in good agreement with experiments and elucidate the pitfalls of using one-dimensional reaction coordinates. The calculations performed here offer the most detailed free energy landscape of solutes embedded in lipid bilayers to date and show that free energy calculations can be used to study complex membrane translocation phenomena.
Implicit inclusion of atomic polarization in modeling of partitioning between water and lipid bilayers
2013. Joakim P. M. Jämbeck, Alexander P. Lyubartsev. Physical Chemistry, Chemical Physics - PCCP 15 (13), 4677-4686Artikel
We propose an effective and straightforward way of including atomic polarization in simulations of the partitioning of small molecules in inhomogenous media based on classical molecular dynamics with non-polarizable force fields. The approach presented here takes advantage of the relatively fast sampling of phase space obtained with additive force fields by adding the polarization effects afterwards. By using pre-polarized charges for the polar and non-polar phases together with a polarization correction term the effects of atomic polarization are effectively taken into account. The results show a clear improvement compared to using the more common setup with one set of charges obtained from gas phase ab initio calculations. It is shown that when proper measures are taken into account computer simulations with non-polarizable force fields are able to accurately determine water-membrane partitioning and preferential location of small molecules in the membrane interior. We believe that the approach presented here can be useful in rational drug design and in investigations of molecular mechanisms of anesthetic or toxic action.
2013. Alexander Mirzoev, Alexander P. Lyubartsev. Journal of Chemical Theory and Computation 9 (3), 1512-1520Artikel
We present software package MagiC, which is designed to perform systematic structure-based coarse graining of molecular models. The effective pairwise potentials between coarse-grained sites of low-resolution molecular models are constructed to reproduce structural distribution functions obtained from the modeling of the system in a high resolution (atomistic) description. The software supports coarse-grained tabulated intramolecular bond and angle interactions, as well as tabulated nonbonded interactions between different site types in the coarse-grained system, with the treatment of long-range electrostatic forces by the Ewald summation. Two methods of effective potential refinement are implemented: iterative Boltzmann inversion and inverse Monte Carlo, the latter accounting for cross-correlations between pair interactions. MagiC uses its own Metropolis Monte Carlo sampling engine, allowing parallel simulation of many copies of the system with subsequent averaging of the properties, which provides fast convergence of the method with nearly linear scaling at parallel execution.
Partial Atomic Charges and Their Impact on the Free Energy of Solvation
2013. Joakim P. M. Jämbeck (et al.). Journal of Computational Chemistry 34 (3), 187-197Artikel
Free energies of solvation (Delta G) in water and n-octanol have been computed for common drug molecules by molecular dynamics simulations with an additive fixed-charge force field. The impact of the electrostatic interactions was investigated by computing the partial atomic charges with four methods that all fit the charges from the quantum mechanically determined electrostatic potential (ESP). Due to the redistribution of electron density that occurs when molecules are transferred from gas phase to condensed phase, the polarization impact was also investigated. By computing the partial atomic charges with the solutes placed in a conductor-like continuum, the charges were effectively polarized to take the polarization effects into account. No polarization correction term or similar was considered, only the partial atomic charges. Results show that free energies are very sensitive to the choice of atomic charges and that Delta G can differ by several k(B)T depending on the charge computing method. Inclusion of polarization effects makes the solutes too hydrophilic with most methods and in vacuo charges make the solutes too hydrophobic. The restrained-ESP methods together with effectively polarized charges perform well in our test set and also when applied to a larger set of molecules. The effect of water models is also highlighted and shows that the conclusions drawn are valid for different three-point models. Partitioning between an aqueous and a hydrophobic phase is also described better if the two environment's polarization is taken into account, but again the results are sensitive to the charge calculation method. Overall, the results presented here show that effectively polarized charges can improve the description of solvating a drug-like molecule in a solvent and that the choice of partial atomic charges is crucial to ensure that molecular simulations produce reliable results.
An extension and further validation of an all atomistic force field for biological membranes
2012. Joakim P. M. Jämbeck, Alexander P. Lyubartsev. Journal of Chemical Theory and Computation 8 (8), 2938-2948Artikel
Biological membranes are versatile in composition and host intriguing molecular processes. In order to be able to study these systems, an accurate model Hamiltonian or force field (FF) is a necessity. Here, we report the results of our extension of earlier developed all-atomistic FF parameters for fully saturated phospholipids that complements an earlier parameter set for saturated phosphatidylcholine lipids (J. Phys. Chem. B, 2012, 116, 3164-3179). The FF, coined Slipids (Stockholm lipids), now also includes parameters for unsaturated phosphatidylcholine and phosphatidylethanolamine lipids, e.g., POPC, DOPC, SOPC, POPE, and DOPE. As the extended set of parameters is derived with the same philosophy as previously applied, the resulting FF has been developed in a fully consistent manner. The capabilities of Slipids are demonstrated by performing long simulations without applying any surface tension and using the correct isothermal-isobaric (NPT) ensemble for a range of temperatures and carefully comparing a number of properties with experimental findings. Results show that several structural properties are very well reproduced, such as scattering form factors, NMR order parameters, thicknesses, and area per lipid. Thermal dependencies of different thicknesses and area per lipid are reproduced as well Lipid diffusion is systematically slightly underestimated, whereas the normalized lipid diffusion follows the experimental trends. This is believed to be due to the lack of collective movement in the relatively small bilayer patches used Furthermore, the compatibility with amino acid FFs from the AMBER family is tested in explicit transmembrane complexes of the WALP23 peptide with DLPC and DOPC bilayers, and this shows that Slipids can be used to study more complex and biologically relevant systems.
Derivation and Systematic Validation of a Refined All-Atom Force Field for Phosphatidylcholine Lipids
2012. Joakim P. M. Jämbeck, Alexander P. Lyubartsev. Journal of Physical Chemistry B 116 (10), 3164-3179Artikel
An all-atomistic force field (FF) has been developed for fully saturated phospholipids. The parametrization has been largely based on high-level ab initio calculations in order to keep the empirical input to a minimum. Parameters for the lipid chains have been developed based on knowledge about bulk alkane liquids, for which thermodynamic and dynamic data are excellently reproduced. The FFs ability to simulate lipid bilayers in the liquid crystalline phase in a tensionless ensemble was tested in simulations of three lipids: 1,2-diauroyl-sn-glycero-3-phospocholine (DLPC), 1,2-dimyristoyl-sn-glycero-3-phosphocholine (DMPC), and 1,2-dipalmitoyl-sn-glycero-3-phospcholine (DPPC). Computed areas and volumes per lipid, and three different kinds of bilayer thicknesses, have been investigated. Most importantly NMR order parameters and scattering form factors agree in an excellent manner with experimental data under a range of temperatures. Further, the compatibility with the AMBER FF for biomolecules as well as the ability to simulate bilayers in gel phase was demonstrated. Overall, the FF presented here provides the important balance between the hydrophilic and hydrophobic forces present in lipid bilayers and therefore can be used for more complicated studies of realistic biological membranes with protein insertions.
M.DynaMix Studies of Solvation, Solubility and Permeability
2012. Aatto Laaksonen, Alexander Lyubartsev, Francesca Mocci. Molecular Dynamics - Studies of Synthetic and Biological Macromolecules, 85-106Kapitel
Modelling chromatin structure and dynamics
2012. Nikolay Korolev (et al.). Current opinion in structural biology 22 (2), 151-159Artikel
The packaging of genomic DNA into chromatin in the eukaryotic cell nucleus demands extensive compaction. This requires attractive nucleosome-nucleosome interactions to overcome repulsion between the negatively charged DNA segments as well as other constraints. At the same time, DNA must be dynamically accessible to the cellular machinery that operates on it. Recent progress in the experimental characterisation of the higher order structure and dynamics of well-defined chromatin fibres has stimulated the attempts at theoretical description of chromatin and the nucleosome. Here we review the present status of chromatin Modelling, with particular emphasis on coarse-grained computer simulation models, the role of electrostatic interactions, and discuss future perspectives in the field.
The polyelectrolyte properties of chromatin
2012. Nikolay Korolev (et al.). Soft Matter 8 (36), 9322-9333Artikel
Double helical DNA is a negatively charged polyelectrolyte and exists in the nucleus of living cells as chromatin, a highly compacted but dynamic complex with histone proteins. The first level of DNA compaction is the linear array of the nucleosome core particles (NCP), which is a well-defined structure of 145-147 bp DNA with the histone octamer, connected by linker DNA. Higher levels of chromatin compaction include two routes which may overlap: intramolecular folding of the nucleosome array resulting in formation of the 30 nm fibre and intermolecular aggregation (self-association) between different arrays (or distant fibres of the same chromosome). This review describes how the polyelectrolyte properties of chromatin are illustrated by experimental results of folding and self-association of well-defined model chromatin, in the form of recombinant nucleosome arrays, and how these properties can be understood from computer modelling. Chromatin compaction shows considerable similarities to DNA condensation. However, the structure of condensed chromatin is sensitive to the detailed molecular features of the nucleosome-nucleosome interactions which include the influence of the histone tails and their modifications.
Calculation of Canonical Properties and Excited States by Path Integral Numerical Methods
2011. P. N. Vorontsov-Velyaminov (et al.). Contributions to Plasma Physics 51 (4), 382-385Artikel
Path integral Monte Carlo method in the expanded ensemble is used for calculation of the ratio of partition functions for different classes of permutations treating the problem of several interacting identical particles (fermions) in an external field. Wang-Landau algorithm is used for adjustment of balancing factors. For systems consisting of greater than two number of particles we propose an advanced variant of our approach which implies calculation of the ratio of positive and negative contributions to the partition function. Densities and energies of the sequence of excited states starting from the ground state for a system of non interacting quantum particles are calculated in turn, one by one, by means of considering systems with artificially excluded lowest energy levels and further obtaining of the ground state of each next system constructed in this way. The idea of evaluation of densities of excited states for a quantum system with interparticle interaction by evaluating density difference between systems of different number of noninteracting Fermi-copies of the system of interest is realized in terms of cyclic expansion formalism for a simple 2D system of two spinless fermions interacting via Coulomb repulsion.
Effective solvent mediated potentials of Na+ and Cl− ions in aqueous solution
2011. Alexander Mirzoev, Alexander P. Lyubartsev. Physical Chemistry, Chemical Physics - PCCP (13), 5722-5727Artikel
The effective solvent-mediated potentials for Na+ and Cl− ions in aqueous solution were calculated in a wide range of temperatures from 0 to 100 °C. The potentials have been determined using the inverse Monte Carlo approach, from the ion–ion radial distribution functions computed in 50 ns molecular dynamics simulations of ions and explicit water molecules. We further separated the effective potentials into a short-range part and an electrostatic long-range part represented by a coulombic potential with some dielectric permittivity. We adjusted the value of the dielectric permittivity to provide the fastest possible decay of the short-range potentials at larger distances. The obtained temperature dependence of the dielectric permittivity follows well the experimental data. We show also that the largest part of the temperature dependence of the effective potentials can be attributed to the temperature-dependent dielectric permittivity.
Molecular Dynamics Simulation Study of Glycerol-Water Liquid Mixtures
2011. Andrei V. Egorov, Alexander P. Lyubartsev, Aatto Laaksonen. Journal of Physical Chemistry B 115 (49), 14572-14581Artikel
To study the effects of water on conformational dynamics of polyalcohols, Molecular Dynamics simulations of glycerol water liquid mixtures have been carried out at different concentrations: 42.9 and 60.0 wt 96 of glycerol, respectively. On the basis of the analysis of backbone conformer distributions, it is found that the surrounding water molecules have a large impact on the populations of the glycerol conformers. While the local structure of water in the liquid mixture is surprisingly close to that in pure liquid water, the behavior of glycerols can be divided into three different categories where roughly 25% of them occur in a structure similar to that in pure liquid of glycerol, ca. 25% of them exist as monomers, solvated by water, and the remaining 50% of glycerols in the mixture form H-bonded strings as. remains of the glycerol H-bond network. The typical glycerol H-bond network still exists even at the lower concentration of 40 wt % of glycerol. The microheterogeneity of water glycerol mixtures is analyzed using time-averaged distributions of the sizes of the water aggregates. At 40 wt % of glycerol, the cluster sizes from 3 to 10 water molecules are observed. The increase of glycerol content causes a depletion of clusters leading to smaller 3-5 molecule clusters domination. Translational diffusion coefficients have been calculated to study the dynamical behavior of both glycerol and water molecules. Rotational-reorientational motion is studied both in overall and in selected substructures on the basis of time correlation functions. Characteristic time scales for different motional modes are deduced on the basis of the calculated correlation times. The general conclusion is that the presence of water increases the overall mobility of glycerol, while glycerol slows the mobility of water.
Recent development in computer simulations of lipid bilayers
2011. Alexander P. Lyubartsev, Alexander L. Rabinovich. Soft Matter 7 (1), 25-39Artikel
Rapid development of computer power during the last decade has made molecular simulations of lipid bilayers feasible for many research groups, which, together with the growing general interest in investigations of these very important biological systems has lead to tremendous increase of the number of research on the computational modeling of lipid bilayers. In this review, we give account of the recent progress in computer simulations of lipid bilayers covering mainly the period of the last 5 years, and covering several selected subjects: development of the force fields for lipid bilayer simulations, studies of the role of lipid unsaturation, the effect of cholesterol and other inclusions on properties of the bilayer, and use of coarse-grained models.
Solute-Solvent Interactions in Aqueous Glycylglycine-CuCl(2) Solutions
2011. Mysore Sridhar Santosh (et al.). Journal of Solution Chemistry 40 (9), 1657-1671Artikel
Acoustical and molecular dynamics studies were carried out to understand the various interactions present in glycylglycine-CuCl(2) aqueous solutions. Amongst these interactions, hydrogen bonding and solute-solvent interactions have been highlighted in this study. The radial distribution function (RDF) was used to investigate solution structure and hydration parameters. Binding of Cu(2+) with various polar peptide atoms reveals the nature and degree of binding. The formation of complex clusters between glycylglycine and water molecules increases the relaxation time. The first hydration shell considerably influences the structure of the second shell, facilitating the formation of an ordered hydrogen bonded network. Both experimental and theoretical results have proved to be efficient in analyzing the behavior of molecules and to give a clear idea on molecular interactions in solutions.
Two-Dimensional Wang–Landau Algorithm for Osmotic Pressure Calculations in a Polyelectrolyte–Membrane System
2011. Nikolai Volkov, Pavel Vorontsov-Velyaminov, Alexander Lyubartsev. Macromolecular Theory and Simulations 20 (7), 496-509Artikel
The Monte Carlo method based on two-dimensional entropic sampling within the Wang–Landau (WL) algorithm is applied to simulation of a continuous model of a polyelectrolyte between membrane surfaces. Membranes are presented by parallel plane surfaces holding either fixed or mobile dipoles (representing lipid headgroups). A strongly charged polyion accompanied by neutralizing counterions is placed between the membranes. Periodic boundary conditions are imposed along X- and Y-axes. The volume of the main cell is varied during the simulation by shifting one of the surfaces along Z-axis. Within two-dimensional WL sampling algorithm we obtain joint density of states as a function of energy and volume in a single run. In order to increase efficiency of our calculations we introduce a number of modifications to the original WL-approach. Various properties of the system over wide temperature and volume or pressure ranges, i.e., conformational energy, heat capacity, and free energy, are obtained from the two-dimensional density of states by simple integration. The osmotic pressure is calculated as a derivative of Helmholtz free energy. Alternatively, properties of the system, including average volume, can be obtained under condition of NPT ensemble. It is shown that the both approaches produce coinciding Posm(V) isotherms. In all considered cases we observe repulsive effective interaction between the membrane surfaces, and repulsion is stronger for the surfaces with fixed dipoles.
Cation-indiced polyelectrolyte-polyelectrolyte attraction in solutions of DNA and nucleosome core particles
2010. Nikolai Korolev, Alexander Lyubartsev, Lars Nordenskiöld. Journal of Colloid and Interface Science 158, 32-47Artikel
The paper reviews our current studies on the experimentally induced cation compaction and aggregation in solutions of DNA and nucleosome core particles and the theoretical modelling of these processes using coarse-grained continuum models with explicit mobile ions and with all-atom molecular dynamics (MD) simulations. Recent experimental results on DNA condensation by cationic oligopeptides and the effects of added salt are presented. The results of MD simulations modelling DNA–DNA attraction due to the presence of multivalent ions including the polyamine spermidine and fragments of histone tails, which exhibit bridging between adjacent DNA molecules, are discussed. Experimental data on NCP aggregation, using recombinantly prepared systems are summarized. Literature data and our results of studying of the NCP solutions are compared with predictions of coarse-grained MD simulations, including the important ion correlation as well as bridging mechanisms. The importance of the results to chromatin folding and aggregation is discussed.
Centroid molecular dynamics
2010. Eugeny A. Polyakov, Alexander P. Lyubartsev, Pavel N. Vorontsov-Velyaminov. Journal of Chemical Physics 133 (19), 194103Artikel
The relation between the accuracy of centroid molecular dynamics correlation functions, and the geometry of the centroid potential is investigated. It is shown that, depending on the temperature, there exist several regimes, and in each of them certain features of the exact Kubo correlation functions are reproduced. The change of regimes is related to the emergence of barriers in the centroid potential. In order to clarify how the above described picture of regimes is modified in real systems when dissipation is important, a methodology is developed to test the accuracy of centroid correlation functions for the model of a particle coupled to a harmonic heat bath. A modification of the centroid molecular dynamics method to include the influence of the heat bath is introduced. Preliminary results of comparison of centroid molecular dynamics with the numerically exact results of filtered propagator functional method are presented.
Electrostatic Origin of Salt-induced Nucleosome Array Compaction
2010. Nikolai Korolev (et al.). Biophysical Journal 99, 1896-1905Artikel
The physical mechanism of the folding and unfolding of chromatin is fundamentally related to transcription but is incompletely characterized and not fully understood. We experimentally and theoretically studied chromatin compaction by investigating the salt-mediated folding of an array made of 12 positioning nucleosomes with 177 bp repeat length. Sedimentation velocity measurements were performed to monitor the folding provoked by addition of cations Na+, K+, Mg2+, Ca2+, spermidine3+, Co(NH3)63+, and spermine4+. We found typical polyelectrolyte behavior, with the critical concentration of cation needed to bring about maximal folding covering a range of almost five orders of magnitude (from 2 μM for spermine4+ to 100 mM for Na+). A coarse-grained model of the nucleosome array based on a continuum dielectric description and including the explicit presence of mobile ions and charged flexible histone tails was used in computer simulations to investigate the cation-mediated compaction. The results of the simulations with explicit ions are in general agreement with the experimental data, whereas simple Debye-Hückel models are intrinsically incapable of describing chromatin array folding by multivalent cations. We conclude that the theoretical description of the salt-induced chromatin folding must incorporate explicit mobile ions that include ion correlation and ion competition effects.
Molecular Dynamics Investigation of Dipeptide - Transition Metal Salts in Aqueous Solutions
2010. Mysore S. Santosh (et al.). Journal of Physical Chemistry B 114 (49), 16632-16640Artikel
Molecular dynamics (MD) simulations of glycylglycine dipeptide with transition metal ions (Mn2+, Fe2+, Co2+, Ni2+, Cu2+, and Zn2+) in aqueous solutions have been carried out to get an insight into the solvation structure, intermolecular interactions, and salt effects in these systems. The solvation structure and hydrogen bonding were described in terms of radial distribution function (RDF) and spatial distribution function (SDF). The dynamical properties of the solvation structure were also analyzed in terms of diffusion and residence times. The simulation results show the presence of a well-defined first hydration shell around the dipeptide, with water molecules forming hydrogen bonds to the polar groups of the dipeptide. This shell is, however, affected by the strong electric field of divalent metal ions, which at higher ion concentrations lead to the shift in the dipeptide−water RDFs. Higher salt concentrations lead also to increased residence times and slower diffusion rates. In general, smaller ions (Cu2+, Zn2+) demonstrate stronger binding to dipeptide than the larger ones (Fe2+, Mn2+). Simulations do not show any stronger association of peptide molecules indicating their dissolution in water. The above results may be of potential interest to future researchers on these molecular interactions.
Molecular Dynamics Simulations of Local Anesthetic Articaine in a Lipid Bilayer
2010. Enamul Mojumdar, Alexander Lyubartsev. Biophysical Chemistry 153 (1), 27-35Artikel
In order to investigate structural and dynamical propertiesof local anesthetic articaine in a model lipid bilayer, a series ofmolecular dynamics simulations have been performed. Simulations werecarried out for neutral and charged (protonated) forms of articaine insertedin fully hydrated dimyristoylphosphatidylcholine (DMPC) lipid bilayer. For comparison purpose, a fully hydrated DMPC bilayer withoutarticaine was also simulated. The length of each simulation was 200 ns. Various properties of the lipid bilayer systems in the presence of both charged and uncharged forms of articaine taken at two different concentrations have been examined: membrane area per lipid, mass density distributions, order parameters, radial distribution functions, head group tilt, diffusion coefficients, electrostatic potential, etc, and compared with results of previous simulations of DMPC bilayer in presence of lidocaine.It was shown that addition of both charged and neutral forms of articainecauses increase of the dipole electrostatic potential in the membraneinterior.
Systematic coarse-graining of molecular models by the Newton inversion method
2010. Alexander Lyubartsev (et al.). Faraday discussions (Online) 144, 43-56Artikel
Systematic construction of coarse-grained molecular models from detailed atomistic simulations, and even from ab initio simulations is discussed. Atomistic simulations are first performed to extract structural information about the system, which is then used to determine effective potentials for a coarse-grained model of the same system. The statistical-mechanical equations expressing the canonical properties in terms of potential parameters can be inverted and solved numerically according to the iterative Newton scheme. In our previous applications, known as the Inverse Monte Carlo, radial distribution functions were inverted to reconstruct pair potential, while in a more general approach the targets can be other canonical averages. We have considered several examples of coarse-graining; for the united atom water model we suggest an easy way to overcome the known problem of high pressure. Further, we have developed coarse-grained models for L- and D-prolines, dissolved here in an organic solvent (dimethylsulfoxide), keeping their enantiomeric properties from the corresponding all-atom proline model. Finally, we have revisited the previously developed coarse-grained lipid model based on an updated all-atomic force field. We use this model in large-scale meso-scale simulations demonstrating spontaneous formation of different structures, such as vesicles, micelles, and multi-lamellar structures, depending on thermodynamical conditions.
Theoretical Approximations to X-ray Absorption Spectroscopy of Liquid Water and Ice
2010. Mikael Leetmaa (et al.). Journal of Electron Spectroscopy and Related Phenomena 177 (2-3), 135-157Artikel
We review methods to compute x-ray absorption spectra (XAS) with special focus on the transition potential approach of Triguero et al. [Phys. Rev. B 58, 8097 (1998)] and its application to calculations on water in condensed phase. We discuss the absolute energy scale, functional dependence, broadening versus sampling of intra- and intermolecular vibrational modes, treatment of the continuum, cluster size convergence as well as compare with periodic calculations and with experiment; periodic and cluster model calculations are found to agree very closely in the relevant near-edge region although neither reproduces the pre-edge and main-edge features in the experimental spectra of thin ice films. The real space grid representation of the wave function in the periodic calculations allows a more extended energy range to be described and we find satisfactory agreement with experiment for higher energy continuum resonances. Two proposed alternative approaches using either the potential from a full core-hole (FCH) or the full core-hole with an excited electron in the lowest state (XCH) are shown to lead to spectra that deviate significantly from experiment.
Assessing the electric-field approximation to IR and Raman spectra of dilute HOD in D2O
2009. Mathias Ljungberg (et al.). Journal of Chemical Physics 131, 034501Artikel
We analyze the validity of the commonly used electric-field (E-field) approximation to vibrational OH stretch Raman spectra of dilute HOD in D2O by computing the OH stretch frequency of all molecules in several different structure models, each containing around 2000 molecules. The calculations are done at the B3LYP level using clusters containing 32 molecules centered around the molecule for which the frequencies are calculated; the large cluster size is required due to significant nonlocal contributions influencing the computed frequencies. The vibrational frequencies are determined using a six-point potential optimized discrete variable representation. Raman and infrared intensities are furthermore computed to generate the spectra. We find that a quadratic fit of E-field versus frequency gives a reasonable representation of the calculated distribution of frequencies. However, the mapping depends significantly on the structural model and is thus not universal. Anharmonic couplings are calculated for several optimized clusters showing a general trend to compress the computed frequency distributions, which is in agreement with dynamical simulations (motional narrowing).
Computer Modeling Reveals that Modifications of the Histone Tail Charges Define Salt-Dependent Interaction of the Nucleosome Core Particles
2009. Ye Yang (et al.). Biophysical Journal 96, 2082-2096Artikel
Coarse-grained Langevin molecular dynamics computer simulations were conducted for systems that mimic solutions of nucleosome core particles (NCPs). The NCP was modeled as a negatively charged spherical particle representing the complex of DNA and the globular part of the histones combined with attached strings of connected charged beads modeling the histone tails. The size, charge, and distribution of the tails relative to the core were built to match real NCPs. Three models of NCPs were constructed to represent different extents of covalent modification on the histone tails: (nonmodified) recombinant (rNCP), acetylated (aNCP), and acetylated and phosphorylated (paNCP). The simulation cell contained 10 NCPs in a dielectric continuum with explicit mobile counterions and added salt. The NCP-NCP interaction is decisively dependent on the modification state of the histone tails and on salt conditions. Increasing the monovalent salt concentration (KCl) from salt-free to physiological concentration leads to NCP aggregation in solution for rNCP, whereas NCP associates are observed only occasionally in the system of aNCPs. In the presence of divalent salt (Mg2+), rNCPs form dense stable aggregates, whereas aNCPs form aggregates less frequently. Aggregates are formed via histone-tail bridging and accumulation of counterions in the regions of NCP-NCP contacts. The paNCPs do not show NCP-NCP interaction upon addition of KCl or in the presence of Mg2+. Simulations for systems with a gradual substitution of K+ for Mg2+, to mimic the Mg2+ titration of an NCP solution, were performed. The rNCP system showed stronger aggregation that occurred at lower concentrations of added Mg2+, compared to the aNCP system. Additional molecular dynamics simulations performed with a single NCP in the simulation cell showed that detachment of the tails from the NCP core was modest under a wide range of salt concentrations. This implies that salt-induced tail dissociation of the histone tails from the globular NCP is not in itself a major factor in NCP-NCP aggregation. The approximation of coarse-graining, with respect to the description of the NCP as a sphere with uniform charge distribution, was tested in control simulations. A more detailed description of the NCP did not change the main features of the results. Overall, the results of this work are in agreement with experimental data reported for NCP solutions and for chromatin arrays.
Conformational characteristics of single flexible polyelectrolyte chain
2009. C.G. Jesudason, Alexander Lyubartsev, Aatto Laaksonen. The European Physical Journal E Soft matter 30, 341-350Artikel
The behavior of a flexible anionic chain of 150 univalent and negatively charged beads connected by harmonic-like potential interactions with each other in the presence of equal number of positive and free counterions, observed in molecular dynamics simulations with Langevin thermostat, is described in a temperature range from 0.1 to 10.0 in reduced units. The total and Coulombic energies, radial distribution functions, radii of gyration, end-to-end distances of the chain are depicted. Our results turned out to be qualitatively similar to the results obtained earlier for a lattice polyelectrolyte model, including temperature maximum of the polyelectrolyte chain and internal phase transition which seems to occur abruptly at low temperatures for all the systems studied, judging from the shape of end-to end distance, gyration radius and energy profiles.
Hierarchical Multiscale Modelling Scheme from First Principles to Mesoscale
2009. Alexander Lyubartsev, Yaoquan Tu, Aatto Laaksonen. Journal of computational and theoretical nanoscience 6 (5), 951-959Artikel
We present a straight-forward implementation of a practical hierarchical multiscale modelling scheme which enables us to start from first-principles atomistic computer simulation and successively coarse-grain the model by leaving out uninteresting degrees of freedom. Using the Car-Parrinello method or our recently developed highly efficient tight-binding-like approximate density-functional quantum mechanical method, we first perform ab initio simulations. From these first-principles simulations we obtain a set of atomistic pair-wise effective interaction potentials to be used as a force field with no empirical data for subsequent classical all-atom simulations while scaling up the system size 2-3 orders of magnitude. The atomistic simulations similarly provide a new set of effective potentials but at a chosen coarse-grain level suitable for large-scale mesoscopic or soft-matter simulations beyond the atomic resolution. Several examples are shown of how this scheme is done based on effective interaction potentials to tie together the various scales of modelling.
One dimensional model for water and aqueous solutions. Part V. Monte Carlo simulation of dilute solutions of hard rod in waterlike particles
2009. Alexander Lyubartsev, Arieh Ben-Naim. Journal of Chemical Physics 131, 204507Artikel
Wehave carried out Monte Carlo simulation on the primitive onedimensional model for water described earlier [A. Ben-Naim, J. Chem.Phys. 128, 024506 (2008)]. We show that by taking intoaccount second nearest neighbor interactions, one can obtain the characteristicanomalous solvation thermodynamic quantities of inert solutes in water. Thismodel clearly demonstrates the molecular origin of the large negativeentropy of solvation of an inert solute in water.
Path-integral-expanded-ensemble Monte Carlo method in treatment of the sign problem for fermions
2009. M. A. Voznesenskiy, P. N. Vorontsov-Velyaminov, Alexander Lyubartsev. PHYS REV E 80 (6), 66702Artikel
Expanded-ensemble Monte Carlo method with Wang-Landau algorithm was used for calculations of the ratio of partition functions for classes of permutations in the problem of several interacting quantum particles (fermions) in an external field. Simulations for systems consisting of 2 up to 7 interacting particles in harmonic or Coulombic field were performed. The presented approach allows one to carry out calculations for low enough temperatures that makes it possible to extract data for the ground-state energy and low-temperature thermodynamics.
Stochastic positive P-representation in problems of quantum statistics. Simulation of one-dimensional Bose-gas with delta-repulsion
2009. Eugeniy Polyakov, Pavel Vorontsov-Velyaminov, Alexander Lyubartsev. Vychislitel'nye Metody i Programirovaniye [Numerical methods and programming] 10, 223-247Artikel
A method of stochastic positive P-representation for the computer simulation of thermal equilibrium and dynamical properties of many-particle quantum systems with interactions is proposed and thoroughly analyzed. The testing procedure of the method includes the evaluation of spatial correlation functions for the one-dimensional Bose-gas with delta-repulsion between particles, both in the state of thermal equilibrium and in the dynamical evolution from a given initial state. This work was supported by the Russian Foundation for Basic Research (project number 08–02–00041) and by the Royal Swedish Academy of Sciences.
Structural Evidence for the Ordered Crystallites of Ionic Liquid in Confined Carbon Nanotubes
2009. Kun Dong (et al.). The Journal of Physical Chemistry C 113 (23), 10013-10020Artikel
Ionic liquids (ILs) are a class of new green materials that have attracted extensive attention in recent decades. Many novel properties not evident under normal conditions may appear when ionic liquids are confined to a nanometer scale. As was observed in the experiment, an anomalous phase behavior from liquid to high melting point perfect crystal occurred when 1-n-butyl-3-methylimidazolium hexafluorophosphate ([bmim][PF6]) ionic liquid was confined in a carbon nanotube. In this work, we performed molecular dynamics (MD) simulations for [bmim][PF6] ionic liquid and provided direct structural evidence that the ionic crystallizes in a carbon nanotube. The ordered ionic arrangement in both the radial and the axial directions can be observed inside the channels of the CNTs to induce the form of crystallites. The ionic stacking and distributing can be determined by the sizes of the CNTs. Hydrogen bonds remain the dominant interactions between cations and anions when the ionic liquid enters into the CNT from the bulk phase. The free energies as the thermal driven forces were calculated, and it is found that it is very difficult for a single anion to enter into the channel of the CNT spontaneously. A more favorable way is through an ion-pair in which a cation “pulls” an anion to enter into the channel of the CNT together. It is predicted that other ionic liquids that possess similar structures, even including the pyridinium-based ionic liquids, can show higher melting points when confined in CNTs.
The Inhomogeneous Structure of Water at Ambient Conditions
2009. Congcong Huang (et al.). Proceedings of the National Academy of Sciences of the United States of America 106, 15214-15218Artikel
Small-angle X-ray scattering (SAXS) is used to demonstrate the presence of density fluctuations in ambient water on a physical length-scale of ≈1 nm; this is retained with decreasing temperature while the magnitude is enhanced. In contrast, the magnitude of fluctuations in a normal liquid, such as CCl4, exhibits no enhancement with decreasing temperature, as is also the case for water from molecular dynamics simulations under ambient conditions. Based on X-ray emission spectroscopy and X-ray Raman scattering data we propose that the density difference contrast in SAXS is due to fluctuations between tetrahedral-like and hydrogen-bond distorted structures related to, respectively, low and high density water. We combine our experimental observations to propose a model of water as a temperature-dependent, fluctuating equilibrium between the two types of local structures driven by incommensurate requirements for minimizing enthalpy (strong near-tetrahedral hydrogen-bonds) and maximizing entropy (nondirectional H-bonds and disorder). The present results provide experimental evidence that the extreme differences anticipated in the hydrogen-bonding environment in the deeply supercooled regime surprisingly remain in bulk water even at conditions ranging from ambient up to close to the boiling point.
2008. Lars Nordenskiöld, Nikolai Korolev, Alexander Lyubartsev. DNA interactions with polymers and surfactants, 209-238Kapitel
Effect of Local Anesthetic Lidocaine on Electrostatic Properties of a Lipid Bilayer
2008. Carl-Johan Högberg, Alexander Lyubartsev. Biophysical Journal 94, 525-531Artikel
The influence of local anesthetic lidocaine on electrostatic properties of a lipid membrane bilayer was studied by molecular dynamics simulations. The electrostatic dipole potential, charge densities, and orientations of the headgroup angle have been examined in presence of different amounts of charged or uncharged forms of lidocaine. Important differences of the membrane properties caused by the presence of the both forms of lidocaine are presented and discussed. Our simulations have shown that both charged and uncharged lidocaine cause almost the same increase of the dipole electrostatic potential in the middle of membrane though for different reasons. The increase, being about 90 mV for 9 mol % of lidocaine and 220 mV for 28 mol% of lidocaine, is of the size which may affect the functioning of voltage-gated ion channels.
Modification of the CHARMM force field for DMPC lipid bilayer
2008. Carl-Johan Högberg, Alexei Nikitin, Alexander Lyubartsev. Journal of Computational Chemistry 29 (14), 2359-2369Artikel
The CHARMM force field for DMPC lipids was modified in order to improve agreement with experiment for a number of important properties of hydrated lipid bilayer. The modification consists in introduction of a scaling factor 0.83 for 1-4 electrostatic interactions (between atoms separated by three covalent bonds), which provides correct transgauche ratio in the alkane tails, and recalculation of the headgroup charges on the basis of HF/6-311(d,p) ab-initio computations. Both rigid TIP3P and flexible SPC water models were used with the new lipid model, showing similar results. The new model in a 75 ns simulation has shown a correct value of the area per lipid at zero surface tension, as well as good agreement with the experiment for the electron density, structure factor, and order parameters, including those in the headgroup part of lipids.
New six-site acetonitrile model for simulations of liquid acetonitrile and its aqueous mixtures
2007. Alexei Nikitin, Alexander Lyubartsev. Journal of Computational Chemistry 28 (12), 2020-2026Artikel
A new six site flexible acetonitrile molecular model is developed. The AMBER force field was used for description of intramolecular parameters, the atomic charges were calculated from a high level ab initio theory and finally the Lennard-Jones parameters were tuned to fit the experimental density and evaporation heat. The obtained in this way model reproduces correctly densities of water-acetonitrile mixtures as well as provides qualitative description of the dielectric permittivity and self-diffusion coefficients.
Bead-Fourier path integral molecular dynamics for identical particles
2005. Sergei Ivanov, Alexander Lyubartsev. Journal of Chemical Physics 123, 034105Artikel
The Bead-Fourier path integral molecular dynamics technique introduced earlier [ S. D. Ivanov, A. P. Lyubartsev, and A. Laaksonen, Phys. Rev. E 67 066710 (2003) ] is applied for simulation of electrons in the simplest molecules: molecular hydrogen, helium atom, and their ions. Special attention is paid to the correct description of electrons in the core region of a nucleus. In an attempt to smooth the Coulomb potential at small distances, a recipe is suggested. The simulation results are in excellent agreement with the analytical solution for the “harmonic helium atom”, as well as with the vibrational potential of the H2 molecule and He ionization energies. It is demonstrated, that the Bead-Fourier path integral molecular dynamics technique is able to provide the accuracy required for the description of electron structure and chemical bonds in cases when electron exchange effects need not be taken into account.
Path-integral-expanded-ensemble Monte Carlo method in treatment of sign problem for fermions
2009. MA Voznesenskiy, PN Vorontsov-Velyaminov, AP Lyubartsev. Physical review. E 80 (6)Artikel
Expanded ensemble Monte Carlo method with Wang – Landau algorithm was used for calculations of the ratio of partition functions for classes of permutations in the problem of several interacting quantum particles (fermions) in an external field. Simulations for systems consisting from 2 and up to 7 interacting particles in harmonic or Coulombic field were performed. The presented approach allows one to carry out calculations for low enough temperatures which makes possible to extract data for the ground state energy and low temperature thermodynamics.
Phase transitions and thermodynamic properties of dense assemblies of truncated nanocubes and cuboctahedra
2012. Nikolai Volkov, Alexander Lyubartsev, Lennart Bergström. Nanoscale 4 (15), 4765-4771Artikel
Inspired by recent advances on the self-assembly of non-spherical nanoparticles, Monte Carlo simulations of the packing and thermodynamic properties of truncated nanocubes and cuboctahedra have been performed. The ergodicity problem was overcome by a modified Wang-Landau entropic sampling algorithm and equilibrium structural and thermodynamic properties were computed over a wide density range for both non-interacting and interacting particles. We found a structural transition from a simple cubic to a rhombohedral order when the degree of truncation exceeds a value of 0.9.
Quantum chemical and molecular dynamics modelling of hydroxylated polybrominated diphenyl ethers
2017. Inna Ermilova, Samuel Stenberg, Alexander P. Lyubartsev. Physical Chemistry, Chemical Physics - PCCP 19 (41), 28263-28274Artikel
A series of 19 hydroxylated polybrominated diphenyl ethers (OH-PBDEs) have been studied using density functional theory (DFT) and molecular dynamics simulations with the purpose of investigating eventual correlations between their physicochemical properties and toxic action. Dissociation constants (pK(a)), solvation free energies and octanol-water partition coefficients (logP) have been computed. Additionally, metadynamics simulations of OH-PBDEs passing through a lipid bilayer have been carried out for four OH-PBDE species. No correlations between computed pKa values and toxicity data have been found. Medium correlations were found between partition coefficients and the ability of OH-PBDEs to alter membrane potential in cell cultures, which is attributed to higher uptake of molecules with larger log P parameters. It was also demonstrated that in lipid bilayers, OH-PBDE molecules differ in their orientational distributions and can adopt different conformations which can affect the uptake of these molecules and influence the pathways of their toxic action.
Diffusion and reaction pathways of water near fully hydrated TiO2 surfaces from ab initio molecular dynamics
2017. Lorenzo Agosta, Erik G. Brandt, Alexander P. Lyubartsev. Journal of Chemical Physics 147 (2)Artikel
Ab initio molecular dynamics simulations are reported forwater-embedded TiO2 surfaces to determine the diffusive and reactive behavior at full hydration. A three-domain model is developed for six surfaces [rutile (110), (100), and (001), and anatase (101), (100), and (001)] which describes waters as hard (irreversibly bound to the surface), soft (with reduced mobility but orientation freedom near the surface), or bulk. The model explains previous experimental data and provides a detailed picture of water diffusion near TiO2 surfaces. Water reactivity is analyzed with a graph-theoretic approach that reveals a number of reaction pathways on TiO2 which occur at full hydration, in addition to direct water splitting. Hydronium (H3O+) is identified to be a key intermediate state, which facilitates water dissociation by proton hopping between intact and dissociated waters near the surfaces. These discoveries significantly improve the understanding of nanoscale water dynamics and reactivity at TiO2 interfaces under ambient conditions.
Curvature sensing by cardiolipin in simulated buckled membranes
Federico Elías-Wolff (et al.).Artikel
Efficient Production of Solar Hydrogen Peroxide Using Piezoelectric Polarization and Photoinduced Charge Transfer of Nanopiezoelectrics Sensitized by Carbon Quantum Dots
2022. Xiaofeng Zhou (et al.). Advanced ScienceArtikel
Piezoelectric semiconductors have emerged as redox catalysts, and challenges include effective conversion of mechanical energy to piezoelectric polarization and achieving high catalytic activity. The catalytic activity can be enhanced by simultaneous irradiation of ultrasound and light, but the existing piezoelectric semiconductors have trouble absorbing visible light. A piezoelectric catalyst is designed and tested for the generation of hydrogen peroxide (H2O2). It is based on Nb-doped tetragonal BaTiO3 (BaTiO3:Nb) and is sensitized by carbon quantum dots (CDs). The photosensitizer injects electrons into the conduction band of the semiconductor, while the piezoelectric polarization directed electrons to the semiconductor surface, allowing for a high-rate generation of H2O2. The piezoelectric polarization field restricts the recombination of photoinduced electron–hole pairs. A production rate of 1360 µmol gcatalyst−1 h−1 of H2O2 is achieved under visible light and ultrasound co-irradiation. Individual piezo- and photocatalysis yielded lower production rates. Furthermore, the CDs enhance the piezocatalytic activity of the BaTiO3:Nb. It is noted that moderating the piezoelectricity of BaTiO3:Nb via microstructure modulation influences the piezophotocatalytic activity. This work shows a new methodology for synthesizing H2O2 by using visible light and mechanical energy.
Exploring High-Pressure Transformations in Low-Z (H2, Ne) Hydrates at Low Temperatures
2022. Paulo H. B. Brant Carvalho (et al.). Crystals 12 (1)Artikel
The high pressure structural behavior of H2 and Ne clathrate hydrates with approximate composition H2/Ne·~4H2O and featuring cubic structure II (CS-II) was investigated by neutron powder diffraction using the deuterated analogues at ~95 K. CS-II hydrogen hydrate transforms gradually to isocompositional C1 phase (filled ice II) at around 1.1 GPa but may be metastably retained up to 2.2 GPa. Above 3 GPa a gradual decomposition into C2 phase (H2·H2O, filled ice Ic) and ice VIII’ takes place. Upon heating to 200 K the CS-II to C1 transition completes instantly whereas C1 decomposition appears sluggish also at 200 K. C1 was observed metastably up to 8 GPa. At 95 K C1 and C2 hydrogen hydrate can be retained below 1 GPa and yield ice II and ice Ic, respectively, upon complete release of pressure. In contrast, CS-II neon hydrate undergoes pressure-induced amorphization at 1.9 GPa, thus following the general trend for noble gas clathrate hydrates. Upon heating to 200 K amorphous Ne hydrate crystallizes as a mixture of previously unreported C2 hydrate and ice VIII’.
Semiconducting piezoelectric heterostructures for piezo- and piezophotocatalysis
2022. Xiaofeng Zhou (et al.). Nano Energy 96Artikel
Piezoelectric semiconductors can be polarized and used in mechanoredox systems and photoredox catalysis. Conventional non-piezoelectric semiconductors have limitations when it comes to charge carrier recombination and slow transport rates in catalytic reactions, which can be overcome by piezoelectric polarization processes in piezoelectric semiconductors. Heterostructures based on semiconducting piezoelectrics often offer enhanced catalytic reactivities that are related to their mechanical, piezoelectric, optical, and electronic characteristics. We review how to use such heterostructures to convert mechanical energy into chemical energy, and how the related piezoelectric polarization tunes the band structures and provides advantages in piezophotocatalysis over regular photocatalysis. We discuss fundamental concepts of piezoelectricity, piezoelectric potential, and examine different piezoelectric heterostructures for piezo- and piezophotocatalysis. A review of dynamic investigations of piezo- and piezophotocatalytic processes is presented. The complementary developments in the understanding of the piezotronic and piezophototronic effects are described, which include the induced charge-transfer mechanisms for piezo- and piezophotocatalytic reactions that can occur with piezoelectric heterostructures. Finally, we derive design principles and suggest future research directions in the emerging field of piezo- and piezophotocatalysis employing semiconductive heterostructures.
Atomistic Molecular Dynamics Simulations of Lipids Near TiO2 Nanosurfaces
2021. Mikhail Ivanov, Alexander P. Lyubartsev. Journal of Physical Chemistry B 125 (29), 8048-8059Artikel
Understanding of interactions between inorganic nanomaterials and biomolecules, and particularly lipid bilayers, is crucial in many biotechnological and biomedical applications, as well as for the evaluation of possible toxic effects caused by nanoparticles. Here, we present a molecular dynamics study of adsorption of two important constituents of the cell membranes, 1,2-dimyristoyl-sn-glycero-3-phosphocholine (DMPC) and 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphoethanolamine (POPE), lipids to a number of titanium dioxide planar surfaces, and a spherical nanoparticle under physiological conditions. By constructing the number density profiles of the lipid headgroup atoms, we have identified several possible binding modes and calculated their relative prevalence in the simulated systems. Our estimates of the adsorption strength, based on the total fraction of adsorbed lipids, show that POPE binds to the selected titanium dioxide surfaces stronger than DMPC, due to the ethanolamine group forming hydrogen bonds with the surface. Moreover, while POPE shows a clear preference toward anatase surfaces over rutile, DMPC has a particularly high affinity to rutile(101) and a lower affinity to other surfaces. Finally, we study how lipid concentration, addition of cholesterol, as well as titanium dioxide surface curvature may affect overall adsorption.
Atomistic Perspective on Biomolecular Adsorption on Functionalized Carbon Nanomaterials under Ambient Conditions
2021. Marzieh Saeedimasine, Erik G. Brandt, Alexander P. Lyubartsev. Journal of Physical Chemistry B 125 (1), 416-430Artikel
The use of carbon-based nanomaterials is tremendously increasing in various areas of technological, bioengineering, and biomedical applications. The functionality of carbon-based nanomaterials can be further broadened via chemical functionalization of carbon nanomaterial surfaces. On the other hand, concern is rising on possible adverse effects when nanomaterials are taken up by biological organisms. In order to contribute into understanding of interactions of carbon-based nanomaterials with biological matter, we have investigated adsorption of small biomolecules on nanomaterials using enhanced sampling molecular dynamics. The biomolecules included amino acid side chain analogues, fragments of lipids, and sugar monomers. The adsorption behavior on unstructured amorphous carbon, pristine graphene and its derivatives (such as few-layer graphene, graphene oxide, and reduced graphene oxide) as well as pristine carbon nanotubes, and those functionalized with OH-, COOH-, COO-, NH2-, and NH3+ groups was investigated with respect to surface concentration. An adsorption profile, that is, the free energy as a function of distance from the nanomaterial surfaces, was determined for each molecule and surface using the Metadynamics approach. The results were analyzed in terms of chemical specificity, surface charge, and surface concentration. It was shown that although morphology of the nanomaterial has a limited effect on the adsorption properties, functionalization of the surface by various molecular groups can drastically change the adsorption behavior that can be used in the design of nanosurfaces with highly selective adsorption properties and safe for human health and environment.
Bottom-Up Coarse-Grained Modeling of DNA
2021. Tiedong Sun (et al.). Frontiers in Molecular Biosciences 8Artikel
Recent advances in methodology enable effective coarse-grained modeling of deoxyribonucleic acid (DNA) based on underlying atomistic force field simulations. The so-called bottom-up coarse-graining practice separates fast and slow dynamic processes in molecular systems by averaging out fast degrees of freedom represented by the underlying fine-grained model. The resulting effective potential of interaction includes the contribution from fast degrees of freedom effectively in the form of potential of mean force. The pair-wise additive potential is usually adopted to construct the coarse-grained Hamiltonian for its efficiency in a computer simulation. In this review, we present a few well-developed bottom-up coarse-graining methods, discussing their application in modeling DNA properties such as DNA flexibility (persistence length), conformation, melting, and DNA condensation.
Computational insight into the hydrogenation of CO2 and carbamic acids to methanol by a ruthenium(II)-based catalyst
2021. Jinqin Yang (et al.). Molecular catalysis 506Artikel
Methanol is a liquid hydrogen carrier and potential platform molecule, and significant efforts are currently devoted to hydrogenate CO2 to methanol. In this work, hydrogenations of CO2 and captured CO2 (as dimethyl-carbamic acid, DMCA, and methylcarbamic acid, MCA) to methanol over Ru-MACHO (pre)catalysts were studied with a Ru-II-catalyst model by Density Functional Theory (DFT) calculations. For the hydrogenation of CO2, the concerted reaction was the rate-determining step involving a synchronous hydride transfer and proton transfer from the Ru-II-catalyst to coordinatively saturated intermediate methanediol. In the hydrogenation cycles of DMCA and MCA, the first hydride transfer reactions were more difficult than the concerted hydride and proton transfer from the Ru-II-catalyst to the aldehyde group of intermediate formaldehyde. These first hydride transfer reactions were identified as the rate-determining steps. The hydrogenations of DMCA and MCA were found much more favourable in methanol production than the direct CO2 hydrogenation, however, formamides could be main intermediate products due to the easier C-O breakage than C-N breakage in gem diols, and during the further hydrogenation of formamides, formaldehyde could be the main intermediate due to the easier C-N breakage than C-O breakage in alkanolamines. In all three hydrogenation cases, the amino (NH) ligand of the Ru-II-catalyst initially remained chemically innocent, and intermediates were stabilized by N-H center dot O hydrogenbonding interactions (HBIs) facilitating the continuation of catalytic hydrogenation cycles, but the NH ligand took part in multi-bond concerted reactions to produce eventually methanol.
First principles characterisation of bio-nano interface
2021. Ian Rouse (et al.). Physical Chemistry, Chemical Physics - PCCP (24)Artikel
Nanomaterials possess a wide range of potential applications due to their novel properties and exceptionally high activity as a result of their large surface to volume ratios compared to bulk matter. The active surface may present both advantage and risk when the nanomaterials interact with living organisms. As the overall biological impact of nanomaterials is triggered and mediated by interactions at the bio-nano interface, an ability to predict those from the atomistic descriptors, especially before the material is produced, can present enormous advantage for the development of nanotechnology. Fast screening of nanomaterials and their variations for specific biological effects can be enabled using computational materials modelling. The challenge lies in the range of scales that needs to be crossed from the material-specific atomistic representation to the relevant length scales covering typical biomolecules (proteins and lipids). In this work, we present a systematic multiscale approach that allows one to evaluate crucial interactions at the bionano interface from the first principles without any prior information about the material and thus establish links between the details of the nanomaterials structure to protein-nanoparticle interactions. As an example, an advanced computational characterization of titanium dioxide nanoparticles (6 different surfaces of rutile and anatase polymorphs) has been performed. We computed characteristics of the titanium dioxide interface with water using density functional theory for electronic density, used these parameters to derive an atomistic force field, and calculated adsorption energies for essential biomolecules on the surface of titania nanoparticles via direct atomistic simulations and coarse-grained molecular dynamics. Hydration energies, as well as adsorption energies for a set of 40 blood proteins are reported.
Pressure-induced amorphization of noble gas clathrate hydrates
2021. Paulo H. B. Brant Carvalho (et al.). Physical Review B 103 (6)Artikel
The high-pressure structural behavior of the noble gas (Ng) clathrate hydrates Ar center dot 6.5 H2O and Xe center dot 7.2 H2O featuring cubic structures II and I, respectively, was investigated by neutron powder diffraction (using the deuterated analogues) at 95 K. Both hydrates undergo pressure-induced amorphization (PIA), indicated by the disappearance of Bragg diffraction peaks, but at rather different pressures, at 1.4 and above 4.0 GPa, respectively. Amorphous Ar hydrate can be recovered to ambient pressure when annealed at >1.5 GPa and 170 K and is thermally stable up to 120 K. In contrast, it was impossible to retain amorphous Xe hydrate at pressures below 3 GPa. Molecular dynamics (MD) simulations were used to obtain general insight into PIA of Ng hydrates, from Ne to Xe. Without a guest species, both cubic clathrate structures amorphize at 1.2 GPa, which is very similar to hexagonal ice. Filling of large-sized H cages does not provide stability toward amorphization for structure II, whereas filled small-sized dodecahedral D cages shift PIA successively to higher pressures with increasing size of the Ng guest. For structure I, filling of both kinds of cages, large-sized T and small-sized D, acts to stabilize in a cooperative fashion. Xe hydrate represents a special case. In MD, disordering of the guest hydration structure is already seen at around 2.5 GPa. However, the different coordination numbers of the two types of guests in the crystalline cage structure are preserved, and the state is shown to produce a Bragg diffraction pattern. The experimentally observed diffraction up to 4 GPa is attributed to this semicrystalline state.
Structural investigation of three distinct amorphous forms of Ar hydrate
2021. Paulo H. B. Brant Carvalho (et al.). RSC Advances 11 (49), 30744-30754Artikel
Three amorphous forms of Ar hydrate were produced using the crystalline clathrate hydrate Ar·6.5H2O (structure II, Fdm, a ≈ 17.1 Å) as a precursor and structurally characterized by a combination of isotope substitution (36Ar) neutron diffraction and molecular dynamics (MD) simulations. The first form followed from the pressure-induced amorphization of the precursor at 1.5 GPa at 95 K and the second from isobaric annealing at 2 GPa and subsequent cooling back to 95 K. In analogy to amorphous ice, these amorphs are termed high-density amorphous (HDA) and very-high-density amorphous (VHDA), respectively. The third amorph (recovered amorphous, RA) was obtained when recovering VHDA to ambient pressure (at 95 K). The three amorphs have distinctly different structures. In HDA the distinction of the original two crystallographically different Ar guests is maintained as differently dense Ar–water hydration structures, which expresses itself in a split first diffraction peak in the neutron structure factor function. Relaxation of the local water structure during annealing produces a homogeneous hydration environment around Ar, which is accompanied with a densification by about 3%. Upon pressure release the homogeneous amorphous structure undergoes expansion by about 21%. Both VHDA and RA can be considered frozen solutions of immiscible Ar and water in which in average 15 and 11 water molecules, respectively, coordinate Ar out to 4 Å. The local water structures of HDA and VHDA Ar hydrates show some analogy to those of the corresponding amorphous ices, featuring H2O molecules in 5- and 6-fold coordination with neighboring molecules. However, they are considerably less dense. Most similarity is seen between RA and low density amorphous ice (LDA), which both feature strictly 4-coordinated H2O networks. It is inferred that, depending on the kind of clathrate structure and occupancy of cages, amorphous states produced from clathrate hydrates display variable local water structures.
To the fast calculation of the solvation free energy. Combining expanded ensembles with L2MC
2021. Alexei Nikitin, Vladislava Milchevskaya, Alexander P. Lyubartsev. Journal of Computational Chemistry 42 (11), 787-792Artikel
A more efficient version of the Expanded Ensembles method for calculation of free energy in molecular-mechanical simulations is proposed. The method is based on the Horowitz L2MC approach to accelerate movement along the alchemical coordinate. It is possible to achieve the same efficiency of the algorithm both with the optimal number of windows and with a larger number of them compared to the original algorithm. Since the optimal number of windows is unknown a priory, the proposed algorithm is more robust than the traditional one. We can choose the number of windows in excess and do not worry about the loss of efficiency. We illustrate the method's efficiency with the computation of the hydration free energies of pyridine and water.
Improved Sampling in Ab Initio Free Energy Calculations of Biomolecules at Solid-Liquid Interfaces
2020. Lorenzo Agosta, Erik G. Brandt, Alexander Lyubartsev. Computation 8 (1)Artikel
Atomistic simulations can complement the scarce experimental data on free energies of molecules at bio-inorganic interfaces. In molecular simulations, adsorption free energy landscapes are efficiently explored with advanced sampling methods, but classical dynamics is unable to capture charge transfer and polarization at the solid-liquid interface. Ab initio simulations do not suffer from this flaw, but only at the expense of an overwhelming computational cost. Here, we introduce a protocol for adsorption free energy calculations that improves sampling on the timescales relevant to ab initio simulations. As a case study, we calculate adsorption free energies of the charged amino acids Lysine and Aspartate on the fully hydrated anatase (101) TiO2 surface using tight-binding forces. We find that the first-principle description of the system significantly contributes to the adsorption free energies, which is overlooked by calculations with previous methods.
Modeling DNA Flexibility
2020. Vishal Minhas (et al.). Journal of Physical Chemistry B 124 (1), 38-49Artikel
Accurate parametrization of force fields (FFs) is of ultimate importance for computer simulations to be reliable and to possess a predictive power. In this work, we analyzed, in multi-microsecond simulations of a 40-base-pair DNA fragment, the performance of four force fields, namely, the two recent major updates of CHARMM and two from the AMBER family. We focused on a description of double-helix DNA flexibility and dynamics both at atomistic and at mesoscale level in coarse-grained (CG) simulations. In addition to the traditional analysis of different base-pair and base-step parameters, we extended our analysis to investigate the ability of the force field to parametrize a CG DNA model by structure-based bottom-up coarse-graining, computing DNA persistence length as a function of ionic strength. Our simulations unambiguously showed that the CHARMM36 force field is unable to preserve DNA's structural stability at over-microsecond time scale. Both versions of the AMBER FF, parmbsc0 and parmbsc1, showed good agreement with experiment, with some bias of parmbsc0 parameters for intermediate A/B form DNA structures. The CHARMM27 force field provides stable atomistic trajectories and overall (among the considered force fields) the best fit to experimentally determined DNA flexibility parameters both at atomistic and at mesoscale level.
Modelling of interactions between Aβ(25-35) peptide and phospholipid bilayers
2020. Inna Ermilova, Alexander Lyubartsev. RSC Advances 10 (7), 3902-3915Artikel
Aggregation of amyloid beta (Aβ) peptides in neuronal membranes is a known promoter of Alzheimer’s disease. To gain insight into the molecular details of Aβ peptide aggregation and its effect on model neuronal membranes, we carried out molecular dynamics simulations of the Aβ(25–35) fragment of the amyloid precursor protein in phospholipid bilayers composed of either fully saturated or highly unsaturated lipids, in the presence or absence of cholesterol. It was found that the peptide does not penetrate through any of the considered membranes, but can reside in the headgroup region and upper part of the lipid tails showing a clear preference to a polyunsaturated cholesterol-free membrane. Due to the ordering and condensing effect upon addition of cholesterol, membranes become more rigid facilitating peptide aggregation on the surface. Except for the case of the cholesterol-free saturated lipid bilayer, the peptides have a small effect on the membrane structure and ordering. It was also found that the most “active” amino-acid for peptide–lipid and peptide–cholesterol interaction is methionine-35, followed by asparagine-27 and serine-26, which form hydrogen bonds between peptides and polar atoms of lipid headgroups. These amino acids are also primarily responsible for peptide aggregation. This work will be relevant for designing strategies to develop drugs to combat Alzheimer’s disease.
Optimization of Slipids Force Field Parameters Describing Headgroups of Phospholipids
2020. Fredrik Grote, Alexander P. Lyubartsev. Journal of Physical Chemistry B 124 (40), 8784-8793Artikel
The molecular mechanics force field Slipids developed in a series of works by Jambeck and Lyubartsev (J. Phys. Chem. B 2012, 116, 3164-3179; J. Chem. Theory Comput. 2012, 8, 2938-2948) generally provides a good description of various lipid bilayer systems. However, it was also found that order parameters of C-H bonds in the glycerol moiety of the phosphatidylcholine headgroup deviate significantly from NMR results. In this work, the dihedral force field parameters have been reparameterized in order to improve the agreement with experiment. For this purpose, we have computed energies for a large amount of lipid headgroup conformations using density functional theory on the B3P86/cc-pvqz level and optimized dihedral angle parameters simultaneously to provide the best fit to the quantum chemical energies. The new parameter set was validated for three lipid bilayer systems against a number of experimental properties including order parameters, area per lipid, scattering form factors, bilayer thickness, area compressibility and lateral diffusion coefficients. In addition, the order parameter dependence on cholesterol content in the POPC bilayer was investigated. It is shown that the new force field significantly improves agreement with the experimental order parameters for the lipid headgroup while keeping good agreement with other experimentally measured properties.
A multiscale analysis of DNA phase separation
2019. Tiedong Sun (et al.). Nucleic Acids Research 47 (11), 5550-5562Artikel
DNA condensation and phase separation is of utmost importance for DNA packing in vivo with important applications in medicine, biotechnology and polymer physics. The presence of hexagonally ordered DNA is observed in virus capsids, sperm heads and in dinoflagellates. Rigorous modelling of this process in all-atom MD simulations is presently difficult to achieve due to size and time scale limitations. We used a hierarchical approach for systematic multiscale coarse-grained (CG) simulations of DNA phase separation induced by the three-valent cobalt(III)-hexammine (CoHex(3+)). Solvent-mediated effective potentials for a CG model of DNA were extracted from all-atom MD simulations. Simulations of several hundred 100-bp-long CG DNA oligonucleotides in the presence of explicit CoHex(3+) ions demonstrated aggregation to a liquid crystalline hexagonally ordered phase. Following further coarse-graining and extraction of effective potentials, we conducted modelling at mesoscale level. In agreement with electron microscopy observations, simulations of an 10.2-kblong DNA molecule showed phase separation to either a toroid or a fibre with distinct hexagonal DNA packing. The mechanism of toroid formation is analysed in detail. The approach used here is based only on the underlying all-atom force field and uses no adjustable parameters and may be generalised to modelling chromatin up to chromosome size.
A multiscale model of protein adsorption on a nanoparticle surface
2019. David Power (et al.). Modelling and Simulation in Materials Science and Engineering 27 (8)Artikel
We present a methodology to quantify the essential interactions at the interface between inorganic solid nanoparticles (NPs) and biological molecules. Our model is based on pre-calculation of the repetitive contributions to the interaction from molecular segments, which allows us to efficiently scan a multitude of molecules and rank them by their adsorption affinity. The interaction between the biomolecular fragments and the nanomaterial are evaluated using a systematic coarse-graining scheme starting from all-atom molecular dynamics simulations. The NPs are modelled using a two-layer representation, where the outer layer is parameterized at the atomistic level and the core is treated at the continuum level using Lifshitz theory of dispersion forces. We demonstrate that the scheme reproduces the experimentally observed features of the NP protein coronas. To illustrate the use of the methodology, we compute the adsorption energies for human blood plasma proteins on gold NPs of different sizes as well as the preferred orientation of the molecules upon adsorption. The computed energies can be used for predicting the composition of the NP-protein corona for the corresponding material.
Cholesterol in phospholipid bilayers
2019. Inna Ermilova, Alexander P. Lyubartsev. Soft Matter 15 (1), 78-93Artikel
Cholesterol is an essential component of all animal cell membranes and plays an important role in maintaining the membrane structure and physical–chemical properties necessary for correct cell functioning. The presence of cholesterol is believed to be responsible for domain formation (lipid rafts) due to different interactions of cholesterol with saturated and unsaturated lipids. In order to get detailed atomistic insight into the behaviour of cholesterol in bilayers composed of lipids with varying degrees of unsaturation, we have carried out a series of molecular dynamics simulations of saturated and polyunsaturated lipid bilayers with different contents of cholesterol, as well as well-tempered metadynamics simulations with a single cholesterol molecule in these bilayers. From these simulations we have determined distributions of cholesterol across the bilayer, its orientational properties, free energy profiles, and specific interactions of molecular groups able to form hydrogen bonds. Both molecular dynamics and metadynamics simulations showed that the most unsaturated bilayer with 22:6 fatty acid chains shows behaviour which is most different from other lipids. In this bilayer, cholesterol is relatively often found in a “flipped” configuration with the hydroxyl group oriented towards the membrane middle plane. This bilayer has also the highest (least negative) binding free energy among liquid phase bilayers, and the lowest reorientation barrier. Furthermore, cholesterol molecules in this bilayer are often found to form head-to-tail contacts which may lead to specific clustering behaviour. Overall, our simulations support ideas that there can be a subtle interconnection between the contents of highly unsaturated fatty acids and cholesterol, deficiency or excess of each of them is related to many human afflictions and diseases.
Curvature sensing by cardiolipin in simulated buckled membranes
2019. Federico Elias-Wolff (et al.). Soft Matter 15 (4), 792-802Artikel
Cardiolipin is a non-bilayer phospholipid with a unique dimeric structure. It localizes to negative curvature regions in bacteria and is believed to stabilize respiratory chain complexes in the highly curved mitochondrial membrane. Cardiolipin's localization mechanism remains unresolved, because important aspects such as the structural basis and strength for lipid curvature preferences are difficult to determine, partly due to the lack of efficient simulation methods. Here, we report a computational approach to study curvature preferences of cardiolipin by simulated membrane buckling and quantitative modeling. We combine coarse-grained molecular dynamics with simulated buckling to determine the curvature preferences in three-component bilayer membranes with varying concentrations of cardiolipin, and extract curvature-dependent concentrations and lipid acyl chain order parameter profiles. Cardiolipin shows a strong preference for negative curvatures, with a highly asymmetric chain order parameter profile. The concentration profiles are consistent with an elastic model for lipid curvature sensing that relates lipid segregation to local curvature via the material constants of the bilayers. These computations constitute new steps to unravel the molecular mechanism by which cardiolipin senses curvature in lipid membranes, and the method can be generalized to other lipids and membrane components as well.
Effects of lipid saturation on amyloid-beta peptides partitioning and aggregation in neuronal membranes
2019. Nikolaos Ntarakas, Inna Ermilova, Alexander Lyubartsev. European Biophysics JournalArtikel
Implicit solvent systematic coarse-graining of dioleoylphosphatidylethanolamine lipids
2019. Saeed Mortezazadeh (et al.). PLOS ONE 14 (4)Artikel
Lamellar and hexagonal lipid structures are of particular importance in the biological processes such as membrane fusion and budding. Atomistic simulations of formation of these phases and transitions between them are computationally prohibitive, hence development of coarse-grained models is an important part of the methodological development in this area. Here we apply systematic bottom-up coarse-graining to model different phase structures formed by 1,2-dioleoylphosphatidylethanolamine (DOPE) lipid molecules. We started from atomistic simulations of DOPE lipids in water carried out at two different water/lipid molar ratio corresponding to the lamellar L-alpha and inverted hexagonal H-II structures at low and high lipid concentrations respectively. The atomistic trajectories were mapped to coarse-grained trajectories, in which each lipid was represented by 14 coarse-grained sites. Then the inverse Monte Carlo method was used to compute the effective coarse-grained potentials which for the coarse-grain model reproduce the same structural properties as the atomistic simulations. The potentials derived from the low concentration atomistic simulation were only able to form a bilayer structure, while both L-alpha and H-II lipid phases were formed in simulations with potentials obtained at high concentration. The typical atomistic configurations of lipids at high concentration combine fragments of both lamellar and non-lamellar structures, that is reflected in the extracted coarse-grained potentials which become transferable and can form a wide range of structures including the inverted hexagonal, bilayer, tubule, vesicle and micellar structures.
2019. Alexander Mirzoev, Lars Nordenskiöld, Alexander Lyubartsev. Computer Physics Communications 237, 263-273Artikel
Molecular simulations of many phenomena related to biomolecular systems, soft matter and nanomaterials require consideration of length scales above 10 nm and time scales longer than 1 mu s, which necessitates the use of coarse-grained (low resolution) models, where each site of the model represents a group of atoms, and where the solvent is often omitted. Our software package MagiC is designed to perform systematic structure-based coarse-graining of molecular models, in which the effective pairwise potentials between coarse-grained sites of low-resolution molecular models are constructed to reproduce structural distribution functions obtained from modeling of systems in a high resolution (atomistic) description. The software takes as input atomistic trajectories generated by an external molecular dynamics package, and produce as an output interaction potentials for coarse-grained models which can be directly used in a coarse-grained simulations package. Here we present a major update (v.3) of the software with substantially improved functionality, compatibility with several major atomistic and coarse-grained simulations packages (GROMACS, LAMMPS, GALAMOST), analysis suite with graphical possibilities, diagnostics, documentation. We describe briefly the coarse-graining methodology, the structure of the software, describe users actions, and illustrate the whole process with two complex examples: cholesterol containing lipid bilayers and condensation of DNA caused by multivalent ions. Program summary Program Title: MagiC Program Files doi : http://dx.doi.org/10.17632/9gnfxyshj8.1 Licensing provisions: GPLv3 Programming language: Fortran, Python Nature of problem: Systematic bottom-up multiscale modeling is a complex multi-stage process, in which results of simulations of a high-resolution (atomistic) model are used to construct a low resolution (coarse-grained) model, providing the same structural properties for the coarse-grained system as for the high-resolution system. Within the approach, structural properties of the high-resolution model are computed in terms of radial distribution functions and distributions of intramolecular degrees of freedom. Then the inverse problem is solved, in which the interaction potentials for the low-resolution model are determined which provide distribution functions coinciding with those obtained in the high-resolution simulations. The low-resolution model can be then used for simulations of the same system on larger length and time-scale. Solution method: The presented software package implements all stages of the systematic structure-based coarse-graining. It works as an integrated pipeline, giving the user ability to easily derive a coarse-grained model for a multicomponent complex molecular system and then use it for large-scale simulations. MagiC implements two approaches to solve the inverse problem: (i) the Inverse Monte Carlo method in which the inverse problem is solved using the Newton-Raphson method, with inversion of the Jacobian for the discretized relationship between interaction potentials and structural distribution functions; and (ii) the Iterative Boltzmann approach in which the inverse problem is solved using approximative relationships neglecting correlations between different degrees of freedom. The inverse solver includes also variational Inverse Monte Carlo approach when some of the coarse-grained potentials are fixed while others vary in order to fit the whole set of reference distribution functions. Additional comments: The MagiC main module can be also used for conventional Monte Carlo simulations of molecular systems described by tabulated pairwise potentials in the canonical ensemble.
A systematic analysis of nucleosome core particle and nucleosome-nucleosome stacking structure
2018. Nikolay Korolev, Alexander P. Lyubartsev, Lars Nordenskiöld. Scientific Reports 8Artikel
Chromatin condensation is driven by the energetically favourable interaction between nucleosome core particles (NCPs). The close NCP-NCP contact, stacking, is a primary structural element of all condensed states of chromatin in vitro and in vivo. However, the molecular structure of stacked nucleosomes as well as the nature of the interactions involved in its formation have not yet been systematically studied. Here we undertake an investigation of both the structural and physico-chemical features of NCP structure and the NCP-NCP stacking. We introduce an NCP-centred set of parameters (NCP-NCP distance, shift, rise, tilt, and others) that allows numerical characterisation of the mutual positions of the NCPs in the stacking and in any other structures formed by the NCP. NCP stacking in more than 140 published NCP crystal structures were analysed. In addition, coarse grained (CG) MD simulations modelling NCP condensation was carried out. The CG model takes into account details of the nucleosome structure and adequately describes the long range electrostatic forces as well as excluded volume effects acting in chromatin. The CG simulations showed good agreement with experimental data and revealed the importance of the H2A and H4 N-terminal tail bridging and screening as well as tail-tail correlations in the stacked nucleosomes.
Computing Curvature Sensitivity of Biomolecules in Membranes by Simulated Buckling
2018. Federico Elías-Wolff (et al.). Journal of Chemical Theory and Computation 14 (3), 1643-1655Artikel
Membrane curvature sensing, where the binding free energies of membrane-associated molecules depend on the local membrane curvature, is a key factor to modulate and maintain the shape and organization of cell membranes. However, the microscopic mechanisms are not well understood, partly due to absence of efficient simulation methods. Here, we describe a method to compute the curvature dependence of the binding free energy of a membrane associated probe molecule that interacts with a buckled membrane, which has been created by lateral compression of a flat bilayer patch. This buckling approach samples a wide range of curvatures in a single simulation, and anisotropic effects can be extracted from the orientation statistics. We develop an efficient and robust algorithm to extract the motion of the probe along the buckled membrane surface, and evaluate its numerical properties by extensive sampling of three coarse-grained model systems: local lipid density in a curved environment for single-component bilayers, curvature preferences of individual lipids in two-component membranes, and curvature sensing by a homotrimeric transmembrane protein. The method can be used to complement experimental data from curvature partition assays and provides additional insight into mesoscopic theories and molecular mechanisms for curvature sensing.
Molecular Dynamics Simulations of Furfural and 5-Hydroxymethylfurfural at Ambient and Hydrothermal Conditions
2018. Fredrik Grote, Inna Ermilova, Alexander P. Lyubartsev. Journal of Physical Chemistry B 122 (35), 8416-8428Artikel
In this work, we present results from molecular dynamics simulations of aqueous solutions of furfural and 5-hydroxymethylfurfural, which are important intermediates in the hydrothermal carbonization processes of biomass conversion. The computations were performed both at ambient and hydrothermal conditions using a two-level factorial design varying concentration, temperature, and pressure. A number of equilibrium and dynamic properties have been computed including enthalpies and free energies of vaporization, free energies of solvation, diffusion coefficients, and rotational/reorientational correlation times. Structural properties of solutions were analyzed using radial and spatial distribution functions. It was shown that the formation of hydrogen bonds among 5-hydroxymethylfurfural molecules is preferred compared to hydrogen bonding between 5-hydroxymethylfurfural and water. In addition, our results suggest that the oxygen atoms in the furan rings of furfural and 5-hydroxymethylfurfural do not participate in hydrogen bonding to the same extent as the oxygen atoms in the hydroxyl and carbonyl groups. It is also observed that furfural molecules aggregate under certain conditions, and we show how this is affected by changes in temperature, pressure, and concentration in agreement with experimental solubility data. The analysis of the computational results provides useful insight into the structure and dynamics of the considered molecules at conditions of hydrothermal carbonization, as well as at ambient conditions.
Stress Relief and Reactivity Loss of Hydrated Anatase (001) Surface
2018. Eugenio Vitale (et al.). The Journal of Physical Chemistry C 122 (39), 22407-22417Artikel
Dissociative and molecular water adsorption on the anatase (001) surface is studied in the context of state-of-the-art density functional theory in large supercells suited for adsorption studies at various water coverage ratios. At low coverage values below 1/4 ML, water adsorption remains dissociative and a network of hydrogen bonds between the so formed hydroxyl groups favors the formation of a ridge surface structure. The hydroxyl patterned (4 X 4) surface thus undergoes a (2 X 4) reconstruction that causes the relief of the large tensile stress measured in the unreconstructed surface along the direction orthogonal to the ridge. This phenomenology is accompanied by the loss of reactivity of the reconstructed surface with respect to the dissociative water adsorption that becomes molecular above 1/4 ML. We also show that the molecular adsorption on the terrace is weaker than the one on the ridge. The present water reconstruction model is discussed and compared to the well-known ADM model of the reconstructed anatase (001) surface in dry environment.
Unperturbed hydrocarbon chains and liquid phase bilayer lipid chains
2018. Alexander L. Rabinovich, Alexander P. Lyubartsev, Dmitrii V. Zhurkin. European Biophysics Journal 47 (2), 109-130Artikel
In this work, the properties of saturated and unsaturated fatty acid acyl chains 16:0, 18:0, 18:1(n-9)cis, 18:2(n-6)cis, 18:3(n-3)cis, 18:4(n-3)cis, 18:5(n-3)cis, 20:4(n-6)cis, 20:5(n-3)cis and 22:6(n-3)cis in a bilayer liquid crystalline state and similar hydrocarbon chains (with CH terminal groups instead of C=O groups) in the unperturbed state characterised by a lack of long-range interaction were investigated. The unperturbed hydrocarbon chains were modelled by Monte Carlo simulations at temperature K; sixteen fully hydrated homogeneous liquid crystalline phosphatidylcholine bilayers containing these chains were studied by molecular dynamics simulations at the same temperature. To eliminate effects of the simulation parameters, the molecular dynamics and Monte Carlo simulations were carried out using the same structural data and force field coefficients. From these computer simulations, the average distances between terminal carbon atoms of the chains (end-to-end distances) were calculated and compared. The trends in the end-to-end distances obtained for the unperturbed chains were found to be qualitatively similar to those obtained for the same lipid chains in the bilayers. So, for understanding of a number of processes in biological membranes (e.g., changes in fatty acid composition caused by environmental changes such as temperature and pressure), it is possible to use, at least as a first approximation, the relationships between the structure and properties for unperturbed or isolated hydrocarbon chains.
All-Atom MD Simulation of DNA Condensation Using Ab lnitio Derived Force Field Parameters of Cobalt(III)-Hexammine
2017. Tiedong Sun (et al.). Journal of Physical Chemistry B 121 (33), 7761-7770Artikel
It is well established that the presence of the trivalent cobalt(III)-hexammine cation (CoHex(3+)) at submillimolar concentrations leads to bundling (condensation) of double stranded DNA molecules, which is caused by DNA DNA attraction induced by the multivalent counterions. However, the detailed mechanism of this process is still not fully understood. Furthermore, in all-atom molecular dynamics (MD) simulations, spontaneous aggregation of several DNA oligonucleotides in the presence of CoHex(3+) has previously not been, demonstrated. In order to obtain a rigorous description of CoHex(3+)-nucleic acid interactions and CoHex(3+)-induced DNA condensation to be used in MD siniulations, we have derived optimized force field parameters of the CoHex(3+) ion. They were obtained from Car Parrinello molecular dynamics simulation of a single CoHex3+ ion in the presence of 126 water molecules. The new set,of force field parameters reproduces the experimentally known transition of DNA from B- to A-form; and qualitatively describes changes of DNA and RNA persistence lengths. We then carried out a 2 mu s long atomistic simulation of four DNA oligomers each consisting of 36 base pairs in the presence of CoHex(3+). We demonstrate that, in this system, DNA molecules display attractive interactions and aggregate into bundle-like structures. This behavior depends critically on the details of the CoHex(3+) interaction with DNA. A control simulation with a similar setup but in the presence of Mg2+ does not induce DNA DNA attraction, which is also in agreement with experiment.
Anesthetics mechanism on a DMPC lipid membrane model
2017. Marzieh Saeedi, Alexander P. Lyubartsev, Seifollah Jalili. Biophysical Chemistry 226, 1-13Artikel
To provide insight into the molecular mechanisms of local anesthetic action, we have carried out an extensive investigation of two amide type local anesthetics, lidocaine and articaine in both charged and uncharged forms, interacting with DMPC lipid membrane. We have applied both standard molecular dynamics simulations and metadynamics simulations to provide a detailed description of the free energy landscape of anesthetics embedded in the lipid bilayer. The global minimum of the free energy surface (equilibrium position of anesthetics in the lipid membrane) occurred around 1nm of the bilayer center. The uncharged anesthetics show more affinity to bind to this region compared to the charged drugs. The binding free energy of uncharged lidocaine in the membrane (-30.3kJ/mol) is higher than uncharged articaine (-24.0kJ/mol), which is in good agreement with higher lipid solubility of lidocaine relative to the articaine. The octanol/water partition coefficient of uncharged drugs was also investigated using expanded ensemble simulations. In addition, complementary standard MD simulations were carried out to study the partitioning behavior of multiple anesthetics inside the lipid bilayer. The results obtained here are in line with previously reported simulations and suggest that the different forms of anesthetics induce different structural modifications in the lipid bilayer, which can provide new insights into their complex membrane translocation phenomena.
Computationally based analysis of the energy efficiency of a CO2 capture process
2017. Bojan Vujic, Alexander P. Lyubartsev. Chemical Engineering Science 174, 174-188Artikel
We propose a model for thermodynamic evaluation of the energy efficiency of a CO2 capture in a temperature-pressure swing adsorption. Major parts of this model are computational prediction of the adsorbed gas loading as a function of temperature and partial CO2 pressure, evaluation of the energy expenses under specified conditions for the working capacity, regenerability of the sorbent and purity of the captured CO2, as well as determination of the most optimal desorption conditions in terms of desorption pressure and temperature. The proposed model can be applied for fast evaluation of the energy costs of the CO2 capture process with the use of both experimental or simulation adsorption data with respect to pressure and temperature. We tested this model analyzing data obtained from Grand Canonical Monte Carlo simulations for more than thousand different zeolite structures. Within our approach it is possible to evaluate a theoretical limit of the energy expenses for each specific material and to use the proposed method in screening different structures for the most efficient sorbent material from the energy efficiency point of view under specified requirements for the working capacity of the process, regenerability and purity of captured CO2. We show that setting realistic from the industrial point of view parameters of the CO2 capture cycle leads to substantial reduction of the number of suitable zeolite structures, and to increase of the energy penalty of the CO2 capture compared to evaluations based on minimization of the parasitic energy only.