Probing Slipids Force Field for Phase Transitions in SOPC Lipid Bilayers with Various Cholesterol Concentrations

: We explore the phase behavior of lipid bilayers containing SOPC (1-stearoyl-2-oleoyl-sn-glycero-3-phosphocholine) with various molar concentrations (0 mol%, 10 mol% and 30 mol%) of cholesterol. To this end, we performed extensive atomistic molecular dynamics simulations in conjunction with the Slipids force field with optimized parameters for the headgroups of phospholipids. We computed thermodynamic and structural quantities describing the ordering of the tails, the mobility of the heads and the arrangement of the lipids in the bilayers. We analyzed the behavior of the named quantities over the temperature range between 271 K and 283 K, where the experimentally determined melting temperature, T m = 279 K, lies, as well as at 400 K, which is used as a reference temperature. The obtained results are compared to available experimental data along with the outcome from molecular dynamics simulations of similar phospholipids containing different amounts of cholesterol. In the temperature interval of interest, we found evidence of the occurrence of a thermal-driven phase transition (melting) in both the pure system and the one with the lower concentration of cholesterol, while in the remaining system, the higher amount of cholesterol in the bilayer smears out the transitional behavior. Thus, we demonstrate the ability of the Slipids force field to predict the phase behavior of bilayers of SOPC and SOPC mixed with cholesterol.


Introduction
The natural boundaries between biological cells are membranes whose main components are phospholipids [1].At different temperatures, pressures and admixtures' concentrations, various phases take place in the lipid bilayers.The occurrence of phase transitions [2,3] significantly affect basic membrane functions, such as cellular organization, nerve impulse propagation [4,5], membrane fusion [6] and others.On the other hand, manipulating the transition of lipid bilayers allows control over lipid fluidity to process the transport and release conditions of drug delivery.At low temperatures, in some cases, the sub-gel (crystalline) or solid-ordered (S O ) phase is realized [7,8].In this state, the lipids are practically immobile, showing a high degree of order and the characteristics of S O [9].By increasing the temperature, a gel phase in two distinct arrangements, depending on whether the carbon tails are tilted (L β ′ ) or upright (L β ), sets in [10].In the gel state, lipids have limited lateral mobility and are coordinatively ordered.A ripple phase (P β ′ ) may also occur in mixed lipid bilayers.In this case, different types of lipid molecules exhibit different phases at the same temperature [10].In the presence of cholesterol, a liquid-ordered phase (L O ) is realized [11].In this phase, there is a simultaneous increase in the tails arrangement (typical of the gel phase) and the heads mobility (typical of a liquid state) [11].Moreover, the lateral packing is looser than that typical to the gel phase, and the lipids' molecules are more tightly packed than the case of liquid phase.Cholesterol intercalation into lipid tails in S O or P β ′ disrupts the crystal lattice and reduces the ordering by favoring the trans-conformation of the lipid tails [12].With a subsequent increase in temperature to a certain temperature, say T m , the lipids melt, and a liquid phase (L α ) is reached.In this state, the arrangement in the tails and the heads vanishes, and the bilayer shows a behavior of a liquid medium [13].
At the main transition temperature T m there is a dramatic change in the tail arrangement, bilayer thickness and lipid area [14].The spontaneous change and the hysteresis indicates that this transition is first-order.While the phase L α is well studied in the gel state, the transitional behavior at T m is not fully understood.It has been suggested that the melting goes through an intermediate hexatic phase [15].Moreover, it is shown [16] that L β ′ is hexatic, and the melting is actually first-order from the hexatic to the liquid bilayer.In the presence of cholesterol above T m there is a disruption of the packing of the acyl tail [17], though their ordering increases, showing an L O phase [14].It is important to note that at high concentrations (about 60 mol%) of cholesterol, the phase transition is hindered [17].
A typical representative of monounsaturated lipids is SOPC (1-stearoyl-2-oleoyl-snglycero-3-phosphocholine).Its molecule, depicted in Figure 1, has a hydrophilic head and two hydrophobic tails -a saturated (sn-1) and an unsaturated (sn-2) chains.It was established that S O is realized at low temperatures, and by increasing the temperature, it turns into L α phase [18].The addition of a small amount of cholesterol into SOPC systems alters their phase behavior into L O .For more details, the reader may the consult the schematic phase diagram proposed in Ref. [19].Molecular dynamics is a powerful tool for the theoretical investigation of the phase behavior of complex lipid membranes at different temperatures.Atomistic simulations give useful insights on the mechanisms underlying their properties and functions.This may be achieved using a force field (FF) that accounts for the different interactions at the atomic level.The most known FFs reproduce satisfactorily well the experimental data obtained by probing membranes.Recent studies [20][21][22], comparing the accuracy of the different FFs, have demonstrated that Slipids FF [23] is the closest to experimentally available data.Recently [24][25][26], we have extended the use of Slipids FF at low temperatures and found that it reproduces the structural characteristics of the cholesterol-free system and its counterparts with concentrations ranging from 10 mol% to 50 mol% cholesterol fairly well.Due to the relatively short MD trajectories that did not capture the main features of the system at a potential phase transition, we were not able to give a conclusive answer as to what extent this FF describes the phase behavior of SOPC systems containing different amounts of cholesterol.
In the present article, lipid bilayers composed of SOPC containing several molar concentrations of cholesterol (0 mol%, 10 mol%, 30 mol%) are probed by extensive atomistic MD simulation over long trajectories with the aid of Slipids FF.Previously [26], we have established that 30 mol% cholesterol is critical for the influence of cholesterol in the SOPC systems, in the sense that beyond this concentration, the phase transition in the system is hampered by cholesterol.Here, the temperature interval of interest is 271 K -283 K, as it spans the temperature range around the experimentally determined T m = 279 K [27] for a pure SOPC bilayer.Moreover, we check whether Slipids FF is able to describe the typical behavior of lipid bilayers across phase transitions and whether it gives meaningful results in the neighborhood of the experimentally determined one.
The rest of the paper is organized as follows: Section 2 provides details on the lipid systems under considerations along with the MD simulation procedure.Section 3 presents the analysis of the simulation data and the main features of the computed thermodynamic and structural quantities and a discussion of the obtained results along with a comparison to the available literature experimental data and simulation outcomes.Finally, Section 4 summarizes our results on molecular modeling of the cholesterol-containing SOPC membranes.

