A machine learning assessment of the two states model for lipid bilayer phase transitions

TL;DR

This study applies machine learning algorithms to classify lipid bilayer phases, providing evidence supporting the two-states model of membrane fluidity through analysis of molecular dynamics data.

Contribution

The paper introduces a machine learning method to classify lipid phases without assuming the order parameter, supporting the two-states model of membrane phase transitions.

Findings

Machine learning accurately classifies lipid configurations as fluid or gel.

Results support the two-states model of lipid bilayer phase transitions.

Method does not rely on prior assumptions about the order parameter.

Abstract

We have adapted a set of classification algorithms, also known as Machine Learning, to the identification of fluid and gel domains close to the main transition of dipalmitoyl-phosphatidylcholine (DPPC) bilayers. Using atomistic molecular dynamics conformations in the low and high temperature phases as learning sets, the algorithm was trained to categorize individual lipid configurations as fluid or gel, in relation with the usual two-states phenomenological description of the lipid melting transition. We demonstrate that our machine can learn and sort lipids according to their most likely state without prior assumption regarding the nature of the order parameter of the transition. Results from our machine learning approach provides strong support in favor of a two-states model approach of membrane fluidity.

Figures

Equations

$$S_{\mathrm{mol}}=\frac{1}{2(N_C-2)}\sum _{j=1}^2\sum _{k=2}^{N_C-1}\frac{1}{2}\left(3\cos ({\theta }_{(j,k)}{)}^2-1\right)$$

$$\ln \left(\frac{P(\mathrm{fluid}|\{D_j\})}{P(\mathrm{gel}|\{D_j\})}\right)$$

$$P(s|\{D_j\})=\frac{P(\{D_j\}|s)P(s)}{P(\{D_j\})},$$

$$P(\{D_j\}|s)\prod _{j=1}^{N_{\mathcal{D}}}\mathrm{d}{D}_j=\prod _{j=1}^{N_{\mathcal{D}}}\frac{\mathrm{d}{D}_j}{\sqrt{2\pi }{\sigma }_{s,j}}\exp \left(-\frac{(D_j-D_{s,j}{)}^2}{2{\sigma }_{s,j}^2}\right)$$

$$\sum _{j=1}^{N_{\mathcal{D}}}\frac{(D_i-D_{\mathrm{fluid},i}{)}^2}{2{\sigma }_{\mathrm{fluid},i}^2}-\sum _{j=1}^{N_{\mathcal{D}}}\frac{(D_i-D_{\mathrm{gel},i})}{2{\sigma }_{\mathrm{gel},i}^2}+\mathrm{Const}$$

$$\mathrm{sign}\left(b+\sum _{j\in \mathcal{J}}{\alpha }_j{y}_j\mathcal{K}({\mathbf{x}}_j,\mathbf{x})\right).$$

$$\mathcal{K}(\mathbf{x},{\mathbf{x}}^{\prime })=\exp \left(-\gamma ||\mathbf{x}-{\mathbf{x}}^{\prime }|{|}^2\right),$$

Full text

Phospholipid molecules play a major structural role in biological membranes Mouritsen_MatterOfFat ; Cevc_Marsh_PhospholipidBilayers where a deep understanding of the physical properties can only be acquired from a detailed knowledge of the lipid assemblies. Thanks to their amphiphilic nature and geometrical characteristics, most phospholipid molecules spontaneously self-assemble in water as bilayers. Supported or free-standing lipid bilayers and vesicles can easily be made, controlled and studied, and have now become standard tools in membrane biophysics studies, referred as model lipid bilayer systems Dimova_Marques_TheGiantVesicleBook . Early studies on pure phospholipid bilayers indicated that lipids were subject to thermodynamic transitions 1976_Mabrey_Sturtevant ; Cevc_Marsh_PhospholipidBilayers ; Heimburg_BiophysicsMembrane ; Marsh_HandbookLipidBilayers2 , with in particular a sharp transition associated to a significant change in enthalpy called main, or melting transition. This transition separates a low temperature well-packed assembly from a high temperature expanded, disordered lipid tail organization. This transition is considered to be of first order, although phase coexistence is difficult to establish experimentally 2012_Armstrong_Rheinstadter . The low temperature phase is commonly referred as the gel state, while the high temperature state is the fluid state. Lipid mixtures also display melting transitions spreading along a finite temperature range, usually accompanied by gel-fluid domain coexistence. It is usually assumed that most biological lipid membranes are found in a fluid state, and many scenarii aiming at explaining the lateral lipid and protein segregation observed in biological membranes involve ordering of the lipid tails.

