Understanding the Interactions between Triolein and Cosolvent Binary Mixtures Using Molecular Dynamics Simulations.

Biodiesel is one of the emerging renewable sources of energy to replace fossil-fuel-based resources. It is produced by a transesterification reaction in which a triglyceride reacts with methanol in the presence of a catalyst. The reaction is slow because of the low solubility of methanol in triglycerides, which results in low concentrations of methanol available to react with triglyceride. To speed up the reaction, cosolvents are added to create a single phase which helps to improve the concentration of methanol in the triglyceride phase. In this study, molecular dynamics simulations are used to help understand the role of cosolvents in the solvation of triglyceride (triolein). Six binary mixtures of triolein/cosolvent were used to study the solvation of triolein at 298.15 K. Results of 100 ns simulations at constant temperature and pressure to simulate mixing experiments show that in the first 10 ns all the binary mixtures remain largely unmixed. However, for the cosolvents that are fully miscible with triolein, the partial densities across the simulation boxes show that the systems are fully mixed in the final 10 ns. Some solvents were found to interact strongly with the polar part of triolein, while others interacted with the aliphatic part. The radial distribution functions and clustering of the solvents around triolein were also used as indicators for solvation. The presence of cosolvents also influenced the conformation of triolein molecules. In the presence of solvents that solubilize it, triolein preferred a propeller conformation but took up a trident conformation when there is less or no solubilization. The results show that tetrahydrofuran is the best solvent at solubilizing triolein, followed by cyclopentyl methyl ether, diethyl ether, and hexane. With 1,4-dioxane, the solubility improves with an increase in temperature. The miscibility of a solvent in triolein is aided by its ability to interact with both the polar and nonpolar parts of triolein.

Introduction

The fast consumption of fossil-fuel-derived energy may lead to global depletion of such sources in the near future. Moreover, the increase in the use of fossil fuels is the main cause of climate change. As such, recent research efforts are focused on renewable and low-cost clean energy technologies to replace dwindling fossil fuels and slow down climate change. Biodiesel has emerged as a promising renewable energy source with potential to replace fossil fuels, owing to its low levels of toxicity, whereas it does not produce gases which are harmful to the environment upon combustion.[1] It has therefore been approved to be used as an alternative to conventional diesel and is blended with petroleum diesel in several countries.[2] Biodiesel is produced via a transesterification reaction in which a short-chain alcohol (methanol, ethanol, propanol, or butanol) is reacted with a triglyceride in the presence of a catalyst. The reaction is usually carried out at a temperature of 60 °C with sodium hydroxide or potassium hydroxide as preferred catalysts. Because the reaction is reversible, in order to promote the formation of products, it is performed in the presence of excess alcohol, with the common ratios of oil to alcohol being 1:6, 1:9, and 1:12.[3−5] The reaction has been found to be slow because of the poor miscibility of methanol with triglycerides,[6,7] owing to the differences in molar masses of the reactants, which causes a delay in reaction because of the resulting mass transfer limitation.[2] The low solubility also explains the slow rate of reaction owing to the low concentration of methanol in the triglyceride phase and successful collisions for the reaction to occur which therefore do not happen frequently. To try to solve this predicament of a two-phase system which causes mass transfer lag, several solutions have been proposed, including agitation,[8] high-temperature reactions, use of longer-chain alcohols,[9] and use of cosolvents.[5,8,10−12] Agitation and high-temperature reactions are energy-consuming processes which will increase the cost of biodiesel production. Ethanol, being one carbon longer than methanol, has a higher solubility in triglycerides than methanol, whereas it is also safer and “greener” as a solvent than methanol. However, the use of ethanol in transesterification is limited by the difficulty of separating the glycerol from the produced biodiesel in the presence of ethanol.[13] The use of a cosolvent is attractive because of the potential to recycle the solvent, which can lower the cost of production.[14]

