Extension of the Slipids Force Field to Polyunsaturated Lipids - The

Subscriber access provided by UNIV OF REGINA

Article

Extension of the Slipids Force Field to Polyunsaturated Lipids Inna Ermilova, and Alexander P. Lyubartsev J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.6b05422 • Publication Date (Web): 30 Nov 2016 Downloaded from http://pubs.acs.org on December 12, 2016

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry B is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Extension of the Slipids Force Field to Polyunsaturated Lipids Inna Ermilova and Alexander P. Lyubartsev∗ Department of Materials and Environmental Chemistry, Stockholm University, SE 106 91, Stockholm, Sweden tel. +46-8161193 E-mail: [email protected]

Abstract 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, 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 centre parallel to the membrane ∗

To whom correspondence should be addressed

1

ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

plane. Furthermore, a tendency of cholesterol molecules to form chain-like clusters in polyunsaturated bilayers was qualitatively observed.

Introduction Lipids are important components of biological systems such as biomembranes. 1,2 A lipid bilayer forming a membrane helps to control the transport of different molecules into or out of a cell. 3,4 Furthermore, the membrane itself is involved in many important processes such as cell division, cell signalling, protein anchoring, cell fusion which are happening on a microscopic level. 5,6 Lipids composed of polyunsaturated fatty acids (PUFAs) are very common in mammalian cells including humans. 7,8 In particular, they are highly abundant in the eye retina and in the brain grey and white matter as one of the side chains of phosphatidylethanolamine (PE) and phosphatidylserine (PS). 7 In addition, linoelic acid is a dominant substituent of cardiolipin, 7 which is found in the heart and in the inner leaflet of mitochondrial membranes. The amount of PUFAs in different cells plays a big role in the human’s health. Some nutritionists suggest consumption of products which are rich with ω-3 and ω-6 fatty acids in order to avoid issues with the health 9 (the terms ω-3 and ω-6 mean polyunsaturated fatty acids which have the first double bond C=C at the third or sixth carbon atom from the end of the chain respectively). It was also discovered that PUFAs are useful for the treatment of mental diseases and memory problems, 10–19 heart diseases and problems with a high blood pressure 20–25 , atherosclerosis, dyslipidemia, inflammation, 26 cancer, 27–30 Parkinson, 31 problems with the vision 32 , diabetes and fatty liver disease. 33,34 The significant positive effects of the consumption of PUFAs were discovered on children’s performance in learning. 35 At the same time some harmful effects of the consumption of PUFAs were discovered: prohemorrhagic effect, negative effects on immune functions, formation of oxidative products by ω-3 polyunsaturated fatty acids. 36 Stillwell et al. 37 have shown that docosahexaenoic acid is not only affecting physical-chemical properties of membranes (such as thickness, lipid pack-

2


Page 2 of 49

Page 3 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


ing, fluidity, elasticity, permeability, etc), but it is activating the important signalling protein kinase C (PKC) as well. The activity of this protein was higher in a failed human heart compared to a healthy one according to Bowling et al. 38 The fact that PUFA lipids affect so wide variety of biological processes is evidently linked to their specific chemical structure with methylene-interrupted cis double bonds, which causes their specific functioning in living organisms. 39 Molecular dynamics (MD) simulations can provide valuable information for understanding of specific features of PUFA lipids, and through the years a number of works on atomistic simulations of bilayers composed of PUFA lipids have been carried out. 40–48 Still, the amount of modelling studies on PUFA is much less than what is available for fully saturated lipids. Partially this can be related to the lack of reliable parametrization of the interactions in the polyunsaturated hydrocarbon chains. Reliability of results obtained by MD simulations depends crucially on the force field (FF) used. In order to get reliable results the force field should accurately describe all interactions in the system of interest. For lipid bilayers, there is a relatively large variety of force fields. Often MD simulations of lipid bilayers are carried out using united atom Gromos FF and its variations: Berger lipids, 49 43A1-S3 FF parameter set 50 or G53A6L parameter set. 51 Use of united-atom force fields is less suitable for modelling of unsaturated lipids because of a weak polarity at the double bonds which is not included into united atom FFs having zero charges of united atom groups. There are also all-atomistic FFs with all hydrogens explicitly included, which were developed for modelling of lipids: several versions of CHARMM, 52–55 AMBER family (GAFF, lipid11, lipid14), 56–59 GLYCAM06, 60 Slipids, 61–63 see also recent review 64 and references therein. All these force fields contain parameters for a double C=C cis bond, however they were derived and validated for an isolated double bond which, for example, is present in such lipids as DOPC or POPC. However application of parameters derived for a single double bond to polyunsaturated lipids having several methylene-separated cis double bonds often leads to inaccurate descriptions of the surface area per lipid, deuterium order parameters, and X-ray