The consensual description of the single component lipid melting transition assumes that dominant molecular conformations evolve from all-trans extended, well oriented hydrocarbon chain conformations in the low temperature phase, to disordered chains melted by rotation isomerism, as proposed in the earliest theoretical proposals 1973_Nagle ; 1974_Marsh ; 1974_Marcelja_2 ; 1979_Pink_Chapman . Lipids in the fluid phase have more configuration entropy and more enthalpy than those in the gel phase, due to lower density and cohesive energy, and higher chain torsion energy. Balance between entropy and enthalpy holds at the melting temperature $T_{m}$. Despite some asymmetry between the two phases, a large number of experimental facts related to melting transition of pure and mixed lipid compositions have been successfully interpreted by means of a phenomenological two-states model, originally proposed by Doniach 1978_Doniach ; 1983_Sugar_Monticelli ; 1996_Heimburg_Biltonen ; 1996_Jerala_Biltonen ; 1999_Sugar_Biltonen ; 2001_Ivanova_Heimburg ; Heimburg_BiophysicsMembrane ; 2011_Wolff_Thalmann ; 2018_Morandi_Marques . This model can be expressed as an Ising model, each lipid taking binary discrete values (say $s=0$ for gel and $s=1$ for fluid) with neighboring lipids being coupled Heimburg_BiophysicsMembrane . In the framework of the two-states model, the Ising variables stand for a coarse-grained description of the lipid tail conformations, assuming that lipids can be classified into two classes, according to their molecular conformations. Within this description, an effective temperature dependent ”magnetic field” $h(T)$ biases the odds in favor of one or the other state, while cooperativity results from nearest neighbor state coupling.

Usually, the determination of the lipid bilayer state relies on a structural scalar order parameter, such as the membrane thickness, given for instance by the head to tail lipid extension or the tail molecular order parameter. We address in this work the validity of the two-states description for the transition of pure DPPC (1,2-dipalmitoyl-sn-glycero-3-phosphatidylcholine) bilayers, using atomistic molecular dynamics (MD) simulations and supervised Machine Learning (ML) classification algorithms. Assessing the two states model requires ones to analyze single lipids and sort them into their respective states. Our ML classification works without reference to an existing or newly defined scalar order parameter. It only relies on a procedure for lifting the orientation degeneracy of each lipid configuration. In addition, the ML approach may serve as testing the relevance of a given scalar order parameter a posteriori.

Machine Learning has already been applied successfully to a number of situations in statistical thermodynamics and phase transitions. For instance, Cubuk et al. used support vector machines to localize plastic flow regions in amorphous structures 2015_Cubuk_Liu , Carrasquilla and Melko revealed the strong aptitude of neural networks models to recognize various spin ordering regimes in condensed matter systems 2017_Carrasquilla_Melko , Le and Tran succeeded in predicting the polymorphism of complex lipid mixtures given a set of structural, chemical and composition parameters, by means of an artificial neural network approach 2019_Le_Tran . We show in the present work that Machine Learning classification tools are also extremely powerful for investigating thermal lipid phase transitions.

DPPC molecules are among the best known phospholipids Marsh_HandbookLipidBilayers2 . They display experimentally a melting transition at 314 K, which is well reproduced by the CHARMM-36/TIP3P force-field for atomistic simulations of lipid bilayers in aqueous solutions 2010_Klauda_Pastor ; best_mackerell_2012 . In our simulations, in agreement with experiments, the bilayer was found in a $L_{\alpha}$ fluid disordered state above the main transition, and in a $P_{\beta^{\prime}}$ rippled phase (Figure 1), i.e a spatially modulated gel phase with peristaltic variations in the bilayer thickness and normal direction 2005_deVries_Marrink ; 2018_Khakbaz_Klauda at lower temperatures. Unlike Khabkaz and Klauda, our simulations did not show evidence of a $L_{\beta^{\prime}}$ tilted gel phase at temperatures below the pretransition, which is expected to take place at $T=307$ K on experimental grounds. This is likely to be due to a difference in size in our system, comprising nearly 3 times as many lipids as the one used in ref 2018_Khakbaz_Klauda . The presence of a ripple phase, whether stable or metastable, renders approaches based on scalar order parameters impractical, while our ML approach was found to work efficiently also in this situation. In what follows, we still name our low temperature state a gel phase, in spite of its ripple character.

The principle of the analysis is as follows. A 212 lipid molecules system was thermalized at low (288 K) and high temperature (358 K) and pressurized with a semi-isotropic barostat. At such temperatures, we assume that lipid conformations are predominantly gel and fluid respectively. Our training set was therefore composed of an equal number of conformations coming from the 288 K MD trajectory (all labelled as gel) and the 358 K trajectory (labelled as fluid). Details regarding the numerical simulation procedure and machine learning implementation details are given as Supplemental Material (SM).

