Assessing molecular simulation for the analysis of lipid monolayer reflectometry

Using molecular simulation to aid in the analysis of neutron reflectometry measurements is commonplace. However, reflectometry is a tool to probe large-scale structures, and therefore the use of all-atom simulation may be irrelevant. This work presents the first direct comparison between the reflectometry profiles obtained from different all-atom and coarse-grained molecular dynamics simulations. These are compared with a traditional model layer structure analysis method to determine the minimum simulation resolution required to accurately reproduce experimental data. We find that systematic limits reduce the efficacy of the MARTINI potential model, while the Berger united-atom and Slipids all-atom potential models agree similarly well with the experimental data. The model layer structure gives the best agreement, however, the higher resolution simulation-dependent methods produce an agreement that is comparable. Finally, we use the atomistic simulation to advise on possible improvements that may be offered to the model layer structures, creating a more realistic monolayer model. Usage: Electronic Supplementary Information (ESI) including all analysis/plotting scripts and figure files, allowing for a fully reproducible, and automated, analysis workflow for the work presented is available at https://github.com/arm61/sim_vs_trad (DOI: 10.5281/zenodo.3254719) under a CC BY-SA 4.0 license. Reduced experimental datasets are available at DOI: 10.15125/BATH-00586, under a CC-BY 4.0 license.