Systems and Methods
We perform atomistic molecular dynamics simulations on three systems of lipid bilayers composed of SOPC with different cholesterol contents (0 mol%, 10 mol% and 30 mol%).The distinct membrane systems will be denoted by SOPC+ n% Chol, where n = 0, 10 and 30 corresponds to the molar concentration of cholesterol in the bilayer.The initial coordinates of the systems are built up with CHARMM-GUI [28].The bilayers contain 256 lipids symmetrically distributed over both layers, i.e., 128 molecules in each monolayer, and are subjected to cubic periodic boundary conditions.The number of lipids is twice that used to parameterize the force field [29].This suggests that the system is large enough to allow saturation of the basic parameters.The water model is TIP3P [30].A separate study showed that the optimal degree of hydration corresponds to 40 water molecules per lipid [31].Thus, the total number of water molecules is slightly different for all studied systems.A detailed analysis of the stability of the ensuing patches was performed in our previous paper [26].The constructed lipid bilayers were studied in the temperature range 271 K-283 K, where the experimentally determined melting temperature T m = 279 K [27,32,33] lies.In the three atomistic structures, a similar heating was carried out gradually by one kelvin from the preceding temperature.The systems at the lowest temperature (271 K) were relaxed.Then, a 1 kelvin heating and subsequent relaxation was undertaken; the procedure is then repeated until 283 K is reached.Moreover, additional simulations were carried out for the three systems at 400 K, that is much higher than T m .At this temperature, we expect all three bilayers to exhibit a fully disordered phase L α .Starting with one and the same equilibrated system to perform simulations to a target temperature has been used in several studies, see, e.g., [34,35], to probe the properties of different lipid systems including SOPC. Figure 2 shows the structures of the patches at the initial and final temperatures at the end of the simulation (3 µs).The final coordinates of all systems at 400 K are shown in Figure S1 in Supplementary Materials.
All MD simulations were performed using GROMACS 2022.4 [36,37] and Slipids 2020 FF [23].Although the field is reasonably correct, the present investigation will be conducted over a wide range of temperatures.The systems are energy-minimized and heated in the NVT ensemble up to the lowest temperature of interest (271 K) for 10 ns.Moreover, for the sake of comparison, all systems are heated up to 400 K, without intermediate steps, using the already described parameters of the simulations.We used the V-rescale thermostat [38] in conjunction with the Berendsen barostat [39].A leap-frog integrator [40] with a time step of 2 fs was used.PME [41] with a cut-off of 1.2 nm and a grid spacing of 0.16 for FFT was applied for the electrostatic potential.Lennard-Jones PME was used for the van der Waals potential with the same cut-off and all bonds constrained.A subsequent relaxation of the systems in the NPT ensemble is performed applying Nose-Hoover's thermostat [42,43] along with Parrinello-Rahman's barostat [44,45] for 3 µs, which was found sufficient to account for the potential occurrence of phase transitions.Time couplings of the barostat and the thermostat is 0.5 ps.The pressure is 1 bar, and the isotropic scaling method is chosen as it gives better results than the semi-isotropic one except for the average area per lipid.The heating of the systems for 200 ns to the next temperature started using as initial coordinates the final ones from the preceding temperature is carried out in the NVT ensemble followed again by a run in the NPT one.All 14 temperatures are reached from the preceding temperature, and the simulation run at each temperature is 3 µs.The analyses were performed over the last 100 ns (from 2.9 µs to 3.0 µs ) of the obtained coordinates; velocities and forces at all temperatures are saved every 1 ps in the production run.To double check the results, two additional simulations with different initial velocities were carried out for each system.Figure S2 shows the total energy profile of the system without cholesterol at 271 K, those of the other two bilayers are analogous.