Raw molecular conformations (spatial coordinates) were processed to remove all translation and rotation degeneracy, and the dimension of the initial conformation space was slightly reduced, resulting respectively in a 100-dimensional reduced coordinate space $\mathcal{X}$, and a 61-dimensional mutual distances space $\mathcal{D}$. The coordinate space $\mathcal{X}$ and the distance space $\mathcal{D}$ provide a very detailed description of individual lipid conformations, and constitute the starting point of the ML approach. They are referred as ”feature spaces” in the context of automated classification (see SM). The processed conformations were then fed to the ML algorithms. Three different algorithms were selected for the purpose of classifying the molecular conformations: Naive Bayes (NB), K-Nearest Neighbors (KNN) and Support Vector Machines (SVM). They were all used as implemented in the Python/Scikit-learn package pedregosa_duchesnay_2011 . The different ML algorithms were tested and found to perform moderately well when used separately, some methods performing better for gel lipids, other methods for fluid lipids. Combining the approaches together, we managed to get a success rate of 88% upon validation, i.e. using 80% of the training set configurations for learning, our best ML algorithm was able to assign 86% of the low temperature configurations to a gel state, and 91% of the high temperature configurations to a fluid state. After training, the ML model was used to analyse simulations at arbitrary temperatures.

One important prediction of the two states model is the presence of minor components in the majority phase, under the form of ”thermal excitations”. A finite fractions of molecules tend to adopt a conformation different from their immediate environment, in spite of the presence of a ”local field” biasing the statistics in favor of the dominant state. In that respect, the training set cannot be considered as containing only pure gel and fluid conformations. The presence of minor components is inherent to the presence of thermal excitations, and it does not seem possible to curate the training set without introducing further unwanted biases into the analysis. Our results show, however, that the training algorithm is not sensitive to the presence of a small fraction of non representative lipid conformations. In other words, the learning procedure was found to be robust so long as the training sets temperatures were chosen far apart from the melting transition.

Let us first analyze how ML predictions differ from those based on standard scalar structural order parameters. Figure (2) shows the distribution of the molecular segmental order parameter $S_{\mathrm{mol}}$, below and above the transition. This scalar observable was defined, for a given lipid and an instantaneous configuration (coming from a MD trajectory frame), by averaging over all the CC bonds in the two aliphatic tails a nematic parameter $(3\cos(\theta)^{2}-1)/2$, $\theta$ being the bond orientation with respect to the bilayer normal direction $z$. Both histograms overlap significantly. Using for instance a threshold value $S_{\mathrm{th}}=0.55$, where the estimated probability density functions (pdf) coincide, it was found that altogether 17% of lipids (1 in 6) were set to be assigned to the opposite state, either gel at 358 K, or fluid at 288 K. The gel $S_{\mathrm{mol}}$ distribution appears to be strongly skewed, as 28% (1 in 3) of the lipid molecules end up in the fluid state at 288 K. The relatively large fraction of lipids with $S_{\mathrm{mol}}\leq 0.55$ found at 288 K is in part due to the ripple structure of the bilayer. This clearly show that using $S_{\mathrm{mol}}$ for categorizing lipid conformations gives, at best, unreliable results. A similar analysis was conducted for two other scalar determinants: the lipid head-to-tail extension $L$, and the area per lipid $A_{l}$ computed from a Voronoi tessellation of the 2d $xy$ projection of lipid center of mass positions (assigning to each lipid a polygonal cell of given area). Both were found to perform worse than the segment orientation $S_{\mathrm{mol}}$ approach (see SM).

To provide further evidence in favor of our approach, we compared the ML predictions to the membrane structure as a function of temperature. Figure (3) shows on the same graph the average area per lipid $A_{l}$, the average volume per lipid $V_{l}$ and the ML prediction for the fraction of lipid assigned to the fluid state. The area per lipid is defined as the projected membrane area divided by the number of lipids per leaflet, neglecting out-of-plane membrane undulations. The volume per lipid results from a Voronoi analysis of the lipid center of masses. The structural values evolve monotonically and reversibly with temperature, from 288 to 358 K. The three curves superimposes very well, which shows that structural data supports the finding of our classification tool. Two snapshots of the membrane upper leaflet are also provided, at low and intermediate temperature. The low temperature membrane shows a small number of isolated small clusters of fluid lipids. A one to one fluid/gel ratio is observed at 318 K (51 $\pm$ 3% of lipids in the fluid phase), a temperature close to the experimental melting temperature (314 K). A snapshot of the upper leaflet at 318 K shows two large distinct domains, separated by a smooth boundary. The fluid state ratio evolves smoothly and reversibly from an estimated value of 14 $\pm$ 2% at 288 K to 95 $\pm$ 1% at 358 K. The very strong correlation between fluid ratio $x_{f}$ and area per lipid $A_{l}$ can be naturally interpreted in the framework of the two-states model by assigning constant values $A_{f}$ and $A_{g}$ to the fluid and gel states, resulting in $A_{l}=x_{f}A_{f}+(1-x_{f})A_{g}$. A similar conclusion could be reached when considering the volume per lipid $V_{l}$.

