Predicting Lipid Eutectics Using Coarse-Grained Molecular Dynamics

Binary triacylglyceride (TAG) systems are known to exhibit nonlinear/eutectic behavior, and predicting this is not straightforward unless the thermodynamic properties of the constituent TAGs are known. While this kind of behavior has been shown previously using molecular dynamics, this was carried out using an atomistic force field and fully saturated symmetric TAGs. In this study, we simulate the changing melting point of binary unsaturated asymmetric TAG systems using a coarse-grained force field. While the simulated melting temperatures of the various TAG ratios are generally found to be less than the empirical values, the nonlinear/eutectic behavior is reproduced very well for the three different binary TAG systems used. Hence, this opens up the possibility of being able to simulate the behavior of different, unknown TAG systems.


■ INTRODUCTION
−4 That is, if two pure materials, be they elements such as lead and tin 3 or ionic liquids, 4 each with their own defined but different melting point, are mixed together, then the melting point of the mixture cannot be predicted by a simple weighted arithmetic mean.Such nonlinear behavior can extend to a eutectic point, i.e., the melting point of the mixture at a specific ratio of the components is lower than the melting point of either of the pure components.This melting point depression may be desired, such as in lead solder 3 ; however, this may not always be the case.
−7 Their behavior, however, is not easy to predict by simply knowing the composition of the mixture's composition.Simplifying fat systems to a binary TAG system, the melting behavior is generally shown in a phase diagram (such as in Figure 1), where the solidus line is the point at which the mixture starts to melt, and the liquidus line is the point at which the mixture is fully melted out, with the melting temperature range changing with TAG ratio.Binary TAG mixtures can show both nonlinear and eutectic behavior, exemplified in Figure 1A,B, respectively.−15 The behavior exhibited by any given TAG combination, or at which TAG ratio the eutectic point will be found, however, is not immediately obvious to predict by a simple analysis of the TAG structures.This means that the melting point for a mixture of TAGs at a given ratio must be either determined empirically or, as in this study, simulated for.Mathematical models, such as the Hildebrand model, 16 to predict the nonlinear behavior have been applied to binary TAG mixtures 10 ; however, these still require knowledge of the melting temperature and enthalpy of fusion of the pure components, which may not always be available.
In a previous study, we have shown that the melting point of various TAGs can be predicted using the using the COarse-Grained Interchangeable Triglyceride-Optimized (COGI-TO) 17 force field (FF). 18In this study, we have extended this and investigated the nonlinear behavior of binary TAG mixtures of 1-palmitoyl-2-oleoyl-3-stearoyl-sn-glycerol (sn-POSt), 1,3-dipalmitoyl-2-oleoyl-sn-glycerol (sn-POP), and 1,3-distearoyl-2-oleoyl-sn-glycerol (sn-StOSt) using the same FF.These are the three main TAGs found in cocoa butter, and  The Journal of Physical Chemistry B their nonlinear behavior has been studied extensively 8,[10][11][12]19 given their importance in the confectionery industry, thus allowing for an easy comparison of the simulated results with empirical data.■ COMPUTATIONAL METHODS Molecular Dynamics (MD) Simulation Settings. All sulations have been carried out using GROMACS 20 2021.3 with the COGITO FF. 17 All equilibrations were done using an NPT ensemble using a time step of 25 fs and a pressure of 1.01325 bar using a v-rescale thermostat and a Berendsen barostat.Anisotropic pressure coupling, with a compressibility of 1 × 10 −5 bar −1 in the x, y, and z directions, was used.Temperature coupling was set at 1 ps, while pressure coupling was set at 10 ps.The cutoff scheme was set to Verlet, with the Coulomb and vdW cutoff distances set to 1.1 nm.The vdW modifier was set to potential shift, and the electrostatics (Coulomb type) was set to particle-mesh Ewald (PME) (PME order = 4). Al systems were heated from 248 to 348 K over 200 ns, giving a linear heating rate of 0.5 K/ns (5 × 10 8 K/s).
Building Binary Mixtures.The coarse-grained representations of sn-POP, sn-POSt, and sn-StOSt were built by mapping the XRD crystalline structures 21 as per Cordina et al. 17 All binary mixtures were built starting from a perfect crystal of the larger molecule (with size increasing as sn-POP < sn-POSt < sn-StOSt).In all cases, an 800-molecule crystal consisting of a single TAG type was built by stacking 10, 2, and 10 unit cells in the a, b, and c directions, respectively.A specified number of molecules was then selected and replaced at random using a custom Python script (see Supporting Information S.1).The number of molecules replaced was chosen to give crystals consisting of 10−90% of the replacing TAG in the binary TAG crystal in steps of 10%.Every binary crystal made from the same two TAGs, and consisting of the same ratio of TAGs, was built individually, i.e., no two binary TAG crystals were identical.
Creating Voids.Using the perfect binary TAG crystals, "crack" voids were generated to more closely represent reality and facilitate the melting onset.The description of these types of voids is given in our previous work. 18In all cases, the ratio of TAGs was kept to within ±1.5% (absolute) of the starting ratio in the perfect crystal, e.g., if the perfect crystal consisted of 10% sn-POSt and 90% sn-StOSt, then all the crystals with a crack void had a ratio of sn-POSt:sn-StOSt of 8.5− 11.5%:91.5−88.5%.Void generation was carried out in an automated fashion using a custom Python script (see Supporting Information S.1).Given that each starting perfect binary TAG crystal was different and the void generation was also randomized, no two void systems were identical.
Melting Point Determination.The determination of the melting point of any binary TAG system was determined using a combination of methodologies described in previous work. 18,22For any given simulation, the near-neighbor occupancy (NNO) was determined for the whole system over the full trajectory. 22The NNO cutoff distance between the middle oleic chain bead (bead 5 of the oleic chain) was set to 1.2 nm for all cases.The NNO was then summed for the whole system at each time point, reducing the dimensionality of the data from 3 to 2 and plotting the values versus temperature (assuming a perfectly linear increase in temperature along the whole trajectory, thus allowing for an easy time-to-temperature conversion).The melting onset temperature for the system was then determined by iterating along the NNO sum vs temperature plot and determining the point at which the number of NNO data points being lower than their respective predicted interval (PI) was more than 95% of the NNO data points remaining. 18For any given TAG 1:TAG 2 ratio and void size combination, the melting point was determined by a simple arithmetic mean of the melting point onset temperature over 50 simulations (see Supporting Information S.2 for the Python script used for melting point determination).