The presence of a cosolvent has been found to help in reducing the mass transfer constraints encountered in the transesterification reaction. It helps in creating a single-phase system which aids in the dissolution of methanol in triglyceride, which not only speeds up the reaction but is also energy efficient because stirring is not required, and it is even possible to perform the reaction at room temperature. Sawangkeaw et al.[15] studied the effect of using hexane and tetrahydrofuran (THF) as cosolvents for the transesterification of vegetable oil, which they found to reduce the viscosity of the palm kernel oil, hence increasing the rate of transesterification. Another important finding was that the addition of the cosolvent produced a single-phase system. However, the reaction was conducted under supercritical conditions, which is a costly process on an industrial scale. Boocook et al.[14,16] used several ethers as cosolvents, including dioxane, THF, and diethyl ether (DEE), in the transesterification of different triglycerides. The solvents were selected depending on how close their boiling points were to the boiling point of methanol, which would make it easier to recover both cosolvents and excess methanol as they can be codistilled and recycled. They were able to carry out the reaction at ambient temperatures without any mechanical stirring, and they found that the separation of glycerol was also four times faster in the presence of cosolvent and that, compared to other ethers, THF was the best cosolvent. These findings were also confirmed by Encicar et al.,[8] who used several cosolvents, including DEE, THF, and acetone, in the transesterification of rapeseed oil. They were also able to achieve high biodiesel yield in short reaction times, even at room temperature, where they found DEE and THF to be the most effective cosolvents, and DEE was easier to recover than THF.

Although THF, DEE, and dioxane have been widely used as cosolvents for transesterification, they have the disadvantage of easily undergoing autoxidation to peroxide, and they are soluble in water. These properties cause a serious safety concern regarding their use as solvents, as over time they tend to form explosive organic peroxides. For this reason, small quantities of peroxide inhibitors are often introduced in ethers to slow down their formation.[17−19] Cyclopentyl methyl ether (CPME) has proven to be useful as a green solvent, as it is stable and, unlike other ethers, does not form peroxides. Additionally, it has low solubility in water, has low toxicity, and is stable under both acidic and basic conditions. It has been used in direct acid-catalyzed transesterification of microalgae as a cosolvent and in purifying the produced biodiesel,[20] although De Jesus et al.[20] only focused on the purity of the product and not the effect of CPME on the reaction rate or glycerol recovery. Even though CPME has a higher boiling point (106 °C) than that of methanol (65 °C), it is nonpolar and can easily penetrate lipid material and completely dissolve it,[21] which makes it a good cosolvent for triglycerides. Gamma valerolactone (GVL) is another green solvent with properties that make it a good cosolvent for biodiesel production. Like CPME, GVL has a boiling point of 205 °C which is higher than that of methanol.[22] It also does not form peroxides, which makes it safer for storage and transportation. GVL can be prepared from waste paper and has been used as an illuminating oil and lighter fluid,[23] where it was found to burn cleanly with low emissions of volatile organic compounds (VOCs). GVL is also a nonpolar solvent which is able to dissolve lipids.[24] Its potential as a blending solvent for biodiesel and petroleum diesel was investigated by Bereczky et al., who found that the presence of GVL reduced the emissions of volatile organic compounds, without affecting engine performance.[25,26] GVL has the potential to be used as a cosolvent, where it has another advantage over other solvents, as it is not necessary to remove GVL at the end of the process because it is itself a fuel already.

Even though there are many experimental studies on the use of cosolvents in the transesterification reaction,[20,27,28] there are still no theoretical investigations on the subject, even though it is important to understand the role of cosolvents in transesterification, which itself is a complicated reversible reaction. There have been no investigations on the molecular-level interactions between the cosolvent and triglycerides, which can, however, be studied using molecular dynamics simulations (MD). MD is a powerful tool that is used to characterize the structural and mechanical properties of systems to gain an in-depth understanding of processes at the atomistic or molecular level. Such direct insight is not always easily obtainable from physical experiments.[29−31]

This paper describes a systematic molecular dynamics study of the physical and structural properties of binary mixtures of the triglyceride triolein with a number of cosolvents. The intermolecular interactions between the triolein and cosolvent can be affected by the polar group and the saturated and unsaturated sites.[32] Triolein is the main constituent of olive oil, and it is obtained from a condensation reaction between glycerol and oleic acid.[32−34] Methyl oleate is the second-most often occurring fatty acid methyl ester (FAMEs) from plant-produced biodiesel.[35] It is interesting to find out how each part of the triglyceride molecule interacts with different cosolvent molecules. Cosolvents act by interacting with both the polar methanol and the nonpolar part of triolein. Since triolein is largely nonpolar due to its aliphatic chains, interaction of the methanol or cosolvent with the aliphatic side chains is essential for the formation of a one-phase triolein/methanol/cosolvent system. The aliphatic part of any cosolvent is therefore critical for this interaction. As the methanol aliphatic part, i.e., just one methyl group is expected to play only an insignificant role, simulations of binary mixtures of triolein and cosolvent molecules can therefore provide important insights into these interactions. This analysis provides information on the role of cosolvent interactions and their relationship to the overall solvation, which will help to provide the basis for the development of a more cost-effective route for the production of biodiesel and give guidance on the properties that make a good cosolvent for transesterification.