The ML predictions agree well with the global structural properties of the membrane and we now consider the local correlation properties. The two-states model combines the internal spontaneous dynamics of each lipid with local interactions promoting cooperativity between neighbors. The interactions, known as $J$ coupling in the Ising model context, are responsible for the sharp structural and thermodynamic changes with temperature, the emergence of local correlations, clustering and domain formation, such as seen in Fig 3. Such couplings create a local field whose effect is to bias the lipid state in favor of the dominant local phase, according to the majority rule. Internal reversible gel $\leftrightarrow$ fluid state transitions occur spontaneously, according to a non conserved parameter, or Glauber dynamics 1963_Glauber . In order to get insight into the local correlations and flip rates, we performed a systematic statistical neighbor analysis, and estimated the conditional conversion states.

Voronoi tessellations provide an operational method for deciding unambiguously which lipid pairs are nearest neighbors. Following a protocol described in the SM, we consider two consecutive MD frames and divide the lipid population into four subsets: (a) lipids categorized as fluid, which stay in the fluid state, (b) lipids categorized as fluid, switching to the gel state, (c) lipids in the gel state, remaining in the gel state and (d) lipids in the gel state switching to a fluid state. At $T=318~{}\mathrm{K}$ the fluid and gel phases compete evenly and the fluid and gel populations are roughly equal. We count for each molecule, the number $n_{g}$ of gel neighbors and $n_{l}$ of fluid neighbors. We find in these conditions that the most typical environment of a lipid in fluid state (a and b) is $n_{g}=2\pm 1$ and $n_{f}=7\pm 1$, while for a lipid in gel state (c and d) $n_{g}=7\pm 1$ and $n_{f}=1\pm 1$. However, the internal lipid state fluctuate spontaneously, and clearly the dynamics is strongly influenced by the local environment, as demonstrated in Figure (4). Indeed, for gel to fluid (case d) and fluid to gel (case b) transitions, we notice that all lipids that are just about to switch are more likely to have an equal number of gel and fluid neighbors (typically $n_{g}=3\pm 1$ and $n_{f}=5\pm 1$). We conclude that lipids subject to internal state transitions are mostly located at the border between domains. The results of our neighbor analysis clearly support the idea that internal state dynamics is under the control of some local field.

As a summary, we trained a Machine Learning algorithm to classify phospholipid molecular conformations obtained by atomistic molecular dynamics simulations. Lipids were sorted into two classes, gel and fluid, according to their similarity with a reference (training) set of conformations originating from low and high temperature trajectories respectively. The efficiency of the ML approach is superior to simple schemes based on scalar order parameters, and deals successfully with the ripple structure at low temperature. The measured fraction of fluid/gel conformations correlates very well with the observed structural changes as temperature evolves. A finite fraction of the minor phase is always present which can be associated to thermal excitations in the framework of a two-states Ising model interpretation. Lipids with similar state tend to cluster into large domains, while spontaneous internal state conversions are more likely to occur at the boundary between domains. The local distribution of neighbors supports the concept of local field and nearest neighbor coupling. We overall conclude that our ML approach provides convincing evidence in favor of the two-states phenomenological model.

We foresee numerous applications of the present approach. A first straightforward extension concerns lipid binary mixtures 1988_Ipsen_Mouritsen and liquid ordered phases caused by the presence of cholesterol 1987_Ipsen_Zuckermann , as well as membranes of more complex composition 2014_Ingolfsson_Marrink . We also anticipate that our ML approach will be useful to study the influence of membrane solutes that are known to influence the thermodynamics of melting in model membranes. This includes hydrophobic pollutants, e.g. pyrene 2014_Franova_Ollila , carbohydrates 1991_Crowe_Crowe ; 2018_Morandi_Marques and synthetic oligomers 2014_Rossi_Monticelli ; 2014_Rossi_Monticelli_2 . A ML approach could then quantify the lipid state alteration induced by these compounds. Importantly, this analysis is also well-suited to study the local lipid environment of membrane proteins, for which the existence of lipid mediated interactions and minor phase nucleation is speculated 1997_Gil_Mouritsen ; 1998_Gil_Zuckermann . For all these cases, significant improvements over approaches relying on scalar order parameters can be anticipated.

V.W. warmly thanks Tiago Espinosa de Oliveira for helping with simulations set-up. The authors gratefully acknowledge support from the high performance cluster (HPC) Equip@Meso from the University of Strasbourg, through grant n*∘* G2018A53.