3



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

form factors. 65 To our knowledge, the only atomistic force field which includes parameters specially developed for PUFA chains is a recent modification of Charmm36 65 which was made by adjustment of the torsion potential between two methylene-separated cis double bonds. Apart from that, there exists a united atom model based on Berger force field 66 and a coarse-grained model based on MARTINI force field 67 which include special parameters for PUFA chains. In this work we extend the Slipids force field to polyunsaturated lipids. The original all-atomistic Slipids FF 61–63 was derived for saturated and monounsaturated lipids with a number of different headgroups and validated by comparisons of computed properties of lipid bilayers with available experimental data: average areas per lipid, bilayer thickness, X-ray and neutron scattering form factors, NMR deuterium order parameters. Slipids force field demonstrates also a good description of the temperature dependencies of experimentally observed bilayer properties including transition to the gel phase, and partitioning of various molecular inclusions across the bilayer. 68–70 In the extension of Slipids to polyunsaturated lipids we adopt the methodology used in the previous works 61–63 based on high-level abinitio computations of partial charges with multiconfigurational averaging, and subsequent parametrization of the torsion potentials. The developed force field is validated by simulations of several phosphatidylcholine (PC) lipids containing different numbers of double bonds in unsaturated tails: 18:1(n-9), 18:2(n-9), 18:3(n-3), 20:4(n-6), 20:5(n-3), 22:5(n-6), 22:6(n-3), and comparisons with available experimental data. Here and below we use numerical notations for the phospholipid chains in the form N : k(n − j) where N denotes the total number of carbon atoms in the chains, k is the number of double bonds and j is the position of the first double bond counted from the methyl terminus of the chain, cis-conformation of double bonds is implied in all cases. Through the text, notation (n − j) in the descriptions of lipids is often omitted for brevity. An example of 18:0-22:6(n-3) PC (phosphatidylcholine) lipid is shown in Fig. 1. Furthermore, the behaviour of bilayers containing highly unsaturated 20:4-20:4 PC and 22:6-22:6

4


Page 4 of 49

Page 5 of 49


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Etors =

6 X

kn (1 + cos(nφ − δn ))

Page 6 of 49

(3)

n=1

The Lennard-Jones and electrostatic interactions are not computed for atoms pairs connected by 1 or 2 covalent bonds. For atom pairs connected by exactly 3 bonds (the so-called 1-4 neighbours), special rules are applied. In the case of Slipids FF, 1-4 electrostatic interactions are scaled by factor 0.83, while special rules (defined by the FF atom types) are applied for 1-4 Lennard-Jones interactions. 61 To optimize FF parameters for atom types related to repeated cis double bonds we used two model molecules: cis-3,cis-6-nonadiene and cis-3,cis-6-dodecadiene (Fig. S1 and S2 of the Supporting Information). The parametrization algorithm involves several steps: 1. We run a MD simulation of model molecules in a liquid phase and extract from it randomly 50 molecular conformations 2. For the chosen conformations we compute partial atomic charges and for each atom type take average over the conformations. 3. We check suitability of the Lennard-Jones parameters by calculating thermodynamic properties of the liquid phase of model molecules and optimize them if necessary 4. After that we optimize the torsional potentials by ab-initio computations of energy while rotating the torsional angle. The procedure can then be reiterated with MD simulations of steps 1)-4) run with the optimized parameters. From the experience of previous works 61,62 it follows that selfconsistency is reached already after the first iteration, so in this work only one iteration was performed. After optimization of the parameters, validation simulations are run for a number of bilayers and computed properties are compared with available experimental data for the respective systems.