Results and Discussion

Partial Densities

Simulated mixing was analyzed by calculating the partial densities of the components across the z-direction of the simulation box using the GROMACS “gmx energy” utility. The results obtained are shown in Figure in which partial densities of triolein (black) and different cosolvents (red) are plotted. Plotted in dotted lines is the average partial density for the first 10 ns of the simulation, while solid lines are partial densities in the last 10 ns. Within the first 10 ns, most of the system is still biphasic, as the three regions along the z direction are still clearly noted. The density at the center of the simulation box is approximately 911 kg/m3, which is similar to the reported literature experimental value for triolein of 907 kg/m3,[10,36] suggesting that the region still consists of triolein only. As the system propagates further, cosolvent molecules move toward the center of the box and exchange with the centrally located triolein molecules and vice versa, causing the partial densities to change. For the solvents that are fully miscible with triolein, the partial density graphs finally flatten out in the last 10 ns. This behavior was observed for CPME, hexane, diethyl ether, and THF. For GVL and dioxane, the graphs do not level off, showing that these two cosolvents are not miscible with triolein or that mixing is slower (Figure (c and f)). To investigate if an increase in the temperature will improve the mixing in triolein/GVL and triolein/dioxane systems, the simulation of these mixtures was repeated at 333.15 K, which is the temperature most commonly used in transesterification reactions.[37,38] The partial densities obtained are shown in Figure S3(a and b). For triolein/dioxane systems, partial densities in the last 10 ns flattened out, especially in the center of the box which originally had only triolein (Figure S3(a)). For the triolein/GVL mixture (Figure S3(b)), there is a change in partial densities, indicating that the higher temperature is improving the movement of both the triolein and GVL molecules. However, the graph has not leveled off, showing that its mixing is still slow even at 333.15 K.

Figure 1

Change in partial density of triolein (black lines) and cosolvents (red lines). The first 10 ns is shown by broken lines, while the last 10 ns is a solid line. (a) THF, (b) CPME, (c) dioxane, (d) DEE, (e) hexane, and (f) GVL.

Solvent Distribution around Triolein Molecules

To understand the interaction of the cosolvent with triolein molecules, densities of cosolvent molecules around different regions of triolein relative to the bulk densities were calculated using TRAVIS.[39] Three regions were selected as shown in Figure , and the geometrical regions are defined as cylinders with 0.600 nm radius. The three regions are defined as the region bounded by the atoms C1–C5 (polar part of triolein where carbonyl oxygen and ester oxygen are located), C6–C15 (region including the C–C double bond), and C13–C19 (entirely saturated part of the hydrocarbon chain). The relative densities reported here are for the entire 100 ns simulation, and running averages are reported. The relative densities are shown in Figure . Apart from dioxane and GVL, all the cosolvent graphs are characterized by three regions: first an initial rapid increase in the relative densities, followed by slowing down, and finally a flat region which represents complete mixing of the cosolvent and the triolein. The equilibration time is approximately 20 ns for DEE, THF, and hexane and 30 ns for CPME. GVL and dioxane do not reach equilibrium, as shown in Figure (e) and Figure (f), in agreement with the results obtained from the partial densities which show that they do not mix with triolein.

Figure 2

Numbering scheme for the triolein R1 chain as used in the text. The R1 chain starts at O1 and ends at C19. R1= R2=R3. For the numbering of all atoms in the triolein molecule, refer to Figure S1.

Figure 3

Density ratio of cosolvents in the geometric region of the cylinder and in the bulk solution. (a) THF, (b) CPME, (c) DEE, (d) hexane, (e) dioxane, and (f) GVL (see Figure for a definition of the geometric region).

Hexane, which is the only alkane used in this study, has been used before as a cosolvent in biodiesel production.[15] As expected, it has higher affinity to the saturated hydrocarbon region, C13–C19 (the local density is approximately equal to bulk density), than to the other two regions (local density is lower than the bulk density in those regions). The flexibility and shape of the hexane molecules also make them easier to pack side by side around the triolein chain which hence leads to a higher density. The van der Waals forces between hexane and the C13–C19 region are stronger, which results in strong attraction.