Introduction
Neutron and x-ray reflectometry techniques are popular in the study of layered structures, such as polyelectrolyte-surfactant mixtures [1], lipid bilayer systems [2], electrodeposited films [3], and dye-sensitised solar cell materials [4]. Unlike other surface-sensitive techniques, such as atomic force microscopy (AFM) or scanning electron microscopy (SEM), reflectometry methods can investigate buried interfaces in addition to the material surface. This is due to the ability of neutrons and x-rays to probe more deeply into a material than an AFM tip or the electron. Additionally, reflectometry techniques can more easily provide information about the average structure over large regions of material, resulting in significantly improved sampling, compared with microscopy techniques [5]. The growth in popularity of reflectometry techniques can be attributed to the significant development of both neutron and x-ray reflectometry instrumentation, such as FIGARO, the horizontal neutron reflectometer at the ILL [6], and the beam deflection system at the I07 beamline of the Diamond Light Source [7].
Typically, the analysis of a neutron or x-ray reflectometry profile is achieved by the application of the Abelès matrix formalism for stratified media [8,9] to a model layer structure. These layer structures are usually defined by the underlying chemistry of the system, for example, the chemically-consistent method that we previously used [10], which accounts for the chemical linkage between the phospholipid head and tail layers. However, there has been growing interest in the use of molecular dynamics (MD) simulations to inform the development of these layer structures. This is due to the fact that the equilibrium structures for soft matter interfaces, that are often of interest in reflectometry studies, are accessible on all-atom simulation timescales [11]. However, to the authors' knowledge, no work has directly compared different levels of simulation coarse-graining in order to assess the required resolution for the accurate reproduction of a given neutron reflectometry profile.
The use of MD-driven analysis of neutron reflectometry usually involves, either the calculation of the scattering length density (SLD) profile from the simulation or the full determination of the reflectometry profile. In the former case, the calculated SLD profile may be compared with the SLD profile determined from the use of a model layer structure analysis method. Bobone et al used such a method to study the antimicrobial peptide trichogin GA-IV within a supported lipid bilayer [12]. A four layer-model consisted of the hydrated SiO 2 layer, an inner lipid head-region, a lipid tail-region, and an outer lipid head-region. The SLD profile from the MD simulations agreed well with that fitted to the reflectometry data from this model layer structure.
The reflectometry profile was calculated explicitly from classical simulation in the works of Miller et al and Anderson and Wilson [13,14]. In these, an amphiphilic polymer at the oil-water interface was simulated by Monte Carlo and MD respectively, and the neutron reflectometry profile found by splitting the simulation cell into a series of small layers and applying the Abelès matrix formalism. There was good agreement between the experimental and calculated reflectometry, for low interface coverages of the polymer. Another study that has made a direct comparison between the atomistic simulation-derived reflectometry and those measured experimentally is that of Darré et al [15]. In that work, NeutronRefTools was developed to produce the neutron reflectometry profile directly from an MD simulation. The particular system studied was a supported 1, 2-dimyristoyl-sn-glycero-3-phosphocholine (DMPC) lipid bilayer, again good agreement was found between the simulation-derived profile and the experimental measurement. However, the nature of the support required a correction for the head-group hydration to be imposed to achieve this agreement.
Koutsioubas used the MARTINI coarse-grained representation of a 1, 2-dipalmitoyl-sn-glycero-3phosphocholine (DPPC) lipid bilayer to compare with experimental reflectometry [16]. This work showed that the parameterisation of the MARTINI water beads was extremely important in the reproduction of the reflectometry data, as the non-polarisable water bead would freeze into crystalline sheets resulting in artefacts in the reflectometry profiles calculated. The work of Hughes et al studied again a DPPC lipid bilayer system [17], albeit an all-atom representation, that was compared with a supported DPPC lipid bilayer system measured with polarised neutron reflectometry. The SLD profile found from MD was varied to better fit the experimental measurement, resulting in good agreement. Additionally, the ability to vary the SLD profile was used to remove artefacts that arose when the MD simulations were merged with the Abelès matrix formalism. This was done to account for regions present in the experiment that were not modelled explicitly.
In all of the examples discussed so far there is no direct comparison between the reflectometry profile determined from simulation and that from the application of a traditional analysis method. Indeed, the only example, to the authors' knowledge where a direct comparison was drawn is the work of Dabkowska et al [18]. That work compares the reflectometry profile from a DPPC monolayer at the air-water interface containing dimethyl sulfoxide molecules with a similar molecular dynamics simulation parameterised with the CHARMM potential model. The use of multimodal analysis allowed the determination of the position of a concentration of DMSO molecules at a particular region within a monolayer and the orientation of such molecules.
The previously mentioned work of Koutsioubas involved the use of the MARTINI coarse-grained force field to simulations the DPPC bilayer system [16]. The use of atomistic simulation for soft matter systems, such as a lipid bilayer, is undesirable as this requires a huge number of atoms to be simulated, due to the large lengthscales involved. The purpose of simulation coarse-graining is to reduce the number of particles over which the forces must be integrated, additionally by removing the higher frequency bond vibrations, the simulation timestep can also be increased [19]. Together, these two factors enable an increase in both simulation size and length. The use of the MARTINI 4-to-1 coarse-grained and the Berger united-atom (where hydrogen atoms are integrated into the heavier atoms to which they are bound) potential models are particularly pertinent for application to lipid simulations as both were developed with this specific application in mind [20,21]. The MARTINI potential model involves integrating the interactions of every four heavy atoms, i.e. those larger than hydrogen, into beads of different chemical nature. This potential model attempts to simplify the interactions of lipid and protein molecules significantly by allowing for only eighteen particle types, defined by their polarity, charge, and hydrogen-bond acceptor/donor character, which are discussed in detail in the work of Marrink et al [20]. Increasing the simulation resolution gives an united-atom potential mode, where all of the hydrogen atoms are integrated into the heavier atoms to which they are bound. One of the most popular unitedatom potential models for lipid simulations is that developed by Berger et al [21], with the original paper being cited 1500 times at the time of writing. Finally, the all-atom Slipid (Stockholm Lipids) lipid potential model was developed in 2012 by Jämbeck and Lyubartsev [22]. All three of these potential models were designed to model lipid bilayer systems.
It is clear that there is substantial interest in the use of classical simulation, and coarse-graining for the analysis of neutron reflectometry data. However, there has been no work to investigate whether the use of atomistic simulations gives more detailed than is required to reproduce the reflectometry profile accurately or to assess whether the application of a coarse-grained representation is suitable to aid in analysis. In this work, three potential models, with different degrees of coarse-graining; namely the Slipid all-atom [22], Berger united-atom [21], and MARTINI coarse-grained potential models [20], are compared in terms of their ability to reproduce neutron reflectometry data. We consider that this work offers a fundamental insight into the potential model resolution that is necessary to accurately reproduce experimental neutron reflectometry measurements. Furthermore, we use the highest resolution simulations to suggest possible adjustments that may be made to the model layer structure analysis methods that are typically used for the rationalisation of neutron reflectometry.