6


Page 7 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Computations of partial charges The initial MD simulations of each of the two dienes were performed using the original Slipids FF 61,62 at temperature 30◦ C, in which parameters for monounsaturated chains were used to describe the polyunsaturated part of the molecules. After 1 ns equilibration, 50 molecular conformations were randomly chosen from the simulated trajectory. For each selected molecular conformation, density functional theory (DFT) single point calculations were carried out with the purpose to extract partial atom charges. In these calculations we used the B3LYP exchange-correlation functional 71–73 with cc-pVTZ basis set. 74 As in the previous works we assumed that the molecules were placed in a polarizable continuum with IEFPCM 75,76 as a solvent model with ǫ = 2.04 which correspond to the typical dielectric permittivity of a hydrophobic media. 61 For calculations of partial charges we used the restrained electrostatic potential (RESP) approach 77 as implemented in the R.E.D. software, 78 in which the difference between the electrostatic potential determined from the electron density and the classical point charges was minimized, while maintaining the total sum of squares of atomic charges possibly low. The data on computed atomic charges for all molecular structures and their conformational averaging are presented in the Supporting Information, Tables S1 - S3. The computed average values of the atomic charges were further averaged over symmetrical atoms and were used as a guide, taking also into account partial charges of the previous version of the Slipid force field, 62 in the final assignment of the partial charges for typical PUFA chains as shown in Figure 2. It is remarkable that charges for atoms at double bonds turned out to be larger for PUFA chains (-0.28/0.14) compared to charges for monounsaturated lipids in the Slipids force field (-0.2/0.11) 61 which were derived using the same methodology. Stronger polarity of PUFA chains, which follows from our quantum-chemical calculations, may result in a weaker hydrophobicity of the interior of PUFA bilayers and affect partitioning of various molecules in the bilayers. After optimization of the partial atom charges, MD simulations of several dienes were carried out in the NPT ensemble at 303 K for 10 ns and the computed density was compared 7




Page 8 of 49

Page 9 of 49




1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 49

2,5-heptadiene. They found that a direct fit to this quantum-mechanical energy surface resulted in too high chain order, especially for SDPC bilayer, and in the final expression for the torsion term they increased the well breadth and reduced the transitional barrier between the minima, to get better agreement with NMR order parameters. In our work, we have followed the approach adopted in derivation of torsion parameters for Slipids FF for saturated 61 and monounsaturated 62 chains, using one-dimensional rotation around torsion angle of interest. First, the geometries of model dienes (cis-3,cis-6-nonadiene and cis-3,cis6-dodecadiene) were optimized using the second order Møller-Plesset perturbation theory (MP2) 82 with correlation-consistent cc-pVDZ basis set, 74 and then molecular conformations were rotated around the torsion angle of interest with the step 10◦ keeping other degrees of freedom frozen. In order to ensure that fixing of the other degrees of freedom while varying torsion angle does not bring major artefacts, we have carried out additional series of computations in which, for each fixed value of the torsion angle of interest, we performed geometry optimization over all other degrees of freedom of the molecule. Results for the quantum energies computed for the fixed and optimized conformations obtained at different values of torsional angle α1 are shown in Fig. S3 of the Supporting Information. The transition barrier over the trans-conformation (180◦ ) shows only a marginal decrease by 0.5 kJ/mol when optimization of molecular structure at each value of α1 angle was made. Somewhat larger difference in energy between the fixed and optimized conformations (within 2 kJ/mol) is seen in the range of α1 angle 110 − 140◦ , and just in this range of α1 the adjacent angle α2 deviates from 110 − 120◦ to 80 − 90◦ (Fig. S4). These results show that fixing of other degrees of freedom while rotating the torsional angle of interest does not bring major deviations from the optimized transition path. The torsion potential was fitted by the form of eq. (3) to the energies obtained by:

