Coarse-grained bond and angle distributions from atomistic simulations: On the systematic parameterisation of lipid models

Accepted Manuscript Title: Coarse-grained bond and angle distributions from atomistic simulations: on the systematic parameterisation of lipid models Author: Samuel Genheden PII: DOI: Reference:

S1093-3263(15)30084-X http://dx.doi.org/doi:10.1016/j.jmgm.2015.11.009 JMG 6633

To appear in:

Journal of Molecular Graphics and Modelling

Received date: Revised date: Accepted date:

23-9-2015 6-11-2015 9-11-2015

Please cite this article as: Samuel Genheden, Coarse-grained bond and angle distributions from atomistic simulations: on the systematic parameterisation of lipid models, Journal of Molecular Graphics and Modelling http://dx.doi.org/10.1016/j.jmgm.2015.11.009 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Coarse-grained bond and angle distributions from atomistic simulations: on the systematic parameterisation of lipid models Samuel Genheden Department of Chemistry and Molecular Biology, University of Gothenburg, 405 30, Göteborg, Sweden [email protected]

1

Graphical abstract

2

Highlights 

Atomistic lipid force fields sample similar bond and angle distributions



Coarse-grained lipid force fields sample distributions that differ from the atomistic distributions



It is possible to re-parameterise the coarse-grained force fields by employing the atomistic distributions



The re-parameterised Elba potential gives a membrane that is closer in structure to the experimental membrane



The re-parameterised Martini potential gives a membrane that is too stiff, most likely because of a inappropriate mapping

3

Abstract Coarse-grained (CG) models are popular alternatives to atomistic (AT) force fields as they enable simulations of larger systems at longer timescale. The bottom-up approach is a systematic parameterisation strategy whereby data from AT simulations are used to determine the CG parameters. This is particular straightforward with the bond and angle parameters as a direct Boltzmann inversion can be used. Still, a reference AT force field has to be chosen. In this study, I compare three common AT force fields (Stockholm lipids, Berger and Gromos) and investigate the sampling of bond and angle distributions in two CG models (Martini and Elba). As a test case, I choose a bilayer of POPC lipids. The AT simulations give distributions that agree to a large extent, especially in the fatty acid tails. However, the AT simulations sample distributions that differ from the distributions observed in CG simulations with respect to both location and width. The bond and angle distributions from the AT simulations are then used to re-parameterise the CG force fields. For the Martini model, this significantly alters the physical behaviour of the membrane, which likely is an effect of the mapping. However, for the Elba model the re-parameterised force field gives a membrane that is in some respects closer to the experimental membrane. Implications for CG parameterisation are discussed. Keywords: coarse-graining, force fields, parameterisation, membrane simulations, Martini, Elba.

4

Introduction Molecular dynamic simulations are attractive alternatives to traditional wet-lab experiments as the simulations provide dynamic insight on (bio)chemical systems at an atomistic (AT) resolution.[1] Typically, the atoms are propagated using classical mechanics guided by an empirical potential, viz. a force field. The most accurate force fields model each atom individually but they become computationally expensive when the aim is to study large systems at biological timescales. Therefore, coarsegrained (CG) models have been developed that group atoms into pseudo-atoms or beads, thereby reducing the number of particles that needs to be propagated.[2,3] The least invasive strategy is simply to remove non-polar hydrogen atoms and merge them with the atom they are bonded to, leading to a united-atom force field.[4] A more drastic reduction in the number of particles can be achieved by merging entire groups of atoms into beads. This has been a successful strategy in popular force fields such as Martini[5-7] TraPPE,[8] and Elba.[9,10] CG models can be parameterised using several strategies.[11] For instance, matching to macroscopic data from experiments is a common approach for biological CG force fields like Martini and Elba.[5,10] Unfortunately, this procedure is typically ad-hoc in nature, rather than systematic, which is unhelpful if new molecules have to be parameterised by researchers not involved in the original parameterisation. This has also traditionally been the case for AT force fields although recent tools to make unbiased parameterisations have been suggested.[12] A more systematic set of approaches sort under the term bottom-up,[13] whereby data from AT simulations are used to derive the CG parameterisation. Several such strategies have been suggested in the literature.[14-18] The most common strategies are force-matching or reconstitution of the pair distribution function.[11,13] The latter can be accomplished with for instance Boltzmann inversion,[19] inverse Monte Carlo,[15] or Iterative Boltzmann Inversion.[16] Although automatic tools to parameterise CG force fields have been suggested,[20-23] they are still less mature and widely used than corresponding tools for AT force fields.[24-26] A critical, and one of the first choices in a systematic bottom-up parameterisation approach is the selection of an AT force field to be used as reference. For biological systems there exists a wide range of force fields [27-30] and it is uncertain what reference to use. A natural question then arises: do the different AT force fields give the same CG parameters? What is the sensitivity of the CG model with respect to the AT parameters? Sensitivity issues have been investigated for some time for AT models[31,32] and more recently in CG modelling.[33] In this paper, I will address this question for the small, albeit illustrative case of the parameterisation of bond and angle parameters in a common lipid, POPC (1-palmitoyl-2-oleoyl-snglycero-3-phosphocholine). First, there are CG parameters available in both Martini and Elba for this lipid and therefore a direct comparison with AT force fields can be made. Second, by re-parameterising the available CG force fields using data from AT simulations, the effect on the physical properties can be observed. I will conclude the paper by discussing future prospects on parameterisation of lipid force fields using a bottom-up approach. Methods System setup. A membrane consisting of 128 POPC and 5120 water molecules were created using the CHARMM Membrane builder web service.[34] This membrane was then simulated with the Stockholm lipids (Slipid),[35] Berger,[36] Gromos,[37] Martini[6] and Elba[10] force fields. The membrane was minimized using 100 steps of steepest

