Comment on
“Capturing Phase Behavior of Ternary Lipid Mixtures with a Reﬁned Martini
Coarse-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 Behavior
of Ternary Lipid Mixtures with a Reﬁned Martini Coarse-Grained Force Field” [Journal of
Chemical Theory and Computation 2018, 14, 11, 6050–6062]. This reﬁned Martini model was reported
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 only
emerged due to an unfortunate choice of simulation parameters, which leads to poor energy conservation.
Speciﬁcally, the constraints on the cholesterol model drained energy out from the membrane,
resulting in two coexisting phases at drastically diﬀerent temperatures. Using the simulation parameters
recommended for the used cholesterol model, this artefact is eliminated, yet so is phase
coexistence, i.e. experimental phase diagrams are no longer reproduced.
It is important to highlight that the present comment was submitted to Chemical Theory and
Computation. However, it was rejected without peer-review by the Editor-in-Chief, who stated that
the journal “rarely publishes such material”.
At room temperature and at certain ratios, the canonical
mixture of dipalmitoylphosphatidylcholine (DPPC)
with two saturated acyl chains, dioleoylphosphatidylcholine
(DOPC) with two monounsaturated acyl chains,
and cholesterol (CHOL) undergoes phase separation into
liquid ordered (Lo) and liquid disordered (Ld) phases
[1]. It is well-documented that the standard implementation
of the coarse-grained (CG) Martini model [2] does
not capture this behavior [3–5]. Therefore, phase separation
and phase coexistence studies performed using
the Martini model replace DOPC by dilinoleoylphosphatidylcholine
(DLiPC) with two diunsaturated acyl
chains as low transition temperature (Tm) lipid [6]. Unfortunately,
no experimental phase diagrams for the
DPPC/DLiPC/CHOL mixture exist to our knowledge,
rendering the comparison between simulations and experiments
qualitative at best.
Recently, Carpenter et al. tackled this outstanding issue
by a careful reﬁnement of the Martini parameters for
DOPC and DPPC. They adjusted the bonded parameters
of each bead separately to reproduce the bond length
and angle distributions extracted from atomistic simulation
data [3]. This approach diverts from the building
block concept of the version 2 family of the Martini force
ﬁeld, where the number of bonded parameters is minimized.
With this extended freedom to ﬁne-tuning of the
interactions, Carpenter et al. achieved almost quantitative
agreement with the experimental phase diagrams for
the mixture of DPPC, DOPC, and CHOL [1, 7]. Considering
the known limitations of the CG models [8], this is
an impressive milestone on the way to accurate modeling
of complex lipid mixtures. In their work, Carpenter
∗ Correspondence email address: matti.javanainen@gmail.com
et al. derived two diﬀerent models; “Extensible parameters”
where DPPC and DOPC head groups share identical
interaction parameters, and “Optimal parametes”,
where they are allowed to diﬀer. The latter parameter
set is of more interest, as it provides the best agreement
with the experimental phase diagram [3].
Here, we demonstrate that the reported liquid–liquid
coexistence in the model by Carpenter et al. (“Optimal
parameters”) is an artifact caused by an unfortunate
choice of simulation parameters for GROMACS. Namely,
the used current cholesterol model by Melo et al. [9] was
parametrized using a very conservative set of parameters
for the LINCS constraint algorithm to properly model the
ring structure of cholesterol. Unfortunately, these more
conservative constraint options were not followed by Carpenter
et al. [3], who instead followed the parameter set
generally recommended for the Martini model. When the
cholesterol model is simulated with these generally used
parameters, i.e. 4th
order expansion of the constraint
coupling matrix (lincs order equal to 4) and 1 iteration
of constraints per integration time step (lincs iter
equal to 1), energy is not properly conserved due to the
presence of virtual sites and corresponding coupled constraints
in this model. This eﬀect, which is also highlighted
in theGROMACS manual [10], manifests itself
especially at large integration time steps.
In the case of the study by Carpenter et al. [3], CHOL
constraints drain out a substantial amount of energy from
their neighborhood, which in the case of the studied
ternary lipid mixtures consists preferentially of DPPC.
Thus, this leads to the cooling down of CHOL and the
surrounding DPPC lipids which together form a sort of
ordered phase. In contrast — to maintain the target
temperature of the thermostat — DOPC-rich region gets
warmer than the target temperature of the thermostat,
thus forming a very disorered phase. This Maxwell de-
arXiv:2009.07767v1[physics.bio-ph]16Sep2020
2
mon causes the temperatures of DPPC (and CHOL) and
DOPC to diverge exponentially as a function of the integration
time step from the target temperature of the
thermostat, therefore resulting in the apparent Lo phase
being dozens of K cooler than the Ld phase. Importantly,
Lo/Ld coexistence regions reported by the authors
and matching the experimental phase diagram in
the DPPC/DOPC/CHOL mixtures are only recovered in
simulations in the presence of this artifact.
We show that this build-up of temperature diﬀerence
can be at least partially prevented by either 1) following
the suggested LINCS parameters for CHOL simulations
by Melo et al., 2) coupling the diﬀerent lipid types
to separate thermostats, or 3) decreasing the simulation
time step. All these options, that improve energy conservation,
cause the Lo/Ld phase coexistence to vanish.
Instead, an almost ideal mixing is instead observed.
Thus, when achieving a proper NPT ensemble, the
model by Carpenter et al. does not present an improvement
over the standard Martini v2.2 lipid model in capturing
the phase coexistence in DOPC/DPPC/CHOL
mixtures [3, 5, 6]. Therefore, further work is required
before a direct comparison of experimental and
Martini-based coarse grained simulations on phaseseparating
lipid mixtures is viable. Meanwhile, the
DPPC/DLiPC/CHOL mixture in the standard Martini
implementation will—despite its obvious limitations—
serve as the main model of choice for studies on Lo/Ld
phase coexistence.
I. RESULTS
To demonstrate the issues in the model by Carpenter
et al., we ﬁrst followed the original simulation
protocol reported in Ref. 3. We performed 15 µs
simulations of 3/3/2 and 2/1/1 compositions of the
DPPC/DOPC/CHOL mixture (see Methods and Set 1
in Table I). Both of these compositions fall in the heart
of the Lo/Ld coexistence region [3]. The simulations were
performed at 298 K using time steps of 10, 15, 20, 25, 30,
35, and 40 fs. Importantly, 35 fs was used in the original
work [3]. We characterized the degree of phase separation
using the mean contact fraction [4] (fmix, see Methods)
extracted from the last 1 µs simulation. For each value of
the integration time step, we extracted fmix and the temperature
of the lipids. With blue markers in Fig. 1A), we
demonstrate the obvious correlation of fmix of the function
of lipid temperature, both of which depend on the
used integration time step. Notably, larger time steps
lead to signiﬁcantly cooler membranes (< 298 K) and
smaller values of fmix, i.e. higher degree of phase separation.
With time steps that maintain the overall temperature
of the lipids reasonably close to the target temperature
of 298 K, both mixtures display nearly ideal
mixing instead. The fact that the target temperature
of the thermostat is not preserved even with small time
steps is in line with a recent report on the imperfect energy
conservation with the New-RF simulation settings
used [11].
We also extracted the temperatures of the diﬀerent
lipid types. As this information not available in the original
trajectories as no velocities were saved, we extended
each simulation by 100 ns, during which we also saved
the instantaneous velocities (see Methods). The temperatures
of each lipid type are shown as a function of the
integration time step in Fig. 1B). It is evident that the
increase in time step causes the low-Tm lipid (DOPC)
to heat up and the high-Tm (DPPC) lipid to cool down.
The behaviors of DPPC and CHOL are similar. For time
steps up to 20 fs, this eﬀect is negligible, yet with the
35 fs time step used by Carpenter et al., the temperature
diﬀerence between the lipid components is already
more than 20 K with DPPC being ∼30 K below its Tm
of 314 K. Thus, it is clear that decreasing the integration
time step improves energy conservation, yet leads to a
smaller level of phase separation.
To pinpoint the culprit of the poor energy conservation,
we performed additional simulations with the input
parameters slightly modiﬁed. First, we used the same
integration time step of 35 fs as Carpenter et al. but
with the LINCS parameters provided by Melo et al. [9]—
namely two LINCS iterations per time step, and an 8th
order expansion of the constraint coupling matrix (Set
3 in Table I). Data for these settings is shown in red in
Fig. 1A). It is evident that this more accurate handling
of the constraints leads to a much smaller degree of cooling
of the lipids. With these settings, phase separation
also essentially vanishes. The temperatures of each lipid
type in the simulations with a time step of 35 fs and
with various simulation settings are shown in Fig. 1C).
The more accurate LINCS settings clearly lead to all the
lipids being at almost the same temperature which also
closely matches the target one. This indicates that energy
conservation can be improved by at least these two
ways—decreasing the time step or increasing the number
of LINCS iterations and the order of the LINCS expansion.
The fact that “LINCS” and “20” overlap in
Fig. 1A) suggests that the eﬀect on phase separation is
independent on the means of improving energy conservation.
However, both approaches lead to a disagreement
between the simulation and the experimental phase diagram,
contrary to the results reported in the paper by
Carpenter et al..
We performed an additional simulation, in which we
attempted to improve energy conservation by using a separate
thermostat for each of the lipid types while still using
an integration time step of 35 fs (Set 2 in Table I). As
shown in Fig. 1A) by green dots, these settings lead to an
almost as cool membrane as in the case when all lipids
are coupled to the same thermostat. Still, the degree
of phase separation is signiﬁcantly decreased. Based on
Fig. 1C), the used temperature coupling does not prevent
CHOL and DPPC from cooling down, yet it does
not allow DOPC to heat up. Therefore, the average lipid
temperature decreases, whereas the temperature diﬀer-
3
0.3
0.4
0.5
Ideal mixing
Together
3/3/2
1020
30
40
Separate
LINCSfmix
Target T
A)
292 294 296 298
0.4
0.5
0.6
0.7
Together
Ideal mixing
2/1/1
10
20
30
40
Separate
LINCS
Lipid temperature (K)
10 20 30 40
260
280
300
320
Time step (fs)
Temperature(K)
DPPC DOPC
CHOL Solvent
B)
Together Separate LINCS
280
290
300
310
Temperature(K)
DPPC CHOL
DOPC Solvent
C)
Figure 1. LINCS settings and temperature coupling aﬀect the degree of phase separation. A) The degree of phase separation
as a function of the simulation temperature of the lipids, both of which depend on the integration time step, indicated by point
labels. Only the labels that are multiples of 10 are shown for clarity. The dotted gray line shows the target temperature of the
thermostat. Data for the two studied mixtures are shown in the two panels, and the fmix values corresponding to ideal mixing
are 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 are
only shown for the 3/3/2 mixture, yet the behavior in the 2/1/1 mixture is essentially identical. C) Temperatures of the lipid
types and the solvent beads as a function of temperature coupling and LINCS options. Data are again shown for the 3/3/2
mixture only.
ence between DPPC and DOPC is minor resulting in no
phase separation.
We also plot the ﬁnal structures of selected simulations
on the top row of Fig. 2. It is evident that proper phase
separation takes place with a time step of 35 fs and following
the simulation options of Carpenter et al.[3] (Set 1
in Table I), whereas decreasing the time step to 10 fs results
in an almost ideal mixing of the lipid types. Similar
behavior is observed with the conservative LINCS parameters
(Set 3 in Table I). However, for the case in which all
lipid types are coupled separately to a thermostat (Set
2 in Table I), some heterogeneity is still observed. The
local temperature maps, calculated from the 100 ns simulation
during which velocities were saved, are shown on
the bottom row of Fig. 2. It is evident that the observed
phase separation with the model by Carpenter et al. goes
hand in hand with the uneven temperature distribution
in the membrane, indicating that the systems is not sampling
the proper equilibrium ensemble. With the simulation
setup of Carpenter et al. [3], the temperatures of the
Lo and Ld phases can diﬀer as much as 100 K between
regions.
II. CONCLUSIONS
From the data presented above, we can draw the following
conclusions: 1) The use of LINCS parameters that
are incompatible with the used CHOL model, combined
with a large integration time step, results in poor energy
conservation. In the case of DPPC/DOPC/CHOL
mixture considered by Carpenter et al. [3], the virtual
site model of CHOL [9] drains out a signiﬁcant amount
of energy from the system, mainly from DPPC which locates
itself close to CHOL. 2) When all lipids are coupled
together to the thermostat, as was done by Carpenter
et al.[3], the cooling of CHOL and DPPC leads to the
heating up of DOPC. This temperature diﬀerence drives
the separation of the lipids into hot (DOPC-rich) and
cool (DPPC- and CHOL-rich) phases, whose tempera-
4
250 K
298 K
350 K
Together, 35 fs Together, 10 fs Separate LINCS
3/3/2
DPPCDPPC
DOPCDOPC
CHOLCHOL
Figure 2. Top row: ﬁnal 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 the
same thermostat (“Together”) and a 35 fs integration time step is used undergoes phase separation. Some heterogeneity is
also seen in other systems with DPPC and CHOL clustering together. Bottom row: the temperature maps calculated from the
100 ns simulations during which velocities are saved (see Methods). For the “Together, 35 fs” system, the two phases show
almost 100 K diﬀerence in temperature. When the lipid types are coupled to separate thermostats (“Separate”), the DPPC-rich
regions are slightly cooler than the remainder of the membrane. For the other systems, the temperature is close to the target
on (298 K) in the entire membrane.
tures diﬀer by signiﬁcantly. This behavior is in obvious
violation of the second law of thermodynamics. 3) The
poor energy conservation can be cured by either decreasing
the integration time step or using the LINCS parameters
that were used in the parametrization of the virtual
site model of CHOL [9]. This better energy conservations
leads to the loss of phase coexistence, and thus to poor
agreement with experimental phase diagrams. Additionally,
by coupling all lipid types to separate thermostats,
the temperature diﬀerence between DPPC and DOPC
is limited, and only a small degree of lipid demixing is
observed.
It is also worth mentioning that the standard implementation
of the current version 2.2 of the Martini
model does not display phase separation of the
DPPC/DOPC/CHOL mixture even when the simulation
settings of Carpenter et al. are used (set 4 in Table I).
This results from the fact that the standard Martini
DOPC model mixes fairly well with cholesterol. Indeed,
a hybrid DPPC/DOPC/CHOL mixture with the parameters
of Carpenter et al. used for DPPC and the standard
Martini parameters used for DOPC did not phaseseparate.
However, with DOPC parameters adapted
from Carpenter et al. and DPPC parameters from the
standard Martini implementation, phase separation was
recovered. This indicates that the unfavorable mixing of
DOPC and CHOL in the model by Carpenter et al. is
central for the co-cooling of DPPC with CHOL, which is
required for the formation of a cool Lo phase.
All in all, while the model of Carpenter et al. seemed
to ﬁnally reproduce the experimental phase diagram of
the well-studied DPPC/DOPC/CHOL lipid mixture, our
results demonstrate that it does due to a simulation artifact.
This artifact relates to a poor energy conservation,
which results from an unfortunate combination of large
integration time steps and not using the the LINCS parameters
required by the CHOL model. These choices
lead to a clear violation of target NPT ensemble. Unfortunately,
this compromises the validity of any ﬁndings
obtained with the model by Carpenter et al. [3], and suggests
that further work is still needed until a Martini-like
coarse-grained model reproduces the phase behavior of
the canonical DPPC/DOPC/CHOL mixture.
III. METHODS
A. Simulations
In the ﬁrst set of simulations, we strictly followed the
protocol of Carpenter et al. [3] and built lipid membranes
with dimensions of 30 × 30 × 15 nm2
, which had a total
of 3040 lipids and ∼77000 solvent beads, out of which
10% were modeled as antifreeze particles. The systems
were set up by the insane tool [12]. We considered two
compositions in the heart of the Lo/Ld coexistence region,
namely DPPC/DOPC/CHOL ratios of 3/3/2 and
2/1/1. We applied restraints with a force constant of
5
2 kJ/(mol×nm2
) to the phosphate beads of the phospholipids
in one leaﬂet. All in all, we used the same equilibration
protocol and simulation parameters as Carpenter
et al. in the commented paper [3]. Namely, the New-RF
simulation parameters [13] were used, unless otherwise
mentioned. To pinpoint the violation of the NPT ensemble
to the used CHOL model, we performed the following
simulations (see also Table I):
First, following Carpenter et al., the lipids and the
solvent were separately coupled to a thermostat with a
target temperature of 298 K, the 15 µs simulations with
time steps of 10, 15, 20, 25, 30, 35, and 40 fs, saving the
trajectories every 1 ns. Secondly, we performed 15 µs
simulations with a 35 fs time step but with all the lipid
types (DPPC, DOPC, and CHOL) coupled separately to
a thermostat. Thirdly, we performed 15 µs simulations
with two iterations of the LINCS algorithm (lincs iter
= 2) and using the 8th
order expansion of the constraint
coupling matrix (lincs order = 8). Here, all lipids were
coupled together to the thermostat.
We also run additional 100 ns simulations, starting
from the ﬁnal structures of the aforementioned simulations,
and saved the velocities (.trr ﬁle) so that the temperatures
of each lipid type could be extracted using gmx
traj. As the used cholesterol model has constraints, we
corrected for the missing degrees of freedom.
B. Analyses
The contact fraction fmix, describing the level of lipid
demixing, was adapted from Ref. 4, and is deﬁned as
fmix =
cUS−S
cUS−S + cUS−US
, (1)
where for example cUS−S refers to contacts between
unsaturated (US, here DOPC) and saturated (S, here
DPPC) lipids. Thus, the smaller the value of fmix, the
sharper the separation is. A contact is registered if the
phosphate 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 separately,
their temperatures were extracted separately and
the averaging was performed over the degrees of freedom.
For CHOL, this was less than 3N due to the constraints.
The temperatures of each lipid type were extracted from
the 100 ns trajectories containing velocities (see Simulations
Methods above) using gmx traj. The values were
corrected for the missing degrees of freedom of CHOL.
The heat maps of temperature distribution were calculated
by projecting the lipid center of mass onto the
macroscopic plane of the membrane. While performing
the binning, each point was weighted by the instantaneous
temperature of the corresponding molecule. For
DPPC and DOPC, the number of degrees of freedom in
each lipid was taken as 3N, while in the case of CHOL
the decrease of degrees of freedom in the presence of constraints
was accounted for.
Table I. Simulated systems. Composition is given as
DPPC/DOPC/CHOL. “∆t” is the integration time step in
fs. “T-coupl.” refers to the way the lipid temperatures are
coupled; the lipids are being coupled either together to one
thermostat or separately to three thermostats.
Composition ∆t T-coupl. lincs order lincs iter
Set 1. Lipids coupled Together
3/3/2 10 Together 4 1
3/3/2 15 Together 4 1
3/3/2 20 Together 4 1
3/3/2 25 Together 4 1
3/3/2 30 Together 4 1
3/3/2 35 Together 4 1
3/3/2 40 Together 4 1
2/1/1 10 Together 4 1
2/1/1 15 Together 4 1
2/1/1 20 Together 4 1
2/1/1 25 Together 4 1
2/1/1 30 Together 4 1
2/1/1 35 Together 4 1
2/1/1 40 Together 4 1
Set 2. Lipids coupled Separately
3/2/2 35 Separately 4 1
2/1/1 35 Separately 4 1
Set 3. Conservative LINCS settings
3/2/2 35 Together 8 2
2/1/1 35 Together 8 2
Set 4. Unmodiﬁed Martini v2.2
3/2/2 35 Together 4 1
ACKNOWLEDGMENTS
M.J., H.M.-S., and B.F. acknowledge support from the
Czech Science Foundation (EXPRO grant 19-26854X).
M.J. thanks the Emil Aaltonen Foundation for funding.
6
[1] S. L. Veatch and S. L. Keller, Biophys. J. 85, 3074 (2003).
[2] S. J. Marrink, H. J. Risselada, S. Yeﬁmov, D. P. Tieleman,
and A. H. De Vries, J. Phys. Chem. B 111, 7812
(2007).
[3] T. S. Carpenter, C. A. L´opez, C. Neale, C. Montour,
H. I. Ing´olfsson, F. Di Natale, F. C. Lightstone, and
S. Gnanakaran, J. Chem. Theory Comput. 14, 6050
(2018).
[4] J. Doma´nski, S. J. Marrink, and L. V. Sch¨afer, BBA –
Biomembranes 1818, 984 (2012).
[5] R. S. Davis, P. Sunil Kumar, M. M. Sperotto, and
M. Laradji, J. Phys. Chem. B 117, 4072 (2013).
[6] H. J. Risselada and S. J. Marrink, Proc. Natl. Acad. Sci.
USA 105, 17367 (2008).
[7] J. H. Davis, J. J. Clair, and J. Juhasz, Biophys. J. 96,
521 (2009).
[8] R. Alessandri, P. C. Souza, S. Thallmair, M. N. Melo,
A. H. De Vries, and S. J. Marrink, J. Chem. Theory
Comput. 15, 5448 (2019).
[9] M. Melo, H. Ing´olfsson, and S. Marrink, J. Chem. Phys.
143, 12B637 1 (2015).
[10] M. J. Abraham, T. Murtola, R. Schulz, S. P´all, J. C.
Smith, B. Hess, and E. Lindahl, SoftwareX 1, 19 (2015).
[11] F. Benedetti and C. Loison, Comput. Phys. Commun.
228, 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.
11, 2144 (2015).
[13] D. H. de Jong, S. Baoukina, H. I. Ing´olfsson, and S. J.
Marrink, Comput. Phys. Commun. 199, 1 (2016).
IV. DATA AVAILABILITY
The inputs and outputs for our simulations with the
model by Carpenter et al. are available as follows:
1. Long simulations with diﬀerent time steps:
• DOI: 10.5281/zenodo.3956709 (3/3/2 mixture)
• DOI: 10.5281/zenodo.3956797 (2/1/1 mixture)
2. Short simulations with diﬀerent time steps and
with 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 with
diﬀerent temperature coupling or LINCS settings:
• DOI: 10.5281/zenodo.3956775 (3/3/2 mixture)
• DOI: 10.5281/zenodo.3956814 (2/1/1 mixture)