In THF/triolein mixtures, the local density of THF around the regions C1–C5 and C13–C19 is slightly higher than the bulk density (Figure (a)) after 20 ns of simulation time, showing that these regions have a high affinity for THF molecules.[40] THF is an aprotic solvent with a dielectric constant of 7.6.[41] It is a slightly polar solvent that can dissolve a wide range of nonpolar and polar compounds. As such, it has an affinity to both polar and nonpolar regions of triolein. The attraction in regions C6–C15 is slightly weaker, which is attributed to steric hindrance from double bonds, which makes packing of THF molecules difficult.

The leveling of the graph takes a bit longer in CPME than in THF (Figure (a and b)), showing that CPME mixes more slowly. The relative densities around C1–C5 and C13–C19 are also higher than at the C6–C15 region.[40] Just as with THF, the presence of the oxygen atoms in CPME molecules causes this attraction, which is also responsible for promoting the mixing with triolein. CPME is also a nonpolar solvent,[17] which explains its high affinity to the C13–C19 region. CPME has a methyl group which might be responsible for slower mixing than seen with THF. Diethyl ether displays similar behavior to THF and CPME (see Figure (d)). The DEE molecule is straight, like hexane, which results in compact packing and a preference for the C13–C19 region.

In GVL/triolein and dioxane/triolein mixtures (Figure (f and e)), the local densities around all the triolein regions are much lower than 1, showing that the triolein/triolein and solvent/solvent interactions are favored compared to triolein/solvent interactions. It should be noted that in the case of dioxane the density was still rising even after 100 ns, suggesting that it might mix with triolein after longer simulation times or at higher temperatures. Results of calculations of relative densities at 333.15 K are presented in Figure S3a and Figure S3b. From these, it can be noted that three distinct regions (rapid initial increase, slowing down, and final constant density ratios) that were observed at 298.15 K in CPME, THF, DEE, and hexane are not clearly defined/visible in Figure S3(a). What can be seen is the rapid increase in the density ratio that only appears to be slowing down at around 60 ns. We can therefore conclude that the dioxane still mixes more slowly than THF, hexane, CPME, and DEE, as it does not reach equilibration before 80 ns even at a higher temperature. Dioxane shows affinity for the polar (C1–C5) and aliphatic (C13–C19) regions in triolein, enabled by the ability of dioxane to exist in a combination of chair conformations, which is nonpolar, and a boat conformation, which is polar.[42] The triolein/GVL local density is still much lower than the bulk density even at 333.15 K, whereas the higher temperature does not seem to have improved the attraction between triolein and GVL.

The miscibility of two liquids depends on the strength of interactions between the systems. The arrangement of solvent molecules around a solute can often give an indication about the types of interactions that occur between the molecules. A snapshot of the distribution of THF around the polar regions of triolein is shown in Figure . In general, the THF oxygen atom seems to point away from the triolein atoms, with the THF/triolein interaction occurring primarily via the carbon atoms. Figure (b) is the snapshot of CPME/triolein, where it can be observed that methyl groups and the cyclopentyl rings are the ones pointing toward the triolein molecule. CPME is also able to cluster at both the polar and nonpolar parts of triolein. In order to gain a better understanding of the solvent orientation and interaction with triolein, site-to-site radial distribution functions (RDFs) have been calculated.

Figure 4

Snapshots of THF molecules around triolein (a) and CPME around triolein (b).

Site on Site Radial Distributions