5

descent followed by 10 ns equilibration at 303 K and 1 bar. Finally, a 250 ns production simulation was run at the same conditions. Two independent repeats were initiated by using different initial coordinates of the system. The simulations are detailed below for each force field. All simulations were run with Gromacs version 5.0.4[38] except the Elba simulations that were run with Lammps.[39] Slipid simulations. The time step was set to 2 fs and the Lincs algorithm[40] was used to constrain all bonds. The non-bonded neighbour list was setup with the Verlet scheme[41] and updated every 10th step. Electrostatic forces were treated using particle-mesh Ewald (PME)[42] and the real-space cut-off was 1.2 nm. The van der Waals forces were switched off from 1.0 to 1.2 nm and long-range corrections were added. The 10 ns equilibration simulation was run employing a weak-coupling thermostat[43] with a 0.5 ps coupling constant and a weak-coupling semi-isotropic barostat with a 10 ps coupling constant. The water and lipid molecules were coupled to independent thermostats. The 250 ns production simulation was run with a Noose– Hover thermostat[44,45] and a Parinello–Rahman barostat[46] both using the same coupling constants as in the equilibration simulation. Berger simulations. Non-polar hydrogen atoms were removed from the lipids using in-house scripts. The time step was set to 2 fs and all bonds were constraint using the Lincs algorithm.[40] The non-bonded neighbour list was setup with the group scheme[41] and updated every 5th step. Electrostatic forces were treated using PME[42] with a 1.0 nm real-space cut-off. The van der Waals forces were cut-off at 1.0 nm and long-range corrections were added. The thermostat and barostat settings were as in the Slipid simulations. Gromos simulations. Non-polar hydrogen atoms were removed from the lipids using in-house scripts. The time step was set to 2 fs and all bonds were constraint using the Lincs algorithm. [40]The non-bonded neighbour list was setup with the group scheme[41] and updated every 2nd step. The electrostatic forces were treated with a generalized reaction field[47] using a 1.4 nm cut-off. The van der Waals forces were cut-off at 1.4 nm and no long-range corrections were added. The thermostat and barostat settings were as in the Slipid simulations. Martini simulations. The POPC molecules were mapped to the CG representation and groups of four water molecules were replaced by a CG water bead using in-house scripts. The time step was set to 40 fs and no bonds were constrained. The nonbonded neighbour list was setup with the Verlet scheme[41] and updated every 10th step. Electrostatic forces were treated using a reaction field[47] and a potential shift from 0 to 1.2 nm. The van der Waals forces were switched off from 0.9 to 1.2 nm. The equilibration simulation was run with a weak-coupling thermostat[43] with a 0.5 ps coupling constant and a weak-coupling semi-isotropic barostat with a 10 ps coupling constant. The production simulation was performed with a velocity rescaling thermostat[48] and a Parinello–Rahman barostat[46] using the same coupling constants as in the equilibration. Elba simulations. The POPC and water molecules were mapped to the CG representation using in-house scripts. The time step was set to 10 fs and the nonbonded neighbour list was updated every 2nd step. The non-bonded forces, which are

6

of a shifted-force kind, were cut-off at 1.2 nm and no long-range corrections were applied.[10] Both the equilibration and production used a Langevin thermostat[49] with a 5 ps coupling constant and a weak-coupling semi-isotropic barostat with a 10 ps coupling constant. AT to CG mapping. The mapping from atoms to CG beads is a non-unique procedure and there was no explicit mapping described in the original publications. Therefore, a mapping recently used by Daily et al.[50] was used for Martini and a novel mapping for Elba is suggested herein. These mappings are illustrated in Figure 1. For Martini, the mapping is based on the 4 to 1 policy,[6] whereas Elba uses a mapping of functional groups combined with a 3 to 1 mapping for the tail beads. To construct CG coordinates from AT coordinates, the following transformation was applied

RCG = TRAT

1)

