Alexander Lyubartsev

Alexander Lyubartsev


Visa sidan på svenska
Works at Department of Materials and Environmental Chemistry
Telephone 08-16 11 93
Visiting address Svante Arrhenius väg 16 C
Room C 548
Postal address Institutionen för material- och miljökemi 106 91 Stockholm


A selection from Stockholm University publication database
  • 2017. Matúš Rebič (et al.). ACS Omega 2 (2), 386-396

    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.

  • 2016. Inna Ermilova, Alexander P. Lyubartsev. Journal of Physical Chemistry B 120 (50), 12826-12842

    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.

  • 2017. I. A. Silanteva (et al.). Nanosystems: Physics, Chemistry, Mathematics 8 (1), 108-120

    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.

  • 2016. Alessio Atzori (et al.). Canadian journal of chemistry (Print) 94 (12), 1181-1188

    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.

  • 2016. Alexander P. Lyubartsev, Alexander L. Rabinovich. Biochimica et Biophysica Acta - Biomembranes 1858 (10), 2483-2497

    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.

  • 2016. Nikolay Korolev, Lars Nordenskiöld, Alexander P. Lyubartsev. Advances in Colloid and Interface Science 232, 36-48

    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.

  • 2016. Erik G. Brandt, Lorenzo Agosta, Alexander P. Lyubartsev. Nanoscale 8 (27), 13385-13398

    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.

  • 2016. Bojan Vujić, Alexander P. Lyubartsev. Modelling and Simulation in Materials Science and Engineering 24 (4)

    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.

  • 2015. Chaitanya A. K. Koppisetty (et al.). Journal of Computer-Aided Molecular Design 29 (1), 13-21

    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.

  • 2015. Erik G. Brandt, Alexander P. Lyubartsev. The Journal of Physical Chemistry C 119 (32), 18126-18139

    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.

  • 2015. Alexander P. Lyubartsev (et al.). Journal of Physics 27 (6)

    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+).

  • 2015. Erik G. Brandt, Alexander P. Lyubartsev. The Journal of Physical Chemistry C 119 (32), 18110-18125

    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.

  • 2015. Alexander P. Lyubartsev (et al.). Journal of Chemical Physics 143 (24)

    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.

  • 2014. Nikolay Korolev (et al.). Polymers 6 (6), 1655-1675

    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.

  • 2014. Alexei M. Nikitin, Yury V. Milchevskiy, Alexander P. Lyubartsev. Journal of Molecular Modeling 20 (3), 2143

    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.

  • 2014. Alexander L. Rabinovich, Alexander P. Lyubartsev. 25th IUPAP Conference on Computational Physics (CCP2013)

    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.

  • 2014. Nikolay Korolev (et al.). Biopolymers 101 (10), 1051-1064

    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.

  • 2014. Joakim P. M. Jämbeck (et al.). Journal of Chemical Theory and Computation 10 (1), 5-13

    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.

  • 2014. Alexander Mirzoev, Alexander P. Lyubartsev. Journal of Computational Chemistry 35 (16), 1208-1218

    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.

  • 2014. Joakim P. M. Jämbeck, Alexander P. Lyubartsev. Journal of Physical Chemistry B 118 (14), 3793-3804

    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.

  • 2013. Bertrand Faure (et al.). Nanoscale 5 (3), 953-960

    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.

  • 2013. Yanping Fan (et al.). PLoS ONE 8 (2), e54228

    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.

  • 2013. Joakim P. M. Jämbeck, Alexander P. Lyubartsev. Journal of Chemical Theory and Computation 9 (1), 774-784

    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.

  • 2013. A. L. Rabinovich, Alexander P. Lyubartsev. POLYMER SCIENCE SERIES C 55 (1), 162-180

    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).

  • 2013. Joakim P. M. Jämbeck, Alexander P. Lyubartsev. Journal of Physical Chemistry Letters 4 (11), 1781-1787

    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.

  • 2013. Joakim P. M. Jämbeck, Alexander P. Lyubartsev. Physical Chemistry, Chemical Physics - PCCP 15 (13), 4677-4686

    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.

  • Article MagiC
    2013. Alexander Mirzoev, Alexander P. Lyubartsev. Journal of Chemical Theory and Computation 9 (3), 1512-1520

    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.

  • 2013. Joakim P. M. Jämbeck (et al.). Journal of Computational Chemistry 34 (3), 187-197

    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.

  • 2012. Joakim P. M. Jämbeck, Alexander P. Lyubartsev. Journal of Chemical Theory and Computation 8 (8), 2938-2948

    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.

  • 2012. Joakim P. M. Jämbeck, Alexander P. Lyubartsev. Journal of Physical Chemistry B 116 (10), 3164-3179

    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.

  • 2012. Aatto Laaksonen, Alexander Lyubartsev, Francesca Mocci. Molecular Dynamics - Studies of Synthetic and Biological Macromolecules, 85-106
  • 2012. Nikolay Korolev (et al.). Current opinion in structural biology 22 (2), 151-159

    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.

  • 2012. Nikolay Korolev (et al.). Soft Matter 8 (36), 9322-9333

    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.

  • 2011. P. N. Vorontsov-Velyaminov (et al.). Contributions to Plasma Physics 51 (4), 382-385

    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.

  • 2011. Alexander Mirzoev, Alexander P. Lyubartsev. Physical Chemistry, Chemical Physics - PCCP (13), 5722-5727

    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.

  • 2011. Andrei V. Egorov, Alexander P. Lyubartsev, Aatto Laaksonen. Journal of Physical Chemistry B 115 (49), 14572-14581

    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.

  • 2011. Alexander P. Lyubartsev, Alexander L. Rabinovich. Soft Matter 7 (1), 25-39

    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.

  • 2011. Mysore Sridhar Santosh (et al.). Journal of Solution Chemistry 40 (9), 1657-1671

    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.

  • 2011. Nikolai Volkov, Pavel Vorontsov-Velyaminov, Alexander Lyubartsev. Macromolecular Theory and Simulations 20 (7), 496-509

    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.

  • 2010. Nikolai Korolev, Alexander Lyubartsev, Lars Nordenskiöld. Journal of Colloid and Interface Science 158, 32-47

    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.

  • 2010. Eugeny A. Polyakov, Alexander P. Lyubartsev, Pavel N. Vorontsov-Velyaminov. Journal of Chemical Physics 133 (19), 194103

    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.

  • 2010. Nikolai Korolev (et al.). Biophysical Journal 99, 1896-1905

    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.

  • 2010. Mysore S. Santosh (et al.). Journal of Physical Chemistry B 114 (49), 16632-16640

    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.

  • 2010. Enamul Mojumdar, Alexander Lyubartsev. Biophysical Chemistry 153 (1), 27-35

    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.

  • 2010. Alexander Lyubartsev (et al.). Faraday discussions (Online) 144, 43-56

    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.

  • 2010. Mikael Leetmaa (et al.). Journal of Electron Spectroscopy and Related Phenomena 177 (2-3), 135-157

    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.

  • 2009. Mathias Ljungberg (et al.). Journal of Chemical Physics 131, 034501

    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).

  • 2009. Ye Yang (et al.). Biophysical Journal 96, 2082-2096

    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.

  • 2009. C.G. Jesudason, Alexander Lyubartsev, Aatto Laaksonen. The European Physical Journal E Soft matter 30, 341-350

    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.

  • 2009. Alexander Lyubartsev, Yaoquan Tu, Aatto Laaksonen. Journal of computational and theoretical nanoscience 6 (5), 951-959

    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.

  • 2009. Alexander Lyubartsev, Arieh Ben-Naim. Journal of Chemical Physics 131, 204507

    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.

  • 2009. M. A. Voznesenskiy, P. N. Vorontsov-Velyaminov, Alexander Lyubartsev. PHYS REV E 80 (6), 66702

    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.

  • 2009. Eugeniy Polyakov, Pavel Vorontsov-Velyaminov, Alexander Lyubartsev. Vychislitel'nye Metody i Programirovaniye [Numerical methods and programming] 10, 223-247

    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.


  • 2009. Kun Dong (et al.). The Journal of Physical Chemistry C 113 (23), 10013-10020

    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.

  • 2009. Congcong Huang (et al.). Proceedings of the National Academy of Sciences of the United States of America 106, 15214-15218

    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-238
  • 2008. Carl-Johan Högberg, Alexander Lyubartsev. Biophysical Journal 94, 525-531

    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.

  • 2008. Carl-Johan Högberg, Alexei Nikitin, Alexander Lyubartsev. Journal of Computational Chemistry 29 (14), 2359-2369

    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.

  • 2007. Alexei Nikitin, Alexander Lyubartsev. Journal of Computational Chemistry 28 (12), 2020-2026

    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.

  • 2005. Sergei Ivanov, Alexander Lyubartsev. Journal of Chemical Physics 123, 034105

    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.

Show all publications by Alexander Lyubartsev at Stockholm University

Last updated: June 30, 2017

Bookmark and share Tell a friend