Radial distribution functions (RDFs) are used to provide the probability of finding an observed particle at a certain distance from the reference molecule. An RDF value greater than 1 shows that the probability of finding a particle at that radius from the reference point is greater than the uniform probability.[43] It is calculated using eq .where Pr is the probability density of finding another atom at any point with the distance r from the central atom and ρ0 is the mean number density. RDFs of the center of mass (COM) of the solvent molecules from the central carbon atom C1 of triolein (C1-COM) were calculated and are presented in Figure . The distribution functions display typical behavior for liquids, which show progressively weaker order at longer distances. The strongest peak is centered at 5–6 Å and a second one at a distance of around 9 Å, followed by a smooth, broad peak at a longer distance of around 14 Å. These three peaks are present in all triolein/cosolvent mixtures and are similar to those observed by Tascini et al.[29] when studying the interaction of tri-cis-6-hexadecenoin with water. The first peak for THF is very sharp and appears at 5 Å, showing a strong interaction between THF and triolein. The CPME peak appears broader and slightly shifted to longer distances compared to THF. Peaks appearing at longer distances include those for GVL, hexane, and diethyl ether. The peak for GVL has the lowest intensity, showing that there is a weaker order in the arrangement of GVL molecules around the C1 of triolein. There is also a small shoulder on the left of the C1-COM RDF, indicating a weak interaction between the polar region of triolein with the acidic hydrogen of GVL. The calculated radial distribution functions for the interaction of the nonpolar part of triolein with cosolvents, C11-COM (cosolvents), are presented in Figure . The three peaks that were observed in the C1-COM are still present, but the positions of these RDF peak maxima and minima are slightly different. The intensities of the C11-COM peaks are weaker in GVL, THF, and dioxane compared to C1-COM, but for hexane and DEE, the C11-COM peaks have higher intensity. This behavior is expected, as the shape of the latter two molecules can easily be arranged in an orderly manner around the hydrocarbon chains. For CPME, there were only minor differences in C1-COM and C11-COM. However, we note that for GVL, which is immiscible with triolein, the C1-COM and C11-COM peak maxima are below 1, as expected, and for miscible systems the peak maxima are above 1. As expected, the increase in temperature for dioxane and GVL mixtures from 298.15 K to 333.15 K also increased the intensity of the peaks of both C1-COM and C11-COM, although the shapes and RDF peak maxima do not change, which shows that the coordination shells of the triolein and dioxane are not affected by temperature changes. Only the probability of the interactions changes, as demonstrated by the increase in peak height showing higher correlations (Figure S5). Even though the g(r) values for C1-COM and C11-COM increase with temperature for both systems, they were still lower than 1 (asymptotic value) for triolein/GVL, showing that the system is not fully mixed.

Figure 5

Radial distribution function (RDF) between C1-triolein and the center of mass of cosolvents.

Figure 6

C1-COM and C11-COM radial distributions functions. (a) THF, (b) CPME, (c) DEE, (d) hexane, (e) dioxane, and (f) GVL (see Figure and S1 for atom numbering).

Since THF is the most widely used cosolvent, its interaction with triolein was further investigated. Selected radial distribution functions of the different sites of THF with the polar and nonpolar regions of triolein are shown in Figure . The closest interaction sites are between O2(triolein)–C1(THF) and O2(triolein)–C4(THF), with both sites peaking at about 3.62 Å. These two peaks are strongly indicative of the importance of these sites in interactions of the triolein with THF. The O2(triolein)–C4(THF) RDF indicates that the interaction is probably dipolar in nature. Since the dipole moment vector of THF bisects the oxygen atom and the C3–C4 bond, if the O2(triolein)–C1(THF) peak was due to the same interaction mode, it would occur at a significantly longer distance than 3.62 Å, indicating that there is another interaction mode primarily via O2(triolein)–C1(THF). Further evidence of multiple THF interactions with triolein is the spacing between peaks 1 and 2 of O2(triolein)–C4(THF). These two peaks are separated by 1 Å, which is too small to be due to a second solvation shell. The existence of multiple favorable interaction modes of triolein with THF is very likely, given that triolein itself occurs in different isomeric forms, as will be discussed later. C19(triolein)–C4(THF) and C19(triolein)–C1(THF) peaks are also strong but occur at longer distances of 3.96 and 3.89 Å, respectively. As expected, the interaction of the THF oxygen atoms with the nonpolar region of triolein is very weak, as shown by the C19(triolein)–O1(THF) RDF. In general, all the interactions of triolein via the THF oxygen atom are weak and occur at longer distances.

Figure 7

Site to site radial distribution functions for triolein–THF interactions (inset is the number of THF atoms, and hydrogen atoms are not shown).

It is instructive here to consider experimental studies of cosolvent-based transesterification. Unfortunately, these studies hardly discuss the physicochemical properties of the ternary system but rather focus on the rate of the reaction and the yield of the biodiesel produced. The yield of biodiesel produced is complex to discuss in terms of solvent interaction with oil because the solvents also have an effect on the ease of separation of biodiesel from a reaction mixture, which may affect yields. A study by Boocock[14,16] showed that less THF is required to form a single phase from a triglyceride and methanol mixture than for dioxane and diethyl ether. Its relatively strong interactions with both the polar and nonpolar regions of triolein may be a reason for this experimental observation. This may make it more effective at solubilizing triolein than other solvents. In light of this, we have also calculated several radial distribution functions for CPME, a greener solvent that is more likely to replace THF as an important solvent for biodiesel production. As shown in Figure , the most prominent peaks are O2(triolein)–C6(CPME) and C19(triolein)–C6(CPME). These indicate that CPME interacts with triolein via the methyl group, and it is able to do so strongly at the hydrocarbon part and at the carbonyl oxygen atom of triolein, showing strong polar interaction in agreement with the results from Figure . The position of the ester O1 (triolein) might hinder interactions, as shown by weaker peaks which also appear at longer distances, as seen in O1(triolein)–C3(CPME), O1(triolein)–C5(CPME), and O1(triolein)–O1(CPME). The O1(triolein)–C6(CPME) appears at 3.8 Å, but it is still weak when compared to other interactions involving C6(CPME).