where R is the Cartesian coordinates (either CG or AT) and T is a transformation matrix with NCG ´ N AT elements (N is the number of beads or atoms). Each element, Tij is equal to zero if atom j does not map to bead i, or 1/m if atom j does map to bead i where m is the total number of atoms mapping to bead i. The weight is not adjusted if an atom is mapped to more than one CG bead. Analysis. The area and volume per lipid were computed from the instantaneous area of the membrane plane and the volume of the simulation box, assuming experimental density of the water molecules. The thickness of the membrane (DHH) was taken as the difference between the peaks in the density of the phosphate atoms/beads.[51] Lipid tail order parameters, SCC were calculated for bead–bead CG bonds according to

SCC 

1 3 cos 2   1 2

2)

where θ is the angle between three consecutive tail beads and the brackets indicate ensemble averages. Below, the average over all SCC in a tail is presented. All analyses were performed on the last 100 ns of a simulation. Re-parameterisation. The bond and angle potentials of the Martini and Elba force fields were re-parameterised using distributions from the simulation with the Slipid force field. A direct Boltzmann inversion was used, assuming that 1) the bonded degrees of freedom are independent and 2) the sampled distribution observed in the AT simulations follow a unimodal Gaussian distribution. The CG potential, U is directly computable form the AT distribution function p[11]

U(x) = -RT ln p(x)

3)

where x is a bond length or the cosine of a valance angle, R is the gas-constant and T is the absolute temperature. The functional form of U is harmonic as given by the force field and p is assumed to be Gaussian distributed. Therefore, the fitting stage can be skipped altogether as it is trivial to see that the equilibrium, veq (either a bond length or the cosine of an angle), is equal to the sample mean

veq = m

4)

7

and the force constant, k, can be calculated from the sample variance, σ2

k = RT / s 2

5)

Simulations with the new CG force field parameters were performed identical to the simulations with the original parameters as discussed above. Results Membrane properties. The AT force fields (Slipids, Berger and Gromos) accurately reproduce the area (AL), volume (VL), and thickness (DHH) of the experiments as can be seen in Table 1. The deviations from the experiments are as expected from earlier simulations with these force fields.[51,52] The CG force fields (Martini and Elba) on the other hand show discrepancies that largely can be attributed to mapping mismatch. Martini, which uses a 4 to 1 mapping, leads to too long lipids and a too thick membrane, whereas Elba, which uses a 3 to 1 mapping, leads to short lipids and a too thin membrane. For Elba, which was parameterised for 18-carbon chain lipids, the 16carbon palmitoyl-chain in POPC cannot be modelled using equal bead sizes. These discrepancies have been discussed previously in the literature.[52] To further quantify the simulated membranes, tail order parameters were calculated based on bead–bead bonds in CG tails. The AT molecules were mapped on either a Martini or Elba lipid model. The average tail order parameters are listed in Table 2. Starting with the Martini mapping, the AT force fields are consistent and show order parameters between 0.40 and 0.43 for the sn-1 chain and between 0.36 to 0.37 for the sn–2 chain. The Martini simulation displays a slightly more flexible lipid. Next, with the Elba mapping, the order parameters are general lower than with the Martini mapping, which probably can be attributed to the difference in bead size between Elba and Martini. Still the AT force fields are consistent and the Elba simulation display lower order. However, the difference between the Elba and AT simulations is considerably larger than the difference between the Martini and AT simulations. Distributions with the Martini mapping. In Figure 2, the distributions of the twelve bonds defined for the POPC molecule in the Martini force field are plotted for simulations with the Slipid, Berger, Gromos or Martini force field. As a reference, the distribution given by the force constant and equilibrium length in the Martini force field is also given (labelled ‘ff’). All of the sampled distributions follow a normal distribution as determined by a Shapiro–Wilk test (see Supplementary Material). There are minor discrepancies between the Martini simulations and the ideal force field distribution, which is expected, because the bonds are not simulated in isolation. A general trend is that Martini samples a much wider distribution than the AT force field; the mean absolute deviation (MAD) of the force constant compared to Martini ranges between 1662 and 2388 kJ / (mol nm2). The deviation in equilibrium length between the Martini force field and the AT models is up to 0.2 nm (for the GL1–C1B bond). The Gromos and Berger distributions are more similar to each other, than to Slipid, which is also expected, as Gromos and Berger are united-atom models. The largest difference between the Gromos and Berger can be seen in the head group region (NC3–PO4 and GL1–GL2 bonds).

8