■ RESULTS AND DISCUSSION
A previous study to determine the melting point of crystals made up of a single type of TAG had shown that the direct heating method was a suitable procedure when using the COGITO FF. 18 Having defects, which in this case are introduced as "cracks" within the system, be it a TAG system or otherwise, has been shown to be very important to avoid superheating of the system and thus avoid overestimation of the melting point.−29 While the proportion of the size of the defect to the crystal size is orders of magnitude larger than that observed in real crystals, this methodology has been shown to be both simple and effective.A more in-depth discussion on this can be found in our previous paper and elsewhere. 18,25Plotting the average melting onset temperature versus the void size in the The Journal of Physical Chemistry B crystal shows a clear decrease in the average melting point with increasing void size until a plateau is reached, followed by a sharp decrease in the melting point at larger void sizes, indicating the mechanical collapse of the system, with the melting point before the drop taken to be the melting point of the system. 18,25This methodology was tested on five different TAGs, including the TAGs being used in this study, with the simulated melting points being in very good agreement with the empirical melting points. 18n this study, a similar methodology was employed to try and determine whether the known nonlinear behavior of binary TAG systems can be reproduced.By creating systems with a mixture of two TAGs, creating a crack void while retaining the ratio of TAGs as close as possible to that in the perfect crystal and heating this system, melting could be observed.In this case, there were two main differences from Cordina et al., 18 namely, the reduced number of simulations used per TAG 1:TAG 2 ratio plus void size system ("TAG ratio/void size combination") and a modification to the melting onset temperature determination for each individual system.
The number of simulations for each TAG ratio/void size combination was reduced from 150 to 50 due to the computational expense involved.Binary combinations of all three TAGs used in this study, ranging from 10 to 90% of a TAG in a system in steps of 10%, with void sizes ranging from 16 to 80 molecules in steps of 8 molecules, result in 12,150 simulations when carrying out 50 simulations for each TAG ratio/void size combination.On inspection of the changing bootstrap error and average melting point with a number of simulation plots, these were determined to be sufficiently converged after 50 runs (see Supporting Information S.3−S.5).In addition, whereas in our previous study, the melting point onset for each simulation was determined from a plot of potential energy versus temperature, 18 in this study, the melting point onset was determined using the system's total NNO instead.We also investigated using the change in volume to determine the melting point onset; however, this metric proved unreliable (see Supporting Information S.6).
The NNO is a quick method to determine whether the same molecules are within a cutoff distance of a given reference molecule in subsequent frames in an MD trajectory.If the same molecules are indeed found to be within this cutoff distance in subsequent frames, this indicates that the reference molecule is in a crystalline state, given the relatively short cutoff distance and the unchanging surrounding molecules, as shown previously. 22The NNO is calculated based on the central oleic bead, whose distance from the same beads on surrounding molecules should only change if the ordered crystalline structure is lost, except for the expected minor fluctuations.In this case, the more robust nature of NNO was judged to be a better metric than the potential energy metric for binary systems.
While the NNO plot (e.g., Figure 2A) gives an easily interpretable visual of the crystalline state or, otherwise, of any given molecule within a system, identifying the exact point at which melting starts requires a rigorous methodology.For this, the methodology using potential energy 18 was adapted by summing up the NNOs of all the molecules in the system for  Orange diamonds indicate empirical liquidus temperatures. 10The solid blue line indicates predicted liquidus. 16Red dots are the simulated melting temperature.
The Journal of Physical Chemistry B each frame, plotting these values versus temperature, and then determining the melting point using the lower predicted interval (PI).
As can be seen in Figure 2, the system is very stable for over half of the simulation, as shown by minimal changes in the NNO plot (Figure 2A), and is also reflected in the stable value for the total NNO (Figure 2B).However, toward the end of the simulation, the system starts to change, observed by a marked decrease in the system's total NNO, which was determined to be the point of melting onset.In the case given in Figure 2, the NNO plot shows that by the end of the simulation, most of the crystal is melted out (darker region of the right-hand side of Figure 2A and a big drop in the total NNO in Figure 2B).This was confirmed by visual inspection of the system at the start (Figure 2C) and after 180 ns (Figure 2D) of the simulation.This gradual melting out of a system was observed previously 18 and was attributed to the much faster heating rate of MD simulations when compared to an empirical measurement.
When plotting the average melting point onset versus the void size for each binary TAG system, a curve can be observed (see Supporting Information S.3−S.5 for all plots); however, the plateau and subsequent drop were not always as clear as for single TAG crystals. 18Taking the 20/80 and 40/60 sn-POSt/ sn-StOSt systems as examples (Figure 3), in the former case, the smooth decrease in the average melting onset temperature with increasing void size is very clear, with a plateau reached at a void size of 56 molecules, followed by a clear drop in the melting point (Figure 3A).In the second case, the plateau (or the melting point onset for that binary TAG system) was determined to be the point at which the difference between the melting points at subsequent void sizes was larger than the previous one.In this case, in the steps from a 16-to 40molecule void, the difference in melting temperature decreases or remains the same, while on going from a 40-to 48-molecule void, the difference in melting temperature increases, giving a discontinuous curve (Figure 3B).This principle was used in all various TAG ratio/void size combinations.
Of all three binary TAG combinations, the average melting onset temperature versus void size curves of the sn-POSt/sn-StOSt (see Supporting Information S.3) were the most difficult to interpret.Taking the 30/70 sn-POSt/sn-StOSt system as an example, the correct void size is not immediately obvious, as there is a marked decrease in melting temperature on going from a 24-to 32-molecule void, which would indicate that the melting point at a 24-molecule void is the melting point for this system.However, when looking at the plot as a whole and taking the other sn-POSt/sn-StOSt systems into consideration, the determined melting temperature at a 32-molecule void size seems to not follow the trend.Going from a 40-to 64molecule void, the curve is very smooth, while when looking at the very smooth 20/80 sn-POSt/sn-StOSt system, the plateau is clearly reached at a 56-molecule void, i.e., at a much larger void size than 24 molecules.Taking this into consideration, the melting point for the 30/70 sn-POSt/sn-StOSt system was thus determined to be at a 64-molecule void.A similar approach was taken for the 80/20 sn-POSt/sn-StOSt system.Given this, the determined melting points for each TAG ratio/void size combination are listed in Table 1.
Plotting these values, and comparing them to empirical measurements 10 and predicted values from the Hildebrand model, 10,16 shows good trend reproduction (Figure 4).The empirical measurements were approximated from Smith et al., 10 and they show the temperature at 5% solid.The predicted liquidus line from the Hildebrand model was calculated using the simulated melting points of the pure TAGs 18 and enthalpies of fusion obtained from the literature. 30,31ll three phase diagrams clearly show the nonlinear mixing behavior of the binary TAG systems, whether from the empirical measurements, the predicted phase diagram, or the simulated temperatures from this study.It is, however, observed that the simulated melting points of the binary systems are generally lower than the empirical or predicted values (Figure 4).However, the general trend is reproduced well.Taking the sn-POP/sn-POSt binary system as an example, a minimum at around 50/50 ratio can be clearly observed (Figure 4A), confirming that the eutectic behavior of this binary system is observed during the simulations.
A possible reason for the lower-than-expected melting points could be extra mini-voids.These are created when, for example, sn-StOSt, a longer chain TAG, is in a crystal with sn-POP, a shorter chain TAG, with the palmitic fatty acid chain represented by one less bead than the stearic chain in the COGITO FF.This is shown in Figure 5, where one can clearly see more space between the ends of the fatty acid chains where sn-POP is present.Given the importance of voids to lower the nucleation melting energy barrier during MD simulations, 18,[23][24][25]27,28 having an artificially created void (the crack void), along with the mini-voids created due to the differing chain lengths, leads to an amplified lowering of the energy barrier, which results in a lower-than-expected melting point.