Figure 8

Site to site radial distribution functions for triolein–CPME interactions.

The strong interaction between triolein O2 and C19 with the methyl group in CPME may suggest that the CH3 groups in hexane will interact even more strongly with the triglyceride. A comparison of the RDFs of O2(triolein)–C6(CPME) and C19(triolein)–C6(CPME) and to O2(triolein)–C6(hexane) and C19(triolein)–C1(hexane) (see inset for numbering of hexane atoms) is shown in Figure , which motivated us to calculate atom–atom correlations to find out which interactions of these two solvents with triolein might help in the dissolution. The RDFs were calculated for C19(triolein)–C6(CPME) and C19(triolein)–C6(hexane), as well as the carbonyl oxygen O2 with C6 (CPME) and C6 (hexane). From Figure it can be noted that the first O2–C6 RDF peaks appear at shorter distances in both CPME and hexane than the C19–C6. In both cases, the CPME peaks are stronger than the hexane peaks for the same interaction sites and appear at shorter distances. The fact that C6 is bonded to O in CPME might be helpful to increase the interaction with triolein. The CPME peak heights do not differ as much compared to the hexane peaks, indicating that CPME distribution around the triolein is more uniform than that of hexane, in agreement with the regional density analysis discussed in Section (see Figure (b and d)). These results show the importance of polarization effects in the interaction between triolein and cosolvent, not just with the polar part of triolein but also with the nonpolar part as well.

Figure 9

Site to site radial distribution functions for triolein–CPME interactions and triolein–hexane (C19–C6) and (O2–C6). Inset is the numbering for both CPME and hexane used in the calculations.

Coordination numbers around the first solvation shell, presented in Table , were calculated from the radial distribution functions C1-COM and C11-COM, according to eq .In both cases, the highest coordination numbers are due to CPME/triolein. With the exception of GVL and 1,4-dioxane, higher coordination numbers occur around C11-COM than C1-COM, in agreement with the region-specific analysis discussed at the beginning of Section . The coordination numbers are consistent with cluster size distribution, which was performed as follows. A total of 5000 triolein reference molecules were randomly selected from the molecular dynamics trajectory frames, and for each reference triolein molecule, the coordinates of the atoms of the reference triolein itself and those of any other triolein or cosolvent molecule found within the first solvation shells were saved.[44] The solvation shells were based on the radial distribution functions of the solvent center of mass (COM) about the triolein C1 atom (see Figure ). The minima in the first solvation shells were found to occur at the following radii of 7.65 Å for CPME, 7.24 Å for DEE, 6.67 Å for dioxane, 7.05 Å for THF, 7.86 Å for hexane, and 7.43 Å for GVL. With this procedure, several clusters were found (Figure S6). Many of the clusters contained a single triolein molecule and several cosolvent molecules, and for the sake of brevity, our analysis is focused on such clusters, as shown in Figure . If dissolution or mixing occurs, one might expect to find high occurrences of a single triolein molecule solvated by several cosolvent molecules. From the plots in Figure , we see that CPME, DEE, hexane, and THF have high occurrences of such clusters, where the number of solvent molecules ranges from 4 to 7. The coordination numbers calculated using RDFs were similar to the most abundant clusters (Table ).

Table 1

Calculated Coordination Numbers and the Most Frequently Observed Cluster, the (TRIO)-(cosolvent) Cluster

			most frequently
	calculated coordination numbers (RDFs),	calculated coordination numbers (RDFs),	observed cluster
cosolvent	C1-COM (cosolvent)	C11-COM (cosolvent)	x	y
THF	5.43	6.29	1	6
CPME	6	6.4	1	6
GVL	2.48	2.2	1	1
dioxane	3.98	3.88	1	3
DEE	4.98	6.00	1	5
hexane	4.12	5.59	1	5

Figure 10

