Comment on "Capturing Phase Behavior of Ternary Lipid Mixtures with a Refined Martini Coarse-Grained Force Field"
CComment on“Capturing Phase Behavior of Ternary Lipid Mixtures with a Refined MartiniCoarse-Grained Force Field”
Matti Javanainen, ∗ Balazs Fabian, and Hector Martinez-Seara
Institute of Organic Chemistry and Biochemistry, Czech Academy of Sciences,Flemingovo n´am. 542/2, 160 00 Prague 6, Czech Republic (Dated: September 17, 2020)We report here on the pitfalls of the simulation model introduced in the “Capturing Phase Be-havior of Ternary Lipid Mixtures with a Refined Martini Coarse-Grained Force Field” [Journal ofChemical Theory and Computation 2018, 14, 11, 6050–6062]. This refined Martini model was re-ported to reproduce experimental phase diagrams for a ternary DOPC/DPPC/cholesterol mixture,including the coexistence of two liquid phases. However, we demonstrate that this coexistence onlyemerged due to an unfortunate choice of simulation parameters, which leads to poor energy conser-vation. Specifically, the constraints on the cholesterol model drained energy out from the membrane,resulting in two coexisting phases at drastically different temperatures. Using the simulation pa-rameters recommended for the used cholesterol model, this artefact is eliminated, yet so is phasecoexistence, i.e. experimental phase diagrams are no longer reproduced.
It is important to highlight that the present comment was submitted to Chemical Theory andComputation. However, it was rejected without peer-review by the Editor-in-Chief, who stated thatthe journal “rarely publishes such material”.
At room temperature and at certain ratios, the canon-ical mixture of dipalmitoylphosphatidylcholine (DPPC)with two saturated acyl chains, dioleoylphosphatidyl-choline (DOPC) with two monounsaturated acyl chains,and cholesterol (CHOL) undergoes phase separation intoliquid ordered (L o ) and liquid disordered (L d ) phases[1]. It is well-documented that the standard implemen-tation of the coarse-grained (CG) Martini model [2] doesnot capture this behavior [3–5]. Therefore, phase sep-aration and phase coexistence studies performed usingthe Martini model replace DOPC by dilinoleoylphos-phatidylcholine (DLiPC) with two diunsaturated acylchains as low transition temperature ( T m ) lipid [6]. Un-fortunately, no experimental phase diagrams for theDPPC/DLiPC/CHOL mixture exist to our knowledge,rendering the comparison between simulations and ex-periments qualitative at best.Recently, Carpenter et al. tackled this outstanding is-sue by a careful refinement of the Martini parameters forDOPC and DPPC. They adjusted the bonded parame-ters of each bead separately to reproduce the bond lengthand angle distributions extracted from atomistic simula-tion data [3]. This approach diverts from the buildingblock concept of the version 2 family of the Martini forcefield, where the number of bonded parameters is mini-mized. With this extended freedom to fine-tuning of theinteractions, Carpenter et al. achieved almost quantita-tive agreement with the experimental phase diagrams forthe mixture of DPPC, DOPC, and CHOL [1, 7]. Consid-ering the known limitations of the CG models [8], this isan impressive milestone on the way to accurate model-ing of complex lipid mixtures. In their work, Carpenter ∗ Correspondence email address: [email protected] et al. derived two different models; “Extensible parame-ters” where DPPC and DOPC head groups share iden-tical interaction parameters, and “Optimal parametes”,where they are allowed to differ. The latter parameterset is of more interest, as it provides the best agreementwith the experimental phase diagram [3].Here, we demonstrate that the reported liquid–liquidcoexistence in the model by Carpenter et al. (“Opti-mal parameters”) is an artifact caused by an unfortunatechoice of simulation parameters for GROMACS. Namely,the used current cholesterol model by Melo et al. [9] wasparametrized using a very conservative set of parametersfor the LINCS constraint algorithm to properly model thering structure of cholesterol. Unfortunately, these moreconservative constraint options were not followed by Car-penter et al. [3], who instead followed the parameter setgenerally recommended for the Martini model. When thecholesterol model is simulated with these generally usedparameters, i.e. th order expansion of the constraintcoupling matrix ( lincs order equal to 4) and 1 itera-tion of constraints per integration time step ( lincs iter equal to 1), energy is not properly conserved due to thepresence of virtual sites and corresponding coupled con-straints in this model. This effect, which is also high-lighted in theGROMACS manual [10], manifests itselfespecially at large integration time steps.In the case of the study by Carpenter et al. [3], CHOLconstraints drain out a substantial amount of energy fromtheir neighborhood, which in the case of the studiedternary lipid mixtures consists preferentially of DPPC.Thus, this leads to the cooling down of CHOL and thesurrounding DPPC lipids which together form a sort ofordered phase. In contrast — to maintain the targettemperature of the thermostat — DOPC-rich region getswarmer than the target temperature of the thermostat,thus forming a very disorered phase. This Maxwell de- a r X i v : . [ phy s i c s . b i o - ph ] S e p mon causes the temperatures of DPPC (and CHOL) andDOPC to diverge exponentially as a function of the in-tegration time step from the target temperature of thethermostat, therefore resulting in the apparent L o phasebeing dozens of K cooler than the L d phase. Impor-tantly, L o /L d coexistence regions reported by the au-thors and matching the experimental phase diagram inthe DPPC/DOPC/CHOL mixtures are only recovered insimulations in the presence of this artifact.We show that this build-up of temperature differencecan be at least partially prevented by either 1) follow-ing the suggested LINCS parameters for CHOL simula-tions by Melo et al. , 2) coupling the different lipid typesto separate thermostats, or 3) decreasing the simulationtime step. All these options, that improve energy con-servation, cause the L o /L d phase coexistence to vanish.Instead, an almost ideal mixing is instead observed.Thus, when achieving a proper NPT ensemble, themodel by Carpenter et al. does not present an improve-ment over the standard Martini v2.2 lipid model in cap-turing the phase coexistence in DOPC/DPPC/CHOLmixtures [3, 5, 6]. Therefore, further work is re-quired before a direct comparison of experimental andMartini-based coarse grained simulations on phase-separating lipid mixtures is viable. Meanwhile, theDPPC/DLiPC/CHOL mixture in the standard Martiniimplementation will—despite its obvious limitations—serve as the main model of choice for studies on L o /L d phase coexistence. I. RESULTS
To demonstrate the issues in the model by Car-penter et al. , we first followed the original simulationprotocol reported in Ref. 3. We performed 15 µ ssimulations of 3/3/2 and 2/1/1 compositions of theDPPC/DOPC/CHOL mixture (see Methods and Set 1in Table I). Both of these compositions fall in the heartof the L o /L d coexistence region [3]. The simulations wereperformed at 298 K using time steps of 10, 15, 20, 25, 30,35, and 40 fs. Importantly, 35 fs was used in the originalwork [3]. We characterized the degree of phase separationusing the mean contact fraction [4] ( f mix , see Methods)extracted from the last 1 µ s simulation. For each value ofthe integration time step, we extracted f mix and the tem-perature of the lipids. With blue markers in Fig. 1A), wedemonstrate the obvious correlation of f mix of the func-tion of lipid temperature, both of which depend on theused integration time step. Notably, larger time stepslead to significantly cooler membranes ( <
298 K) andsmaller values of f mix , i.e. higher degree of phase sepa-ration. With time steps that maintain the overall tem-perature of the lipids reasonably close to the target tem-perature of 298 K, both mixtures display nearly idealmixing instead. The fact that the target temperatureof the thermostat is not preserved even with small timesteps is in line with a recent report on the imperfect en- ergy conservation with the New-RF simulation settingsused [11].We also extracted the temperatures of the differentlipid types. As this information not available in the orig-inal trajectories as no velocities were saved, we extendedeach simulation by 100 ns, during which we also savedthe instantaneous velocities (see Methods). The temper-atures of each lipid type are shown as a function of theintegration time step in Fig. 1B). It is evident that theincrease in time step causes the low- T m lipid (DOPC)to heat up and the high- T m (DPPC) lipid to cool down.The behaviors of DPPC and CHOL are similar. For timesteps up to 20 fs, this effect is negligible, yet with the35 fs time step used by Carpenter et al. , the tempera-ture difference between the lipid components is alreadymore than 20 K with DPPC being ∼
30 K below its T m of 314 K. Thus, it is clear that decreasing the integrationtime step improves energy conservation, yet leads to asmaller level of phase separation.To pinpoint the culprit of the poor energy conserva-tion, we performed additional simulations with the inputparameters slightly modified. First, we used the sameintegration time step of 35 fs as Carpenter et al. butwith the LINCS parameters provided by Melo et al. [9]—namely two LINCS iterations per time step, and an 8 th order expansion of the constraint coupling matrix (Set3 in Table I). Data for these settings is shown in red inFig. 1A). It is evident that this more accurate handlingof the constraints leads to a much smaller degree of cool-ing of the lipids. With these settings, phase separationalso essentially vanishes. The temperatures of each lipidtype in the simulations with a time step of 35 fs andwith various simulation settings are shown in Fig. 1C).The more accurate LINCS settings clearly lead to all thelipids being at almost the same temperature which alsoclosely matches the target one. This indicates that en-ergy conservation can be improved by at least these twoways—decreasing the time step or increasing the num-ber of LINCS iterations and the order of the LINCS ex-pansion. The fact that “LINCS” and “20” overlap inFig. 1A) suggests that the effect on phase separation isindependent on the means of improving energy conserva-tion. However, both approaches lead to a disagreementbetween the simulation and the experimental phase di-agram, contrary to the results reported in the paper byCarpenter et al. .We performed an additional simulation, in which weattempted to improve energy conservation by using a sep-arate thermostat for each of the lipid types while still us-ing an integration time step of 35 fs (Set 2 in Table I). Asshown in Fig. 1A) by green dots, these settings lead to analmost as cool membrane as in the case when all lipidsare coupled to the same thermostat. Still, the degreeof phase separation is significantly decreased. Based onFig. 1C), the used temperature coupling does not pre-vent CHOL and DPPC from cooling down, yet it doesnot allow DOPC to heat up. Therefore, the average lipidtemperature decreases, whereas the temperature differ- . . . f m i x Target T A)
292 294 296 2980 . . . . T e m p e r a t u r e ( K ) DPPC DOPCCHOL Solvent B) Together Separate LINCS280290300310 T e m p e r a t u r e ( K ) DPPC CHOLDOPC Solvent C) Figure 1. LINCS settings and temperature coupling affect the degree of phase separation. A) The degree of phase separationas a function of the simulation temperature of the lipids, both of which depend on the integration time step, indicated by pointlabels. Only the labels that are multiples of 10 are shown for clarity. The dotted gray line shows the target temperature of thethermostat. Data for the two studied mixtures are shown in the two panels, and the f mix values corresponding to ideal mixingare highlighted by dashed gray lines. The blue markers show data for simulations in which all lipids are coupled Together to the thermostat (Set 1 in Table I). Green markers show data for simulations in which lipid types are coupled to
Separate thermostats (Set 2 in Table I). Red markers show data for simulations in which more conservative
LINCS parameters are used(Set 3 in Table I). B) Temperatures of the lipid types and the solvent beads as a function of the integration time step. Data areonly shown for the 3/3/2 mixture, yet the behavior in the 2/1/1 mixture is essentially identical. C) Temperatures of the lipidtypes and the solvent beads as a function of temperature coupling and LINCS options. Data are again shown for the 3/3/2mixture only. ence between DPPC and DOPC is minor resulting in nophase separation.We also plot the final structures of selected simulationson the top row of Fig. 2. It is evident that proper phaseseparation takes place with a time step of 35 fs and fol-lowing the simulation options of Carpenter et al. [3] (Set 1in Table I), whereas decreasing the time step to 10 fs re-sults in an almost ideal mixing of the lipid types. Similarbehavior is observed with the conservative LINCS param-eters (Set 3 in Table I). However, for the case in which alllipid types are coupled separately to a thermostat (Set2 in Table I), some heterogeneity is still observed. Thelocal temperature maps, calculated from the 100 ns sim-ulation during which velocities were saved, are shown onthe bottom row of Fig. 2. It is evident that the observedphase separation with the model by Carpenter et al. goeshand in hand with the uneven temperature distributionin the membrane, indicating that the systems is not sam-pling the proper equilibrium ensemble. With the simula-tion setup of Carpenter et al. [3], the temperatures of the L o and L d phases can differ as much as 100 K betweenregions. II. CONCLUSIONS
From the data presented above, we can draw the fol-lowing conclusions: The use of LINCS parameters thatare incompatible with the used CHOL model, combinedwith a large integration time step, results in poor en-ergy conservation. In the case of DPPC/DOPC/CHOLmixture considered by Carpenter et al. [3], the virtualsite model of CHOL [9] drains out a significant amountof energy from the system, mainly from DPPC which lo-cates itself close to CHOL. When all lipids are coupledtogether to the thermostat, as was done by Carpenter et al. [3], the cooling of CHOL and DPPC leads to theheating up of DOPC. This temperature difference drivesthe separation of the lipids into hot (DOPC-rich) andcool (DPPC- and CHOL-rich) phases, whose tempera-
250 K
298 K
350 KTogether, 35 fs Together, 10 fs Separate LINCS
DPPCDPPCDOPCDOPCCHOLCHOL
Figure 2. Top row: final structures of chosen simulations with the 3/3/2 DPPC/DOPC/CHOL mixture after 15 µ s. DPPC,DOPC, and CHOL are shown in red, green, and yellow, respectively. Only the system where all lipids are coupled to thesame thermostat (“Together”) and a 35 fs integration time step is used undergoes phase separation. Some heterogeneity isalso seen in other systems with DPPC and CHOL clustering together. Bottom row: the temperature maps calculated from the100 ns simulations during which velocities are saved (see Methods). For the “Together, 35 fs” system, the two phases showalmost 100 K difference in temperature. When the lipid types are coupled to separate thermostats (“Separate”), the DPPC-richregions are slightly cooler than the remainder of the membrane. For the other systems, the temperature is close to the targeton (298 K) in the entire membrane. tures differ by significantly. This behavior is in obviousviolation of the second law of thermodynamics. Thepoor energy conservation can be cured by either decreas-ing the integration time step or using the LINCS param-eters that were used in the parametrization of the virtualsite model of CHOL [9]. This better energy conservationsleads to the loss of phase coexistence, and thus to pooragreement with experimental phase diagrams. Addition-ally, by coupling all lipid types to separate thermostats,the temperature difference between DPPC and DOPCis limited, and only a small degree of lipid demixing isobserved.It is also worth mentioning that the standard im-plementation of the current version 2.2 of the Mar-tini model does not display phase separation of theDPPC/DOPC/CHOL mixture even when the simulationsettings of Carpenter et al. are used (set 4 in Table I).This results from the fact that the standard MartiniDOPC model mixes fairly well with cholesterol. Indeed,a hybrid DPPC/DOPC/CHOL mixture with the param-eters of Carpenter et al. used for DPPC and the stan-dard Martini parameters used for DOPC did not phase-separate. However, with DOPC parameters adaptedfrom Carpenter et al. and DPPC parameters from thestandard Martini implementation, phase separation wasrecovered. This indicates that the unfavorable mixing ofDOPC and CHOL in the model by Carpenter et al. iscentral for the co-cooling of DPPC with CHOL, which isrequired for the formation of a cool L o phase. All in all, while the model of Carpenter et al. seemedto finally reproduce the experimental phase diagram ofthe well-studied DPPC/DOPC/CHOL lipid mixture, ourresults demonstrate that it does due to a simulation arti-fact. This artifact relates to a poor energy conservation,which results from an unfortunate combination of largeintegration time steps and not using the the LINCS pa-rameters required by the CHOL model. These choiceslead to a clear violation of target NPT ensemble. Un-fortunately, this compromises the validity of any findingsobtained with the model by Carpenter et al. [3], and sug-gests that further work is still needed until a Martini-likecoarse-grained model reproduces the phase behavior ofthe canonical DPPC/DOPC/CHOL mixture. III. METHODSA. Simulations
In the first set of simulations, we strictly followed theprotocol of Carpenter et al. [3] and built lipid membraneswith dimensions of 30 × ×
15 nm , which had a totalof 3040 lipids and ∼ insane tool [12]. We considered twocompositions in the heart of the L o /L d coexistence re-gion, namely DPPC/DOPC/CHOL ratios of 3/3/2 and2/1/1. We applied restraints with a force constant of2 kJ/(mol × nm ) to the phosphate beads of the phospho-lipids in one leaflet. All in all, we used the same equili-bration protocol and simulation parameters as Carpenter et al. in the commented paper [3]. Namely, the New-RFsimulation parameters [13] were used, unless otherwisementioned. To pinpoint the violation of the NPT ensem-ble to the used CHOL model, we performed the followingsimulations (see also Table I):First, following Carpenter et al. , the lipids and thesolvent were separately coupled to a thermostat with atarget temperature of 298 K, the 15 µ s simulations withtime steps of 10, 15, 20, 25, 30, 35, and 40 fs, saving thetrajectories every 1 ns. Secondly, we performed 15 µ ssimulations with a 35 fs time step but with all the lipidtypes (DPPC, DOPC, and CHOL) coupled separately toa thermostat. Thirdly, we performed 15 µ s simulationswith two iterations of the LINCS algorithm ( lincs iter= 2 ) and using the 8 th order expansion of the constraintcoupling matrix ( lincs order = 8 ). Here, all lipids werecoupled together to the thermostat.We also run additional 100 ns simulations, startingfrom the final structures of the aforementioned simula-tions, and saved the velocities (.trr file) so that the tem-peratures of each lipid type could be extracted using gmxtraj . As the used cholesterol model has constraints, wecorrected for the missing degrees of freedom. B. Analyses
The contact fraction f mix , describing the level of lipiddemixing, was adapted from Ref. 4, and is defined as f mix = c US − S c US − S + c US − US , (1)where for example c US − S refers to contacts betweenunsaturated (US, here DOPC) and saturated (S, hereDPPC) lipids. Thus, the smaller the value of f mix , thesharper the separation is. A contact is registered if thephosphate beads of the lipids were within 1.1 nm.The temperatures of the lipids were extracted using gmx energy . When the lipid types were coupled sepa-rately, their temperatures were extracted separately andthe averaging was performed over the degrees of freedom.For CHOL, this was less than 3 N due to the constraints.The temperatures of each lipid type were extracted fromthe 100 ns trajectories containing velocities (see Simula-tions Methods above) using gmx traj . The values werecorrected for the missing degrees of freedom of CHOL.The heat maps of temperature distribution were cal-culated by projecting the lipid center of mass onto themacroscopic plane of the membrane. While performingthe binning, each point was weighted by the instanta-neous temperature of the corresponding molecule. ForDPPC and DOPC, the number of degrees of freedom ineach lipid was taken as 3 N , while in the case of CHOLthe decrease of degrees of freedom in the presence of con-straints was accounted for. Table I. Simulated systems. Composition is given asDPPC/DOPC/CHOL. “∆ t ” is the integration time step infs. “ T -coupl.” refers to the way the lipid temperatures arecoupled; the lipids are being coupled either together to onethermostat or separately to three thermostats.Composition ∆ t T -coupl. lincs order lincs iter Set 1. Lipids coupled
Together
Separately
LINCS settings3/2/2 35 Together 8 22/1/1 35 Together 8 2Set 4.
Unmodified
Martini v2.23/2/2 35 Together 4 1
ACKNOWLEDGMENTS
M.J., H.M.-S., and B.F. acknowledge support from theCzech Science Foundation (EXPRO grant 19-26854X).M.J. thanks the Emil Aaltonen Foundation for funding. [1] S. L. Veatch and S. L. Keller, Biophys. J. , 3074 (2003).[2] S. J. Marrink, H. J. Risselada, S. Yefimov, D. P. Tiele-man, and A. H. De Vries, J. Phys. Chem. B , 7812(2007).[3] T. S. Carpenter, C. A. L´opez, C. Neale, C. Montour,H. I. Ing´olfsson, F. Di Natale, F. C. Lightstone, andS. Gnanakaran, J. Chem. Theory Comput. , 6050(2018).[4] J. Doma´nski, S. J. Marrink, and L. V. Sch¨afer, BBA –Biomembranes , 984 (2012).[5] R. S. Davis, P. Sunil Kumar, M. M. Sperotto, andM. Laradji, J. Phys. Chem. B , 4072 (2013).[6] H. J. Risselada and S. J. Marrink, Proc. Natl. Acad. Sci.USA , 17367 (2008).[7] J. H. Davis, J. J. Clair, and J. Juhasz, Biophys. J. ,521 (2009).[8] R. Alessandri, P. C. Souza, S. Thallmair, M. N. Melo,A. H. De Vries, and S. J. Marrink, J. Chem. TheoryComput. , 5448 (2019).[9] M. Melo, H. Ing´olfsson, and S. Marrink, J. Chem. Phys. , 12B637 1 (2015).[10] M. J. Abraham, T. Murtola, R. Schulz, S. P´all, J. C.Smith, B. Hess, and E. Lindahl, SoftwareX , 19 (2015).[11] F. Benedetti and C. Loison, Comput. Phys. Commun. , 146 (2018).[12] T. A. Wassenaar, H. I. Ing´olfsson, R. A. B¨ockmann, D. P.Tieleman, and S. J. Marrink, J. Chem. Theory Comput. , 2144 (2015). [13] D. H. de Jong, S. Baoukina, H. I. Ing´olfsson, and S. J.Marrink, Comput. Phys. Commun. , 1 (2016). IV. DATA AVAILABILITY
The inputs and outputs for our simulations with themodel by Carpenter et al. are available as follows:1. Long simulations with different time steps: • DOI: 10.5281/zenodo.3956709 (3/3/2 mixture) • DOI: 10.5281/zenodo.3956797 (2/1/1 mixture)2. Short simulations with different time steps andwith velocities saved: • DOI: 10.5281/zenodo.3956761 (3/3/2 mixture) • DOI: 10.5281/zenodo.3956812 (2/1/1 mixture)3. Additional simulations of the 3/3/2 mixture withdifferent temperature coupling or LINCS settings: • DOI: 10.5281/zenodo.3956775 (3/3/2 mixture) ••