Results and Discussion
A visual analysis of the bilayers in Figure 2 hints that there is a change in the structure of the pure bilayer.At lower temperatures, there is a wave-like pattern that suggests a state similar to P β ′ .As the temperature increases, the wavy shapes gradually vanishes.A less pronounced behavior is observed at low temperatures for the system containing 10 mol% cholesterol.For the system having 30 mol% cholesterol, the bilayer remains flat regardless of the temperature.
To examine the potential occurrence of a phase transition in cholesterol containing SOPC lipids, the main features of a number of thermodynamic and structural quantities were extracted from the simulation data at different temperatures in an interval centered around the experimental melting temperature.The transitional behavior of the systems under consideration is investigated through temperature dependence of the enthalpy and the heat capacities at constant pressure.The changes in the tails ordering is measured via the order parameter and tilt angles.While the movement of the heads and their arrangement is determined through the area per lipid, lateral diffusion coefficient and radial distribution function.

Thermodynamic Properties
Lipid bilayers composed of SOPC undergo a first-order phase transition from L β to L α [46].A stepwise increase is expected to show up in the enthalpy, say H, at the phase transition temperature, while the heat capacity profile C p exhibits a sharp peak at T m .
Following GROMACS manual, the heat capacity of a membrane is determined by the enthalpy fluctuations [47,48] where R is the gas constant, and N is the number of molecules.
For the studied systems, the enthalpy and the heat capacity at constant pressure were calculated using GROMACS tools.The results are shown on Figure 3, where block averages were used to obtain the statistical errors.