Triolein–cosolvent cluster distribution within the first coordination shell. Triolein is the reference molecule.

Orientation of Triolein Molecules in Pure and Solvated Triolein

Next, we investigated the effect of solvent molecules on the conformation of triglyceride chains. Based on the relative orientations of the triglyceride chain, the triglyceride conformers can be classified broadly as being shaped like a propeller, a tuning fork, a chair, or a trident (see illustrations in Figure ).[29,45] Following the classification by Tascini et al.,[29] the relative amounts of these conformers in both pure and solvated triolein were calculated, as shown in Table . In bulk triolein, the trident and propeller are the most common conformations with the percentage of trident only slightly higher at 28%. Triolein has longer chains than sebum triglyceride, as investigated by Tascini et al., which might be responsible for differences in the orientation of the molecules in the bulk. For example, in pure sebum triglyceride, the probability of finding propeller conformers was 27.73%, tuning fork 39.06%, chair 17.58%, and trident 15.63%.[29] The addition of cosolvent appears to affect the distribution of these conformers. When no solvation occurs (GVL and dioxane), the trident is the predominant conformer, where it appears that the solvent approaches the molecule at the C1, as shown by stronger peaks in the C1-COM RDFs (Figure ). This causes all the hydrocarbon chains to face the opposite side from the solvent, thereby causing trident conformation. This conformation allows the hydrocarbon chains to be closer to each other and harder for solvent molecules to approach and solubilize them. Solvents that are soluble in triolein have high propeller as well as trident orientations (THF, CPME, DEE, and hexane), while those that are only slightly soluble or insoluble have a high composition of trident conformations, as reported in Table .

Figure 11

Illustrations of conformation of triolein molecules. (a) Propeller, (b) tuning fork, (c) chair, and (d) trident.

Table 2

Probability of Finding Different Conformers for Molecules in Bulk-TG and a Cosolvent/Triolein Mixture

	triolein	GVL	hexane	dioxane	THF	CPME	DEE
propeller	26.70	21.46	28.18	27.42	27.89	28.84	28.18
tuning fork	23.04	25.54	21.77	20.7	23.45	21.93	22.44
chair	22.48	21.91	21.80	21.22	21.80	21.66	20.43
trident	27.78	31.09	28.25	30.66	26.86	27.57	28.95

Conclusions

We have used molecular dynamics simulations to study the interaction of cosolvents with triolein. The density of pure triolein was calculated first and was found to be 911 kg/m3, which is in agreement with experimental data in the literature. The radial distribution functions of C1-COM and C11-COM showed the presence of three peaks in decreasing intensity, which is typical of liquids. From these results, the miscibility of the solvents and triolein can be related to the intensity of the first peak, with THF being the most miscible followed by CPME > DEE > hexane, whereas the miscibility of 1,4-dioxane improves with an increase in temperature. The coordination numbers for all systems have been obtained by integrating the respective RDFs. The results show that CPME can cluster most around one triolein molecule, followed by THF > DEE > HEX > dioxane and GVL.

The changes in partial densities across the simulation boxes were also used to investigate the mixing of triolein with cosolvents. It was shown that in the first 10 ns all the systems were biphasic. For miscible systems, the partial densities of solvent and triolein densities gradually become uniform, demonstrating a change from a biphasic system to a single-phase system in the last 10 ns. When mixing is slow, the partial density does not become uniform in the last 10 ns. Preferential solvation was noticed in different regions of triolein. It was found that some solvents were more attracted to the polar part of triolein, while others were more attracted to the aliphatic part of triolein. The presence of cosolvents also influenced the conformation of triolein molecules. In the presence of solvents that dissolved it, triolein prefers to exist in the propeller conformation, but when there is less or no solvation, it prefers the trident conformation.

Based on the simulation results alone, one may conclude that THF is the best cosolvent for transesterification because it has strong interactions with the polar and nonpolar regions of triolein. Experimentally, however, the main reason for the use of THF has been that it dissolves triglycerides and that its boiling point is very close to that of methanol, making solvent recovery easy at the end of the process. Consideration of the interaction of cosolvent with different regions of triglycerides as done in this study may be an important factor in the search for new green cosolvents for transesterification. Thus, a potential alternative to THF is CPME which is both safer and greener, as it interacts strongly with polar and nonpolar regions of triglycerides. The nonpolar hexane may be better at solubilizing saturated long-chain triglycerides, albeit being too hazardous.