Neutron reflectometry measurements
The neutron reflectometry measurements analysed in this work have been previously published by Hollinshead et al [23] and full details of the experimental methods can be found in this previous publication. These measurements concern the study of a monolayer of 1, 2-distearoyl-sn-phosphatidylcholine (DSPC) at the airwater interface. The neutron reflectometry measurements were conducted on seven isotopic contrasts of the lipid and water. These contrasts were made up from four lipid types; fully-hydrogenated lipid (h-DSPC), headdeuterated lipid (d 13 -DSPC), tail-deuterated lipid (d 70 -DSPC), and fully-deutered lipid (d 83 -DSPC), were paired with two water contrasts; fully-deuterated water D 2 O and air-contrast matched water (ACMW), where D 2 O and H 2 O are mixed such that the SLD is zero. The pairing of the fully-hydrogenated lipid with ACMW was not used due to the lack of scattering available from such a system. Measurements were conducted at four different surface pressures; 20mN m −1 , 30 mN m −1 , 40 mN m −1 and 50 mN m −1 . Table 1 outlines the shorthands used to refer to the different contrast pairings in this work.

Molecular dynamics simulations
The DSPC monolayer simulations were made up of lipid molecules modelled with three potential models, each of a different particle grain-size. The Slipids potential model is an all-atom representation of the lipid molecules [22], which was used alongside the single point charge (SPC) water model [24], with a timestep of 0.5 fs, the SHAKE, RATTLE, and PLINCS methods were used to constrain the C-H bonds [25,26]. The Berger potential model is obtained by the integration of the hydrogen atoms into the heavy atoms to which they are bound, producing a united-atom potential model [21]; again the SPC water model was used. This potential model was simulated with an increased timestep of 1 fs. It is noted that these timesteps are shorter than those typically used for both forcefields, and timesteps of up to 2 fs have been applied previously [21,22]. Finally, the lowest resolution potential model used was the MARTINI [20] alongside the polarisable MARTINI water model [27], to avoid the freezing issues observed previously [16]. The MARTINI 4-to-1 heavy atom beading allows for the use of a 20 fs timestep. For the Slipids and Berger potential model a short-range cut-off of 10 Å was used, while for the MARTINI potential model the cut-off was extended to 15 Å. All simulations were conducted with temperature coupling to a heat bath at 300 K and a leap-frog integrator, and run using GROMACS 5.0.5 [28][29][30][31] on 32 cores of the STFC Scientific Computing resource SCARF. The simulation was of a monolayer, therefore the Ewald 3DC correction was applied to allow for the use of x/y-only periodic boundary conditions [32]. A close-packed 'wall' of non-interacting dummy atoms was placed at each side of the simulation cell in the zdirection to ensure that the atoms could not leave the simulation cell.
The starting simulation structure was generated using the molecular packing software Packmol [33]. This was used to produce a monolayer of 100 DSPC molecules, with the head groups oriented to the bottom of the simulation cell. A 6 Å layer of water was then added such that it overlapped the head groups, this was achieved using the solvate functionality in GROMACS 5.0.5. Examples of a dry and a wet monolayer can be seen in figure 1 for the Berger potential model representation. A general protocol was then used to relax the system at the desired surface coverage, reproducing the effects of a Langmuir trough in silico. This involved subjecting the system to a semi-isotropic barostat, with a compressibility of 4.5×10 −5 bar −1 for the Slipids and Berger simulations and 3.0×10 −4 bar −1 for the MARTINI simulations. The pressure in the z-dimension was kept constant at 1 bar, while it was increased in the x-and y-dimensions isotropically. This allowed for the surface area of the interface to reduce, as the lipid molecules have a preference to stay at the interface, while the total volume of the system stayed relatively constant, as the water molecules move down to relax the pressure in the z-dimension. When the xy-surface area is reached that is associated with the area per molecule (APM) for each surface pressure, described by the experimental surface pressure-isotherm (figure 2), given in table 2, the coordinates were saved and used as the starting structure for the equilibration simulation. This equilibration simulation involved continuing the use of the semi-isotropic barostat, with the xy-area of the box fixed, allowing the system to relax at a pressure of 1 bar in the z-dimension. Following the application of the pair of semiisotropic barostats, the thickness of the water layer was typically in the region of 30 Å. The equilibration period was 1 ns, following which the 50 ns NVT ensemble production simulations were run, on which all analyses were conducted.