I Simulations

All simulations were performed using GROMACS 2016.4 berendsen_vandrunen_1995 ; abraham_lindahl_2015 along with the CHARMM-36 all-atom force-field best_mackerell_2012 (June 2015 version). A lipid bilayer made of 106 lipid molecules per leaflet, each containing 130 explicit atoms, was created using CHARMM-GUI jo_im_2008 ; wu_im_2014 ; jo_im_2009 ; jo_im_2007 . It was hydrated with two 8 nm thick water layers on each side (connected through periodic boundary conditions), using the TIP3P water model, for a total of 29826 solvent molecules. The force field parameters for DPPC molecules were provided directly by CHARMM-GUI brooks_karplus_2009 ; lee_im_2016 .

The above system was first subject to energy relaxation using steepest descent energy minimization, followed by a 10 ps NVT thermalization stage at 288 K. Then, the bilayer was subject to a 1 ns NPT equilibration run coupled to a semi-isotropic barostat (1 bar in all directions). The system was then further equilibrated at the desired temperature with the same semi-isotropic barostat during a second NPT equilibration step of 10 ns. Molecular dynamics production runs of 50 ns were finally generated at the same temperature and with the same semi-isotropic barostat. The analysis were performed on the last 25 ns of simulations. All time steps were set to 2 fs.

All the molecular dynamics simulations used the leap-frog integration algorithm hockney_eastwood_1974 . Temperature and pressure were kept constant using respectively a Nosé-Hoover thermostat nose_1984 ; hoover_1985 (correlation time $\tau_{T}=0.4$ ps) and a Parrinello-Rahman semi-isotropic barostat nose_klein_1983 ; parrinello_rahman_1998 (correlation time $\tau_{P}=2.0$ ps, compressibility $4.5\times 10^{-5}$ bar*-1*).

Lipid and water molecules were separately coupled to the thermostat. Following GROMACS recommendations for the CHARMM-36 all-atom force field, 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) were also set to 1.2 nm. Fast smooth Particle-Mesh Ewald electrostatics was selected 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 was applied to the Van der Waals interactions. We did not account for long range energy and pressure corrections, and constrained all the hydrogen bonds of the system using the LINCS algorithm.

II System Analysis

II.1 Determination of structural parameters

Values of the average area per lipid $A_{l}$ and order parameter $S_{\mathrm{mol}}$ of the bilayer were respectively obtained using the GROMACS built-in commands gmx energy and gmx order.

For measurements of individual lipid properties (area, order parameter, elongation), atom positions were collected from trajectories using the Python 3 MDAnalysis module gowers_beckstein_2016 ; michaud_beckstein_2011 . The individual areas per lipid were obtained from Voronoi tessellations of the two-dimensional projections of the lipid center of masses, computed using the Voro++ library rycroft_2009 . The individual volumes per lipid were derived from three-dimensional tessellations of the lipid centers of mass, again using the Voro++ library. Note that the bilayer geometry requires a specific tessellation procedure: this was done by introducing ghosts lipids in the water regions. Without these ghost lipids, the tessellation cells cannot be correctly defined and are unbounded across the membrane-water interface, thus overestimating significantly the individual volume per lipid. Ghosts lipids are mirror images of bilayer lipids across the local lipid-water interface (cf. Fig. S1). After the tessellation was made, ghost lipids described in the previous section and their corresponding cells were removed the lists, and only the volumes of physical lipids were collected and analyzed.

The molecular order parameter $S_{\mathrm{mol}}$ of individual lipids was calculated by measuring, for every $N_{C}-2=14$ non-terminal carbon atoms $k=2\ldots 15$ of the 2 tails of the lipids, the angle formed between the $z$-axis of the system directed along $\vec{u}_{z}$ and the vector $\overrightarrow{C_{(j,k-1)}C_{(j,k+1)}}$ defined by the carbon atoms surrounding atom $k$ within the same tail $j=1,2$. The order parameter $S_{\mathrm{mol},(j,k)}$ of the atom $k$ is obtained from the 2nd Legendre polynomial $P_{2}$ using $\cos(\theta_{(j,k)})=\vec{u\vphantom{c}}_{z}\cdot\overrightarrow{C_{(j,k-1)}C_{(j,k+1)}}$, and averaging over $j$ and $k$:

[TABLE]

II.2 Next-nearest neighbors statistics

After completion of the 3d Voronoi tessellation using Voro++, a list of next-nearest neighbors was established for each lipid center of mass. We also collected the areas of the polygonal surfaces separating each pair of neighboring Voronoi cells. The ghost lipids and their corresponding faces were removed from the lists. The neighbor lists were further curated by removing all the neighbor pairs for which the corresponding face area accounted for less than 1% of the total surface area of each Voronoi cell in contact. The number of next-nearest neighbors were finally counted to build the coordination statistics $(n_{g},n_{f})$, where each lipid molecule has $n_{g}$ gel and $n_{f}$ fluid neighbors.

