Aisling C.
Stewart
,
Martin J.
Paterson
and
Stuart J.
Greaves
*
Institute of Chemical Sciences, Heriot-Watt University, Edinburgh EH14 4AS, UK. E-mail: s.j.greaves@hw.ac.uk
First published on 9th September 2021
The CC double bonds in unsaturated fatty acids found within atmospheric aerosols play an important role in tropospheric pollution cycles. An understanding of the affinity of such double bonds for the surface of these aerosols is thus vital to the accurate modelling of reactions occurring at these interfaces. To simulate the structures of the liquid surfaces of a series of C18 fatty acids with increasing degrees of unsaturation molecular dynamics calculations have been carried out using the GROMACS suite of programs, employing the general Amber force field. The relative surface coverage of HCCH, carboxylic acid, methyl and CH2 groups has been probed and found to be significantly different from that predicted based on a purely statistical distribution of groups or assumed based on X-ray diffraction data. HCCH, methyl, and CH2 groups dominate the surface, with the methyl groups being found disproportionally at the surface relative to their bulk concentration. The HCCH surface exposure does not increase linearly as the acids become more unsaturated. The distinct structure found at the interface extends only a few nanometers into the sample, suggesting that surface order only affects the first layer of molecules at the interface, but is key to the exposure of highly reactive chemical functional groups to atmospheric oxidants.
Environmental significanceUnsaturated fatty acids are key components of organic atmospheric aerosols, with oleic acid also being used as a mimic for aerosol droplets in both experimental and computational work. The presence of alkene groups within aerosols plays an important role in the chemistry of these aerosols; the oxidation of unsaturated groups leads to the formation of secondary aerosols. Models of aerosol structure and behaviour often assume that the distribution of functional groups at surfaces will be stoichiometric. This molecular dynamics study of the liquid interfaces of fatty acids highlights the enhanced surface presence of particular functional groups as compared to bulk liquid samples. The implications of this surface preference on the reactivity of atmospheric aerosols are discussed. |
The importance of aerosols as systems of high chemical reactivity within the atmosphere has been well documented.6,8,10–12 The surfaces of aerosols provide sites for the adsorption and reaction of gaseous species, which can drive atmospheric processes such as those involved in the evolution of greenhouse gases and other pollutants13–15 and in the destruction of Antarctic ozone.16–18 In addition to this aerosols are in themselves important atmospheric pollutants, with enhanced concentrations known to have a detrimental effect on human health.15,19,20 Aerosols also have an important role to play as cloud condensation nuclei and for direct and indirect radiative forcing.12,18,21–23 The exact effect of aerosols within atmospheric systems depends on their composition and structure. This is in turn affected by the reactions of these aerosols with atmospheric radicals, which can oxidise the organic components of primary aerosols, emitted from both natural and anthropogenic sources, and secondary organic aerosols (SOAs).18,24
Long chain unsaturated fatty acids, such as oleic acid, are chemically significant components of atmospheric aerosols,25–31 and have been used in laboratory studies of proxies for aerosols.32–35 The presence of carbon–carbon double bonds within such molecules is of particular importance to the chemistry that occurs at the surface of these aerosols. Alkene groups can be oxidised by radicals, such as odd oxygen and odd hydrogen species, chemically ageing the aerosols and generating SOAs.27,30,33,35–38 This altering of the chemical nature of aerosols naturally brings with it a change in their activities in chemical processes within the atmosphere, and being able to model the nature and rates of formation of secondary aerosols is thus of great interest to atmospheric chemists.
As a component in marine, arboreal and anthropogenic aerosols, the oleic acid family has also previously been used in experiments on model systems for aerosols.30,33,35,39–41 Given the highly reactive nature of CC groups,39,42–44 it is of particular interest to carry out comparative studies using oleic acid and other unbranched C18 fatty acids, as such studies allow the investigation of the importance of degree of unsaturation on the rate of reactions without this being significantly affected by differences in mass or steric bulk. Such experiments have shown that the uptake of OH and Cl2 by a wetted surface increases linearly (up to n = 3) with the number of double bonds in the hydrocarbon chain of the fatty acid present.42,43 The importance of carbon–carbon double bonds at the surface of a liquid is thus central to determining the rate at which oxidation can occur. It is therefore of particular relevance to atmospheric chemists to gain an insight into the preference of alkene groups for the surface, as compared to other groups within fatty acids, such as the carboxylic acid group and alkane groups. Studies of model systems for aerosols, however, sometimes make assumptions about fatty acids on the surfaces of these based on work on bulk liquid samples45,46 or solid salts of these acids,47 without explicit knowledge about the distribution and orientation of fatty acid functional groups at the gas–liquid interface. Interfaces are, however, known to be energetically, and as result often chemically, distinct from the bulk of a system7,8,48 and an understanding of the surface preferences of different molecules and functional groups is vital in order to improve the modelling of surfaces, including those of atmospheric relevance.
While interactions at the gas–liquid interface show many similar features to those at the gas–solid interface,7,49,50 the disordered nature of liquids often makes their surface structure complex, meaning predictions of the nature and distribution of functional groups at these interfaces can be hard to predict using simple models.7,51 Molecular dynamics (MD) simulations are a popular way of modelling the interactions of systems containing large numbers of atoms. Indeed, a range of hydrocarbon and ionic liquid surfaces have already been studied using this method.52–56
This paper presents the findings of MD simulations using the general Amber force field (GAFF),57 an all-atom force field, to simulate the structure of surfaces of four members of the oleic acid family with increasing degrees of saturation: oleic acid, linoleic acid, linolenic acid and stearidonic acid (Fig. 1). The concentration of different functional groups at the surface of the liquid, as compared to the bulk, is discussed. Simulations have been carried out at 273 K, 298 K and 333 K, to investigate the effect of increasing energy on the surface preferences of the functional groups within these molecules.
The resulting samples of pure fatty acids underwent 100–200 cycles of energy minimization using the steepest-descent algorithm, until the energy of the system had converged. A constant volume (NVT) equilibration stage was then carried out for 5 ns at 298 K using a velocity Verlet algorithm, with a time step of 1 fs. Following on from this constant pressure (NPT) equilibration was carried out in two stages. The first stage was carried out at 298 K and employed a leap-frog algorithm combined with a modified Berendsen thermostat with a time constant of 0.1 ps for the temperature coupling.68 This was carried out at a time step of 1 fs until the total density of the system had stabilised, with this typically requiring a total length of 5–10 ns. The next pressure equilibration stage was carried out at the temperature of interest using the velocity Verlet algorithm69 and Martyna-Tuckerman–Tobias–Klein pressure coupling.70 This equilibration stage was also carried out for between 5 and 10 ns, until the total density of the system had stabilised (for final bulk density values see Table S1 of the ESI S1†).
Having equilibrated a box of a bulk sample of the mixture this was converted into a ‘slab’, by extending the simulation by a factor of 3 along the direction of the longest side of the box (the ‘z’ direction). This generated interfaces between the liquid and the regions of vacuum above and below it. The slab of molecules was twice as thick as it was wide with volumes of vacuum above and below the slab that each had the same dimensions as the slab itself (see ESI Section S3†). The number of molecules used within the simulation was selected based on prior tests of the interface versus bulk properties of slabs of different thicknesses (see ESI Section S2†). The size of the slab was chosen so as to ensure that this was thick enough for there to be a clear distinction between the surface and the bulk while taking into account the higher computational cost of simulating larger numbers of molecules.
The slab was then subjected to a series of temperature equilibration and annealing cycles, which proceeded as follows. Temperature equilibration was carried out at the temperature of interest. As with the previous temperature equilibration process a velocity Verlet algorithm was used, with a time step of 1 fs. This equilibration process was carried out in 5 ns steps until the distribution of functional groups throughout the slab was found to be the same for two consecutive steps. Typically it took a total of 15–20 ns for this to be the case. The slab was then annealed in order to destroy any preferential positioning of the molecules with respect to either each other or to the surface. This was achieved by carrying out an annealing cycle at a higher temperature, typically around 600 K for a short length of time (≈20 ps). The exact temperature and length of time were chosen on a slab-to-slab basis, as longer times and higher temperatures had a greater effectiveness at destroying the internal structure of the slab, but also risked causing large numbers of molecules to evaporate off the surface of the slab. The distribution of molecules throughout the slab was checked at the end of the annealing step in order to ensure that this had been successful. Following annealing the slab was cooled back to the temperature of interest and then underwent a further temperature equilibration cycle. This was carried out as before, until two of the consecutive 5 ns equilibration steps gave the same distribution of functional groups throughout the slab. Annealing and equilibration cycles repeated until the distributions reached by the slabs at the end of two equilibration stages were the same, thus indicating that an equilibrium distribution had been reached.
Following this equilibration process two ‘production runs’ were carried out, with the slab undergoing two temperature equilibration runs for 20 ns at the temperature of interest, with a final annealing step, as described above, in between the two runs. All calculations employed 3D periodic boundary conditions and particle-mesh Ewald long-range electrostatics. A timeline summarizing the stages in the generation, equilibration and production runs of the slabs is shown in Fig. S4 of the ESI.†
Fig. 2 presents MD snapshots looking down onto the final surfaces of slabs of the four acids at 298 K, as shown using the Visual Molecular Dynamics (VMD, Version 1.9.3) software.72 These were taken at the end of the final 20 ns production runs and show typical equilibrated positions of the molecules with relation to the surface. Colour-coding has been used to highlight the distributions of the functional groups present at the surface. It can be seen that the surface coverage of carboxylic acid groups (blue) is very low in all of these slabs. Visual comparisons of the HCCH and CH2 distributions across the different species are, however, made more complicated, as on moving from oleic acid to stearidonic acid the degree of unsaturation, and therefore the ratio of alkene groups to CH2 groups, increases. As a result the surface of oleic acid seems dominated by CH2 groups (green) and that of stearidonic acid by HCCH groups (black), but it is not clear by inspection whether this is purely based on a statistical distribution of these differing numbers of groups, or if there is an underlying preference for a particular moiety to be found at the interface.
Information about the surface area coverage of different groups can be obtained in a more quantitative manner by the use of solvent accessible surface area (SASA) analysis,73,74 included as gmx sasa in the GROMACS software package. This calculates the total surface area of the sample that would be accessible to a solvent molecule of a given size, with this being simulated as a hard sphere probe of a chosen radius approaching from the vacuum side of the interface. This probe records the atoms that it comes into contact with and the total surface area contact that it has with each of these. This then allows information to be gained about both the number of atoms of each functional group at the surface and the total surface area coverage by this functional group.
Fig. 3 shows the results of carrying out SASA analysis on the four samples, as measured by a probe of radius 0.15 nm. This corresponds to the van der Waals radius of an OH radical,75 one of the most common atmospherically relevant species that are likely to react with aerosols75–77 (for a discussion of how changing the probe size affects the outcome of the SASA analysis see ESI S6†). For each slab, analysis was carried out at 100 ps intervals over the course of the two 20 ns production runs, following the methodology of McKendrick and coworkers.54,78 The graphs shown in Fig. 3 represent the averaged values across each of these intervals and production runs. Additional details of the analysis procedure are given in ESI Sections S6–S8.† Carrying out analysis at 100 ps intervals allowed the surface to be observed over time, giving evidence that the slabs were well equilibrated during the production runs, as the distributions of the groups did not vary significantly over the course of these. The time interval was chosen so as to ensure that enough frames had been taken so as to be representative of the entire run, without sampling the same configuration multiple times (see ESI, Section S7† for details). Analysis was only carried out at times >2.5 ns, to ensure that the atoms had fully relaxed to their equilibrium positions at the temperature of interest prior to probing. When counting those atoms that were considered to be present at the surface a minimum surface area of contact was required in order for the atoms of interest to be counted, with this threshold being set in order to avoid spurious signals from the bulk (ESI S8† describes how this threshold was chosen). The absolute values for Fig. 3 are shown in Tables S3–S5 of the ESI (Section S9†) along with their corresponding error values, which are typically around 3%.
The graphs in Fig. 3 show the SASA coverage of different functional groups within each of the acid samples. The data is ordered so that the degree of unsaturation of the acids increases on moving from left to right across each of the plots. These plots represent the surface area of each of the groups that can come into contact with probe that is similar in size to an OH radical, and thus are related to the chance of those groups reacting in the atmosphere. It should be noted that the probe used in SASA calculations is inert, and thus the analysis here cannot account for any changes to the surface structure caused the polarising abilities of the incoming species, however such interactions are beyond the scope of MD simulations and as such we present here a description of the native structure of the interface.
It can be seen from these graphs that there is only a very small surface area coverage from the COOH groups, with significantly larger amounts coming from the methyl and CH2 groups. The comparison between the COOH and methyl group coverages is particularly interesting, as in each of the molecules there are the same number of methyl and acid groups (one) and these are present at the termini of the molecules, so if the chains were orientated at random, one would expect a similar surface area coverage for the two groups. Indeed, as the van der Waals radius of the COOH group is greater than that of a methyl group, a higher surface area coverage of COOH would be predicted based purely on statistical coverage.
This observed non-statistical distribution of groups at the surface is likely to be based on the energetic preferences of different groups for the surface or the bulk. Positioning the non-polar methyl groups at the surface, with the polar acid groups oriented inwards increases the strength of the interactions within the sample as a whole. Polar groups are able to form stronger intermolecular interactions both with each other and with other non-polar molecules via induction effects. Orienting COOH groups so that they are directed outwards towards an area of vacuum would mean that these could be involved in fewer interactions with other molecules and it is thus favourable to have the methyl groups at the surface, with a corresponding below-statistical presence of each of the other groups at the surface. This effect would be expected to be more prominent at lower temperatures, as higher slab energies would mean that intermolecular interactions within the slab would become less important, and surfaces with more equal coverages of different groups, and thus higher entropies, more favourable. The surface coverage of groups presented in Fig. 3, however, seems to be temperature independent, and thus it appears that this effect is not significant over the temperature range studied.
The surface area coverage of HCCH groups within an aerosol is of interest to atmospheric chemists, as reactions involving these groups play a particularly significant role in the conversion of primary to secondary aerosols. It can be seen in Fig. 3 that the surface presence of the alkene groups varies significantly from acid to acid; this is to be expected, as the degree of unsaturation increases on moving across the series of the oleic acids. Dividing the data for each functional group by the number of those functional groups in the molecule allows for a comparison of trends between the different acids. Fig. 4 shows the results of this; here the percentage coverages from Fig. 3 have been divided by the number of functional group units present within the acid of interest, before the results were renormalized to 100.
Fig. 4 Weighted surface area coverages of different groups, generated by dividing the SASA areas (Fig. 3) by the number of groups present in the molecule of interest and normalizing the results to 100. |
It can be seen from Fig. 4 that the surface preference of COOH is higher in oleic and linoleic acid than in their more saturated counterparts: there is hardly any COOH present at the surface of the more unsaturated acids. The bars for the alkene group coverage are a similar height for each of the acids, suggesting that the coverage of these groups at the surface increases roughly in proportion with their number in the samples as a whole, although surface HCCH coverage is always below statistical, as a result of the enhanced CH3 presence there.
For each of the acids there is a much lower than expected proportional CH2 presence at the surface. This can already be seen in Fig. 3 – although there are more CH2 than any other group in each of the samples, the surface area coverage of these is not overwhelming, with the surface area coverage of CH2 and CH3 often similar, although each of the acid molecules has only one CH3 and between 8 and 14 CH2 groups depending on the species. Fig. 4 shows this contrast more starkly, by dividing the surface area of the CH2 groups by the presence of these in each molecule. It can be seen that there a much lower than expected CH2 presence at the surface. This suggests that there is a preference for the termini of the molecules to be present at, or even protruding from, the interfaces, a result previously seen in united atom MD studies of alkanes.53,79
An additional method to investigate ordering of molecules relative to the interface is to calculate the concentrations of functional groups at different distances from the centre of the slab (z = 0). The gmx density function of GROMACS was used for this. This splits the slab along the z direction into equal slices of a given thickness and calculates the densities of the functional groups of interest within these. Plots of distance from the centre of the slab versus concentration for each of the functional groups are shown in Fig. 5. Here each slab was split into 100 slices, with this number being a compromise between the greater resolution afforded by smaller slice sizes and the need to have enough atoms within each slice for reliable statistical analysis. The graphs shown have been averaged across slices an equal distance from the centre of the slab in the +z and −z directions in order to gain a clearer picture of the preferences of different functional groups for the interface of the slab. Just as for the SASA analysis the results have been symmetrised as the upper and lower surfaces of the slab have identical energetic properties. For each sample the analysis was carried out on the two production runs, with the results shown as the average across these two runs. The error bars represent one standard error of the mean. As for the SASA analysis only times > 2.5 ns were included.
Fig. 5 Partial density analysis for each of the slabs studied, showing the distributions of each of the functional groups studied at different distances from the centre of the slab (z = 0). Each curve is the average of results from two 20 ns production runs, separated by an annealing step, with the error bars corresponding to one standard error about this mean. Each of the curves has been area normalised to 0.5, to take into account the different numbers of each of the functional groups present within the samples. The curves represent the averages over the positive and negative z directions. Colour scheme matches Fig. 2–4. |
On examination of the graphs in Fig. 5 it is clear that the surfaces of the slabs (i.e. z ⪆ 3.5) have a distinct distribution of functional groups which differs from that of the bulk (z = 0 → 3.5). This distribution can be understood by starting at the extreme positive z values and moving in towards the centre of the slab; at the extreme outermost parts of the slabs it is almost exclusively CH3 groups present, with regions of enhanced HCCH and CH2 concentration appearing consistently slightly closer to the centre of the slab, and then the COOH inside this. This supports the conclusions drawn from the SASA analysis (highlighted in Fig. 3 and 4), in which the greatly enhanced presence of methyl groups, and diminished presence of acid groups, at the surface was observed. However, whereas SASA analysis only probes those atoms that are directly present at the interface, partial density analysis allows for information to be gained on the ordering of molecules further into the slab. From Fig. 5 it can be seen that while there is a lower presence of COOH groups at the surface, there is a region of enhanced concentration in the sub-interfacial region at z ≈ 4.3 nm. This corresponds to approximately a single molecule length below the surface. The enhanced sub-surface concentration of COOH coincides with a diminished concentration of methyl in all molecules, and of HCCH in linoleic, linolenic and stearidonic acids. This gives further evidence of the preferential configuration of those molecules directly at the surface, so that the CH3 groups are positioned pointing outwards, as discussed with regards to the SASA analysis. There is, however, only a single peak for each of the functional groups in the positive and negative z directions with the distributions settling down to their bulk values in the centre of the slab. It seems that this preferential ordering of groups is only present in those molecules directly at the interface, a long-distance order extending inwards towards the centre of the slab is not observed. Other than the preference of methyl groups to be at the surface there was no observed ordering of the chains with relation to each other, as is discussed in the ESI S10.† There was also no ordering observed in any directions parallel to the surfaces, showing that the effects discussed above where purely a result of the presence of the interfaces (see ESI S11†).
Solvent accessible surface area analysis was used to probe the coverage of different functional groups present directly at the interface. We find that there was a much higher than statistical number of methyl groups observed at the surface for all acids and at all temperatures studied compared to that predicted based on the stoichiometry of the sample. Conversely the proportions of CH2, HCCH and COOH groups present at the interface was surprisingly low (and essentially zero for COOH in linolenic and stearidonic acids) for all samples. This was postulated to be due to the favourable energetics of having the non-polar CH3 at the interface, with the more polar groups directed towards the centre of the slab. The exposure of HCCH at the interface was seen to increase linearly with the increase in the number of double bonds present in the molecules up to linolenic acid (n = 3) as found experimentally by Li et al.42 and Nah et al.43 However, this trend did not continue to stearidonic acid (n = 4), which displayed a similar level of exposure to linolenic acid.
Partial density analysis has been used in order to probe the concentrations of each of the functional groups at different distances from the interface. This confirms the higher concentration of methyl groups directly at the interface, with enhanced concentrations of each of the other groups in the regions immediately inwards of this. For each of the functional groups only a single peak is observed in the partial density plot, with values for all groups returning to their bulk values a single molecule length in from the interface. This suggests that any non-statistical distribution of functional groups as a result of the presence of the surface is only present in the first layer of molecules at this interface, with no evidence of long-range ordering in the direction perpendicular to the interface.
The presence of non-statistical coverages of groups at the interface is relevant when modelling the reactions that take place at the surface of atmospheric aerosols. The greater affinity of less reactive methyl groups for the surface could lead to slower rates of reaction for these processes than would be calculated using models which assume that the surface functional group concentrations are the same as in the bulk. The effect of the presence of the surface on the orientation of molecules at an interface must therefore be explicitly taken into account when modelling the rates of reactions that occur there. Whilst the presence of water and other species within an aerosol could alter the orientations of fatty acid molecules present within these, as compared to in the pure samples simulated here, the presence of such groups is likely to only enhance this surface preference for methyl groups, as further species that can form hydrogen bonds within the slab would increasingly favour COOH groups being directed inwards. The use of pure fatty acids as model systems for aerosols is indeed common,40,45,80,81 and the functional group surface preferences found here using MD simulations could be confirmed using highly surface sensitive experimental techniques such as reactive atom scattering.7,52,78,82 A technique that can provide surface atom (as opposed to surface layer) resolution using reactive OH probes. This would allow this MD methodology to be used to make predictions about the surfaces of more complex systems, such as those that would be present on the surface of atmospheric aerosols, which are currently not experimentally tractable.
Footnote |
† Electronic supplementary information (ESI) available. See DOI: 10.1039/d1ea00043h |
This journal is © The Royal Society of Chemistry 2021 |