Abelès matrix formalism
To compare with the simulation-derived reflectometry profiles, a modified version of the chemically-consistent surfactant monolayer model previously used in the group was applied [10,36]. This model is implemented as a  class that is compatible with the Python package refnx [37,38] and is made up of two layers; the head-layer at the interface with the solvent and the tail-layer at the interface with the air. The head components have a calculated scattering length, b h , (found as a summation of the neutron scattering lengths of the individual atoms, see table S1 of the ESI) and a component volume, V h . These make up a head-layer with a given thickness, d h , and interfacial roughness, σ h , and within this layer, some volume fraction of solvent may intercalate, f h . The tail components also have a similarly calculated scattering length, b t , and component volume, V t . This tail-layer also has a given thickness, d t , and interfacial roughness, σ t . A maximum value for the thickness of the tail-layer was imposed, this value was taken from the Tanford equation [39], where n is the number of carbon atoms in the chain, and so for DSPC t 24.3 t = Å. The SLD of the tail and head layers used in the Abelès matrix formalism can, therefore, be found as, where, SLD s is the scattering length density of the subphase (water), and i indicates either the tail-or head-layer; it is assumed that the tail layer contains no solvent or air, i.e. f t =0 in agreement with the work of Campbell et al [40]. To ensure that the number density of the head components and pairs of tail components is the same, the following constraint was included in the model [41], A single value for the interfacial roughness was fitted for all interfaces, which was limited to be no less than 3 Å, as there is only a single lipid type in each monolayer [40]. Therefore, any roughness at the air-water interface is carried equally through all the layers, in a conformal fashion [42]. The modifications over the previous implementation were that the tail component volume was constrained, based on the APM (taken from the surface pressure isotherm), resulting in the monolayer model and simulation-derived models being equally constrained by the calculated surface coverage. Additionally, the head component volume was constrained to a value of 339.5 3 Å , in agreement with the work of Kučerka et al [43] and Balgavý et al [44]. A uniform background, limited to lie within 10 % of the highest q-value reflected intensity, and a scale factor were then determined using refnx to offer the best agreement between the calculated reflectometry profile and that measured experimentally.
In this work, the experimental data from all seven contrasts were co-refined to a single monolayer model, where the head thickness, tail thickness, and interfacial roughness were allowed to vary. The values of the head and tail scattering lengths, along with the super and subphase SLDs are given in table S1. For each co-refinement of seven neutron reflectometry measurements, there were in total five degrees of freedom in the fitting process, and the fitting was performed using a differential evolution algorithm, which has been shown to be particularly useful in the analysis of reflectometry data [45,46]. To obtain uncertainties on the fitted model, Markov chain Monte Carlo sampling, enabled by the emcee package [47] was used to assess the probability distribution function for each parameter. In the MCMC sampling, 200 walkers were used over 1000 iterations, following equilibration of 200 iterations. The use of MCMC sampling allowed for Bayesian inference of the PDF for each of the variables and their respective interactions and the Shapiro test to be used to assess if each PDF was normally distributed. Parameters that were shown to be normally distributed are given with symmetric confidence intervals, while those that failed the Shapiro test are given with asymmetric confidence intervals (95 % confidence intervals in both cases). The Abelès matrix formalism was used to calculate the reflectometry profiles as described in the ESI.

