# Assessing molecular simulation for the analysis of lipid monolayer
reflectometry

A. R. McCluskey<EMAIL_ADDRESS><EMAIL_ADDRESS>Department of Chemistry, University of Bath, Claverton Down, Bath, BA2 7AY, UK
Diamond Light Source, Harwell Campus, Didcot, OX11 0DE, UK  J. Grant
Computing Services, University of Bath, Claverton Down, Bath, BA2 7AY, UK  A.
J. Smith  Diamond Light Source, Harwell Campus, Didcot, OX11 0DE, UK  J. L.
Rawle  Diamond Light Source, Harwell Campus, Didcot, OX11 0DE, UK  D. J.
Barlow  Institute of Pharmaceutical Science, King’s College London, London,
SE1 9NH, UK  M. J. Lawrence  Division of Pharmacy and Optometry, University of
Manchester, Manchester, M13 9PT, UK  S. C. Parker  Department of Chemistry,
University of Bath, Claverton Down, Bath, BA2 7AY, UK  K. J. Edler
<EMAIL_ADDRESS>Department of Chemistry, University of Bath, Claverton
Down, Bath, BA2 7AY, UK

###### Abstract

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) available: 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.2600729) under a CC
BY-SA 4.0 license. Reduced experimental datasets are available at
https://researchdata.bath.ac.uk/id/eprint/586, under a CC-BY 4.0 license.

##  I Introduction

Neutron and X-ray reflectometry techniques are popular in the study of layered
structures, such as polyelectrolyte-surfactant mixtures Llamas  _et al._
(2018), lipid bilayer systems Waldie  _et al._ (2018), electrodeposited films
Beebee  _et al._ (2019), and dye-sensitised solar cell materials McCree-Grey
_et al._ (2015). 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 Renaud  _et al._ (2009). 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 Campbell  _et al._
(2011), and the beam deflection system at the I07 beamline of the Diamond
Light Source Arnold  _et al._ (2012).

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 Abelès (1948); Parratt (1954) 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 McCluskey
_et al._ (2019a), 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 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 Scoppola and Schneck (2018).
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 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 Bobone  _et al._ (2013). A four layer-model consisted of the hydrated
SiO2 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 Miller  _et al._
(2003); Anderson and Wilson (2004). 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._ Darré  _et al._ (2015).
Darré _et al._ , 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-3-phosphocholine (DPPC) lipid bilayer to compare
with experimental reflectometry Koutsioubas (2016). 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 Hughes  _et al._ (2016), 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._ Dabkowska  _et al._ (2014). This 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 Koutsioubas
(2016). 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 lengths scales 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 Pluhackova and
Böckmann (2015). 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 Marrink  _et al._ (2007); Berger  _et al._ (1997).

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._ Marrink  _et al._ (2007). 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
united-atom potential models for lipid simulations is that developed by Berger
_et al._ Berger  _et al._ (1997), 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 Jämbeck
and Lyubartsev (2012). 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 Jämbeck and Lyubartsev (2012), Berger united-atom Berger  _et
al._ (1997), and MARTINI coarse-grained potential models Marrink  _et al._
(2007), 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.

##  II Methodology

###  II.1 Neutron reflectometry measurements

The neutron reflectometry measurements analysed in this work have been
previously published by Hollinshead _et al._ Hollinshead  _et al._ (2009) 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 air-water 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), head-deuterated 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 D2O and air-contrast
matched water (ACMW), where D2O and H2O 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;
$20\text{\,}\mathrm{mN}\text{\,}{\mathrm{m}}^{-1}30\text{\,}\mathrm{mN}\text{\,}{\mathrm{m}}^{-1}40\text{\,}\mathrm{mN}\text{\,}{\mathrm{m}}^{-1}50\text{\,}\mathrm{mN}\text{\,}{\mathrm{m}}^{-1}$.
Table 1 outlines the shorthands used to refer to the different contrast
pairings in this work.

Table 1:  The different contrasts of lipid and water investigated in this work. Shorthand | Lipid contrast | Water contrast
---|---|---
h-D2O | h-DSPC |  D2O
d_13-ACMW |  d_13-DSPC | ACMW
d_13-D2O |  d_13-DSPC |  D2O
d_70-ACMW |  d_70-DSPC | ACMW
d_70-D2O |  d_70-DSPC |  D2O
d_83-ACMW |  d_83-DSPC | ACMW
d_83-D2O |  d_83-DSPC |  D2O