■ CONCLUSIONS
In this study, we have shown the suitability of the COGITO FF to simulate the nonlinear behavior of binary TAG systems.While this study was carried out on known binary systems, being able to reproduce the nonlinear trends of such systems increases the confidence in the prediction of the behavior of unknown binary TAG systems.
The results in this study show that the predicted melting points of the binary systems are lower than expected, possibly as a result of the increased number of mini-voids created by the The Journal of Physical Chemistry B different TAG chain lengths.Given that the COGITO FF was not specifically parametrized to be able to reproduce melting points of mixed TAG crystals, the ability to correctly identify the binary composition at which the eutectic point is observed is an indication of the broad applicability of this FF.However, this methodology provides a powerful and simple way to determine the nonlinear/eutectic behavior of any binary crystalline TAG system, especially for those where the thermodynamic properties (melting temperature and enthalpy of fusion) are not known and, hence, where mathematical predictive models cannot be used.This study thus provides a methodology to determine binary TAG behavior without any empirical measurements.

■ ASSOCIATED CONTENT
* sı Supporting Information The Supporting Information is available free of charge at https://pubs.acs.org/doi/10.1021/acs.jpcb.3c06297.S.1: Python script to generate perfect binary TAG crystals with random distribution of the TAGs, followed by the creation of crack voids; S.2: Python script to determine the melting point of a single system and any given TAG ratio−void size combination; S.