In Figure 3, the distributions of the nine valance angles defined for POPC are shown. As with the bond distributions, the sampled angle distributions follow a normal distribution. In general, Martini samples valance angles that are smaller than the ideal force field behaviour, this is especially pronounced in the tail regions. As with the bond distributions, the AT force fields tend to display a narrower distribution than the Martini force field, although the effect is much less than what was observed for the bonds; the MAD of the force constant is between 8 and 18 kJ/mol. The largest deviations between the differences force fields are again seen in the head group region, but Berger and Gromos tend to sample very similar distributions. Overall, the agreement in equilibrium angles is fairly consistent between the different force fields. Distributions with the Elba mapping. In Figure 4 the distributions of the 13 bonds defined for the POPC molecule in the Elba force field are plotted, for simulations with either Elba or AT force fields. A Shapiro-Wilk test show that the sampled distributions follow a normal distribution (see Supplementary Material). The Elba simulations show in general shorter equilibrium bond lengths compared to the ideal force field distribution, except for the n–p bond. As with the Martini mapping, the CG simulation show a much wider distribution than the AT force field; the MAD for the force constant is between 2389 and 3705 kJ / (mol nm2) with individual differences up to 10 000 kJ / (mol nm2) for the t2a–t3a bond. The difference in equilibrium length between the AT force field and Elba is up to 0.1 nm, for the n–p bond but many of the individual differences are smaller. The AT force fields are fairly consistent; the largest difference is seen in the head groups. For instance, there is a 0.05 nm difference in equilibrium length between Slipid and Gromos/Berger for the g–e1 bond. In Figure 5 the distributions of the twelve angles are plotted, for simulation with either Elba or AT force fields. As with the bonds, the sampled distributions are normal-distributed. The actual simulation with Elba displays different equilibrium angles than the ideal force field as also seen with the Martini model. Many of the equilibrium angles are 180 degrees, but none of the simulated distributions are centred at this angle. There is also a tendency for many of the simulated distributions to be less wide than the ideal distributions, especially in the tail region. Overall, the agreement for the angle distributions when comparing to AT force fields is better than the equivalent comparisons for the bond distributions, the MAD for the force constant is between 32 and 66 kJ/mol. The largest variations among the force fields are seen around the glycerol beads, viz. the p–g–e1 and p–g–e2 angles. Berger and Gromos are generally more similar to each other than to Slipid. Re-parameterised force fields. Because the AT simulations sample different distributions than the CG dittos, it is of interest to see if the AT simulations can be used to re-parameterise the CG force fields. Naturally, one has to make a choice on the reference AT force field. Because the AT simulations largely agree except in the head groups, I here choose the Slipid simulations as the reference. Furthermore, as the bond and angle potentials in both Martini and Elba are harmonic, it is trivial to extract the force field parameters from the sample distribution as outlined in the Methods section. The new force field parameters are listed in Tables S1 and S2 in the supplementary material. In Table 3, various structural membrane properties are listed as observed in the simulations with re-parameterised CG force fields.

9

It is clear that the Martini membrane is greatly affected by the changes in the force fields. Both the AL and VL change drastically with 0.2 nm2 and 0.2 nm3, respectively, whereas DHH is unaffected. Furthermore, the lipid tail ordering is much stiffer, with increases of up to 0.44 All these effects are a consequence of the much shorter distances between the tail beads, which leads to a more ordered tail and a smaller membrane area. However, this does not affect the thickness of the membrane. This is also clear from the membrane properties if the bond or angle parameters are changed individually. Changing only the angle parameters have very little effect on any of the membrane properties, whereas changing only the bond parameters has a large effect. The effects on the Elba simulations are different. Overall the order parameters are increased slightly with ≈ 0.1, which leads to slightly less flexible lipids. This leads to a membrane with a 0.1 nm2 smaller area and a thicker membrane by 0.2 nm. Interestingly, this effect stems from both the re-parameterisation of the bond and angle parameters. Re-parameterisation of either bonds or angles gives a decrease in both AL and VL. However, the bond changes leads to a 0.1 nm thinner membrane, whereas the angle changes leads to a 0.4 nm thicker membrane. Discussion Developing systematic parameterisation methods is an advantageous pursuit: by using a systematic approach, automatic tools can be developed that ensures that the parameterisation of new molecules are consistent regardless of whoever the researcher is, hopefully leading to a situation where different molecules can be combined effortlessly in new and interesting systems, i.e. that the parameterisation is transferable. For the bottom-up parameterisation of new CG models, two critical choices have to be made: first, a mapping from AT to CG needs to be determined and second, a reference AT force field needs to be chosen[11]. In this paper, I am concerned with the latter choice, because the first one has to some extent already been discussed in the literature.[50,53] An investigation of the sensitivity of the parameterisation on the choice of AT force field requires a good test case. Here, I have chosen a simple POPC membrane. The POPC molecule is a common lipid to simulate due to its abundance in mammalian cell membranes,[54] and it is also interesting from a modelling perspective. Its two fatty acid tails, palmitoyl and oleoyl are naturally of different length, which is a challenge for CG models such a Martini or Elba that utilises a fixed bead size.[55] Furthermore, the two tails have different saturation, requiring different CG parameters, especially angle parameters. Therefore, despite its simplicity this choice of test case can illustrate subtle behaviours. Furthermore, because it is simple in characteristics, it should also be straightforward to parameterise. If the bottom-up approach despite this leads to unwanted effects, the validity of this approach on more complex molecules has to be questioned. As examples of two different CG mappings, I have chosen the Martini and Elba models, which are CG force fields used to simulate large biomolecular systems. They were both initially parameterised in a more or less ad hoc fashion by matching to macroscopic data from experiments, such as heat of vaporisation or area per lipid.[5,10] Whereas Martini includes parameters for an impressive range of lipid molecules, the Elba force field is less developed. Therefore, it is of interest to investigate if a systematic approach can shorten this gap, because the Elba model has several advantages over the Martini model, such as a more physical treatment of electrostatics and seamless mixing with AT force field in a hybrid scheme.[53]