Enthalpy of the Phase Transition
The enthalpy profile of the pure system composed of SOPC shows two cusps at 273 K and 277 K, corresponding to −36.66 kJ•mol −1 and −36.00 kJ•mol −1 , respectively.The behavior at 273 K may be traced back to the thermal fluctuations at the transition temperature of the water surrounding the bilayers, while the cusp at the higher temperature may be attributed to the melting of SOPC at T m = 277 K from a gel to the L α state.This result compares fairly well with the melting temperature T m = 279 K found from experimental data [27,32,33].We would like to point out that the hydrated system studied here is similar to the one probed experimentally in Ref. [49].The SOPC system containing 10 mol% cholesterol seems to exhibit an L O phase at 276 K corresponding to an enthalpy of −35.38 kJ•mol −1 .Just like the cholesterol-free bilayer at 272 K, a transition related to water is again observed.The changes in the enthalpy profile are significantly smaller and shifted by one Kelvin down from the SOPC + 0% Chol system, and the presence of cholesterol in the membranes lowers their melting temperature.In the SOPC + 30% Chol system, there is no pronounced change in the slope of the enthalpy, and there is a high probability that the system shows a L O state throughout the considered temperature interval.
In all three lipid bilayers, there is a clear tendency for the enthalpy to increase with temperature.This thermodynamic property increases with increasing cholesterol concentration, as well.Experimental data were obtained via differential scanning calorimetry indicating a transition enthalpy of 36-38 kJ•mol −1 [50][51][52] in fair agreement with our simulation results.The enthalpy and molar heat capacity values for the three systems at 400 K are significantly higher (18-16 kJ•mol −1 ), thus, showing that the systems are in a completely liquid state, while in the investigated temperature interval depending on the cholesterol concentration, the systems may exhibit a transition from S O phase to L α state of the lipids.

Molar Heat Capacity
The heat capacity profiles decrease as the temperature increases.The presence of cholesterol in the systems leads to a systemically higher C p values, following the same trend as the enthalpy, in good agreement with experimental findings [27].For membranes composed of DPPC, DPPE, DPPA, a C p of about 1.20 kJ•mol −1 •K −1 was reported at 283 K [53].In a DOPC/DPPC mixed bilayer, a heat capacity of 0.490 kJ•mol −1 •K −1 was obtained at 273 K [54].
The cusp in the heat capacity profile is yet another evidence for a possible phase transition in the studied bilayers.In the system without cholesterol there are two pronounced peaks of C p at 273 K and 277 K with values 0.967 kJ•mol −1 •K −1 and 0.865 kJ•mol −1 •K −1 , respectively.In the presence of 10 mol% cholesterol in the SOPC membranes, a second sharp peak takes place at 276 K, as noted in the enthalpy, corresponding to 1.075 kJ•mol It is worth noting that the increase in cholesterol may lead to the liquefaction of the lipid membranes, and we may expect a decrease in the temperature of the phase transition from The results obtained for the SOPC + 30% Chol system shows that the heat capacity is significantly higher compared to the other two systems.With increasing temperature there is a slight decrease in the C p profile.The computed heat capacity reproduces to a large extent the behavior of the enthalpy.There is a shift in the values due to an increase in the concentration of cholesterol in the model structures, and peaks show up for the pure membrane.Their location coincides with the above mentioned temperatures, with smaller amplitudes.Accordingly, a phase transition was reported within the simulation with Slipids FF.

Effect of Simulation Time
The results presented here are obtained for simulation runs that lasted 3 µs.To study the effect of the simulation duration on the profile of C p , we have recorded the outcome at four different simulation times, namely 1 µs, 2 µs and 3 µs.Analyses were performed on a trajectory from the final 100 ns of each period.The behavior of C p for the three membranes is shown on Figure 4.
It can be seen that the potential transition temperature is attained at the long simulation runs.In the simulation of the cholesterol-free system in the 1 µs long trajectory, there is no peak, while in SOPC + 30% Chol, there is a significant increase in C p at 272 K.In longer trajectories of 2 µs, peaks in the systems SOPC + 0% Chol and SOPC + 10% Chol set in.Even at 3 µs, there are obvious indications of peaks in these systems while in the SOPC + 30% Chol, the profile starts smoothing.Thus, short or insufficiently long simulations could lead to misleading results.
Long relaxation times provide a more accurate description of the behavior of the systems and leads to the equilibration of the bilayers at the targeted temperatures.In the studied structures, clear changes are visible only beyond 2 µs trajectories.The simulation time should take into account the size of the structures (the number of atoms), though it is not reasonable to spend excessive computational time.So far, it can be concluded that atomistic simulations of systems of similar size need simulation runs of about 3 µs.The transitional effects are already clearly pronounced, and the changes should be reflected in the other parameters of the system.Hence, it can be concluded that in the SOPC + 0% Chol system, there are two peaks that might correspond to phase transitions at 273 K and 277 K. Additionally, in a membrane containing 10 mol% cholesterol, there are two possible transition temperatures at 272 K and 276 K, while in the SOPC + 30% Chol system, no phase transition takes place within the temperature interval of interest.For all membranes, the higher values of C p around 273 K may be ascribed to the strong effect of the aqueous environment.