Simulation-derived analysis
The ESI also includes a Python class that is compatible with refnx [37,38] allowing for simulation-derived reflectometry profiles to be obtained, using a similar method to that employed in previous work, such as Dabkowska et al [18]. The Abelès matrix formalism is applied to layers, the SLD of which is drawn directly from the simulation, and the thickness of which is defined. The layer thickness used was 1 Å for the Slipid and Berger potential model simulations, with an interfacial roughness between these layers is defined as 0 Å. For the MARTINI potential model, a layer thickness of 4 Å was used, with an interfacial roughness of 0.4 Å, a detailed discussion for the rationale behind this is available in the ESI. Each of the 50 ns production simulations were analysed with a frequency of ns 10 1 -, and the SLD profiles were determined by summing the scattering lengths, b j , for each of the atoms in a given layer. where, V n is the volume of the layer n, obtained from the simulation cell parameters in the plane of the interface and the defined layer thickness. Again a uniform background, limited to lie within 10 % of the highest q-value reflected intensity, and a scale factor were then determined using refnx.

Comparison between monolayer model and simulation-derived analysis
The agreement between the models from each method was assessed using the following goodness-of-fit metric, following the transformation of the data into Rq 4 space, where q i is a given q-vector, which depends on the neutron wavelength and reflected angle, R exp (q i ) is the experimental reflected intensity, R sim (q i ) is the simulation-derived reflected intensity, and δ R exp (q i ) is the resolution function of the data. The number of water molecules per head group, wph, was also compared between the different methods. This was obtained from the chemically-consistent model by considering the solvent fraction in the head-layer, f h , the volume of the head group, V h , and taking the volume of a single water molecule to be 29.9 Å 3 (from the density of water as 997 kg m −3 ), V wph 29.9 29.9 . 7 h h The number densities, in the z-dimension, for each of the three components (lipids heads, tails, and water) may be obtained directly from the MD simulation trajectory. In order to determine the number of water molecules per headgroup from the MD simulations, a head-layer region was defined as that which contained 60 % of the lipid head number density. The ratio between the water density and the lipid head density was then found within this head-layer region.

Simulation trajectory analysis
In order to use the MD trajectory to guide the future development of the chemically-consistent layer model, it was necessary to investigate the solvent penetration into the head group region of the lipids, the roughness of each interface and the lipid tail length. The solvent penetration was determined using the intrinsic surface approach, as detailed by Allen et al [48,49]. The intrinsic surface approach enables the calculation of the solvent penetration without the effect of the monolayer roughness. This involves taking the z-dimension position of each water molecule with respect to an anchor point, in this work the anchor point was the phosphorus atom of the lipid head that was closest to the water molecule in the xy-plane. The roughness was probed by investigating the variation in positions for the start, middle, and end of each of the head and tail groups. The start of the lipid head was defined as the nitrogen atom, the middle the phosphorus and the end the tertiary carbon, while the start of the lipid tail was defined as the carbonyl carbon atom, the middle the ninth carbon in the tail and the end the final carbon atom in the tail. The distribution of each of these atom types was determined by finding the 95 % quantile for the position in the z-dimension and comparing the spread of the mean and the upper quantile. Finally, the tail length distance, t t was found as the distance from the carbonyl carbon atom to the final primary carbon atom of the lipid tail. All of these analyses used MDAnalysis package [50,51] and the scripts that were used can be found in the ESI. Figure 3 presents the reflectometry and SLD profiles from each of the different methods, both the traditional layer model and the three potential model simulations, at an APM associated with a surface pressure of 30 mN m −1 . This work will focus discussions on the data at this surface pressure, however other surface pressures showed similar trends and can be found in the ESI. In addition, the χ 2 for each contrast, average χ 2 , and standard deviation for each method are given in table 3 for each contrast.