Figure 1 .
Figure1.Examples of (A) nonlinear mixing behavior and (B) eutectic behavior of two TAGs with different pure melting points.The top line is the liquidus (point at which the mixture is fully melted), and the bottom line is the solidus (point at which the mixture starts to melt).Adapted from Machridachis-Gonzaĺez et al.12Copyright 2020, MDPI.

Figure 2 .
Figure 2. (A) Near-neighbor occupancy plot.The sidebar shows the number of molecules around a given reference molecule.(B) Plot of sum of NNO vs temperature (°C).(C) Minimized (fully crystalline) system (red beads = sn-POSt molecules; green beads = sn-StOSt molecules).(D) System after 180 ns (most of the system is melted, while parts are still crystalline).All plots are for 50% sn-POP/50% sn-StOSt with an 80-molecule crack void (run 10).

Figure 3 .
Figure 3. Plots of average melting onset temperature vs void size for (A) 20/80 sn-POSt/sn-StOSt and (B) 40/60 sn-POSt/sn-StOSt.Orange dots/ line represents the average melting onset temperature.Blue error bars are the 95% confidence interval.The red line indicates the void size/melting point onset for the system.

Figure 5 .
Figure 5. Cross section of a 50% sn-POP (red beads)/50% sn-StOSt (green beads) crystal.The yellow circle shows the area with the ends of the stearic chains in close proximity.The blue circle shows the area with mini-voids created by the presence of sn-POP in a crystal with sn-StOSt.
3: results for different void sizes for different sn-POSt:sn-StOSt ratios; S.4: results for different void sizes for different sn-POP:sn-POSt ratios; S.5: results for different void sizes for different sn-POP:sn-StOSt ratios; S.6: plots of average melting point onset temperature vs void size (determined using volume change) for different sn-POP:sn-POSt ratios; S.7: plots of average melting point onset temperature vs void size (determined using volume change) for different sn-POSt:sn-StOSt ratios; S.8: plots of average melting point onset temperature vs void size (determined using volume change) for different sn-POP:sn-StOSt ratios (PDF) ■ AUTHOR INFORMATION Corresponding Author Tell Tuttle − Department of Pure and Applied Chemistry, University of Strathclyde, Glasgow G1 1XL, U.K.; orcid.org/0000-0003-2300-8921;Phone: +44 141 548 2290; Email: tell.tuttle@strath.ac.uk

Table 1 .
18mulated Melting Point (°C) for All TAG Ratios for All Binary Systems a Simulated temperatures for 0/100 and 100/0 ratio systems, i.e., for the pure TAGs, obtained from Cordina et al.18 a