The Behavior of Lipid Tails
To determine the structural changes in the lipid molecules due to a phase transition, the deuterium order parameter S CD [55] and the tilt angles of the tails relative to the membrane surface were calculated.Results were obtained for the saturated (sn-1) and unsaturated (sn-2) tails.The temperature-dependent data for the order parameter for sn-1 are shown on Figure 5 and for sn-2 tails in Figure S3 in Supplementary Materials.The tilt angles for the two SOPC tails are shown on Figure 6.

Order Parameter
In the SOPC lipid bilayer without cholesterol, two more ordered S CD profiles (see e.g., Figure 5) are visible at 271 K and 273 K for sn-1 (saturated tail), while the one at 272 K is slightly more disordered.At low temperatures, membrane folding is observed, and the lipid tails are highly stretched (Figure 2); however, the highest S CD value is reported at 273 K. Strong ordering in the tails is also present at 277 K, which is the possible melting temperature in the studied system.At the other temperatures, the values of S CD are relatively close.The data for sn-2 (unsaturated tail) show the same behavior with values slightly lower than those for the saturated tail due to the presence of a double bond.
In the system containing 10 mol% cholesterol, the profile at 272 K shows that the system is highly ordered.The obtained value for sn-1 is a bit higher than both profiles at 0 mol% cholesterol.The presence of cholesterol leads to a higher degree of arrangement of the tails, but on the other hand, there is also a liquefaction of the bilayer.These are characteristics of an L O phase.The data for the unsaturated tail also repeat this trend but of a lower height.
For SOPC + 30% Chol, both tails show no ordering, and the values at the investigated temperatures are close to each other.The S CD values are between the other two systems.The tails are more disordered than with 10 mol% cholesterol, but the bilayer is the most fluid compared to the other two membranes.In this system, there is no signature of a phase transition.
For the cholesterol-free SOPC membrane at 303 K, a value S CD = 0.33 was obtained [56].The change in the order parameter with cholesterol concentration has been well studied [57,58].The absence of a phase transition in the system containing 30 mol% concentration of cholesterol may be explained by the larger amount of cholesterol that hinders the mobility of tails in the system, and the L O phase is always reached [19].
Let us mention that the obtained behavior does not straightforwardly imply the presence of a phase transition.Moreover, the calculated average parameters of the interrow are greatly underestimated and show rather liquid states of the lipids even in the presence of cholesterol.The Slipids FF, even its improved version [23], has a systemic issue.The profile of the three systems at 400 K show that these are completely disordered.This is evidence that, despite the high S CD values, all the membranes are initially in an ordered state (S O ) at low temperatures and under a transition into a less ordered states at higher temperatures, respectively, L α and L O phases.