Traditional analysis
The chemically-consistent model was used to determine the structure of the lipid monolayer,   pressure is increased, as expected (and as found previously [52,53]), the overall thickness of the monolayer increases. The thickness increase for the lipid tails may be associated with the straightening of the tails with respect to the interface normal, while the thickness increase of the head groups has been noted previously for DSPC [23].
It would be anticipated that as the surface pressure increases, there would be a corresponding decrease in the volume fraction of solvent in the head group [54]. However, for DSPC, the volume fraction of the solvent appears to be constant (or even increase slightly) with increasing surface pressure. We believe that this is due to the decision to constrain the volume of the lipid head, which may decrease with increasing surface pressure. It has been noted previously that the interfacial roughness will increase with increasing surface pressure [55], this can be observed with the slight increase between 20 mN m −1 to 50 mN m −1 .
Hollinshead et al [23] suggest a tail volume of 972 Å 3 from the density data. However, the values found in this work are substantially lower, at∼850 Å 3 . This reduction, of∼12 %, agrees well with the work of Campbell et al [40] and Small [56], which suggest that under the surface pressure investigated in this work a reduction of the tail volume of up to15 % may be observed. We believe that the model layer structure from the chemicallyconsistent method provides a satisfactory description of the monolayer structure. However, the use of an MDdriven analysis method may provide greater insight into the chemical nature of the monolayer.

MARTINI
It is clear from figure 3 and table 3, that the MARTINI potential model simulations do not effectively reproduce the reflectometry profile, with a clear difference between the model and data. The SLD profiles derived from the MARTINI simulations contain significant dislocations, which lead to artefacts in the resulting reflectometry profile, and therefore the poor agreement with the data.
It is noted that the agreement with the contrasts containing D 2 O is particularly poor. This is most likely an artefact of the structuring effect from the wall at the bottom of the simulation cell on the polarisable MARTINI water, this may be reduced through the use of a less-ordered wall structure [16]. Alternatively, it may be possible to completely remove the presence of this structuring through the inclusion of∼10 % of antifreeze MARTINI beads alongside the normal MARTINI water. However, this method has been noted to also give structuring effects in the presence of ordered walls [57].
Another artefact present in the MARTINI potential model simulations, particularly notable in the d 83 -ACMW and d 70 -ACMW contrasts where the reflectometry fringe at low-q is substantially broader than represented in the data, is that the length of the hydrocarbon tail in the simulation was found to be 16.60 1.88 This is significantly less than the 24.3 Å estimated by the Tanford equation. The reduction in the tail length is due to the nature of the MARTINI's 4-to-1 beading process, as DSPC has a hydrocarbon tail consisting of 18 carbon atoms, and it is not possible to bead such a chain accurately with the MARTINI potential model. In this work, a MARTINI lipid molecule was used with 4 MARTINI beads making up the chain; corresponding to an all-atom hydrocarbon chain of 16 atoms. Applying the Tanford equation to a hydrocarbon chain of such a length results in an anticipated length of 18.7 Å, which agrees better with that found from the simulation.
The requirement for a 4-to-1 beading structure of the MARTINI potential model is a significant weakness in the utility of this potential model in this work. A better method may be limiting experiments to systems that can be modelled exactly or the use of a 2-to-1 beading model. However, we are not aware of an off-the-shelf 2-to-1 coarse-grained potential model that is commonly applied to lipid molecules. Table 3 shows that both the Slipid and Berger potential models agree well with the experimental data, with small values for the χ 2 . While figure 3 shows that the SLD profiles both appear qualitatively similar to those from the model layer structure method. Furthermore, the quality of agreement between these higher-resolution potential models and the model layer structure is relatively similar. However, the model layer structure still offers a better fit to the experimental data than those determined from MD simulation. The result that the model layer structure offers better agreement with the data than those from even all-atom simulation is to be expected, simply by considering the level of constraint present implicitly when determining the reflectometry profile directly from a simulation. While the model layer structure constrains the layer model to be chemicallyconsistent, those from MD simulation have real chemical constraints present in the simulation; e.g. the bonding of atoms, and the non-bonded potentials. The quality of the agreement from this multi-modal analysis technique is sufficient for such a method to be applied regularly in the analysis of neutron reflectometry. Both the Slipid and Berger simulations produced values for the tail length that were in better agreement with the Tanford equation than the MARTINI simulation. For the Slipid simulation, the tail length was found to be 20.17 7 respectively. These are in good agreement with the 7.69±0.76 found from the monolayer model method in conjunction with equation (7).