III Machine Learning discrimination of lipid states

III.1 Lipid classification procedure

The lipid classification process involves three steps:

• •

The molecular conformation of a lipid is recorded as a list of atom positions. Lipids are then shifted and rotated in order to remove all rotation and translation degeneracy.

• •

The lipid molecular conformations are further simplified, reducing each single lipid conformation to a set of 100 $(r,z)$ coordinates (configuration space $\mathcal{X}$), or 61 mutual distances (configuration space $\mathcal{D}$, see below for details).

• •

424 conformations were selected with the purpose of training the algorithms. A number of Machine Learning procedures were tested and compared. A combination of 4 algorithms was found to maximize the training success rate. We therefore combine the 4 algorithms in our final classification procedure (see below for details).

III.2 Definition of the configuration spaces

Machine Learning algorithms share the ability to discriminate well multidimensional data. Our purpose is to feed the algorithms with single lipid conformations and have them sorted into two classes. We observe that the raw chemical formula of DPPC (CAS 2644-64-6) is PNO8C40H80. Disregarding hydrogens, each lipid molecule comprises 50 ”heavy atoms”, and therefore the associated single lipid configurations belong to a 150 dimensions vector space. As some configurations can be mapped onto each other by means of a spatial displacement (translation and/or rotation), the set of configurations has only 144 independent degrees of freedom.

For simplicity, we decided to work with slightly smaller configuration spaces, by projecting further the 144 dimensional molecular conformations onto two configuration spaces, $\mathcal{X}$ and $\mathcal{D}$, that we now define. Within the first approach, one determines the proper inertial frame of each lipid configuration, locating the center of mass, fitting the position of all atoms with a 3D line to find the longest axis. Lipid configurations can be recast into the inertial frame, using the center of mass as origin, and the longest axis (smaller moment of inertia) as vertical axis. This almost always results in having the lipid directed along the bilayer normal, with two possible orientations. When necessary, the lipid is flipped in order to place by convention the phosphocholine group into the positive $z$ upper half space. This combination of reorientations lifts entirely the translation and orientation degeneracy of the original coordinate space. New cylindrical coordinates $\{r_{i},\theta_{i},z_{i}\}$, $i=1\ldots 50$ can be associated to the resulting lipid conformation. We decided then to disregard entirely the angles $\theta_{i}$, keeping only the coordinate subset $\{X_{i}=r_{i},z_{i}\}$, $i=1\ldots 50$. This defines a $N_{\mathcal{X}}=100$ dimensional space $\mathcal{X}$ that will subsequently be referred as ”Coordinate space”.

For the second approach, starting from the original 150 lipid coordinates, we calculate a set of mutual euclidean distances between pairs of atoms. Two atoms participate in a pair when they are separated by 6 positions along the chemical graph defining the molecule (see Fig S2). Enumerating the possibilities of the tree-like graph, one finds 61 non equivalent pairs of atoms, associated to 61 distances $\{d_{j}\},j=1\ldots 61$. This defines a $N_{\mathcal{D}}=61$ dimensional ”Distance space” $\mathcal{D}$, without translation and rotation degeneracy.

In both cases, the original dimension of the problem is reduced (from 144 to 100 and 61 respectively) but yet the configuration spaces $\mathcal{X}$ and $\mathcal{D}$ are large and preserve to a large extent the complexity of the original conformations. As our comparisons show, the classification of lipid states is efficient, whether one starts from $\mathcal{X}$ or $\mathcal{D}$.

Three different algorithms commonly used for Machine Learning classifications were used in this Letter. A short description of these algorithms is given below.

III.3 Naive Bayes

Given a configuration space $\mathcal{D}=\{D_{j}\},j=1\ldots N_{\mathcal{D}}$, and two classes $s=\{\mathrm{gel},\mathrm{fluid}\}$, the Bayesian approach assumes that it exists a joint probability distribution $P(s,\{D_{j}\})$, that weighs the respective likelihood of each state $s$ once a configuration $\{D_{j}\}$ is provided. The Bayesian model makes a decision regarding the state $s$ by comparing the conditional probability densities $P(\mathrm{gel}|\{D_{j}\})$ and $P(\mathrm{fluid}|\{D_{j}\})$. To be precise, the Bayesian approach attributes a fluid label to a configuration if the sign of

[TABLE]

is positive, and a gel label otherwise.

Each Bayesian approach provides a mathematical model $P(\{D_{j}\}|s)$ describing the expected configuration distribution for a given class $s$: fluid and gel. There is in principle entire freedom in choosing the model, but the efficiency and the optimization requirements limit such choices in practice. Training a Bayesian model therefore amounts to finding the most realistic function $P(\{D_{j}\}|s)$ as far as classifying a given training set of data is concerned.