Tilt Angles
The tilt angles are calculated relative to a vector connecting the first and last carbon atoms of the lipid chain and the z axis perpendicular to the membrane surface.The vector through cholesterol connects the oxygen atom of the hydroxyl group and the last carbon atom of the molecule.
The data for the tilt angles, shown on Figure 6, reflect the trends reported for S CD , which is more pronounced at the sn-1 tails.In the SOPC + 0% Chol system, the profile of both tails largely replicates the results obtained for the molar heat capacity.There are two maxima (at temperatures 273 K, 277 K), and the angles for the remaining systems are in a narrow range.In SOPC + 10% Chol, there is a slight decrease in the values of the angle due to the intercalation of cholesterol between the lipid tails.Again, there is an overlap with the C p profile, and the two peaks at 272 K and 276 K are clear.As expected, the 30 mol% cholesterol coagulant system has no pronounced peak, and the angle data overlap with that of SOPC + 10% Chol, again due to the placement of the cholesterol molecule.The differences in the systems containing cholesterol at some temperatures fall within the margin of error, and it is not possible to distinguish the behaviors of both systems.
The data for the saturated chain tilt angles of the cholesterol-containing systems overlap.The angles for the unsaturated tail are much larger due to the double bond.Here too, the system SOPC + 0% Chol has noticeable maxima at the named temperatures.For the bilayer containing 10 mol% cholesterol, the height of the maximum is smaller but outside the standard deviation.The data for the SOPC + 30% Chol system are analogous to those for the saturated tails.Similarly, at some temperatures, the profiles of both systems containing cholesterol overlap within the standard errors.
Notice that in contrast to the order parameter, the tilt angle clearly outlines the possible melting temperatures also suggested by the profile of the molar heat capacity.The obtained values for the tilt angles for all systems fall within the experimentally determined interval of tilt angles for a wide range of lipids [59,60].

The Behavior of the Lipid Heads at the Phase Transition
Another important criterion to quantify the degree of order in the membranes is the average area per lipid A L and the mobility of the lipid heads in the bilayer, which is determined by the lateral diffusion coefficient D L .Both quantities are calculated for the phosphorus atom in the hydrophilic part of the molecule.The mean area and its fluctuations give the isothermal area compressibility K A .The data for average area per lipid and isothermal area compressibility are presented in Figure 7. Diffusion coefficient values are presented in Figure 8.

Area per Lipid and Isothermal Area Compressibility
The mean areas of SOPC lipids were calculated using the Qhull package of VMD [61].The calculated average area per lipid (Figure 7) are typical of the liquid state rather than the gel phase.It is worth noting that the obtained values for the areas at 400 K in all three systems are about 65-66 Å 2 and are significantly higher than their counterparts in the investigated temperature interval.There is a rather systematic increase in the mean area per lipids due to a force field and not so much to the arrangement of the lipids.In the systems SOPC + 0% Chol and SOPC + 10% Chol slight decreases are observed for the areas per lipid at the relevant temperatures (277 K and 276 K).This shows that the first reported peak in the C p profile for both membranes is associated with water molecules, while the second one reports the behavior of lipid ones.
With the increase in the concentration of cholesterol, there is a decrease in the areas for SOPC lipid, and the smallest reported areas are in SOPC + 30% Chol.The effect of cholesterol on the bilayers is conclusive, as the reported area changes are outside the calculated standard errors.The effect of cholesterol has also been reported at physiological temperature [62].From X-ray lamellar diffraction for a membrane containing SOPC with 10 mol% cholesterol, an area of 62 Å 2 [63] was obtained.Values for the cholesterol-free system are also close to those reported in the literature [64,65].
The compressibility properties of lipid bilayers can be determined via [9] where k B is the Boltzmann constant, n L the number of lipids in a leaflet and A L is the standard deviation of the mean area per lipid.The effect of cholesterol in all investigated systems is depicted in Figure 7 right.Notice that as the concentration of cholesterol increases, the ability of the bilayer to compress decreases.This is evident from the lower values of SOPC + 30% Chol compared to the other two systems.In SOPC + 0% Chol, we also have the largest strands associated with bilayer shrinkage (see Figure 2).By increasing the temperature of SOPC + 0% Chol, there are fluctuations in the modulus, while for the system with 30 mol% concentration of cholesterol, the values remain more or less similar throughout the temperature interval.The calculated values for the cholesterol-free bilayer are in agreement with those found in the literature for similar systems.For a bilayer composed only of POPC, a lipid also with an unsaturated ring at 303 K with Slipids FF, a compressibility modulus of about 271.1 mN•m −1 was reported [66].Experimentally determined values for area compressibility modulus were obtained for DPPC at 323 K K A = 231 mN•m −1 and for the POCC K A = 180-330 mN•m −1 at 298 K [67].
The obtained area per lipid and the isothermal area compressibility account for the influence of cholesterol, but there is a signature of a phase-transition-like behavior.The visual changes in the cholesterol-free bilayer (Figure 2) and the calculated data varying within the studied temperature range indicate that Slipids FF is able to describe the phase transitions in these systems.