Etorsional = (ECCSD(T )/BBS − EM D ) − (ECCSD(T )/BBS − EM D )min

10


(5)

Page 11 of 49




1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 49

region may deteriorate fitting in the most populated low energy area, which is much more important for the correct reproduction of the conformational behaviour. Table 1: Parameters for dihedrals CEL1=CEL1-CTL2-CEL1 (α) and CEL1=CEL1-CTL2CTL2 (β) nα 1 2 3 4 5 6

δα , degrees 0.00 0.00 0.00 0.00 0.00 0.00

kα . kJ/mol 4.331 5.161 -0.415 3.135 1.59 0.137

nβ 1 2 3 4 5 6

δβ , degrees 0.00 0.00 0.00 0.00 0.00 0.00

kβ , kJ/mol 3.967 3.89 -0.138 1.892 0.421 0.421

We have computed distributions of dihedrals α and β for the model molecules by MD simulations at T=303K, the results are displayed in Figure 5. Two probability maxima are seen at values of dihedrals at around -120◦ and 120◦ for both types of dihedrals. Similar results for the distribution of dihedrals in polyunsaturated chains were achieved in works of Saiz and Klein 83 and Klauda et al. 65 Figure 5 shows that dihedrals which were excluded from the optimization of the potential for torsional angles (around cis-conformation) never appear during simulations. The region of conformations around 180◦ is weakly populated showing that these dihedrals can change conformations from ”right” to ”left” and vise versa on the time scale of the simulations.

Simulation Details All molecular dynamics simulations were done with Gromacs-4.6.7 84 software using the leapfrog integrator with the time step 2 fs. Covalent bonds were constrained by the LINCS algorithm 85 with exception for water for which bonds were constrained by the SETTLE method. 86 Model molecules simulations were carried out in a cubic periodic box containing 340 molecules in the isotropic NPT ensemble at 1 bar pressure. Bilayer simulations were carried out in the N P T ensemble with semianizotropic Berendsen barostat 87 keeping the constant pressure at 1 bar separately in Z-direction and in the XY-plane. The barostat 12


Page 13 of 49




1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 49

Table 2: Simulated systems details Lipid, PC 18:0-18:1(n-9) 16:0-18:2(n-9) 16:0-22:6(n-3) 18:0-18:3(n-3) 18:0-20:4(n-6) 18:0-20:5(n-3) 18:0-22:5(n-6) 18:0-22:6(n-3) 20:4(n-6)-20:4(n-6) 22:6(n-3)-22:6(n-3) 20:4(n-6)-20:4(n-6) 22:6(n-3)-22:6(n-3)

Name SOPC PiLPC PDPC SAPC

SDPC DAPC DDPC DAPC-CHOL DDPC-CHOL

Nlip 128 128 128 128 128 128 128 128 128 128 200 200

Nw 5120 3840 3840 3840 3840 3840 3840 3840 3840 3840 10000 10000

Nchol 0 0 0 0 0 0 0 0 0 0 100 100

T (K) 303 303, 293, 303, 313 303 303 303 303 303 303 293 293 293 293

Time, µs 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 2.4 2.4

All quantum chemical computations were carried out using Gaussian09 package 92

Results and discussion Volumetric parameters of the bilayers Average area per lipid is a commonly used fast test to evaluate quality of the force field. The problem of such validation is that different experimental methods use different model assumptions which may lead for a large spread of experimental average areas for the same lipid. For example for 18:0-18:1 PC experimental values of the average area per lipid were found in different works in the range between 61.4 and 71 ˚ A2 . 64 Neutron and X-ray scattering data are considered as the most reliable sources of experimental determination of the average areas, however they are not available for many lipid species. Still, due to its simplicity a comparison of the simulated and experimental average lipid area becomes a necessary element of the analysis, aiming also to determine whether equilibration of the bilayer was reached or not. The average areas per lipids were computed for all bilayers simulated in this work and they are given in Table 3 in comparison with experimental results. Evolution of the area

14


Page 15 of 49




1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 49

Table 3: Average areas per lipids (T = 303K) Lipid, PC 18:0-18:1