10

The first question that needs to be addressed is whether the different AT force fields lead to different bond and angle distributions. If this were the case, it would be difficult to rationalise the selection of one AT force field before the other. Encouragingly, the different AT force fields agree to a large extent. For the bond distributions, the AT force fields agree both in location and width at least for the saturated fatty acid bonds (see Figures 2 and 4). There are minor differences around the double bond in the oleoyl chain, e.g. the C2–D3B bond in the Martini mapping and the t2a–t3a bond in the Elba mapping. The situation for the angle distribution is similar. However, the largest differences are seen in the head group and glycerol regions. For instance, with the Martini mapping, there is a difference of 0.1 nm in location between the Slipid and Berger distributions for the GL1–GL2 and PO4–GL1 bonds. Similar behaviour is seen with the Elba mapping. Naturally, this reflects the difficulties of mapping the glycerol backbone. Contrary to the fatty acid tails, the glycerol region is physically and chemically more complex. This is also seen in the different mapping approaches: Martini maps the glycerol region with two beads, whereas Elba uses three. More surprising is the observation that the distribution for the bond between the choline and phosphate bead differ between the AT force fields, since this is a more straightforward mapping. Because the AT force fields agree at least in the tails, the Slipid force fields were the reference in the subsequent reparameterisation, but it should be emphasized that either Berger or Gromos could have been chosen. Naturally, other AT force fields could have been used as references, e.g. Charmm[56] or GaffLipid.[57] However, because Slipid, Berger and Gromos agree to such a large extent, it seems unlikely that other AT force fields would alter the conclusions herein. In addition, Slipid borrows a majority of the bonded terms from the Charmm force field[35] and such is expected to sample similar bond and angle distributions. In a novel parameterisation campaign, one would directly use the AT distributions to parameterise the CG model and then perform a CG simulation that hopefully reproduces properties of interest (such as AL). Herein, a direct comparison between the AT and the CG bond and angle distributions can be made because POPC is parameterised in both Elba and Martini. Because the Elba and Martini models of the POPC membrane are fairly accurate in reproducing experimental observables (see Table 1), this comparison is in some respect like looking at the correct solution. However, it should be noted that these observables are macroscopic in nature and can be produced by many microscopic models. Looking at the distributions, a clear trend is that the bond distributions are much narrower in the AT simulations compared to the CG simulations (see Figures 2 and 4). This is probably an effect of the very stiff bonds between atoms in the AT simulations. The CG bonds on the other hand are weak and are present more as a mean to hold beads together more than to model a physical interaction. Furthermore, the CG force constant is usually equal for most of the bonds in the molecule. Still the relatively larger discrepancy in force constant could lead to a problem when using a bottom-up approach. The location of the distributions also displays some interesting results. For the Martini mapping, the location of the distribution for many of the tail beads differ significantly when comparing AT and CG simulations, which most likely is an effect of a mapping where an atom is allowed to be mapped to more than one bead (see Figure 1). For the Elba mapping, the differences between the AT and CG simulations are smaller. Looking instead at the angle distributions it is clear that they agree to a larger extent (see Figures 3 and 5). This is true for both the Martini and Elba mapping, with the exception of some angles in the head group or glycerol backbone. The actual CG

11