Lateral Diffusion Coefficient
From Einstein's relation [68], the lateral diffusion coefficient, say D L , is calculated.Only the movement of the phosphorus atoms in the hydrophilic head in horizontal directions parallel to the surface of the membrane was taken into account, while that in the vertical direction was excluded.Two minima are observed at temperatures 273 K and 277 K.It is noteworthy that the calculated lateral diffusion coefficient at both temperatures read, respectively, 0.162 × 10 −8 cm 2 s −1 and 0.332 × 10 −8 cm 2 s −1 .As the temperature increases, the usual increase in the lateral diffusion takes place.In the other two systems containing cholesterol, there are no clearly pronounced minima.Increasing the temperature, we observe a slight increase in the diffusion coefficient which indicates a tendency towards liquefaction of the bilayer.Moreover, at high temperatures, there is an increase in the fluctuations of D L .In SOPC + 30% Chol, there is a jump at 283 K.In the calculated data, there is no clear dependence of the lateral diffusion coefficient on cholesterol.Similar conclusions were reported experimentally with NMR spectroscopy [69].The obtained values for the lateral diffusion coefficient are in the experimentally determined interval for the L O phase of lipids [7].For other unsaturated lipids (DOPC and POPC), a 0.09 × 10 −8 cm 2 s −1 coefficient was obtained by NMR [70].
Expectedly, for SOPC + 0% Chol, there are changes in the diffusion coefficient profile that prove the manifestation of a phase transition.The description of the heads' motion using the Slipids FF has a weak dependence on temperature.On the other hand, the simulation run is long enough for the changes to occur as reported in the behavior of the heat capacity.This suggests certain difficulties for Slipids FF to describe quite accurately the parameters related to the heads' dynamics of the SOPC lipid, yet for the cholesterol-free system, there are clear indications of a phase transition.

Radial Distribution Function
Additional insights on the arrangement of the molecules in the bilayer can be gained via the behavior of the radial distribution function (RDF).To this end, the most probable distance between the phosphorus atoms in the lipid heads was calculated.From the RDFs of all systems, we show the height of the first and second peaks in Figure 9 and the RDF profile for SOPC + 10% Chol system in Figure S4 in Supplementary materials.RDF profiles for all systems in the entire temperature range are not significantly different.One can distinguish three peaks, which correspond to a long-range arrangement of the lipids.In all three bilayers, there is a high degree of order of the lipid molecules, which is consistent with the gel and L O phase.The first two peaks indicate two preferred distances for the phosphorus atoms.The most probable distance (position of the peaks) is not affected neither by the temperature nor by the amount of cholesterol in the system.For all systems in the studied temperature interval, the position of the first peak is at 0.642 nm and the second at 0.894 nm.
The height of the peaks shows a slight dependence on temperature.The values for SOPC + 0% Chol slightly decrease with temperature, and there is a smaller fluctuation in the height.A similar behavior is found for SOPC + 10% Chol where the magnitude of height shows a weak dependence on temperature.It is noteworthy that at 282 K, there is a pronounced peak (Figure S4 in Supplementary materials).An arrangement of the heads has been noticed even at their small mobility in the profile of D L .For the other parameters, no such effect has been reported so far, and it is probably due to local changes in the membrane and not to temperature dependence.The values for SOPC + 30% Chol are slightly lower compared to the other two systems, and again, the peak magnitude is not strongly affected by temperature.
It can be concluded that the arrangement of the heads is independent of temperature and the amount of cholesterol.This fact is also confirmed by the obtained data for the lateral diffusion coefficient, again determined for the phosphorus atom.The lipid bilayers are in a gel or in an ordered state characterized by high ordering of the molecules in the membranes.