˚2 Computed, A 66.7 ± 0.5

16:0-18:2 18:0-18:3 18:0-20:4 18:0-20:5 18:0-22:5 16:0-22:6 18:0-22:6

67.2 ± 0.5 68.9 ± 0.6 69.6 ± 0.6 70. ± 0.6 67.1 ± 0.6 68.8 ± 0.6 68.6 ± 0.7

˚2 Experimental area/lipid, A 67 (low angle X-ray scattering) 93 65.5 ± 1.3 (small angle neutron scattering) 94 66 (Langmuir film balance) 95 66.6 (NMR) 96 70.6 (NMR) 96 69.1 (NMR) 96 68.7 ± 0.4 (X-ray) 40 70. (Langmuir film balance) 95 68.2 ± 0.4 (X-ray) 40 71.6 (NMR) 96

Table 4: Distance between head-groups DHH , hydrophobic DC and Luzzati DB thicknesses Lipid, PC 18:0-18:1 16:0-22:6 18:0-22:5 18:0-22:6

A A DHH exp , ˚ DHH comp , ˚ 93 38.5± 0.8 39.2 37.5± 0.8 38.1 98,99 38.1± 0.8 37.9 ± 1 40 38.1± 1 37.9 ± 1 40

˚ ˚ 2DC exp , A 2DC comp , A 97 28.50 29.9 29.00 28 99 31.00 30.5 ± 1 40 33.00 30.5 ± 1 40

A DB comp , ˚ 41.00 44.00 46.00 46.00

A DB exp , ˚ 97 40.00 no data 48.5 ± 1 40 48.5 ± 1 40

can be determined from the electron density profiles. 97,100 We collected simulated and available experimental data for 18:0-18:1 PC, 16:0-22:6 PC, 18:0-22:5 PC and 18:0-22:6 PC lipids in Table 4. All volumetric data for the simulated bilayers were found in agreement with experimental results within computational and experimental uncertainty, in cases when such data were available.

Deuterium order parameters A profile of deuterium order parameters measurable in NMR experiments is a property widely used in validation of force fields for lipid bilayer simulation. The order parameter for a specific CH (or CD) bond is determined by 1 |SCD | = | h3cos2 θ − 1i| 2 where θ is the angle between direction of the bond and the bilayer normal.

16


(6)

Page 17 of 49




1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

We have computed order parameters for saturated and unsaturated lipid tails of the studied bilayers using formula (6). Figure 7 shows comparison of the results for 18:0-22:6(n-3) PC, 18:0-22:5(n-6) PC, and 16:0-22:6(n-3) PC obtained in our simulations with experimental data from several sources. 40,101,103 For bilayers consisting of 18:0-18:3(n-9) PC, 18:0-20:4(n-6) PC and 18:0-20:5(n-3) PC simulated order parameters are given in the Supporting Information, Figures S8 - S10, together with experimental data for sn-1 chains from Holte et al. 104 Statistical uncertainty of the simulated order parameters is always within 0.01. In all cases simulations show overall good agreement with experimental order parameters. In all simulated bilayers, order parameters for saturated sn-1 chains are about 0.2 in the beginning of the chain and gradually decreases toward the chain end. One can note slightly higher order parameters in the beginning of sn-1 chain for 18:0-22:5 PC than for other lipids both in simulation and experimental data. The difference between order parameters for the sn-1 chain of the 18:0-22:5 PC and 18:0-22:6 PC bilayers, displayed in Figure S11 of the Supporting Information, demonstrates remarkable similarity between computed and experimental data of Eldho et al, 40 with reproduction of the maximum of the difference in the middle of the chain. For unsaturated tails the order parameters are generally lower. For most of carbons (except C4 and C5 in 18:0-22.5(n-6) PC), the difference between simulations and experimental results is within the same interval (about 0.02) as the difference between two experiments in case of 18:0-22:6(n-3) bilayer. Simulated order parameters for C2 carbon of sn-2 chains show splitting for the two hydrogens, more noticeable for the most unsaturated PUFA (16:0-22:6, 18:0-22:6, 18:0-20:4 and 18:0-20:5). Splitting of deuterium order parameters at C2 carbon was also observed in some experiments for other lipids. 105 For 18:0-22:5(n-6) PC, simulations overestimate the order parameter for C4 and C5 carbons of the unsaturated chain by 0.05 - 0.07, but also reproduce the experimentally observed minimum of the order parameter at around 15-th carbon and its increase toward the end of the chain. For 16:0-18:2(n-9) PC bilayer we have carried out simulations in the range of temperatures