simulations seems in many instances to sample narrower angle distributions than dictated by the force field as in clear from both Figures 3 and 5, which indicate that the bonded terms are affected by non-bonded interactions. Because of the significant difference when comparing bond distributions with the Martini mapping, it is perhaps not surprising that using Slipid to re-parameterise the CG force field leads to drastic changes in the membrane. Clearly, the bond lengths observed in the AT simulations are too short to be used, leading to very stiff bonds (see Table 3) and a too low area per lipid. This is likely an effect of the mapping, so if a bottom-up approach should be used to re-parameterise Martini, a different mapping needs to be devised, although this is a non-trivial problem.[11] However, using the mapping to re-parameterise the angle parameters is successful as already demonstrated in the literature.[50] An alternative to changing the mapping is to adjust the non-bonded parameters. The Elba force field can, on the contrary be reparameterised using the Slipid force field without changing the membrane too much. In the original parameterisation, the Elba model gives a membrane with a too high AL and a too thin DHH (see Table 1), whereas the re-parameterised model gives for instance a too low AL and a too thin DHH (but closer to the experimental value). In fact, it can be argued that the re-parameterised Elba model is a better model than the original. Furthermore, the tail bonds become stiffer as shown by the although they are still more flexible than in the AT simulations. The relative success of the reparameterised Elba model also shows that the larger force constant observed in AT simulations compared to CG simulations has a minor effect on the ability to reproduce macroscopic properties. By design the current investigation is limited in nature. It has been shown that the choice of AT force field probably only has a minor effect on the CG parameterisation and that instead the mapping is the critical choice. The generality of these conclusions cannot be addressed herein but requires a more quantitative approach, which has been discussed previously[33]. Another issue to keep in mind is the potential danger in re-parameterising the bonded parameters without reparameterising the non-bonded ones. Generally, all the parameters in a force field are a detailed balance between different forces and changing only some parts could potentially disrupt this balance. For instance, it is possible that a proper POPC membrane with the re-parameterised Martini force field could be recovered by making the lipid beads smaller (i.e. changing the Lennard-Jones parameters). Therefore, the re-parameterised Martini model herein, can be seen as the starting point for further optimization of the non-bonded parameters. Still, this study has established that it is possible to parameterise CG lipid models using a bottom-up approach granted a good mapping is provided as with the Elba model. A discussion on how such a mapping should be designed is however outside the scope of this paper. Conclusions It has been shown that different AT simulations to a large extent sample the same bond and angle distributions, except in regions where CG modelling is ambiguous and challenging. Thus, the choice of the reference force field seems to be of minor importance. The major factor influencing the usefulness of a bottom-up parameterisation strategy seems instead to be the mapping. The mapping used for the Martini model leads to a re-parameterisation that displays too short and stiff bonds in the lipid tails, which drastically affects the physical properties of the membrane. On the other hand, the mapping used for the Elba model leads to a re-parameterisation

12

that is on a par or is in fact slightly better than the original parameterisation when it comes to reproducing structural experimental observables. Thus, it has to be concluded that a bottom-up approach holds great promise in parameterising bonded potentials in popular CG models for biomolecular applications. Acknowledgements I acknowledge the Wenner–Gren foundations for financial support and C3SE at Chalmers Technical University and iSolutions at the University of Southampton for computational resources. Sarah Witzke is acknowledged for fruitful discussions and Patricia Saenz Méndez for proof-reading.

13

References [1] J.L. Klepeis, K. Lindorff-Larsen, R.O. Dror, D.E. Shaw, Curr. Opin. Struct. Biol. 19 (2009) 120–7. [2] P.J. Bond, J. Holyoake, A. Ivetac, S. Khalid, M.S.P. Sansom, J. Struct. Biol. 157 (2007) 593–605. [3] M.G. Saunders, G.A. Voth, Annu. Rev. Biophys. 42 (2013) 73–93. [4] D.P. Tieleman, J.L. Maccallum, W.L. Ash, C. Kandt, Z. Xu, L. Monticelli, J. Phys. Condens. Matter 18 (2006) S1221–34. [5] S.J. Marrink, A.H. de Vries, A.E. Mark, J. Phys. Chem. B 108 (2004) 750–760. [6] S.J. Marrink, H.J. Risselada, S. Yefimov, D.P. Tieleman, A.H. de Vries, J. Phys. Chem. B 111 (2007) 7812–24. [7] D.H. de Jong, G. Singh, W.F.D. Bennett, C. Arnarez, T.A. Wassenaar, L. V. Schäfer, X. Periole, D.P. Tieleman, S.J. Marrink, J. Chem. Theory Comput. 9 (2013) 687–697. [8] K.A. Maerzke, J.I. Siepmann, J. Phys. Chem. B 115 (2011) 3452–65. [9] M. Orsi, D.Y. Haubertin, W.E. Sanderson, J.W. Essex, J. Phys. Chem. B 112 (2008) 802–15. [10] M. Orsi, J.W. Essex, PLoS One 6 (2011) e28637. [11] W.G. Noid, J. Chem. Phys. 139 (2013) 090901. [12] L.-P. Wang, T.J. Martinez, V.S. Pande, J. Phys. Chem. Lett. 5 (2014) 1885– 1891. [13] E. Brini, E.A. Algaer, P. Ganguly, C. Li, F. Rodríguez-Ropero, N.F.A. van der Vegt, Soft Matter 9 (2013) 2108–2119. [14] F. Ercolessi, J.B. Adams, Europhys. Lett. 26 (1994) 583–588. [15] A.P. Lyubartsev, A. Laaksonen, Phys. Rev. E 52 (1995) 3730–3737. [16] D. Reith, M. Pütz, F. Müller-Plathe, J. Comput. Chem. 24 (2003) 1624–36. [17] S. Izvekov, M. Parrinello, C.J. Burnham, G.A. Voth, J. Chem. Phys. 120 (2004) 10896–913. [18] A. Chaimovich, M.S. Shell, J. Chem. Phys. 134 (2011) 094112. [19] W. Tschöp, K. Kremer, J. Batoulis, T. Bürger, O. Hahn, Acta Polym. 49 (1998) 61–74. [20] V. Rühle, C. Junghans, A. Lukyanov, K. Kremer, D. Andrienko, J. Chem. Theory Comput. 5 (2009) 3211–3223. [21] H.A. Karimi-Varzaneh, H.-J. Qian, X. Chen, P. Carbone, F. Müller-Plathe, J. Comput. Chem. 32 (2011) 1475–1487. [22] A. Mirzoev, A.P. Lyubartsev, J. Chem. Theory Comput. 9 (2013) 1512–1520. [23] T. Bereau, K. Kremer, J. Chem. Theory Comput. 11 (2015) 2783–2791.. [24] J. Wang, R.M. Wolf, J.W. Caldwell, P.A. Kollman, D.A. Case, J. Comput. Chem. 25 (2004) 1157–74. [25] A.K. Malde, L. Zuo, M. Breeze, M. Stroet, D. Poger, P.C. Nair, C. Oostenbrink, A.E. Mark, J. Chem. Theory Comput. 7 (2011) 4026–4037. [26] K. Vanommeslaeghe, E. Hatcher, C. Acharya, S. Kundu, S. Zhong, J. Shim, E. Darian, O. Guvench, P. Lopes, I. Vorobyov, A.D. Mackerell, J. Comput. Chem. 31 (2010) 671–90. [27] W.D. Cornell, P. Cieplak, C.I. Bayly, I.R. Gould, K.M. Merz, D.M. Ferguson, D.C. Spellmeyer, T. Fox, J.W. Caldwell, P.A. Kollman, J. Am. Chem. Soc. 117 (1995) 5179–5197. [28] W.L. Jorgensen, D.S. Maxwell, J. Tirado-Rives, J. Am. Chem. Soc. 118 (1996) 11225–11236.