Conclusions
We investigated the transitional behavior of SOPC bilayers with added cholesterol molecular dynamics with the aid of the Slipids force field.In the quest of a probable temperature-driven phase transition in cholesterol containing SOPC lipid bilayers, we performed a thorough analysis of the thermodynamics and structural properties of three distinct cholesterol concentration systems, namely, pure SOPC and SOPC bilayers containing 10 and 30 mol% cholesterol.We examined the behavior of thermodynamic quantities, such as the enthalpy and molar heat capacity, as well as the structural properties describing the ordering behavior of the lipids' heads (area per lipid, lateral diffusion coefficient, radial distribution function) and tails (order parameter, tilt angle) in a temperature range around the experimental melting temperature.
The behavior of the thermodynamic quantities of the pure system and that with 10 mol% cholesterol shows changes with temperature that are characteristic to the occurrence of temperature-driven phase transitions between two distinct phases.The bilayer with the highest concentration has smooth temperature-dependent thermodynamic quantities.In the calculated enthalpy and heat capacity, two cusps were found at different temperatures.The low temperature (273 K for SOPC + 0% Chol and 272 K for SOPC + 10% Chol) most probably corresponds to a transition of the water environment, while the second one (respectively, 277 K and 276 K) is related to changes in the phase of the corresponding lipid bilayers.A detailed inspection of the structural quantities confirms that phase transition takes place within the temperature range examined here.The calculated order parameter of the tails show singularities at the named temperatures for the pure system, and the tilt angle confirms changes in the arrangement of the molecules.Their behavior is rather characteristic of a liquid phase at all temperatures; on the other hand, the calculations at even higher temperatures show that this effect may be traced back to the very definition of the force field.The tilt angles of both cholesterol-containing systems overlap.In the average areas per lipid, a slight decrease was reported for the pure system and SOPC + 10% Chol, which would point to a probable phase transition from a solid-ordered to a liquid-ordered phase.The other quantities describing the properties related to the lipid heads show no clear trend with the increase in temperature, while the effect of cholesterol is more pronounced.
The results of this study show that the Slipids force field for most investigated quantities is accurate in the description of phase transitions in SOPC lipid bilayers with and without cholesterol in the investigated temperature range that covers the vicinity of the experimental melting temperature.The structural quantities related to the hydrophobic tails show a change in the behavior of the lipids at some specific temperatures.On the other hand, the behavior of the heads shows a weak dependence upon cholesterol.To conclude, it is most likely that a phase transition takes place within the considered interval close to the experimental value.

Supplementary Materials:
The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/chemistry6040031/s1. Figure S1: Side view of SOPC lipid bilayers for the studied concentrations at temperature 400 K; Figure S2: The total energy of the pure SOPC system at temperature 271 K for the entire simulation duration of 3 µs; Figure S3: The order parameter S CD for the unsaturated sn-2 tails for all systems.

Figure 1 .
Figure 1.Schematic structure of the molecules of SOPC and cholesterol.

Figure 2 .
Figure 2. Side view of SOPC lipid bilayers at the studied concentrations for the initial and final temperatures.

Figure 3 .
Figure 3. Enthalpy H (left) and molar heat capacity C p (right) versus temperature for the studied systems.The error bars are smaller than the point size.

Figure 4 .
Figure 4. Changes in the heat capacity profile C p during the molecular dynamics simulations in three time windows of 100 ns duration each.

Figure 5 .
Figure 5.The order parameter S CD for the saturated sn-1 tails for all systems.

Figure 9 .
Figure 9.The magnitude of the first and second peaks in the radial distribution function versus temperature.
Figure S4: Radial distribution function RDF g P−P of the distance between the phosphorus atoms at the initial and final temperatures for SOPC + 10% Chol.Author Contributions: Conceptualization, N.I. and H.C.; Methodology, N.I. and H.C.; Software, N.I. and H.C.; Writing-original draft, N.I.; Writing-review and editing, N.I. and H.C.; Funding acquisition, H.C. All authors have read and agreed to the published version of the manuscript.

Funding:
This work was supported by Grant No BG05M2OP001-1.001-0008funded by the Operational Programme Science and Education for Smart Growth, co-financed by the European Union through the European Regional Development Fund.Institutional Review Board Statement: Not applicable.Informed Consent Statement: Not applicable.Data Availability Statement: The data generated within this research are available upon request.