Comparison of other simulations
The 50 ns production run for the Slipids potential model simulation required 13 days of using 32 cores of the SCARF computing resource. This is non-trivial and therefore not necessarily applicable to all neutron reflectometry experiments. However, we note that the use of a 2 fs simulation timestep could reduce this time significantly. Additionally, figure 4 shows the results from the first 5 ns of the Slipid potential model, at an APM associated with a surface pressure of 30 mN m −1 , and already good agreement with the data is apparent. It is important to keep in mind that this length of simulation required may be extremely system specific. Furthermore, recent developments of molecular dynamics simulations on graphical processing units (GPUs) may allow for significant speed up of the simulations. The nearly as accurate Berger potential model simulations (which are only marginally less accurate) took approximately 2 days, on the same compute resource. This suggests that by using a larger timestep, shorter simulations, and the power of GPU-based molecular dynamics engines it may be possible to run these simulations alongside experiments at large facilities to aid interpretation and analysis.

Using the Slipid simulations to improve the monolayer model
Despite the model layer structure offering a small improvement in agreement over the Slipid potential model simulation, we believe that it is possible to use these chemically constrained MD simulations to improve the existing monolayer model. For example, figure 5 considers the solvent penetration of the lipid heads, using the intrinsic surface approach to remove the effect of the interfacial roughness. It is clear that the plot is not stepwise as is obtained from the uniform solvation model that is commonly used in traditional layer models. Nor is the distribution sigmoidal, as there is a small deviation in the region of the ester group of the lipid heads. This is either due to the hydrophilic interaction of the carbonyl moiety or from pockets of water forming at the air-water interface. Regardless of the mechanism, this suggests that a different solvation model should be considered for a realistic description of the solvent penetration. Figure 5 also shows that, without the presence of the roughness, the distribution of the head groups is relatively normal. This agrees well with the method used previously to fit the experimental data by Hollinshead et al [23], where Gaussian functions were used to describe the lipids head and tail groups. However, the tail group distribution is not distributed in a Gaussian fashion, and this previous method failed to account for any roughness in the interface.
Previous work has suggested that when only a single lipid type is present, the roughness between the layers should be conformal in nature, that is it should be carried uniformly through the layers [40,42]. However, from the investigation of the SLD profiles in figure 3(b) it appears that the roughness between the lipid tails and the air is dramatically different from that at the lipid head-water interface. In an effort to quantify the interfacial roughness in the simulations, we have used the method outlined in section 2.6. The values for the mean, 95 % quantile, and the spread between these for the z-dimension position for atoms representative of the start, middle, and end of each of the lipid head and tails are given in table 5, for an APM associated with a surface pressure of 30 mN m −1 with the other surface pressures available in the ESI. From this table, it is clear that at the very start of the lipid molecule (at the head) the roughness is very large with a value of ∼10 Å for the nitrogen atom. However this decreases slightly within the lipid head, reaching a value of 8.6Å for the end of the head group. There is then a substantial decrease noted in the lipid tail, going from ∼8.5 Å at the start of the tail to ∼1.5 Å at the end. We believe that this indicates the presence of a highly non-conformal roughness in the lipid monolayer of a single lipid type and therefore in future, it is important to consider this possibility in the use of model layer structure method.

Conclusions
This work presents, for the first time, a direct comparison between a traditional method for analysis of neutron reflectometry measurement with analysis derived from a range of all-atom and coarse-grained molecular dynamics simulations; using the all-atom Slipid, the united-atom Berger, and the coarse-grained MARTINI potential models. It was found that the MARTINI potential model did not accurately model the lipid monolayer system, likely, due to the limitations of the 4-to-1 beading system when applied to a carbon tail containing 18 atoms. The Berger and Slipid potential models both showed good agreement with the experimental data, however, the best agreement was obtained by the traditional monolayer model. This would be expected given that the monolayer model contains many more 'degrees of freedom' than the simulations which are severely chemically constrained by the potential model.
Finally, some points from the highest resolution, Slipid, simulations were noted that may be used to improve the traditional monolayer model. For example, it is desirable to model non-uniform solvation of the head group region which would enable a more accurate modelling of the lipid monolayer and the use of a conformal roughness may not be the best constraint to apply.