###  II.2 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 Jämbeck
and Lyubartsev (2012), which was used alongside the single point charge (SPC)
water model Berendsen  _et al._ (1987), with a timestep of
$0.5\text{\,}\mathrm{fs}$, the SHAKE, RATTLE, and PLINCS methods were used to
constrain the C-H bonds Miyamoto and Kollman (1992); Hess (2008). 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
Berger  _et al._ (1997); again the SPC water model was used. This potential
model was simulated with an increased timestep of $1\text{\,}\mathrm{fs}$. It
is noted that these timesteps are shorter than those typically used for both
forcefields, and timesteps of up to $2\text{\,}\mathrm{fs}$ have been applied
previously Berger  _et al._ (1997); Jämbeck and Lyubartsev (2012). Finally,
the lowest resolution potential model used was the MARTINI Marrink  _et al._
(2007) alongside the polarisable MARTINI water model Yesylevskyy  _et al._
(2010), to avoid the freezing issues observed previously Koutsioubas (2016).
The MARTINI 4-to-1 heavy atom beading allows for the use of a
$20\text{\,}\mathrm{fs}$ timestep. For the Slipids and Berger potential model
a short-range cut-off of $10\text{\,}\mathrm{\SIUnitSymbolAngstrom}$ was used,
while for the MARTINI potential model the cut-off was extended to
$15\text{\,}\mathrm{\SIUnitSymbolAngstrom}$. All simulations were conducted
with temperature coupling to a heat bath at $300\text{\,}\mathrm{K}$ and a
leap-frog integrator, and run using GROMACS 5.0.5 Berendsen  _et al._ (1995);
Lindahl  _et al._ (2001); van der Spoel  _et al._ (2005); Hess  _et al._
(2008) 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 Yeh and
Berkowitz (1999). A close-packed “wall” of non-interacting dummy atoms was
placed at each side of the simulation cell in the _z_ -direction to ensure
that the atoms could not leave the simulation cell.