14

[29] C. Oostenbrink, A. Villa, A.E. Mark, W.F. van Gunsteren, J. Comput. Chem. 25 (2004) 1656–76. [30] R.B. Best, X. Zhu, J. Shim, P.E.M. Lopes, J. Mittal, M. Feig, A.D. Mackerell, J. Chem. Theory Comput. 8 (2012) 3257–3273. [31] C.F. Wong, J. Am. Chem. Soc. 113 (1991) 3208–3209. [32] G.J. Rocklin, D.L. Mobley, K.A. Dill, J. Chem. Theory Comput. 9 (2013) 3072– 3083. [33] J.W. Wagner, J.F. Dama, G.A. Voth, J. Chem. Theory Comput. 11 (2015) 3547– 3560. [34] E.L. Wu, X. Cheng, S. Jo, H. Rui, K.C. Song, E.M. Dávila-Contreras, Y. Qi, J. Lee, V. Monje-Galvan, R.M. Venable, J.B. Klauda, W. Im, J. Comput. Chem. 35 (2014) 1997–2004. [35] J.P.M. Jämbeck, A.P. Lyubartsev, J. Phys. Chem. B 116 (2012) 3164–79. [36] O. Berger, O. Edholm, F. Jähnig, Biophys. J. 72 (1997) 2002–13. [37] D. Poger, W.F. Van Gunsteren, A.E. Mark, J. Comput. Chem. 31 (2010) 1117– 1125. [38] B. Hess, C. Kutzner, D. van der Spoel, E. Lindahl, J. Chem. Theory Comput. 4 (2008) 435–447.435–447. [39] S. Plimpton, J. Comput. Phys. 117 (1995) 1–19. [40] B. Hess, H. Bekker, H.J.C. Berendsen, J.G.E.M. Fraaije, J. Comput. Chem. 18 (1997) 1463–1472. [41] S. Páll, B. Hess, Comput. Phys. Commun. 184 (2013) 2641–2650. [42] T. Darden, D. York, L. Pedersen, J. Chem. Phys. 98 (1993) 10089. [43] H.J.C. Berendsen, J.P.M. Postma, W.F. van Gunsteren, A. DiNola, J.R. Haak, J. Chem. Phys. 81 (1984) 3684. [44] S. Nosé, J. Chem. Phys. 81 (1984) 511. [45] W. Hoover, Phys. Rev. A 31 (1985) 1695–1697. [46] M. Parrinello, A. Rahman, J. Appl. Phys. 52 (1981) 7182. [47] I.G. Tironi, R. Sperb, P.E. Smith, W.F. van Gunsteren, J. Chem. Phys. 102 (1995) 5451. [48] G. Bussi, D. Donadio, M. Parrinello, J. Chem. Phys. 126 (2007) 014101. [49] P.H. Hünenberger, Adv. Polym. Sci. 173 (2005) 105–147. [50] M.D. Daily, B.N. Olsen, P.H. Schlesinger, D.S. Ory, N.A. Baker, J. Chem. Theory Comput. 10 (2014) 2137–2150. [51] D. Poger, A.E. Mark, J. Chem. Theory Comput. 6 (2010) 325–336. [52] S. Genheden, J. W. Essex, J. Chem. Theory Comput. 11 (2015) 4749-4759 [53] V.A. Harmandaris, D. Reith, N.F.A. van der Vegt, K. Kremer, Macromol. Chem. Phys. 208 (2007) 2109–2120. [54] N.H. Tattrie, J.R. Bennett, R. Cyr, Can. J. Biochem. 46 (1968) 819–24. [55] A.P. Eichenberger, W. Huang, S. Riniker, W.F. van Gunsteren, J. Chem. Theory Comput. 11 (2015) 2925–2937. [56] Klauda, J.B., Venable, R.M., Freites, J.A., O’Connor, J.W., Tobias, D.J., Mondragon-Ramirez, C., Vorobyov, I., MacKerell, A.D., Pastor, R.W., J. Phys. Chem. B 114 (2010) 7830–43. [57] Dickson, C.J., Rosso, L., Betz, R.M., Walker, R.C., Gould, I.R., Soft Matter 8 (2012) 9617–9627. [58] J.L. Klepeis, K. Lindorff-Larsen, R.O. Dror, D.E. Shaw, Curr. Opin. Struct. Biol. 19 (2009) 120–7.