18


Page 18 of 49

Page 19 of 49




1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 49

shown previously that the low order parameters in the area of double bonds are mainly related to the specific orientations of the CH bonds at the angle to the bilayer normal close to the magic angle of 54◦ , 45,48 and such behaviour was also observed in our simulations.

X-ray form factors X-ray scattering form factors are considered as the most direct test of the quality of the FF for lipid bilayer simulation since experimental form factors are not subjected to model-based interpretations. 113 In simulations, X-ray scattering factor can be straightforwardly computed by the Fourier transform of the electron density: 114

S(q) =

Z

D/2

(ρ(z) − ρw ) cos(qz)dz

(7)

−D/2

where ρ(z) is the electron density along z-coordinate (computed in simulations taking into account partial charges of atoms), ρw is the electron density of bulk water and D is size of the simulation box in Z-direction, which corresponds to the D-spacing (the distance between parallelly oriented bilayers in the sample) of bilayers in scattering experiments. The scattering factor for 18:0-18:1 PC lipid bilayer, computed using eq. (7) with integration step 0.2 ˚ A, is displayed in Figure 9 together with experimental form factor from paper by Kuˇcerka et al. 97 This result (red line in the figure) essentially repeats the simulation for SOPC lipid with the previous version of Slipids FF (blue line), 62 and taken together with results for the average lipid area, volumetric properties and order parameters for this lipid shows that update of the torsion potential for dihedral next to a double bond (β in Figure 3) does not affect the computed properties for monounsaturated lipids. For polyunsaturated lipids, we found in the literature scattering data only for 18:0-22:6 PC bilayer in the work by Eldho et el. 40 We plot these data, together with our results obtained by the Fourier transform of the electron density, in Figure 10. The experimental results 40 combines data from three different experiments which were performed at three

20


Page 21 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 9: Scattering form factors for 18:0-18:1 PC bilayer. ”comp (1)”: computed with the FF of this work; ”comp (2)’: computed with the previous version of Slipids FF; 62 ”exp”: experimental X-ray form factor 97 different values of D-spacing. Experimental points are placed close to the simulated ones that gives additional arguments that our model provides realistic description of PUFA lipids.

Figure 10: Computed (red line) and experimental X-ray scattering form factors for 18:022:6 PC bilayer. The experimental points are taken from paper 40 for different D-spacings. D-spacing in the simulations is 66.5 ˚ A.

21



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 49

Reorientational dynamics Experimentally, dynamical properties of lipid bilayers can be probed by NMR relaxation experiments. Particularly,

13

C spin-lattice relaxation time T1 can provide information on

the reorientational dynamics of CH-bonds, which happens on pico- and nano-second time scale and which is also accessible by molecular dynamics simulations. 40 T1 time is related to the spectral density function J(ω) 40 : 1 1 = N T1 10

h 1 γC γH 3 2π rCH

2

× (J(ωC − ωH ) + 3J(ωC ) + 6J(ωC + ωH ))

(8)

where N is the number of protons bound to the considered carbon atom, ωC and ωH are resonance frequencies of C and H nuclei, γC and γH are gyromagnetic ratios respectively, and rCH is the length of the CH bond. In computer simulations, the spectral density J(ω) can be calculated from the reorienational time correlation function C(t) = hP2 (~µ(0)~µ(t))i of the CH unit vector µ ~ (t) (P2 being the second Legendre polynomial):

J(ω) =

Z

∞

C(t)cos(ωt)dt

(9)

0

In the limit of motional narrowing (suggesting that τ ω

Extension of the Slipids Force Field to Polyunsaturated Lipids - The

Recommend Documents