The starting simulation structure was generated using the molecular packing
software Packmol Martínez  _et al._ (2009). 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\text{\,}\mathrm{\SIUnitSymbolAngstrom}$ 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\text{\times}{10}^{-5}\text{\,}{\mathrm{bar}}^{-1}$ for the Slipids and
Berger simulations and $3.0\text{\times}{10}^{-4}\text{\,}{\mathrm{bar}}^{-1}$
for the MARTINI simulations. The pressure in the _z_ -dimension was kept
constant at $1\text{\,}\mathrm{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 moved 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\text{\,}\mathrm{bar}$ in the _z_ -dimension.
Following the application of the pair of semi-isotropic barostats, the
thickness of the water layer was typically in the region of
$30\text{\,}\mathrm{\SIUnitSymbolAngstrom}$. The equilibration period was
$1\text{\,}\mathrm{ns}$, following which the $50\text{\,}\mathrm{ns}$ NVT
ensemble production simulations were run, on which all analyses were
conducted.

Figure 1: The DSPC monolayer (a) without water layer and (b) with water layer, visuallised using VMDHumphrey  _et al._ (1996). Figure 2: The experimental surface pressure isotherm for DSPC on water, taken from the work of Kubo _et al._ Kubo  _et al._ (2001). Table 2:  The areas per molecule (APM) associated with particular surface pressures and the size of the _x_ \- and _y_ -cell dimension for a simulation of 100 lipid molecules. $\pi$/$\mathrm{mN}\text{\,}{\mathrm{m}}^{-1}$ |  APM/${\mathrm{\SIUnitSymbolAngstrom}}^{2}$ |  _xy_ -cell length/$\mathrm{\SIUnitSymbolAngstrom}$
---|---|---
20 | 47.9 | 69.1
30 | 46.4 | 68.1
40 | 45.0 | 67.1
50 | 44.6 | 66.0

###  II.3 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 McCluskey  _et al._ (2019a, b). This model is
implemented as a class that is compatible with the Python package refnx Nelson
and Prescott (2019); Nelson  _et al._ (2018) 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, $\sigma_{h}$, and within this layer, some volume
fraction of solvent may intercalate, $\phi_{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, $\sigma_{t}$. A maximum value for the thickness of the tail-layer
was imposed, this value was taken from the Tanford Equation Tanford (1980),

$t_{t}=1.54+1.265n,$ (1)

where $n$ is the number of carbon atoms in the chain, and so for DSPC
$t_{t}=$24.3\text{\,}\mathrm{\SIUnitSymbolAngstrom}$$. The SLD of the tail and
head layers used in the Abelès matrix formalism can, therefore, be found as,

$\text{SLD}_{i}=\frac{b_{i}}{V_{i}}(1-\phi_{i})+\text{SLD}_{s}(\phi_{i}),$ (2)

where, $\text{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. $\phi_{t}=0$ in agreement with
the work of Campbell _et al._ Campbell  _et al._ (2018). 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 Braun  _et al._
(2017),

$\phi_{h}=1-\bigg{(}\frac{d_{t}V_{h}}{V_{t}d_{h}}\bigg{)}.$ (3)

A single value for the interfacial roughness was fitted for all interfaces,
which was limited to be no less than
$3\text{\,}\mathrm{\SIUnitSymbolAngstrom}$, as there is only a single lipid
type in each monolayer Campbell  _et al._ (2018). Therefore, any roughness at
the air-water interface is carried equally through all the layers, in a
conformal fashion Kozhevnikov (2012). The modifications over the previous
implementation were that the tail component volume was constrained, based on
the APM (taken from the surface pressure isotherm),

$V_{t}=d_{t}\text{APM},$ (4)

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\text{\,}{\mathrm{\SIUnitSymbolAngstrom}}^{3}$, in agreement with the
work of Kučerka _et al._ Kučerka  _et al._ (2004) and Balgavý _et al._ Balgavý
_et al._ (2001). A uniform background, limited to lie within
$10\text{\,}\mathrm{\char 37\relax}$ 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 Wormington
_et al._ (1999); Björck (2011). To obtain uncertainties on the fitted model,
Markov chain Monte Carlo sampling, enabled by the emcee package Foreman-Mackey
_et al._ (2013) 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\text{\,}\mathrm{\char 37\relax}$ confidence intervals in both cases). The
Abelès matrix formalism was used to calculate the reflectometry profiles as
described in the ESI.

###  II.4 Simulation-derived analysis

The ESI also includes a Python class that is compatible with refnx Nelson and
Prescott (2019); Nelson  _et al._ (2018) allowing for simulation-derived
reflectometry profiles to be obtained, using a similar method to that employed
in previous work, such as Dabkowska _et al._ Dabkowska  _et al._ (2014). 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\text{\,}\mathrm{\SIUnitSymbolAngstrom}$ for the Slipid
and Berger potential model simulations, with an interfacial roughness between
these layers is defined as $0\text{\,}\mathrm{\SIUnitSymbolAngstrom}$. For the
MARTINI potential model, a layer thickness of
$4\text{\,}\mathrm{\SIUnitSymbolAngstrom}$ was used, with an interfacial
roughness of $0.4\text{\,}\mathrm{\SIUnitSymbolAngstrom}$, a detailed
discussion for the rationale behind this is available in the ESI. Each of the
$50\text{\,}\mathrm{ns}$ production simulations were analysed with a frequency
of $10\text{\,}{\mathrm{ns}}^{-1}$, and the SLD profiles were determined by
summing the scattering lengths, $b_{j}$, for each of the atoms in a given
layer.

$\text{SLD}_{n}=\frac{\sum_{j}{b_{j}}}{V_{n}},$ (5)

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\text{\,}\mathrm{\char
37\relax}$ of the highest $q$-value reflected intensity, and a scale factor
were then determined using refnx.

###  II.5 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,

$\chi^{2}=\sum_{i=1}^{N_{\text{data}}}\frac{[R_{\text{exp}}(q_{i})-R_{\text{sim}}(q_{i})]^{2}}{[\delta
R_{\text{exp}}(q_{i})]^{2}},$ (6)

where $q_{i}$ is a given $q$-vector, which depends on the neutron wavelength
and reflected angle, $R_{\text{exp}}(q_{i})$ is the experimental reflected
intensity, $R_{\text{sim}}(q_{i})$ is the simulation-derived reflected
intensity, and $\delta R_{\text{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, $\phi_{h}$, the volume
of the head group, $V_{h}$, and taking the volume of a single water molecule
to be $29.9\text{\,}{\mathrm{\SIUnitSymbolAngstrom}}^{3}$ (from the density of
water as $997\text{\,}\mathrm{kg}\text{\,}{\mathrm{m}}^{-3}$),

$\text{wph}=\frac{\phi_{h}V_{h}}{29.9-29.9\phi_{h}}.$ (7)

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\text{\,}\mathrm{\char 37\relax}$ 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.

###  II.6 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._ Allen  _et
al._ (2016); Pandit  _et al._ (2003). 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\text{\,}\mathrm{\char
37\relax}$ 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 Gowers  _et al._ (2016); Michaud-Agrawal  _et al._ (2011) and the
scripts that were used can be found in the ESI.

##  III Results & Discussion

Figure 3: A comparison of the reflectometry and SLD profiles obtained from (a)
the chemically-consistent layer model, (b) the Slipid simulation, (c) the
Berger simulation, and (d) the MARTINI simulation, at an APM associated with a
surface pressure of $30\text{\,}\mathrm{mN}\text{\,}{\mathrm{m}}^{-1}$. From
top-to-bottom the contrasts are as follows; d_83-D2O, d_83-ACMW, d_70-D2O,
d_70-ACMW, h-D2O, d_13-D2O, d_13-ACMW. The different contrast reflectometry
profiles have been offset in the _y_ -axis by an order of magnitude and the
SLD profiles offset in the _y_ -axis by
$10\text{\times}{10}^{-6}\text{\,}{\mathrm{\SIUnitSymbolAngstrom}}^{-2}$, for
clarity.

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\text{\,}\mathrm{mN}\text{\,}{\mathrm{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
$\chi^{2}$ for each contrast, average $\chi^{2}$, and standard deviation for
each method are given in Table III for each contrast.

Table 3:  The goodness-of-fit between the calculated and experimental
reflectometry profile at a surface pressure of
$30\text{\,}\mathrm{mN}\text{\,}{\mathrm{m}}^{-1}$.

Contrast | Monolayer model | Slipid | Berger | MARTINI
---|---|---|---|---
h-D2O |  |  |  |

& 154.91 107.69 1427.77 d_13-ACMW 74.39 79.71 72.39 124.33 d_13-D2O 36.10
225.26 87.18 1987.96 d_70-ACMW 112.91 74.30 91.28 345.40 d_70-D2O 183.04
622.06 549.92 1873.57 d_83-ACMW 78.36 134.18 315.34 706.97 d_83-D2O 273.28
331.01 418.31 3128.58 Average$\pm$Standard deviation $113.70\pm 80.04$
$231.63\pm 240.29$ $234.59\pm 232.66$ $1370.66\pm 1046.56$

###  III.1 Traditional analysis

The chemically-consistent model was used to determine the structure of the
lipid monolayer, Table III.1 gives the optimum values for the parameters that
were varied in the model. It is clear from this Table, that as the surface
pressure is increased, as expected (and as found previously Mohwald (1990);
Vaknin  _et al._ (1991)), 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
Hollinshead  _et al._ (2009).

Table 4:  The values for the parameters allowed to vary in the fitting of the
chemically-consistent model, at each surface pressure measured.

& $18.20\pm{0.23}$ $3.02^{+0.08}_{-0.02}$ $35.57\pm{2.99}$
$871.67^{+11.27}_{-11.02}$ 30 $12.27\pm{0.49}$ $18.33\pm{0.24}$
$3.01^{+0.07}_{-0.01}$ $40.38^{+2.29}_{-2.47}$ $850.38\pm{10.95}$ 40
$13.54\pm{0.49}$ $18.60^{+0.22}_{-0.22}$ $3.03^{+0.12}_{-0.03}$
$44.30\pm{2.10}$ $836.95\pm{9.88}$ 50 $14.27\pm{0.46}$ $19.20^{+0.22}_{-0.27}$
$3.10^{+0.22}_{-0.10}$ $46.68\pm{1.77}$ $856.25\pm{12.08}$

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
Bayerl  _et al._ (1990). 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 Lu  _et al._ (1994), this can be observed with the
slight increase between
$20\text{\,}\mathrm{mN}\text{\,}{\mathrm{m}}^{-1}50\text{\,}\mathrm{mN}\text{\,}{\mathrm{m}}^{-1}$.

Hollinshead _et al._ Hollinshead  _et al._ (2009) suggest a tail volume of
$972\text{\,}{\mathrm{\SIUnitSymbolAngstrom}}^{3}$ from the density data.
However, the values found in this work are substantially lower, at $\sim
850\text{\,}{\mathrm{\SIUnitSymbolAngstrom}}^{3}$. This reduction, of $\sim
12\text{\,}\mathrm{\char 37\relax}$, agrees well with the work of Campbell _et
al._ Campbell  _et al._ (2018) and Small Small (1984), which suggest that
under the surface pressure investigated in this work a reduction of the tail
volume of up to $15\text{\,}\mathrm{\char 37\relax}$ may be observed. We
believe that the model layer structure from the chemically-consistent method
provides a satisfactory description of the monolayer structure. However, the
use of an MD-driven analysis method may provide greater insight into the
chemical nature of the monolayer.

###  III.2 MARTINI

It is clear from Figure 3 and Table III, 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 D2O 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. It is noted that this may be reduced through the use of a less-ordered
wall structure Koutsioubas (2016). Alternatively, it may be possible to
completely remove the presence of this structuring through the inclusion of
$\sim 10\text{\,}\mathrm{\char 37\relax}$ 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 Marrink  _et
al._ (2010).

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.65}_{-1.88}$ $\mathrm{\SIUnitSymbolAngstrom}$. This is
significantly less than the $24.3\text{\,}\mathrm{\SIUnitSymbolAngstrom}$
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\text{\,}\mathrm{\SIUnitSymbolAngstrom}$, 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.

###  III.3 Comparison of other simulations

Table III shows that both the Slipid and Berger potential models agree well
with the experimental data, with small values for the $\chi^{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 chemically-consistent, 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^{+1.41}_{-7.39}$ $\mathrm{\SIUnitSymbolAngstrom}$, while for the Berger
simulations a value of $19.80^{+1.59}_{-8.17}$
$\mathrm{\SIUnitSymbolAngstrom}$ was obtained. Neither is quite as large as
the $24.3\text{\,}\mathrm{\SIUnitSymbolAngstrom}$ from the Tanford equation,
however, it should be noted that this value is considered a maximum for the
_fully extended_ carbon tail.

Using the molecular dynamics simulations, and the model layer structure it is
possible to compare the number of water molecules per head group. From the
Slipids and Berger simulations, the number of water molecules per head group
was found to be $6.41^{+1.63}_{-0.76}$ and $5.49^{+0.68}_{-0.53}$
respectively. These are in good agreement with the $7.69\pm{0.76}$ found from
the monolayer model method in conjunction with Equation 7.

The $50\text{\,}\mathrm{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\text{\,}\mathrm{fs}$ simulation timestep could reduce this time
significantly. Additionally, Figure 4 shows the results from the first
$5\text{\,}\mathrm{ns}$ of the Slipid potential model, at an APM associated
with a surface pressure of $30\text{\,}\mathrm{mN}\text{\,}{\mathrm{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.

Figure 4: The reflectometry and SLD profiles obtained from the first
$5\text{\,}\mathrm{ns}$ of the Slipid potential model simulation, at an APM
associated with a surface pressure of
$30\text{\,}\mathrm{mN}\text{\,}{\mathrm{m}}^{-1}$. From top-to-bottom the
contrasts are as follows; d_83-D2O, d_83-ACMW, d_70-D2O, d_70-ACMW, h-D2O,
d_13-D2O, d_13-ACMW. The different contrast reflectometry profiles have been
offset in the _y_ -axis by an order of magnitude and the SLD profiles offset
in the _y_ -axis by
$10\text{\times}{10}^{-6}\text{\,}{\mathrm{\SIUnitSymbolAngstrom}}^{-2}$, for
clarity.

###  III.4 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: The simulation time-averaged intrinsic density profile of the water
molecules (blue dots) and lipid components (head groups: green dots, tail
groups: red dots), where the phosphorus atoms of the lipid heads create the
intrinsic surface at $z=$$0\text{\,}\mathrm{\SIUnitSymbolAngstrom}$, at an APM
associated with a surface pressure of
$30\text{\,}\mathrm{mN}\text{\,}{\mathrm{m}}^{-1}$ and the equivalent
scattering length density from the chemically-consistent model (orange line);
similar data for the other surface pressures can be found in the ESI.

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._ Hollinshead  _et al._ (2009), 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 Kozhevnikov (2012); Campbell  _et al._
(2018). 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 II.6. The values for the mean, $95\text{\,}\mathrm{\char
37\relax}$ 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 III.4, for an APM associated with a
surface pressure of $30\text{\,}\mathrm{mN}\text{\,}{\mathrm{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 $\sim 10\text{\,}\mathrm{\SIUnitSymbolAngstrom}$
for the nitrogen atom. However this decreases slightly within the lipid head,
reaching a value of 8.6 $\mathrm{\SIUnitSymbolAngstrom}$ for the end of the
head group. There is then a substantial decrease noted in the lipid tail,
going from $\sim 8.5\text{\,}\mathrm{\SIUnitSymbolAngstrom}$ at the start of
the tail to $\sim 1.5\text{\,}\mathrm{\SIUnitSymbolAngstrom}$ 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.

Table 5:  The mean, $95\text{\,}\mathrm{\char 37\relax}$ quantile, and their
spread for the _z_ -dimension position of atoms representative of difference
parts of the lipid, at an APM associated with a surface pressure of
$30\text{\,}\mathrm{mN}\text{\,}{\mathrm{m}}^{-1}$.

& 76.6 10.1 Mid-Head 67.7 76.6 9.0 End-Head 70.8 79.3 8.6 Start-Tail 1 72.2
80.3 8.1 Start-Tail 2 73.0 81.7 8.6 Mid-Tail 1 80.9 87.1 6.2 Mid-Tail 2 82.3
87.9 5.6 End-Tail 1 91.1 93.3 2.2 End-Tail 2 92.4 93.5 1.1

##  IV 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.

##  V Author Contributions

The initial experiments were conducted by D.J.B. and M.J.L.; the analysis
methodology was developed by A.R.M. with input from J.G., A.J.S., J.L.R.,
S.C.P. and K.J.E.; A.R.M. wrote the manuscript, with input from all authors.

###### Acknowledgements.

ARM is grateful to the University of Bath and Diamond Light Source for co-
funding a studentship (Studentship Number STU0149). This work benefited from
the computing resources provided by STFC Scientific Computing Department’s
SCARF cluster. We thank Robert D. Barker for insightful discussion.

## References

  * Llamas  _et al._ (2018) S. Llamas, E. Guzmán, A. Akanno, L. Fernández-Peña, F. Ortega, R. A. Campbell, R. Miller, and R. G. Rubio, J. Phys. Chem. C 122, 4419 (2018).
  * Waldie  _et al._ (2018) S. Waldie, T. K. Lind, K. Browning, M. Moulin, M. Haertlein, V. T. Forsyth, A. Luchini, G. A. Strohmeier, H. Pichler, S. Maric, and M. Cárdenas, Langmuir 34, 472 (2018).
  * Beebee  _et al._ (2019) C. Beebee, E. B. Watkins, R. M. Sapstead, V. C. Ferreira, K. S. Ryder, E. L. Smith, and A. R. Hillman, Electrochim. Acta 295, 978 (2019).
  * McCree-Grey  _et al._ (2015) J. McCree-Grey, J. M. Cole, and P. J. Evans, ACS Appl. Mater. Interfaces 7, 16404 (2015).
  * Renaud  _et al._ (2009) G. Renaud, R. Lazzari, and F. Leroy, Surf. Sci. Rep. 64, 255 (2009).
  * Campbell  _et al._ (2011) R. A. Campbell, H. P. Wacklin, I. Sutton, R. Cubitt, and G. Fragneto, Eur. Phys. J. Plus 126, 107 (2011).
  * Arnold  _et al._ (2012) T. Arnold, C. Nicklin, J. Rawle, J. Sutter, T. Bates, B. Nutter, G. McIntyre, and M. Burt, J. Synchotron Radiat. 19, 408 (2012).
  * Abelès (1948) F. Abelès, Ann. Phys. 12, 504 (1948).
  * Parratt (1954) L. G. Parratt, Phys. Rev. 95, 359 (1954).
  * McCluskey  _et al._ (2019a) A. R. McCluskey, A. Sanchez-Fernandez, K. J. Edler, S. C. Parker, A. J. Jackson, R. A. Campbell, and T. Arnold, Phys. Chem. Chem. Phys. xx, xxx (2019a).
  * Scoppola and Schneck (2018) E. Scoppola and E. Schneck, Curr. Opin. Colloid Interface Sci. 37, 88 (2018).
  * Bobone  _et al._ (2013) S. Bobone, Y. Gerelli, M. De Zotti, G. Bocchinfuso, A. Farrotti, B. Orioni, F. Sebastiani, E. Latter, J. Penfold, R. Senesi, F. Formaggio, A. Palleschi, C. Toniolo, G. Fragneto, and L. Stella, (BBA) - Biomembranes 1828, 1013 (2013).
  * Miller  _et al._ (2003) A. F. Miller, M. R. Wilson, M. J. Cook, and R. W. Richards, Mol. Phys. 101, 1131 (2003).
  * Anderson and Wilson (2004) P. M. Anderson and M. R. Wilson, J. Chem. Phys. 121, 8503 (2004).
  * Darré  _et al._ (2015) L. Darré, J. Iglesias-Fernandez, A. Kohlmeyer, H. Wacklin, and C. Domene, J. Chem. Theory Comput. 11, 4875 (2015).
  * Koutsioubas (2016) A. Koutsioubas, J. Phys. Chem. B 120, 11474 (2016).
  * Hughes  _et al._ (2016) A. V. Hughes, F. Ciesielski, A. C. Kalli, L. A. Clifton, T. R. Charlton, M. S. P. Sansom, and J. R. P. Webster, Acta Crystallogr. D 72, 1227 (2016).
  * Dabkowska  _et al._ (2014) A. P. Dabkowska, L. E. Collins, D. J. Barlow, R. Barker, S. E. McLain, M. J. Lawrence, and C. D. Lorenz, Langmuir 30, 8803 (2014).
  * Pluhackova and Böckmann (2015) K. Pluhackova and R. A. Böckmann, J. Phys. Condens. Matter 27, 323103 (2015).
  * Marrink  _et al._ (2007) S. J. Marrink, H. J. Risselada, S. Yefimov, D. P. Tieleman, and A. H. de Vries, J. Phys. Chem. B 111, 7812 (2007).
  * Berger  _et al._ (1997) O. Berger, O. Edholm, and F. Jähnig, Biophys. J. 72, 2002 (1997).
  * Jämbeck and Lyubartsev (2012) J. P. M. Jämbeck and A. P. Lyubartsev, J. Phys. Chem. B 116, 3164 (2012).
  * Hollinshead  _et al._ (2009) C. M. Hollinshead, R. D. Harvey, D. J. Barlow, J. R. P. Webster, A. V. Hughes, A. Weston, and M. J. Lawrence, Langmuir 25, 4070 (2009).
  * Berendsen  _et al._ (1987) H. J. C. Berendsen, J. R. Grigera, and T. P. Straatsma, J. Phys. Chem. 91, 6269 (1987).
  * Miyamoto and Kollman (1992) S. Miyamoto and P. A. Kollman, J. Comput. Chem. 13, 952 (1992).
  * Hess (2008) B. Hess, J. Chem. Theory Comput. 4, 116 (2008).
  * Yesylevskyy  _et al._ (2010) S. O. Yesylevskyy, L. V. Schäfer, D. Sengupta, and S. J. Marrink, PLoS Comput. Biol. 6, e1000810 (2010).
  * Berendsen  _et al._ (1995) H. J. C. Berendsen, D. van der Spoel, and R. van Drunen, Comput. Phys. Commun. 91, 43 (1995).
  * Lindahl  _et al._ (2001) E. Lindahl, B. Hess, and D. van der Spoel, J. Mol. Model. 7, 306 (2001).
  * van der Spoel  _et al._ (2005) D. van der Spoel, E. Lindahl, B. Hess, G. Groenhof, A. E. Mark, and H. J. C. Berendsen, J. Comput. Chem. 26, 1701 (2005).
  * Hess  _et al._ (2008) B. Hess, C. Kutzner, D. van der Spoel, and E. Lindahl, J. Chem. Theory Comput. 4, 435 (2008).
  * Yeh and Berkowitz (1999) I.-C. Yeh and M. L. Berkowitz, J. Chem. Phys. 111, 3155 (1999).
  * Martínez  _et al._ (2009) L. Martínez, R. Andrade, E. G. Birgin, and J. M. Martínez, J. Comput. Chem. 30, 2157 (2009).
  * Humphrey  _et al._ (1996) W. Humphrey, A. Dalke, and K. Schulten, J. Mol. Graph. 14, 33 (1996).
  * Kubo  _et al._ (2001) I. Kubo, S. Adachi, H. Maeda, and A. Seki, Thin Solid Films 393, 80 (2001).
  * McCluskey  _et al._ (2019b) A. R. McCluskey, A. Sanchez-Fernandez, K. J. Edler, S. C. Parker, A. J. Jackson, R. A. Campbell, and T. Arnold, “Lipids_at_airdes (Version 1.0),” http://doi.org/10.5281/zenodo.2577796 (2019b).
  * Nelson and Prescott (2019) A. R. J. Nelson and S. W. Prescott, J. Appl. Crystallogr. 52, 193 (2019).
  * Nelson  _et al._ (2018) A. R. J. Nelson, S. W. Prescott, I. Gresham, and A. R. McCluskey, “Refnx v0.1.1,” http://doi.org/zenodo.2222048 (2018).
  * Tanford (1980) C. Tanford,  _The Hydrophobic Effect: Formation of Micelles and Biological Membranes_ , 2nd ed. (John Wiley & Sons, New York, 1980).
  * Campbell  _et al._ (2018) R. A. Campbell, Y. Saaka, Y. Shao, Y. Gerelli, R. Cubitt, E. Nazaruk, D. Matyszewska, and M. J. Lawrence, J. Colloid Interface Sci. 531, 98 (2018).
  * Braun  _et al._ (2017) L. Braun, M. Uhlig, R. von Klitzing, and R. A. Campbell, Adv. Colloid Interface Sci. 247, 130 (2017).
  * Kozhevnikov (2012) I. V. Kozhevnikov, Crystallogr. Rep. 57, 490 (2012).
  * Kučerka  _et al._ (2004) N. Kučerka, M. A. Kiselev, and P. Balgavý, Eur. Biophys. J. 33, 328 (2004).
  * Balgavý  _et al._ (2001) P. Balgavý, N. Kučerka, V. I. Gordeliy, and V. G. Cherezov, Acta. Phys. Slovaca 51, 53 (2001).
  * Wormington  _et al._ (1999) M. Wormington, C. Panaccione, K. M. Matney, and D. K. Bowen, Philos. Trans. R. Soc. London Ser. A 357, 2827 (1999).
  * Björck (2011) M. Björck, J. Appl. Crystallogr. 44, 1198 (2011).
  * Foreman-Mackey  _et al._ (2013) D. Foreman-Mackey, D. W. Hogg, D. Lang, and J. Goodman, Publ. Astron. Soc. PAc. 125, 306 (2013).
  * Allen  _et al._ (2016) D. T. Allen, Y. Saaka, L. C. Pardo, M. J. Lawrence, and C. D. Lorenz, Phys. Chem. Chem. Phys. 18, 30394 (2016).
  * Pandit  _et al._ (2003) S. A. Pandit, D. Bostick, and M. L. Berkowitz, J. Chem. Phys. 119, 2199 (2003).
  * Gowers  _et al._ (2016) R. Gowers, M. Linke, J. Barnoud, T. Reddy, M. Melo, S. Seyler, J. Domański, D. Dotson, S. Buchoux, I. Kenney, and O. Beckstein, in  _Python in Science Conference_ (Austin, Texas, 2016) pp. 98–105.
  * Michaud-Agrawal  _et al._ (2011) N. Michaud-Agrawal, E. J. Denning, T. B. Woolf, and O. Beckstein, J. Comput. Chem. 32, 2319 (2011).
  * Mohwald (1990) H. Mohwald, Annu. Rev. Phys. Chem. 41, 441 (1990).
  * Vaknin  _et al._ (1991) D. Vaknin, K. Kjaer, J. Als-Nielsen, and M. Lösche, Biophys. J. 59, 1325 (1991).
  * Bayerl  _et al._ (1990) T. M. Bayerl, R. K. Thomas, J. Penfold, A. Rennie, and E. Sackmann, Biophys. J. 57, 1095 (1990).
  * Lu  _et al._ (1994) J. R. Lu, E. A. Simister, R. K. Thomas, and J. Penfold, J. Phys. Condens. Matter 6, A403 (1994).
  * Small (1984) D. Small, J. Lipids Res. 25, 1490 (1984).
  * Marrink  _et al._ (2010) S. J. Marrink, X. Periole, D. P. Tieleman, and A. H. de Vries, Phys. Chem. Chem. Phys. 12, 2254 (2010).