15

Figure 1 – Atoms to CG bead mapping for Martini (above) and Elba (below). The shape of the region covering the atoms is only for illustrative purposes, all CG beads are spherical. The name of the beads is shown next to each region.

16

Figure 2 – Bond distributions for the Martini mapping from simulations with AT and Martini force fields. The distributions given by the force field parameters are denoted as “ff”.

17

Figure 3 – Angle distributions for the Martini mapping from simulations with AT and Martini force fields. The distributions given by the force field parameters are denoted as “ff”.

18

Figure 4 – Bond distributions for the Elba mapping from simulations with AT and Elba force fields. The distributions given by the force field parameters are denoted as “ff”.

19

Figure 5 – Angle distributions for the Elba mapping from simulations with AT and Elba force fields. The distributions given by the force field parameters are denoted as “ff”.

20

Table 1 – Properties of the POPC membrane with different force fields Type A [nm2] V [nm3] D L

AT

Slipids

AT

Berger Gromos

AT CG

Martini

CG

Elba Experiment[51]

L

HH [nm]

0.649 ±0.001

1.206 ±0.000

3.76 ±0.01

0.636 ±0.002

1.212 ±0.001

3.79 ±0.01

0.636 ±0.004

1.246 ±0.001

3.76 ±0.04

0.646 ±0.001

1.378 ±0.000

4.26 ±0.01

0.706 ±0.001 0.62 – 0.68

1.221 ±0.001

3.23 ±0.03

1.256

3.70

Uncertainties are the standard errors of the mean over the two repeats

21

Table 2 – Average order parameters, for POPC with different force fields sn-1 chain

sn-2 chain Martini mapping

Slipid

0.41

0.36

Berger

0.40

0.36

Gromos

0.43

0.37

Martini

0.37

0.31

Elba mapping Slipid

0.39

0.30

Berger

0.38

0.32

Gromos

0.41

0.33

0.23

0.17

Elba

The uncertainty is equal to or less than 0.01

22

Table 3 – Membrane properties in simulations with re-parameterised Martini and Elba force fields.1 Martini Bond + angle Bond

Angle

Bond + angle

Elba Bond

Angle

AL [nm2] Observed Difference

0.50 0.15

0.49 0.16

0.64 0.01

0.62 0.09

0.69 0.02

0.63 0.08

1.15 0.07

1.15 0.07

1.21 0.01

VL [nm3] Observed Difference

1.14 0.24

1.02 0.35

1.38 0.00

DHH [nm] Observed Difference

4.26 0.00

4.30 4.26 3.46 3.09 3.64 -0.04 0.00 -0.24 0.14 -0.41 sn-1 chain Observed 0.79 0.79 0.40 0.30 0.25 0.28 Difference -0.42 -0.42 -0.03 -0.07 -0.02 -0.05 sn-2 chain Observed 0.75 0.75 0.34 0.25 0.19 0.24 Difference -0.44 -0.44 -0.04 -0.09 -0.03 -0.07 1 Both the observed value and the difference compared to simulation with the original parameterisation are shown. The parameterisation is done on bonds and angles (bond + angle), or either bonds or angles individually.

23