The Bayes theorem provides the connection between the conditional probabilities entering in the choice function (2) and the model:

[TABLE]

where $P(s)$ represents the prior distribution of $s$, i.e. the statistical distribution of $s$ in the absence of any configuration related information, and $P(\{D_{j}\})=P(\mathrm{gel},\{D_{j}\})+P(\mathrm{fluid},\{D_{j}\})$ a normalization factor which cancels out in eq (2). The Naive Bayes (NB) gaussian model assumes that $P(\{D_{j}\}|s)$ factorizes as a product of independent gaussian distributions of $D_{j}$.

[TABLE]

Training the algorithm means finding the best mean value $D_{s,j}$ and standard deviation $\sigma_{s,j}$ for every parameter $D_{j}$ in the distance space $\mathcal{D}$ and each class $s$. The number of parameters to compute turns out to be equal to twice the dimension of $\mathcal{D}$. In our case, the training set contains an equal number of gel and fluid conformations, and there is no a priori bias between classes, meaning that $P(\mathrm{gel})=P(\mathrm{fluid})=1/2$. Therefore, the number of parameters to determine during training is $2\times 61=122$.

To sum up, the Naive Bayes approach classifies a lipid conformation by computing a quadratic score function in the conformation space,

[TABLE]

and deciding whether a $\{D_{i}\}$ lies closer to $\{D_{\mathrm{fluid},i}\}$ or to $\{D_{\mathrm{gel},i}\}$ according to this generalized distance.

III.4 K-Nearest Neighbors

The K-Nearest Neighbors (KNN) classification method is based on defining a distance between any arbitrary pairs of objects to discriminate. A natural choice is the Euclidean norm of the feature space, here the $N_{\mathcal{X}}$ dimensional coordinate space $\mathcal{X}$.

The KNN algorithm finds the $K$ nearest neighbors within the training set, of each new configuration to classify. Decision is taken based on majority rule, i.e. the most abundant class found among the $K$ closest neighbors. The optimal $K$ is determined during the training and validation process, and was set equal to 5, a typical value for this method.

III.5 Support Vector Machines

Support Vector Machines (SVM) classify data by means of linear separation in high dimensional representation spaces. Denoting $\bm{\phi}$ an arbitrary element in a given representation space $\mathcal{R}$, the binary decision is given by the sign of the affine expression $\bm{w}\cdot\bm{\phi}+b$, with $b$ a numerical constant and $\bm{w}$ the hyperplane of separation normal vector. In a few favorable cases, it is possible to use directly the data definition space as representation space. However, in many practical situations, efficient classification can only be achieved by mapping the data (e.g. $\mathcal{X}$ or $\mathcal{D}$) onto a larger representation space, namely $\bm{x}\mapsto\bm{\phi}(\bm{x})$. Training a SVM corresponds to choosing a suitable representation space $\mathcal{R}$, and finding the optimal $b$ and $\bm{w}$. As shown in boser_vapnik_1992 ; cortes_vapnik_1995 , given a training set $\{\bm{x}_{i}\},i=1\ldots n$, the optimal $\bm{w}$ can always be expressed as a linear combination $\bm{w}=\sum_{i}\alpha_{i}y_{i}\bm{\phi}(\bm{x}_{i})$ with either positive or vanishing $\alpha_{i}$ coefficients, and $y_{i}=\pm 1$, depending on the class (fluid 1, gel -1) of the corresponding vector data $\bm{x}_{i}$.

The subset of vectors $\{\bm{x}_{i}\}$ participating in the definition of $\bm{w}$ with non vanishing coefficients $\alpha_{i}>0$ forms the so-called support vectors. Denoting $\mathcal{J}$ the sequence of indices of support vectors, and $\mathcal{K}(\bm{x},\bm{x}^{\prime})$ the product $\bm{\phi}(\bm{x})\cdot\bm{\phi}(\bm{x}^{\prime})$ in $\mathcal{R}$, called kernel function, the SVM decision function for any vector data $\bm{x}$ reads:

[TABLE]

The SVM training optimization problem can therefore be formulated without any explicit reference to the representation space $\mathcal{R}$, nor the mapping $\bm{\phi}(x)$. It only requires an explicit positive kernel function $\mathcal{K}(\bm{x},\bm{x}^{\prime})$. As explained in boser_vapnik_1992 ; cortes_vapnik_1995 , there are efficient quadratic optimization algorithms for determining the support vectors $\bm{x}_{j}$, the non-vanishing coefficients $\alpha_{j}$ and the shift constant $b$.

In this study, we used the standard radial basis kernel function