With regard to the miscibility behavior with other triglycerides other than triolein, we note that since triglycerides differ only on the aliphatic part by chain length and degree of saturation it is expected that triglycerides with shorter aliphatic chains will mix better with more polar solvents. More saturated triglycerides usually have stronger interactions between chains of different triglycerides due to improved packing efficiency. Consequently, nonpolar solvents like hexane are needed to dissolve such triglycerides. Finally, this study has considered only the miscibility of cosolvents with triolein. In transesterification, the single phase is formed between the cosolvent, oil, and methanol. Further studies of these ternary systems are planned to better understand how cosolvents solubilize the triglycerides.

Methods

System Preparation

Molecular dynamics simulations were employed to investigate the interaction between triglycerides and selected cosolvents. Six systems of triolein/solvent systems were used in this study, i.e., triolein/GVL, triolein/CPME, triolein/hexane, triolein/THF, triolein/1,4-dioxane, and triolein/diethyl ether. The structures of all cosolvent molecules are shown in Figure , and triolein is in Figure S1. All the interactions between the molecules were described by the GROMOS54a7 force field.[46,47] The force field parameters, including the charges, were obtained by employing the Automated Topology Builder and repository (ATB) engine.[48] The construction of the configurations, the molecular dynamics simulations, and some of the analyses were performed using the GROMACS software.[49−51] The initial configuration of the pure triolein system was constructed as follows: 300 molecules of triolein were placed in a cube with dimensions of approximately 8 × 8 × 8 nm. The triolein molecules were then energy minimized using the steepest descent algorithm, followed by a short isothermal isochoric dynamics (NVT) simulation to relax the system to 298.15 K. The molecules were then subjected to isothermal isobaric dynamics (NPT) with T and P set to 298.15 K and 1 bar, respectively, for 20 ns in order to equilibrate the system. The pressure was controlled through a Parrinello–Rahman barostat[52] with a compressibility of 4.5 × 10–5 bar–1 and a coupling constant of 2 ps, and the temperature was controlled with a modified Berendsen thermostat[53,54] with a temperature coupling constant of 0.1 ps. The time step for the simulations was set to 1 fs. The Verlet cutoff scheme[55] was used for the van der Waals and short-range electrostatic interactions and was kept at a distance of 1.2 nm. Long-range electrostatic interactions were computed using the Particle Mesh Ewald algorithm.[56,57] After this equilibration, the density of the system was calculated and was found to be 911 kg/m3, which is comparable to the experimental value in the literature of 907 kg/m3[10] (see Figure S2). The equilibrated system was then elongated in the z axis, keeping the triolein molecules at the center of the box. 3600 cosolvent molecules were then added to the empty space within the extended box, i.e., 1800 molecules to each side of the triolein, as shown in Figure , to form the initial configurations.

Figure 13

Structures of cosolvents used in this study.

Figure 12

Snapshot of the initial configuration of triolein (gray) and X solvent (green).

Simulation of Mixtures

Next, the prepared system configurations were energy minimized using the steepest descent algorithm with convergence assumed when the maximum force is smaller than 100 kJ mol–1 nm–1. The minimization was followed by a short NVT simulation of 100 ps to relax the system to 298.15 K. NPT semi-isotropic runs were then performed for 100 ns using a time step of 1.0 fs and temperature of 298.15 K and pressure of 1 bar. The pressure of the system was controlled with the Parrinello–Rahman barostat,[52] with a coupling constant of 2 ps. The compressibility of the simulation box was 4.5 × 10–5 bar–1 along the z direction and was kept incompressible in the x and y directions. The temperature was kept constant using a modified Berendsen thermostat[53,54] with a temperature coupling constant of 0.1 ps. The equations of motions were integrated by means of the leapfrog algorithm.[58] The trajectories and energies of the systems were recorded every 1 ps. The LINCS algorithm is employed to constrain bond lengths.[59] A cutoff radius of 1.2 nm is used for the short-range Lennard-Jones (LJ) interactions, whereas long-range electrostatics were treated with the particle mesh Ewald (PME) method with a truncation at the same distance as the LJ cutoff and a spacing for the PME grid of 0.16 nm. Analytical tail corrections in the potential energy are used to compensate for the truncation in LJ interactions. The semi-isotropic pressure coupling simulation procedure allows the box to change only along the z-axis.

All the images of the systems were visualized with the VMD program,[60] and analysis of the trajectories was achieved using TRAVIS software[39,43] and Gromacs.[50,51]