Ripple-like instability in the simulated gel phase of finite size phosphocholine bilayers
Vivien Walter, Céline Ruscher, Adrien Gola, Carlos M. Marques, Olivier Benzerara, Fabrice Thalmann
RRipple-like instability in the simulated gel phase of finitesize phosphocholine bilayers
Vivien Walter a, ∗ , C´eline Ruscher b , Adrien Gola b , Carlos M. Marques b ,Olivier Benzerara b , Fabrice Thalmann b, ∗∗ a Department of Chemistry, King’s College London, Britannia House, 7 Trinity Street,SE1 1DB, London, United Kingdom b Institut Charles Sadron, CNRS and University of Strasbourg, 23 rue du Loess, F-67034Strasbourg, France
Abstract
Atomistic molecular dynamics simulations have reached a degree of ma-turity that makes it possible to investigate the lipid polymorphism of modelbilayers over a wide range of temperatures. However if both the fluid L α and tilted gel L β states are routinely obtained, the P β ripple phase ofphosphatidylcholine lipid bilayers is still unsatifactorily described. Perform-ing simulations of lipid bilayers made of different numbers of DPPC (1,2-dipalmitoylphosphatidylcholine) molecules ranging from 32 to 512, we demon-strate that the tilted gel phase L β expected below the pre-transition cannotbe obtained for large systems ( >
94 DPPC molecules) through common sim-ulations settings or temperature treatments. Large systems are instead foundin a disordered gel phase which display configurations, topography and ener-gies reminiscent from the ripple phase P β observed between the pretransitionand the main melting transition. We show how the state of the bilayers be-low the pretransition can be controlled and depends on thermal history andconditions of preparations. A mechanism for the observed topographic in-stability is suggested. Keywords:
Lipid bilayers, Molecular dynamics simulations, Phasetransition ∗∗∗
Email addresses: [email protected] (Fabrice Thalmann), [email protected] (Fabrice Thalmann)
Preprint submitted to Elsevier February 9, 2021 a r X i v : . [ c ond - m a t . s o f t ] F e b . Introduction Lipid membranes are fundamental components of living organisms, fortheir pivotal role in the structure and the biochemistry of the cells. Theirproperties are for a large part the consequence of the organisation of thephospholipid molecules that compose them to a large extent. This organi-sation is best revealed by experimenting with artificial lipid bilayers [1, 2].Through a perfect control of the molecular composition, an extended rangeof geometries and the possibility to insert membrane proteins and all sortof molecules, artificial bilayers have become a standard tool in modern bio-physics [3, 4].A striking property of pure phospholipids bilayers is to exhibit a num-ber of thermodynamic transitions upon temperature changes, meaning thatthey can be found in several phases including a dense structured L β tiltedgel phase at low temperature or a disordered fluid phase at higher tempera-ture L α [3, 4, 5, 6]. Decades ago, a new phase have been identified in the mostcommon sort of phospholipid, the phosphocholines (PC lipids) [7, 8, 3]. Thisnew phase, specific to the PC lipids and called the ripple phase P β , is char-acterised by important corrugation along the bilayer and the alternation be-tween interdigitated leaflets and non-interdigitated ones. The ripple phase isobserved experimentally in pure bilayers of 1,2-dipalmitoylphosphatidylcholine(DMPC) between 16 and 24 ◦ C, 1,2-dipalmitoylphosphatidylcholine (DPPC)between 34 and 41 ◦ C and 1,2-distearoylphosphatidylcholine (DSPC) between51 and 55 ◦ C [9, 3].This structure is sometimes presented as resulting from an alternation ofgel and fluid lipid configurations [10]. There is no consensus regarding thecause[11, 12, 13, 14, 15], and the experimental structure is still subject todetailed investigations [16]. The interest in the ripple phase has grown overthe years, and intense research in its nature are made jointly in experimentaland numerical sciences [17, 18, 19].In a recent paper, Khakbaz and Klauda investigated the phase transitionof DMPC and DPPC [20]. They reported formation of a structure resemblingthe P β phase for DMPC bilayers in a range of temperature, while for DPPConly a transition from the tilted to the fluid phase was observed in bilayerscomposed of 70 lipids. More recently, we investigated the phase transitionof DPPC bilayers using Machine Learning algorithms [21]. In membranes2omposed of 212 lipids, we observed a transition at 315 K from a fluid phaseto a condensed disordered gel phase similar to the ripple phase that persistswell below the pre-transition temperature of 307K and the L β tilted phasewas never obtained directly from a quench of the L α state.In this article, we provide a detailed analysis of the nature of the phaseof DPPC simulated with the CHARMM36 force field below its melting tem-perature. We investigate the state of DPPC bilayers at 288 K, for differentsystem sizes and different thermal routes. While the L β phase was observedfor small systems, we found that a disordered gel phase, reminiscent from the P β phase, occurs in larger systems. In particular, we noticed the formationof corrugations whose amplitude increases with the system sizes investigated.By applying different thermal treatments, we characterize the metastabilityof DPPC and show that disordered gel is the preferred phase of simulatedDPPC at 288K, in the thermodynamic limit.
2. Material & Methods
DPPC bilayers were obtained using the CHARMM-GUI online MembraneBuilder [22, 23, 24, 25]. The size of the membrane was controlled via thenumber of lipids, ranging from 32 to 512 and with equal amounts of lipidsin both leaflets. The bilayers were hydrated with water blocks of 10 nmon each side to prevent any interaction of the leaflets through the PBCsalong the Z-axis. The exact composition of each system used is given in theSupplementary Informations.The constructed system were minimized in energy using a steepest de-scent algorithm and equilibrated by running two NVT simulations at 288 Kfor 10 ps, with respectively a 0.001 ps and a 0.002 ps step; and two NPTsimulations at 288 K and 1 bar, with a 0.002 ps step for respectively 100 psand 1 ns.
Unless specified, the following conditions and parameters have been usedfor all simulations. All simulations were performed using GROMACS 2016.4 [26,27] along with the CHARMM-36 all-atom force-field [28] (June 2015 version).The force field parameters for the lipid molecules were provided directly byCHARMM-GUI [29, 30]. 3ll the molecular dynamics simulations used the leap-frog integration al-gorithm [31] with a time step set to 2 fs. The temperature was controlledduring the simulation using a Nos´e-Hoover thermostat [32, 33] with a cor-relation time of τ T = 0 . τ P = 2 . . × − bar − ).Lipid and water molecules were separately coupled to the thermostat.Following GROMACS recommendations for the CHARMM-36 all-atom forcefield, a Verlet cut-off scheme on grid cells was used with a distance of 1.2 nm,and non-bonded interactions cut-offs (Van der Waals and Coulombic) werealso set to 1.2 nm. Fast smooth Particle-Mesh Ewald electrostatics wasselected for handling the Coulombic interactions, with a grid spacing of 4 nm.A standard cut-off scheme with a force-switch smooth modifier at 1.0 nm wasapplied to the Van der Waals interactions. We did not account for long rangeenergy and pressure corrections, and constrained all the hydrogen bonds ofthe system using the LINCS algorithm.Molecular dynamics production were run for a minimum of 50 ns. Whentemperature treatments occurred, the systems were simulated for another50 ns to allow for the bilayer to reach equilibrium. Regardless of the lengthof the simulation, the analysis were always performed on the last 25 ns ofthe simulations. The areas per lipid of the systems were measured using two differentmethods: (i) measure by projection of the bilayer in the XY plane of thebox, noted A pH , and (ii) measure by meshing of the water-lipid interfaces ofeach leaflet of the bilayer, written A mH . The projected area per lipid wasmeasured directly in Gromacs by measuring the area of the simulation boxin the XY plane and by dividing it by the number of lipids per leaflet. Themeshed area per lipid was computed using Ovito 2.9. To do so, the surface ofthe water blocks at the water-lipid interface was meshed using a probe sphereradius of 6 and a smoothing level of 30 after removing the lipid molecules.The area of the meshing was then divided by the total number of lipids inthe bilayer.The amplitude and the period of the corrugations were measured usingthe meshing of the water-lipid interface collected for the measurement of thearea per lipid. The position of the vertices of the meshing were interpolatedusing SciPy [36] on a (200,200) uniform grid of the size of the simulation4ox, and the height of these 200 ×
200 points were analysed. This grid wasdirectly used to generate the height maps. The amplitude of the corrugations,represented by the root mean square (RMS) height of the points on the grid,were calculated by subtracting the mean and then measuring the standarddeviation of the height σ ( h ). The period of the corrugations were calculatedby computing the azimuthal power spectrum density (PSD) of the points ofthe grid and determining the frequency of the highest peak of the PSD.To measure the local chain tilt direction, the positions of the atoms ofthe tails of the lipids were extracted from the simulation files using MDAnal-ysis [37, 38]. The vectors from the first carbon atom of each tail to all thecarbon atom of their respective tails were computed, before calculating themean vector. The vector field was then created by binning the membrane in n × m squares and averaging over time the tilt direction vectors found in eachsquare ( X i , Y i ). As a consequence, the vector field account for the averagedirection of the tails found at this location on the membrane over time andnot to the average direction of a lipid over time.The enthalpy of the systems were extracted from the simulation usingdirectly the tool provided with Gromacs, gmx energy . The systems usedstrictly had the exact same number of atoms, for both lipids and watermolecules, to prevent changes due to the system composition. The (tilted)gel - fluid phase transition enthalpy was measured using the method from [39].Briefly, the total enthalpy of a system simulated at different temperaturesover a wide range, here 283 to 358 K, was collected. The effects of thetemperature increase on the enthalpy, besides any phase transition, wereremoved from the measurement by subtracting the affine baseline measuredin each phase. The values were then divided by the total number of lipidsin the system, and the transition enthalpy was calculated by integratingthe difference between the fitted gel and fluid baseline over the range oftemperature at which the transition occurs (here 308 to 318 K, accountingfor Gromacs accuracy in setting a temperature). For the tilted to disorderedgel phase transition, the enthalpy was measured by bringing the two systemsat the exact same temperature. The difference in enthalpy measured betweenthese two systems was then divided by the total number of lipids in thesystem. 5 . Results We considered DPPC bilayers prepared with CHARMM-GUI at temper-ature T = 288K whose size is ranging from 32 to 512 lipids. For smallbilayers made of 64 lipids or less, we observed the formation of a smooth L β tilted gel phase, in agreement with the numerous experiments and simula-tions reported in the literature [40]. The geometry of these bilayers remainedstable over time, even after a 50 ns simulation. When larger systems ( > cf Figure 1).
Figure 1: DPPC bilayer made of 212 lipids simulated at 288 K right after construction andequilibration. 2.5 nm thick slices of the system along each axis are shown on the bottom,highlighting the corrugations of the membrane. Top and bottom lipid leaflets are colouredrespectively in cyan and orange while the phosphorus atoms are shown as silver beads.Hydrogen atoms were removed from the screenshot to improve readability. The PBC boxis displayed using the blue plain lines.
The corrugations takes place along both simulations box axis ( x and y ).As already observed and investigated in details for DMPC by Khakbaz andKlauda [20], the molecular configurations of the lipids do not appear uniformalong the corrugations, and the typical stretched tails of the gel phase seemto turn into the typical disordered tail configurations of the fluid phase in the6hin portions of the corrugations where interdigitation happens. These varia-tions in configuration can be highlighted by computing and mapping the localsegment order parameter S mol of each atom of each lipid ( cf Figure 2(a)). -0.50.5 1.00.5 (a) − − − Relative height (Å) -15 -10 -5 0 5 10 (b)
Figure 2: (Left) Topography of the upper leaflet of (Top) 256 and (Bottom) 64 lipidbilayers, highlighting the corrugations observed in the membrane. Both contour plotshares the same color code and contour line scale, in ˚A. The white dashed line show thesystem before replication. (Right) Representation of slices of the lipid membranes on theleft, with a color code on the tail atoms corresponding to their order parameters.
We systematically characterized the corrugation and probed its evolutionwith system size. To quantify it in a more specific way, we investigated thetopographic elevation function for each leaflet (see Figure 2(b)) from whichwe defined the amplitude of the corrugation as the root mean square (RMS)height h RMS of the two water-lipid interfaces of the bilayers. The heights of7he interfaces were obtained by meshing the water surface and removing themean height of each leaflet. In these circumstances, h RMS is equal to thestandard deviation of the heights σ ( h ). The leaflet corrugation amplitudesare identical for both leaflets even though cross section pictures might suggestotherwise (Figure 3(a)). The corrugation amplitude increases with systemsize but saturates for large systems (Figure 3(b)). The longitudinal periodof the corrugations always coincide with the periodic boundary conditions(PBC, see Figure 3(c)). The power spectrum densities (PSD) used to measurethe periods are given in the Supplementary Materials. (a) R M S A m p li t ud e ( Å ) Leaflet:
TopBottom (b) P e r i od ( Å ) Leaflet:
TopLeaflet L x ( Å ) (c) Figure 3: Effect of the size of the simulation box on the geometry of bilayers. (a) 2.5 nmthick slices of the systems made with different amount of lipid molecules and simulatedat 288 K. Measurement of the (b) amplitude, or h RMS and (c) period of the corrugationsobserved in different systems versus the number of lipids. The evolution of the size of thebox is shown on top of the period. Dashed lines are respectively the cubic root and squareroot fits of the data.
To better quantify the nature of the corrugation, we investigated the areaper lipid using two methods: i) projection in the ( XY ) plane of the box, i.ethe area A pH of the box divided by the number of lipid; ii) meshing of the8
200 400 600 α ( % ) A r ea p e r li p i d ( Å ²) Projected A pH Meshed A mH Figure 4: (Top) Differences between the area per lipid measured by surface meshing and byprojection in the XY plane α (in %) and (Bottom) evolution of each type of area per lipidfor different system sizes. Blue dashed line on the bottom graph is the average projectedarea per lipid. water lipid interface that enables the determination of the interfacial area A mH . We found that the difference between interfacial and projected area perlipid α = A mH − A pH A pH (1)provides an excellent characterization of the corrugated phase (Figure 4).Indeed, α is three times larger for corrugated systems ( α = 14 ± α = 5 ± cf Supplementary materials).The local chain tilt direction in the XY plane was also determined andmapped (Figure 5(a)). As expected, tilts are uniform in the L β state andpresent a random short-range correlation in the fluid L α state. This wasconfirmed by the distribution of the angles with respect to the X-axis, as afunction of the height on the leaflet. In the corrugated systems, the localchain tilts displays long-range variations. The angle distribution of the cor-rugated systems shows that small heights have random orientation similar tothe fluid, while small heights are subject to less variations.We also probed the sensitivity of the corrugation to modifications of spacegroup. A simple change from cubic to rectangular simulation boxes with L x = 2 L y shows commensurate corrugations periods. Moreover we carryout simulations in hexagonal and monoclinic simulation boxes counting 128lipids. Both of these systems were found to be corrugated, with modulationvectors directed along the PBC/crystallographic directions (Figure 6).To release residual stress on the bilayer configuration that could havebeen brought by the semi-isotropic barostat, an anisotropic barostat was ap-plied to the systems, with 1 bar and a 4 . × − bar − compressibility setalong each axes, and the pressure crossed terms set to 0 for both pressure andcompressibility. We found that releasing residual anisotropic stress did notmodify the structure of the bilayers. After this simulation run, the systemwas still found in a corrugated state, as shown in Figure 7, with a relativearea increase α of 11.9 ± L dβ , to make a clear distinction betweenthese configurations and the ripple phase P β that has only be experimentallyreported above the pre-transition temperature.Finally one can wonder whether the appearance of the corrugation isrestricted to DPPC or DMPC. To answer this question, we probed the effectof the tail and the head groups by considering the longer-tailed DSPC and theethanolamine-based DPPE. As shown in Figure 8, the large DSPC systemswere found in the disordered gel state while the DPPE systems remained in10 − − Relative height (Å) -15 -10 -5 0 5 10 (a) A ng l e ( r a d ) A ng l e ( r a d ) −10 0 10 Relative height (Å) A ng l e ( r a d ) (b) Figure 5: (a) Orientation in the XY plane of the tails of each of the 256 lipids forminga DPPC bilayer (Top) in the disordered gel phase, (Middle) in the tilted gel phase and(Bottom) in the fluid phase. Only one leaflet of the membrane is shown for each system.(b) Respective scatter distribution of the angles between the tails of the lipids and theX-axis of the membrane, as a function of the height on the leaflet. The measurement wereperformed on both top and bottom leaflets, respectively colored in cyan and in orange. the expected homogeneous tilted gel state. The respective area differences ofthese systems are 12 ± ± a)(b) Figure 6: Topography of bilayers (a) made of 188 DPPC molecules in a simulation boxwith L x = 2 L y and (b) made of 128 DPPC molecules in an hexagonal simulation box,both simulated at 288 K.Figure 7: Slice of the 212 DPPC molecule system obtained after being simulated at 288 Kfor 50 ns with an anisotropic barostat. E H E A D P C H E A D D S P C D PP C D PP E
64 LIPIDS 212 LIPIDS P T A I L SS T A I L S Figure 8: Profiles of bilayers made of 64 or 212 lipids, either with DPPC, DSPC or DPPE,after a 50 ns simulation at 288 K. Only the DPPE membrane remained in the tilted gelphase when constructed with a large number of lipids.
The previous observations suggest that disordered L dβ and tilted L β gelstates are two competing states whose appearance seem to be correlated tothe system size for freshly thermalised systems generated with CHARMM-GUI. In what follows we probe the sensitivity of the stability of these phasesto different routes of thermal treatments.We focus our analysis on two system sizes: 64 lipids, which has beenfound in the L β phase, and 256 lipids which has been found in the L dβ .Both systems were subjected to the following thermal treatment: startingfrom T = 288 K , systems were annealed at T = 358K in the fluid phase.These systems are respectively named pc64-A and pc256-A, for annealing.They were then cooled down to T = 288K in two different ways: eitherwith a brutal fast cooling, named here quenched, or with a slow gradualcooling of 1K/ns that we denote gentle cooling. The quenched systems arenoted pc64-AQ and pc256-AQ (annealing-quenching), and the cooled systemscalled pc64-AC and pc256-AC (annealing-cooling).We first notice that whatever the thermal history, large systems werealways found in the disordered gel state as shown in Figure 9. Gentle coolingof small systems allowed them to recover the L β phase while quenching leadto L dβ phase. As in this latter case, either phases could be obtained depending13n the thermal history, we conclude that the disordered gel is metastable withrespect to the tilted gel for small systems. Furthermore, we can note thatboth the gently cooled and the quenched systems have a final projected areaper lipid A pH (respectively 50.1 ± ± ) close to the averageprojected area of 50 ˚A found for systems of all sizes, hence correcting theodd value of 52.3 ˚A found before temperature treatment ( cf Figure 4).As large systems of 256 lipids were never spontaneously found in the L β state, we decided to force them into this state by duplicating along both x and y directions the 64 L β system. The flat tilted duplicated system was foundto be stable on the simulation time scale, suggesting that L β could be eitherstable or metastable in large systems too. The nature of the relevant stablethermodynamic phase for large systems remains an open question, while oursimulations clearly favor the disordered gel state L dβ . Measurements of thedifference in area α can be found in the Supplementary Materials. Since the small systems made of 64 lipids can be prepared and controlledto reach all the observed phases, we used them to compare the energeticproperties of the respective gel phase. Having in mind the idea of a com-plex underlying potential energy landscape composed of several minima lo-cated at different energy levels, we probe how far energetically the disor-dered gel phase stand from the tilted gel phase minimum. Therefore weperformed an energy minimization using conjugate gradient to remove ther-mal fluctuations on the tilted and disordered gel configurations. We found∆ E P = E P (disordered) − E P (tilted) = 64kJ/mol, meaning that if both statesare metastable, the tilted system is the most stable one. In addition, we alsocomputed the difference in enthalpy between the disordered and the tiltedthermalised gel phases. To this aim, we used the same initial states, afterquenching or gentle cooling but without energy minimisation. Since the twosystems share the same atomic compositions (same number of lipid, watermolecules and constraints), the difference in enthalpy (kinetic and potentialenergy according to the force field, plus the P V contribution) at a giventemperature should be characteristic of the enthalpy difference between thetwo states. We found an enthalpy difference of 12 ± C H A R MM - G U I CoolingQuenchingAnnealingThermalisationConstruction64
Lipids
Lipids
Lipids pc64-ACpc64-AQpc64-Tpc256d-Tpc256-AQ / pc256d-AQpc256-Tpc256-AC / pc256d-AC
Figure 9: Synopsis of the different thermal trajectories, showing the systems obtainedfor different cooling rates as well as the systems they originated from. Only the resultsfrom the 256 DPPC system obtained by replication of the 64 lipid system are shown forthe quenching and cooling experiments. Results for the 256 DPPC system constructedby CHARMM-GUI, as well as the intermediate systems, are shown in the SupplementaryMaterials. found a difference in enthalpy of 3 . ± . L dβ phase (see Supplementary materials).The enthalpy of the L dβ −→ L α transition was also measured in our smallsystems. This was performed by simulating the 64 DPPC system at temper-atures ranging from 283 to 358 K and by removing the change in enthalpydue to the change in temperature (data shown in Supplementary Materials).15he transition enthalpy was found equal to 27.3 kJ/mol, which is also com-parable with the experimental reported values circa
4. Discussion
A careful inspection of the corrugations shows that the topographic mod-ulation is imposed by the periodicity of the simulation box. The corrugationshows finite system size dependence, as its amplitude increases up to sizes ofthe order of L x = 8 nm where the modulation saturates to a value of 7 ± ± α , introduced in equation (1) can be taken as a relevantorder parameter for the transition between the L β and the disordered gelphases L dβ . The latter being reminiscent from P β ripple phase, α could beseen as a critical parameter to investigate the existence of the ripple phaseand discriminate it from the tilted gel or from the fluid phase.However, unlike experiments, the numerical instability occurs along twoorthogonal directions, or along the hexagonal axes. Non-square boxes fail toselect only one modulation direction. We nevertheless think that the numer-ical corrugation instability is related to the experimental ripple instability,as also suggested by the dependence of the presence of corrugation to thechemical nature of the heads and tails of the lipids. Indeed, experimentshave shown that the ripple phase is specific from the phosphocholine (PC)lipids [40, 42].Another striking observation is the insensitivity of large systems (256lipids) to thermal treatment which have been systematically ended in the dis-ordered gel state. On the opposite small systems can alternate between both16hases. However, when thermalised at low temperatures from CHARMM-GUI or slowly quenched, they end up preferentially into the tilted gel state,suggesting that in the range of temperatures investigated the L β state isthermodynamically favored. By contrast the L β state is never selected spon-taneously by larger systems. For the latter, MD suggests that the disorderedgel state L dβ is the most preferred state for all temperatures below melting.Metastability and kinetic effects are certainly significant and may hide totrue nature of the stable phase.Assuming that the difference between tilted and disordered gels has some-thing to do with the pretransition, we found an enthalpy difference 3 timeslarger than the premelting latent heat at 288K, and of the same order of mag-nitude at 305K, close to the observed experimental transition. This pointstowards the relevance of disordered gel state L dβ as a ripple state analogue.Moreover the higher enthalpy of the disordered gel with respect to the tiltedgel is consistent with a L β −→ L dβ −→ L α sequence of transitions as tem-perature is increased.The reason of the outcome of a ripple instability below the melting tem-perature is nothing but obvious. Our simulations point out to a competitionbetween an homogeneous tilted state and an inhomogeneous corrugated state.The transition between these states is discontinuous. The corrugated state isnot very tilted, and partially melted, or disordered, and interdigitated. Wesuggest now a possible mechanism explaining the observed situation. Thetilted phase can be understood as the result of a frustration between lipidheadgroups which try to increase their exposure to water in the interfaceregion, lipid chains which try to reach an optimal packing density as a resultof cohesive forces, and chain stiffness for which the introduction of gauche dihedral angles is unfavorable at low temperatures. Tilt allows lipid to op-timize simultaneously those three constraints. On the other hand, melting alipid chain enables the release of the constraint acting on the chain stiffness,and makes it possible to increase the hydration free-energy by reducing themembrane thickness. Below melting temperature, thermodynamics makes itunfavorable to melt all the lipid molecules. However, some local disorderingof the lipids may still be favorable, increasing the hydration of the headgroupswithout need of spending too much energy in melting the chains.Based upon those considerations, we designed a simple one dimensionallipid chain model that supports the idea that in an temperature range justbelow melting, the homogeneous tilted state energetically unstable with re-spect to a local corrugation of the bilayer, see Figure 10. Details of the17 igure 10: Schematic representation of the mechanical instability occurring in a simple onedimensional lipid chain model, where the tilted state has a higher energy than a locallydisordered state. The gain in energy originates mostly from the hydration term. Bluesolid rods stand for gel state, red dashed rods for fluid state. parameterisation of the model and further results are presented in the Sup-plementary Information section. We conclude that the thickness modulationmight indeed be caused by a subtle interplay between headgroup hydration,hydrophobic chain packing, trans-gauche isomerisation and tilt elasticity en-ergy terms.
5. Conclusion
We have successfully demonstrated in this work how the size of the simu-lation box influences the ripple-like instability in a PC membrane simulatedwith the Charmm36 force field at low temperature, where it is usually ex-pected to be in the tilted gel L β phase. This unexpected organisation, whichwe called the disordered gel phase L dβ , does not appear in small systems,which is consistent with the results from Khakbaz and Klauda [20] as wellas with our results in a previously published paper [21]. The energy andgeometry analyses demonstrated that this disordered gel phase has a lot incommon with the P β phase. Furthermore, this instability was not observedwith PE lipids, in agreement with experimental findings. For small systems,we found ways of preparing the system in either tilted gel or disordered gelstates by acting on the thermal treatment. More work is needed to deter-mine whether the Charmm36 force-field can describe a one-dimensional P β spatial thickness modulation, with the right periodicity, and at which tem-peratures. Our work suggests that simulations will have to be guided to thedesired structure. Finally, we conclude that the ripple instability looks likea generic mechanism adopted by the phosphatidylcholines lipids to increasethe headgroup hydration while still satisfying the packing constraints, at theexpense of a mild cost in disordering/melting a small fraction of the chains.18 cknowledgements V. W. warmly thanks T.E. de Oliveira for helping with simulations set-up. A.G. acknowledges partial support from the Investissements d’Avenirprogram “D´eveloppement de l’ ´Economie Num´erique” through the SMICEproject. The authors gratefully acknowledge support from the high per-formance cluster (HPC) Equip@Meso from the University of Strasbourg,through grants no. G2019A131C and G2020A140C.
References [1] O. G. Mouritsen, Life-as a matter of fat: the emerging science oflipidomics, Springer, 2005.[2] R. Dimova, C. M. Marques (Eds.), The Giant Vesicle Book, 1st ed.,CRC Press, Taylor and Francis, 2019.[3] D. Marsh, Handbook of Lipid Bilayers, 2nd ed., CRC Press, Boca Raton,2013.[4] G. Cevc, D. Marsh, Phospholipid Bilayers. Physical Principles and Mod-els, John Wiley & Sons, New-York, 1987.[5] S. Mabrey, J. M. Sturtevant, Investigation of phase transitions in lipidsand lipid mixtures by high sensitivity differential scanning calorimetry,Proceedings of the Natural Academy of Sciences USA 73 (1976) 3862–3866.[6] T. Heimburg, Thermal Biophysics of Membranes, Wiley-VCH, 2007.[7] A. Tardieu, V. Luzzati, F. C. Reman, Structure and polymorphismof the hydrocarbon chains of lipids: A study of lecithin-water phases,Journal of Molecular Biology 75 (1973) 719–733.[8] W. J. Sun, S. Tristram-Nagle, R. M. Suter, J. F. Nagle, Structure of theripple phase in lecithin bilayers., Proceedings of the National Academyof Sciences 93 (1996) 7008–7012. doi: .[9] R. N. A. H. Lewis, N. Mak, R. N. McElhaney, A differential scan-ning calorimetric study of the thermotropic phase behavior of model19embranes composed of phosphatidylcholines containing linear satu-rated fatty acyl chains, Biochemistry 26 (1987) 6118–6126. doi: .[10] T. Heimburg, A model for the lipid pretransition: coupling of rippleformation with the chain-melting transition, Biophysical Journal 78(2000) 1154–1165.[11] S. Doniach, A thermodynamic model for the monoclinic (ripple) phaseof hydrated phospholipid bilayers, The Journal of Chemical Physics 70(1979) 4587–4596. doi: .[12] T. C. Lubensky, F. C. MacKintosh, Theory of “ripple” phases of lipidbilayers, Physical Review Letters 71 (1993) 1565–1568. doi: .[13] J.-B. Fournier, Coupling between membrane tilt-difference and di-lation: A new “ripple” instability and multiple crystalline inclusionsphases, Europhysics Letters (EPL) 43 (1998) 725–730. doi: .[14] C. Misbah, J. Duplat, B. Houchmandzadeh, Transition to ripple phasesin hydrated amphiphiles, Physical Review Letters 80 (1998) 4598–4601.doi: .[15] K. Sengupta, V. A. Raghunathan, Y. Hatwalne, Role of tilt order inthe asymmetric ripple phase of phospholipid bilayers, Physical ReviewLetters 87 (2001). doi: .[16] K. Akabori, J. F. Nagle, Structure of the dmpc lipid bilayer ripple phase,Soft Matter 11 (2015) 918–926.[17] A. H. de Vries, S. Yefimov, A. E. Mark, S. J. Marrink, Molecular struc-ture of the lecithin ripple phase, Proceedings of the National Academyof Sciences 102 (2005) 5392–5396. doi: .[18] O. Lenz, F. Schmid, Structure of symmetric and asymmetric ”ripple”phases in lipid bilayers, Physical Review Letters 98 (2007).2019] A. Debnath, F. M. Thakkar, V. Kumaran, P. K. Maiti, K. G. Ayappa,Laterally structured ripple and square phases with one and two dimen-sional thickness modulations in a model bilayer system, Soft Matter 10(2014) 7630–7637.[20] P. Khakbaz, J. B. Klauda, Investigation of phase transitions of satu-rated phosphocholine lipid bilayers via molecular dynamics simulations,Biochimica et Biophysica Acta (BBA) - Biomembranes 1860 (2018) 1489– 1501.[21] V. Walter, C. Ruscher, C. M. Marques, O. Benzerara, F. Thalmann, Amachine learning study of the two states model for lipid bilayer phasetransitions, Phys. Chem. Chem. Phys. (2020) –.[22] S. Jo, T. Kim, W. Im, Automated builder and database of pro-tein/membrane complexes for molecular dynamics simulations, PLoSONE 2 (2007) 880.[23] S. Jo, T. Kim, V. G. Iyer, W. Im, Charmm-gui: A web-based graphicaluser interface for charmm, Journal of Computational Chemistry 29(2008) 1859–1865.[24] S. Jo, J. B. Lim, J. B. Klauda, W. Im, Charmm-gui membrane builderfor mixed bilayers and its application to yeast membranes, BiophysicalJournal 97 (2009) 50–58.[25] E. L. Wu, X. Cheng, S. Jo, H. Rui, H. K. Song, E. M. Davila-Contreras,Y. Qi, J. Lee, V. Monje-Galvan, R. M. Venable, J. B. Klauda, W. Im,Charmm-gui membrane builder toward realistic biological membranesimulations, Journal of Chemical Theory and Computation 35 (2014)1997–2004.[26] H. Berendsen, D. van der Spoel, R. van Drunen, Gromacs: A message-passing parallel molecular dynamics implementation, Computer PhysicsCommunications (1995).[27] M. J. Abraham, T. Murtola, R. Schulz, S. P´all, J. C. Smith, B. Hess,E. Lindahl, Gromacs: High performance molecular simulations throughmulti-level parallelism from laptops to supercomputers, SoftwareX(2015). 2128] R. B. Best, X. Zhu, J. Shim, P. E. M. Lopes, J. Mittal, M. Feig, A. D.MacKerell Jr., Optimization of the additive charmm all-atom proteinforce field targeting improved sampling of the backbone phi, psi andside-chain khi1 and khi2 dihedral angles, Journal of Chemical Theoryand Computation 8 (2012) 3257–3273.[29] B. R. Brooks, C. L. Brooks III, A. D. MacKerell Jr, L. Nilsson, R. J.Petrella, B. Roux, Y. Won, G. Archontis, C. Bartels, S. Boresch,A. Caflisch, L. Caves, Q. Cui, A. R. Dinner, M. Feig, S. Fischer, J. Gao,M. Hodoscek, W. Im, K. Kuczera, T. Lazaridis, J. Ma, V. Ovchinnikov,E. Paci, R. W. Pastor, C. B. Post, J. Z. Pu, M. Schaefer, B. Tidor, R. M.Venable, H. L. Woodcock, X. Wu, W. Yang, D. M. York, M. Karplus,Charmm: The biomolecular simulation program, Journal of Computa-tional Chemistry 30 (2009) 1545–1614.[30] J. Lee, X. Cheng, J. M. Swails, M. S. Yeom, P. K. Eastman, J. A.Lemkul, S. Wei, J. Buckner, J. C. Jeong, Y. Qi, S. Jo, V. S. Pande,D. A. Case, C. L. Brooks III, A. D. MacKerell Jr, J. B. Klauda, W. Im,Charmm-gui input generator for namd, gromacs, amber, openmm, andcharmm/openmm simulations using the charmm36 additive force field,Journal of Chemical Theory and Computation 12 (2016) 405–413.[31] R. W. Hockney, S. P. Goel, J. W. Eastwood, Quiet high-resolutioncomputer models of a plasma, Journal of Computational Physics 14(1974) 148–158.[32] S. Nos´e, A molecular dynamics method for simulations in the canonicalensemble, Molecular Physics 52 (1984) 255–268.[33] W. G. Hoover, Canonical dynamics: Equilibrium phase-space distribu-tions, Physical Review A 31 (1985) 1695.[34] S. Nose, M. L. Klein, Constant pressure molecular dynamics for molec-ular systems, Molecular Physics 50 (1983) 1055–1076.[35] M. Parinello, A. Rahman, Polymorphic transitions in single crystals: Anew molecular dynamics method, Journal of Applied Physics 52 (1998)7182.[36] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy,D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright,22. J. van der Walt, M. Brett, J. Wilson, K. Jarrod Millman, N. May-orov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. Carey, ˙I. Polat,Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cim-rman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald,A. H. Ribeiro, F. Pedregosa, P. van Mulbregt, S. . . Contributors, Scipy1.0: Fundamental algorithms for scientific computing in python, NatureMethods 17 (2020) 261–272.[37] N. Michaud-Agrawal, E. J. Denning, T. B. Woolf, O. Beckstein, Md-analysis: A toolkit for the analysis of molecular dynamics simulations,Journal of Computational Chemistry 32 (2011) 2319–2327.[38] R. J. Gowers, M. Linke, J. Barnoud, T. J. E. Reddy, M. N. Melo, S. L.Seyler, D. L. Dotson, J. Domanski, S. Buchoux, I. M. Kenney, O. Beck-stein, Mdanalysis: A python package for the rapid analysis of moleculardynamics simulations, Proceedings of the 15th Python in Science Con-ference (2016) 98–105.[39] T. E. de Oliveira, F. Leonforte, L. Nicolas-Morgantini, A.-L. Fameau,B. Querleux, F. Thalmann, C. M. Marques, Fluid bilayer phase inaqueous mixtures of fatty alcohol and cationic surfactant, Phys. Rev.Research 2 (2020) 013075. doi: .[40] B. A. Cunningham, A.-D. Brown, D. H. Wolfe, W. P. Williams, A. Brain,Ripple phase formation in phosphatidylcholine: Effect of acyl chain rel-ative length, position, and unsaturation, Physical Review E 58 (1998).[41] A. H. de Vries, S. Yefimov, A. E. Mark, S. J. Marrink, Molecular struc-ture of the lecithin ripple phase, PNAS 102 (2005) 5392–5396.[42] J. Katsaras, S. Tristram-Nagle, Y. Liu, R. L. Headrick, E. Fontes, P. C.Mason, J. F. Nagle, Clarification of the ripple phase of lecithin bilayersusing fully hydrated, aligned samples, Physical Review E 61 (2000).23 ipple-like instability in the simulated gel phase of finitesize phosphocholine bilayers-Supplementary Materials-
Vivien Walter a , C´eline Ruscher b , Adrien Gola b , Carlos M. Marques b ,Olivier Benzerara b , Fabrice Thalmann b a Department of Chemistry, King’s College London, Britannia House, 7 Trinity Street,SE1 1DB, London, United Kingdom b Institut Charles Sadron, CNRS and University of Strasbourg, 23 rue du Loess, F-67034Strasbourg, FrancePreprint submitted to Elsevier February 9, 2021 a r X i v : . [ c ond - m a t . s o f t ] F e b ppendix A. System Composition The atomic composition of all systems used in this work is given in theTable A.1 below. Each DPPC molecule is made of 130 atoms, DSPC of 142atoms and DPPE of 121. The average number of water molecules per lipidis 163 ± Table A.1:
Atomic composition and size of the different simulation systems usedin this work.
Type System x , L ya (nm) L za (nm)DPPC
32 Lipids 4,160 15,570 2.91 22.5564 Lipids 7,020 31,869 4.24 21.5794 Lipids 12,220 46,590 5.12 21.65128 Lipids 16,640 63,897 5.98 21.71170 Lipids 22,100 84,315 6.87 21.77188 Lipids b x ) 11.09 19.74(L y ) 5.55212 Lipids 27,560 89,478 7.69 19.01256 Lipids 33,280 127,926 8,46 21.71256 Lipids c DSPC
64 Lipids 9,088 32,526 4.29 21.74212 Lipids 30,104 66,447 7.82 15.18
DPPE
64 Lipids 7,744 30,582 4.01 23.07212 Lipids 25,652 62,925 7.44 15.64 a Values measured on the last frame of the equilibrated system beforesimulation. b System constructed by replication of the system made of 94 lipids alongthe X axis for the L x = 2 L y system. c System constructed by replication of the system made of 64 lipids alongboth X and Y axes. 2 ppendix B. Corrugation formation and characterization
The topography of all the DPPC bilayers simulated at 288 K right afterconstruction is shown below. The figures are separated into top (Figure B.1)and bottom leaflets (Figure B.2). The complete sets of leaflets for the non-square systems, both rectangular (with L x = 2 L y , Figure B.3) and hexagonalor frustrated (Figure B.4). To facilitate the comparison between systems, allthe color codes and contour lines were scaled with respect to the system withthe largest amplitude. The color scale of the height is set as the distance fromthe mean height of the leaflet (0) toward the center of the bilayer (positivevalues, yellow/red). Negative height values are height further from the centerof the bilayer than the mean height (purple/blue).Figure B.7 is a rendering generated in Ovito 2.9 to illustrate the surfacemeshing performed to measure the meshed area per lipid A mH of the mem-brane. The mesh was done on the water molecules at the interface with thelipid bilayer, with a probe sphere radius of 6 and a smoothing level of 30.The meshed surface was fed to a two-dimensional Fast-Fourier Trasnform(FFT) routine, and the squared moduli of the Fourier coefficients (powerspectrum density PSD) plot as a function of the x and y components of thewave numbers k x , k y . Figures B.5, B.6 shows that no significant peak existexcept but the periodicity of the simulation box. The projected area per lipid A pH was simply calculated by multiplying together the size of the box alongthe X and Y axes and dividing by the number of lipid per leaflet. Theseareas were measured for all the square systems and gathered in Table B.2.The relative increase of area per lipid α from A pH to A mH was also measuredfor all these and shown in Figure B.8. For this Figure, the nomenclatureused is { tail }{ head }{ number of lipids } - { temperature treatment } . For exam-ple, pe64-T is a 64 DPPE (pe) bilayer thermalised (T) after construction ofthe system. The different temperature treatments used are thermalisation(T), annealing (A), quenching (Q) and gentle cooling (C) - AC indicates thatthe system was obtained first by annealing followed by a quenching.The efficiency of the relative increase of area per lipid α to be used as anorder parameter distinguishing between the tilted and disordered phase hasbeen compared to the usual order parameters of the carbon atom n of the tail k of each lipid, noted S mol ( k, n ), used to discriminate between the gel andfluid phases. The parameter S mol ( k, n ) was calculated using the followingformula: 3 mol ( k, n ) = (cid:28) θ k,n − (cid:29) , (B.1)with θ k,n the angle formed between the vector generated by the atoms n − n + 1 along the tail k chain and the vector orthogonal to the plane ofthe bilayer, here Z. The results are illustrated in Figure B.9. While the tailorder parameters could discriminate most of the disordered systems (94, 128,212 and 256) from the tilted ones (32, 64), the system made of 170 DPPCmolecules could not be undoubtedly told apart from neither the tilted northe disordered phases. Same results was observed by measuring the averagetail order parameter h S mol i of all atoms of both tails (cf. Figure B.10). Incontrast, the Figure B.8 shows that the relative increase of area α could easilyclassify the same system as a system in the disordered phase.4 a) 32 Lipids (b) 64 Lipids (c) 94 Lipids(d) 128 Lipids (e) 170 Lipids (f) 212 Lipids(g) 256 Lipids (h) 256 Lipids (replicated) (i) 512 Lipids Figure B.1: Top leaflets. a) 32 Lipids (b) 64 Lipids (c) 94 Lipids(d) 128 Lipids (e) 170 Lipids (f) 212 Lipids(g) 256 Lipids (h) 256 Lipids (replicated) (i) 512 Lipids Figure B.2: Bottom leaflets. a)(b) Figure B.3: Corrugation in a system where L x = 2 L y . (a) Top and (b) bottom leafletsare both shown. Color code shows in blue the thick portions of the membrane while redare the thin ones. igure B.4: Corrugation in a hexagonal system. (a) Top and (b) bottom leaflets are bothshown. Color code shows in blue the thick portions of the membrane while red are thethin ones. .00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 Frequency (Å −1 ) P o w e r D e n s i t y ( Å ) Leaflet:
TopBottom (a) 32 Lipids
Frequency (Å −1 ) −1 P o w e r D e n s i t y ( Å ) Leaflet:
TopBottom (b) 64 Lipids
Frequency (Å −1 ) P o w e r D e n s i t y ( Å ) Leaflet:
TopBottom (c) 94 Lipids
Frequency (Å −1 ) P o w e r D e n s i t y ( Å ) Leaflet:
TopBottom (d) 128 Lipids
Figure B.5: Power spectrum densities of the topography used to calculate the period ofthe corrugations. The frequency of the mean background value (0 ˚A − ) has been removedfrom the spectrum to improve the clarity. Plain lines are the spectrum along the X axis,while dashed lines are along the Y axis. Plain black lines are the average PSD of themembrane. .000 0.025 0.050 0.075 0.100 0.125 0.150 Frequency (Å −1 ) −1 P o w e r D e n s i t y ( Å ) Leaflet:
TopBottom (a) 170 Lipids
Frequency (Å −1 ) P o w e r D e n s i t y ( Å ) Leaflet:
TopBottom (b) 212 Lipids
Frequency (Å −1 ) −1 P o w e r D e n s i t y ( Å ) Leaflet:
TopBottom (c) 256 Lipids
Frequency (Å −1 ) P o w e r D e n s i t y ( Å ) Leaflet:
TopBottom (d) 512 Lipids
Figure B.6: (Continuation of Figure B.5). a)(b) Figure B.7: Rendering of DPPC bilayers made with 256 lipids (a) generated by CHARMM-GUI and simulated at 288 K or (b) produced by replication of a 64 DPPC bilayer found inthe tilted gel phase; with the surface meshing generated on the water-membrane interface(using water molecules as the reference) to calculate the meshed area per lipid. Therespective area difference of the systems were calculated at 13 ± ± able B.2: All mean values of area per lipid measured in this work on thesystems simulated at 288 K, either by projection of the area in the XY planeof the system box ( A pH ) or by meshing of the surface of the bilayer using Ovito( A mH ). Error bars are the standard deviation of the measured values. T h e r m a li s e d ( T ) C oo li n g ( A C ) Q u e n c h i n g ( A Q ) T y p e S y s t e m A H p A H m A H p A H m A H p A H m D PP C L i p i d s . ± . . ± . ---- ( p c ) L i p i d s . ± . . ± . . ± . . ± . . ± . . ± . L i p i d s . ± . . ± . ---- L i p i d s . ± . . ± . ---- L i p i d s . ± . . ± . ---- L i p i d s . ± . . ± . ---- L i p i d s . ± . . ± . . ± . . ± . . ± . . ± . L i p i d s a . ± . . ± . . ± . . ± . . ± . . ± . L i p i d s . ± . . ± . ---- D S P C L i p i d s . ± . . ± . ---- ( s c ) L i p i d s . ± . . ± . ---- D PP E L i p i d s . ± . . ± . ---- ( p e ) L i p i d s . ± . . ± . ---- a System constructed by replication of the system made of 64 lipids.12 igure B.8: Distribution of the mean α of all the square systems simulated at 288 K. Theblue bands highlights the two ranges of values of α in which the tilted and disorderedphases were observed. All systems with a α <
10 % were visually classified as tilted gelphase, while all systems with α ≥
10 % were classified as disordered gel. Error bars arethe standard deviations of each distribution. Atom number O r d e r p a r a m e t e r Figure B.9: Tail order parameters of the carbon atoms from the DPPC tails measured fordifferent systems made of increasing amount of lipids per bilayers. Plain and dashed linesare respectively the tails sn1 and sn2 of the lipids. While generally the system visuallyclassified as disordered have a very order parameter, occasionally some systems were foundwith an order parameter close to the one of the tilted gel phase - here 170 DPPC. < S m o l > Figure B.10: Mean order parameters of the carbon atoms from both DPPC tails measuredin Figure B.9 for different systems made of increasing amount of lipids per bilayers. Errorbars are the standard deviation over all the atom of both chains. ppendix C. Influence of the thermal history In this section are presented the screenshots of the systems that were notused to prepare the figures in the
Influence of the thermal history section ofthe main text: the systems made of 256 DPPC molecules but obtained viaCHARMM-GUI instead of replication of the small systems (Figure C.11),all the systems before the thermalisation (Figure C.12) or in the fluid phaseafter annealing (Figure C.13).The increase of area per lipid α of all the systems obtained during thedifferent temperature treatments were gathered in Figure B.8. Similarly tothe previous section, we measured the tail order parameters S mol ( k, n ) forthe different systems made of 64 DPPC molecules but obtained via differenttemperature treatment (Figure C.14). The tail order parameters confirmedthat the gently cooled system was in a state similar to the one before an-nealing, noted gel here, but drastically different from the one obtained afterquenching. The latter was found different from both gel and fluid phases,which highlight its nature of being of a new phase. (a) Gentle cooling(b) Quenching Figure C.11: Screenshots of the bilayers with 256 DPPC molecules constructed usingCHARMM-GUI and simulated at 288 K, followed by an annealing over the gel-fluid melt-ing temperature before cooled down to 288 K, either with a (a) gentle cooling or a (b)quenching. a) DPPC 64 - CHARMM-GUI (b) DPPC 64 - Equilibration(c) DPPC 256 - CHARMM-GUI(d) DPPC 256 - Equilibration(e) DPPC 256 - Replication(f) DPPC 256 - Equilibration after replication Figure C.12: Screenshots of the DPPC bilayers constructed and used in the experimentson temperature treatments, either using CHARMM-GUI or by replication of a 64 DPPCbilayer simulated at 288 K. All systems are shown before and after a NVT and NPT shortequilibration at 288 K, but always before the 50 ns thermalisation at 288 K. a) DPPC 64(b) DPPC 256(c) DPPC 256 - Replication Figure C.13: Screenshots of the DPPC bilayers used in the experiments on temperaturetreatments after annealing at 358 K. In this work, all these systems were classified as beingin the fluid phase. Atom number k S m o l ( k ) System: gelfluid quenchingcooling
Figure C.14: Average order parameters of the carbon atoms from the DPPC tails measuredin four 64 lipids bilayer systems: the initial system at 288 K in the gel phase (blue), thesame system in the fluid phase after annealing at 358 K (red), and the result of bothgentle cooling (gray) or quenching (orange) of the annealed system back to 288 K. Whilethe system gently cooled goes back the configuration of the initial system in the gel phase,the quenched system reaches another phase: the disordered gel phase. Plain and dashedlines correspond respectively to the sn1 and sn2 tails of the lipids. ppendix D. Thermodynamics of tilted and disordered states In this section, we present the detailed results obtained and gathered towrite the
Thermodynamics of tilted and disordered states section of the maintext.As mentioned in the main text, we compute the energy minimizationusing the conjugate gradient method implemented in Gromacs. Startingfrom tilted and disordered configurations, we monitor the energy evolutionwith time as shown in Figure D.16 and we find that V d − V t = 64kJ/mol.The corresponding inherent structures are depicted in Figure D.15.In order to compare the enthalpy of the systems in the tilted and dis-ordered states, it was obviously critical to force a system in both phasesusing different temperature treatments. We picked the system at 64 DPPCmolecules, where all phases were observed. However, while trying to obtain adisordered system at 305 K to compare the enthalpy of the disordered phaseto the one of the real ripple phase, a quenching from 358 to 305 K (∆ T =53) lead to a tilted gel phase. We made several attempts, and found that, re-gardless of the final temperature after quenching, a difference of temperature∆ T of at least 70 K is required to force the system in the disordered phase.The difference attempts and their results are summarised in Figure D.17.The enthalpy of the (tilted) gel to fluid phase transition was measuredby collecting the total enthalpy of a 64 DPPC molecules system after simu-lations at temperature ranging from 283 to 358 K (cf. Figure D.18). Aftersubtraction of the baseline (here set to the enthalpy of the fluid systems)and dividing by the total number of lipids, the transition enthalpy was mea-sured by integrating the enthalpy between the gel line and the fluid line, overthe range where the transition occurs (here 308 to 318 K). By doing so, wemeasured a difference of enthalpy ∆ H = 27.3 kJ/mol.20 a) From gently cooled system(b) From quenched system Figure D.15: Screenshots of the DPPC 64 bilayers obtained after energy minimisationusing a conjugate gradient algorithm (a) from a gently cooling and found initially in atilted gel phase, or (b) from a quenching and found initially in a disordered gel phase.Despite the energy minimisation, both systems remain in their initial phase.
500 1000 1500 2000
Time (ps) V d - V t ( kJ / m o l )
100 1100 2100−36000−35500−35000−34500 V ( kJ / m o l ) DisorderedTilted
Figure D.16: Evolution of potential energies of the disordered system V d and of the tiltedsystem V t during an energy minimisation simulation using a conjugate gradient integrator.The inset shows the evolution of the difference in potential energy between both systemover time. igure D.17: Visualisation of the different DPPC 64 systems obtained via quenching fordifferent starting and final temperatures. To force the small system into a disorderedphase, a ∆ T of at least 70 K is required. igure D.18: Measurement of the enthalpy of a 64 DPPC system at different temperaturesranging from 288 to 358 K to compute the ∆ H of the gel-fluid transition. The baselinecalculated on the systems from 318 to 358 K was subtracted from all points, and the ∆ H was calculated from the gray zone corresponding to the jump from one phase to another.Here ∆ H = 27.3 kJ/mol. ppendix E. A simple mechanical model that could explain theobserved instability We briefly introduce a one-dimensional minimal model displaying an in-teresting thickness modulation instability, that supports the view that thecorrugation observed in the simulations could be explained by a competitionbetween hydration energy, tilt elasticity and partial melting of the lipids.System parameterization
The system consists in a chain comprising a number N ∼
10 rods, eachone representing a single lipid. Each rod is given by its two end points (Fig-ure E.19). The lower end-point has coordinates ( U i ,
0) and remains at themid-plane position. Only one leaflet is considered in the model, the otherleaflet being assumed to follow the same behavior by symmetry. In addition,each rod possesses an internal state reminiscent from the Doniach two-statesapproach of the melting transition [1, 2]. Our variables are • the rods upper positions defined as ( X i = Lx i , Z i ), i ∈ [0 , N − x i ∈ [0 , L -periodicity assumed for X i , 1-periodicity for x i . • the rods lower positions ( U i = Lu i , u i ∈ [0 , L -periodicity assumedfor U i , 1-periodicity for u i . • the rod state S i = 0 ,
1, where 0 is gel , 1 is fluid , melted or disordered .We also define the periodic difference a (cid:5) b = ( a − b − b a − b + c ) casu-ally obtained by the python a-b -np.floor(a-b+0.5) macro. The periodicdifference is an implementation of the periodic boundary conditions minimalimage convention. Shift indices modulo N are denoted with brackets, e.g. [ i + 1] ≡ ( i + 1)[ N ], with [ i + 1] ∈ [0 , N − rules • the system tries to optimize its exposure to the solvent. Each beadattempts to reach separately the optimal ”area” per lipid (here lengthper lipid) A . The cost associated with it is related to the contourlength of the rod upper beads H solvent = C N X i =1 (cid:2) L ( x [ i +1] (cid:5) x i ) + ( Z [ i +1] − Z i ) − A (cid:3) (E.1)25 igure E.19: Sketch of the mechanical model. Blue solid rod: gel state, red dashed rod:fluid or disordered state. Also represented are the four main contributions to the energycost (hydration, packing, contour length and tilt). the ”volume” (here area) of a lipid rod cannot depart too much froma target value V , whichever state the lipid occupies. This is mainlyan incompressibility constraint imposed by the hydrophobic tails. Ifone takes the 1D lipid area the same as a trapezoid (approximate buthopefully accurate expression), the cost is14 LZ i ( x [ i +1] (cid:5) x [ i − + u [ i +1] (cid:5) u [ i − ) (E.2)then the corresponding cost function reads H packing = C N X i =1 (cid:2) L ( x [ i +1] (cid:5) x [ i − + u [ i +1] (cid:5) u [ i − ) Z i − V (cid:3) (E.3) • the lipids in the gel state behave as stiff rods with well-defined contourlength L c , mimicking an all- trans chain conformation. Therefore a statedependent energy term is introduced H all − trans = C N X i =1 S i (cid:2) L ( x i (cid:5) u i ) + Z i − L c (cid:3) (E.4) • the lipids in the fluid state have a restricted tilt, with an elastic coef-ficient of the order of magnitude of the bending modulus. Meanwhile,tilt restriction is necessary to frustrate the lipids in the gel state. With-out it, the area per lipid in the gel state could reach the optimal value A at the expense of an excessive tilt, obviously not seen experimen-tally nor numerically. Therefore, we introduce irrespective of the state,a tilt elasticity term H tilt restriction = C N X i =1 (cid:2) ( x i (cid:5) u i ) (cid:3) (E.5) • finally, though not strictly required, a repulsion term between end-tailsimproves the look of the resulting solutions. H cosmetic = C N X i =1 ( u [ i +1] (cid:5) u i ) (E.6)27 each disordered, or melted lipid, brings an additional thermodynamicfree-energy contribution, irrespective of the position and orientation. H melting = ε ( T ) N X i =1 (1 − S i ) (E.7)State mismatch energy contributions could also easily be introduced. • a global chain tension term can be introduced as usual H tension = − F e L (E.8)The control parameters below can be set arbitrarily. Our choice hereis • A optimal interfacial area per lipid, close its fluid phase value, here0.8 • L c optimal contour length in all-trans configuration, only in the gelphase, here 2. • V optimal chain packing volume, cohesion forces, here 0.8*1.5 = 1.2 • F e optional traction or tension force, here 0 • ε ( T ) Gibbs free-energy of melting, sends the fluid phase above the gelphase, depends linearly on temperature, vanishing around T m , here 0.08per bead. • Elastic parameters – C conjugated to the real area, here 1. – C conjugated to the packing volume; here 1. – C conjugated to the contour length, only in gel state, here 1. – C conjugated to tilt, related to Deserno-Terzi [3] or Brown [4]coupling, here 0.05 – C conjugated to the separation between the lipid end-tails, op-tional , improves the chain structure appearance, here 0.0528nce the model is set, one seeks for a minimum of the total energy. Min-imization is done by gradient descent. Analytical gradients H ,x i , H ,u i , H ,Z i can be obtained and evaluated directly. A repeated iteration ( x i ; u i ; Z i ) → ( x i − δt · H ,x i ; u i − δt · H ,u i ; Z i − δt · H ,Z i ) leads to a minimum of energy H after a number of ca × iterations, δt ∼ − , a process similar to therelaxation stage in Molecular Dynamics.The gradient with respect to the chain scaling parameter L , provides aIrving-Kirkwood-Buff kind of virial expression, that matches the externalapplied force F e . Local mechanical equilibrium is achieved automatically solong as H ,x i and H ,u i vanish separately. Forces oriented parallel to the z direction are disregarded. For completeness, H ,L reads H ,L = N X i =1 (cid:8) C L (cid:2) L ( x [ i +1] (cid:5) x i ) + ( Z [ i +1] − Z i ) − A ) (cid:3) ( x [ i +1] (cid:5) x i ) + 2 C Z i (cid:2) LZ i ( x [ i +1] (cid:5) x [ i − + u [ i +1] (cid:5) u [ i − ) − V (cid:3) × ( x [ i +1] (cid:5) x [ i − + u [ i +1] (cid:5) u [ i − )+4 C LS i (cid:2) L ( x i (cid:5) u i ) − L c (cid:3) ( x i (cid:5) u i ) + 2 C L ( x i (cid:5) u i ) (cid:9) (E.9)Let us summarize the main results of the 1d model. We investigated pairs( N g , N f ) of N g consecutive gel lipids and N f fluid lipids (Figure E.21). Thestructure of the (12 ,
0) all-gel minimum was found to be an homogeneoustilted phase . The area per lipid was 0 . < A = 0 .
8, see Figure E.20. Thestructure of the (0 ,
12) all-fluid minimum was a thin non tilted phase . Thearea per lipid was 0 . ’ A . The energy was on purpose higher than the all-gel phase. Figure E.20 shows that a minimum of energy is obtained for (10 , , L dβ state. We note however that parameters have not yet been really optimized29 igure E.20: Top, energy of a 12 beads chain as a function of the number of contiguousmelted lipids. Bottom, projected (circles) and contour (square) length as a function of thenumber of melted lipids with respect to experimental data (area per lipid in fluid and gel phase, tilt,membrane thickness in fluid and gel phase, enthalpy of melting per lipid. Asit is, this current model does not predict any lower critical lateral size for theapparition of the instability. It seems that the conformation (4 , , ,
2) haslower energy than (10 ,
2) which suggests that the ripple period could be veryshort. A gel-fluid neighbour mismatch penalty might help with enlarging thelateral corrugation size. 30 igure E.21: From left to right and top to bottom, minimal energy configurations forstates ( N g , N f ) = (12,0),(10,2),(8,4),(6,6),(4,8), (2,10),(0,12). eferences [1] S. Doniach, Thermodynamic fluctuations in phospholipid bilayers, TheJournal of Chemical Physics 68 (1978) 4912–4916. URL: http://link.aip.org/link/?JCP/68/4912/1 .[2] T. Heimburg, Thermal Biophysics of Membranes, Wiley-VCH, 2007.[3] M. M. Terzi, M. Deserno, Novel tilt-curvature coupling in lipid mem-branes, The Journal of Chemical Physics 147 (2017) 084702. URL: https://doi.org/10.1063/1.4990404 .[4] M. C. Watson, E. S. Penev, P. M. Welch, F. L. H. Brown, Thermal fluctu-ations in shape, thickness, and molecular orientation in lipid bilayers, TheJournal of Chemical Physics 135 (2011) 244701. doi:10.1063/1.3660673