[TABLE]

with a default value $\gamma$ equal to the inverse of the dimension of the configuration space. This choice assumes that each component of $\bm{x}$ is of order 1. When using the coordinate space $\mathcal{X}$, the trained SVM ends up using 297 non vanishing support vectors and coefficients $\alpha_{j}$, $\gamma=1/100$, $b=-0.416$ (167 fluid, 130 gel support vectors). When considering the distance space $\mathcal{D}$, the trained SVM used 161 non vanishing coefficients $\alpha_{j}$, $\gamma=1/61$ and $b=0.134$ (100 fluid, 61 gel support vectors).

IV Training

The Machine Learning analysis performed in this Letter were conducted using the Scikit-Learn module for Python 3 pedregosa_duchesnay_2011 . An unbiased selection of lipid conformations extracted from a trajectory at low temperature (288 K) was part of the training set, with a label gel. Similarly, an unbiased selection of conformations from a trajectory at high temperature (358 K) was added to the training set with a label fluid. The number of gel and fluid conformations were in equal number in the training set. As customary, 20 % of conformations were removed from the training sets, and used for verification and scoring purposes. The training set therefore consists of a sequence of 370 vectors (elements of the configuration spaces $\mathcal{X}$ or $\mathcal{D}$), used for building the prediction model.

V Asserting the predictive capacity of each model

Once defined the training set, with properly labelled gel and fluid states, 20% of the lipids in each phase were taken apart for forming a validation set. After training the models (i.e. optimizing the parameters with respect to the 80% remaining conformations) a prediction score was separately calculated for the gel and fluid conformations in the validation set. The overall procedure was repeated 10 times, each time with the same training set, but independently drawn validation subsets.

The resulting scores shown in Fig. S3 shows that the NB predictive capacity is independent from the phase of the lipid considered. KNN performs better for fluids than for gels. At the opposite, ”coordinates” $\mathcal{X}$-SVM and ”distances” $\mathcal{D}$-SVM seems to perform better in the gel phase.

We therefore decided to combine the predictions of the above models, retaining those who perform best in each phase. The following decision chain was implemented:

If the 4 models agree on the same prediction for the phase, this prediction is retained; 2. 2.

if the $\mathcal{D}$-SVM algorithm predicts a gel phase, the lipid configuration is assumed to be gel; 3. 3.

if the $\mathcal{X}$-SVM algorithm predicts a fluid phase, the lipid configuration is assumed to be fluid; 4. 4.

if the NB algorithm predicts a gel phase, the lipid configuration is assumed to be gel; 5. 5.

if none of the above conditions have been met, the KNN algorithm makes the final decision on the configuration classification.

VI Comparison between Machine Learning decisions and structural characterizations of the lipid configurations

Machine Learning predictions were compared to two typical lipid structural properties: the carbon carbon (CC) order parameter $S_{\mathrm{mol}}$ along the chains and the area per lipid $A$ in the 2d Voronoi tessellation of the lipid projected centers of mass. The corresponding results are presented in Fig. S4(a) and (b).

The order parameter curves (Fig. S4(a)) clearly discriminate among the low temperature (288 K) and high temperature (358 K) fluid phases, in agreement with published results on these systems klauda_pastor_2010 .

The order parameter of the atoms in the lipid tails at low (288 K) and high (358 K) temperature is characteristic from membranes in the gel and fluid phases respectively Cevc_Marsh_PhospholipidBilayers . The phase transition can be clearly seen in Fig. S4(b) as a significant variation in the evolution of the area per lipid $A_{l}$ around 321 K.

Experimental structural values are available at 323 K nagle_suter_1996 ; nagle_tristram_2000 ; kucerka_nagle_2006 ; kucerka_katsaras_2008 . Nagle et al. obtained for DPPC an area per lipid equals to 64 $\pm$ 1 Å2 significantly close to the value circa 60 Å2 we obtained in our simulations. Using the average DPPC bilayer thickness reported by Nagle et al., we could estimate an experimental volume per lipid of 1220 $\pm$ 50 Å3, which agrees fairly with our Voronoi value of 1300 Å3

VII Naive classifications

The distributions of the areas per lipid and molecular elongations at low and high temperatures are shown in Fig. S5. Using a naive classification scheme based on a single threshold value for either of the two previous scalar parameters would at best result in a prediction accuracy of respectively 69% and 67%.

Fig. S6 compares the histogram of molecular order parameters $S_{\mathrm{mol}}$ as a function of the temperature of the lipid bilayer from which the configurations are extracted (288 K or 358 K), and as a function of the result of the Machine Learning classification procedure. The difference between the distribution at 288 K and the distribution in the a posteriori gel state ensemble indicates that a small fraction of lipids in the fluid state are already present at 288 K.