# A Gaussian model of fluctuating membrane and its scattering properties

Cedric J. Gommes<EMAIL_ADDRESS>Department of Chemical Engineering,
University of Liège B6A, Allée du Six Août 3, 4000 Liège, Belgium  Purushottam
S. Dubey  Forschungszentrum Jülich GmbH, Jülich Center for Neutron Science at
the Heinz Maier Leibnitz Zentrum, Lichtenbergstrasse 1, 85747 Garching,
Germany  Andreas M. Stadler  Forschungszentrum Jülich GmbH, Jülich Center for
Neutron Science, 52425 Jülich, Germany  Institute of Physical Chemistry, RWTH
Aachen University, Landoltweg 2, 52056 Aachen, Germany  Baohu Wu
Forschungszentrum Jülich GmbH, Jülich Center for Neutron Science at the Heinz
Maier Leibnitz Zentrum, Lichtenbergstrasse 1, 85747 Garching, Germany  Orsolya
Czakkel  Institut Laue-Langevin, 71 avenue des Martyrs, 38042 Grenoble Cedex
9, France  Lionel Porcar  Institut Laue-Langevin, 71 avenue des Martyrs, 38042
Grenoble Cedex 9, France  Sebastian Jaksch  Forschungszentrum Jülich GmbH,
Jülich Center for Neutron Science at the Heinz Maier Leibnitz Zentrum,
Lichtenbergstrasse 1, 85747 Garching, Germany  European Spallation Source
(ESS) ERIC, Partikelgatan 2, 224 84 Lund, Sweden  Henrich Frielinghaus
Forschungszentrum Jülich GmbH, Jülich Center for Neutron Science at the Heinz
Maier Leibnitz Zentrum, Lichtenbergstrasse 1, 85747 Garching, Germany  Olaf
Holderer<EMAIL_ADDRESS>Forschungszentrum Jülich GmbH, Jülich
Center for Neutron Science at the Heinz Maier Leibnitz Zentrum,
Lichtenbergstrasse 1, 85747 Garching, Germany

###### Abstract

A mathematical model is developed, to jointly analyze elastic and inelastic
scattering data of fluctuating membranes within a single theoretical
framework. The model builds on a non-homogeneously clipped time-dependent
Gaussian random field. This specific approach provides one with general
analytical expressions for the intermediate scattering function, for any
number of sublayers in the membrane and arbitrary contrasts. The model is
illustrated with the analysis of small-angle x-ray and neutron scattering as
well as with neutron spin-echo data measured on unilamellar vesicles prepared
from phospholipids extracted from porcine brain tissues. The parameters fitted
on the entire dataset are the lengths of the chain and head of the molecules
that make up the membrane, the amplitude and lateral sizes of the bending
deformations, the thickness fluctuation, and a single parameter characterizing
the dynamics.

##  I Introduction

Membranes are ubiquitous in living systems, where they play a central role in
controlling the interactions and exchanges of cells and organelles with their
environment.Watson (2015) They are also common in synthetic systems, such as
vesicles, lipo- or polymer-somes used as drug carriers or reaction
compartments.Discher and Eisenberg (2002); Guinart  _et al._ (2023) Developing
analytical tools to investigate the nanometer-scale structure of membranes and
their dynamics, is key to understanding their formation mechanisms, their
physicochemical properties, and in the case of biological or neuronal
membranes how they fulfill their function.Krugmann  _et al._ (2020, 2021) In
that general context, scattering methods (of either x-rays or neutrons) play a
unique role because they enable one to characterize structures in-situ in
their natural environment with nanometer resolution.Büldt  _et al._ (1978);
Pabst  _et al._ (2010)

Scattering methods provide invaluable yet indirect structural and dynamical
information in the form of correlation functions.Van Hove (1954); Sivia
(2011); Squires (2012); Glatter (2018); Luo  _et al._ (2023) A central aspect
of any scattering investigation is, therefore, the development of suitable
data analysis methods, by which the reciprocal-space data is converted to
structurally significant information in real space. In that process, models
are almost unavoidable.Pedersen (1997); Svergun (1999); Gommes (2018)

The simplest mathematical models of membranes describe them as a succession of
flat layers with specific scattering-length densities. In that spirit, the
main scattering characteristics of, say, a phospholipid membrane can be
captured assuming a central hydrophobic layer squeezed between two hydrophilic
layers. The description can be refined to account for more complex scattering-
length density profiles, in relation e.g. to uneven water dissolution in the
various segments,Kiselev  _et al._ (2008) or to the curvature of the
membrane.Kiselev  _et al._ (2002); Pencer  _et al._ (2006); Chappa  _et al._
(2021) These types of models, however, do not capture the deformation of the
actual membrane shape away from the ideal flat-layer structure.Monzel and
Sengupta (2016) Such deformations contribute to a background scattering,Hamley
(2022) and capturing them in a model would provide physical insights into the
propensity of the membrane to bending and/or compression. Both deformation
modes are central for a variety of physiological processes, including
endo/exocytosis McMahon and Gallop (2005) and cellular adhesion.Zidovska and
Sackmann (2006)

More realistic membrane models can be produced through molecular dynamics
simulations.Venable  _et al._ (1993); Sodt  _et al._ (2014); Chakraborty  _et
al._ (2020) This offers the possibility of building scattering data analyses
on real physical interactions but the focus of this type of simulations is
generally on molecular-scale processes, while the deformations of membranes
are collective modes that involve thousands of molecules. Because of
computational constraints, simulations are often restricted to small volumes,
which is expected to overestimate the stiffness - i.e. underestimate the
amplitude - of large-scale deformation modes. The accessible timescales are
also in the picosecond range, which makes these models suitable for analyzing
the diffusion of individual molecules. However, molecular-dynamic simulations
of nanometer-scale deformations of membranes occurring over nanoseconds remain
a challenge.Hayward  _et al._ (2022)

Interestingly, the models classically used to analyze the dynamics of membrane
deformation through inelastic neutron scattering are distinctly different from
those used to analyze the structure of the membrane through elastic
scattering. When analyzing the data of neutron spin-echo experiments,Luo  _et
al._ (2023) a central role is played by the theoretical results of Zilman and
Granek concerning the bending fluctuations of two-dimensional films,Zilman and
Granek (1996, 2002) whereby the time-dependence of the intermediate scattering
function is described as stretched exponentials. These developments have been
since generalized and adapted to other contexts,Nagao and Seto (2023); Granek
_et al._ (2024) also to characterize thickness fluctuations.Nagao (2009);
Woodka  _et al._ (2012) With such models, however, the inelastic and elastic
scattering are analyzed independently of one another, and this prevents one
from building on the known structure of the membrane to understand its
dynamics.

Here, we have developed a mathematical model of a fluctuating membrane, to
jointly analyze elastic and inelastic scattering data within a single
theoretical framework, with arbitrary sets of contrasts between the various
phases that make up the membrane. The fluctuating parts of the model are
captured through the statistics of a time-dependent Gaussian random field.
Models built on clipped Gaussian fields have been used for many years to
analyze scattering data from emulsions,Berk (1987, 1991); Teubner (1991); Chen
_et al._ (1996) gelsRoberts (1997); Gommes and Roberts (2008) and porous
materials.Levitz (1998); Prehal  _et al._ (2017); Gommes (2018) The present
membrane modelling builds on three generalizations of the classical modelling
approach. First, it uses a non-homogeneous clipping procedure to generate a
structure with a more controlled morphology.Gommes and Pirard (2009); Gommes
_et al._ (2020) It also allows for the presence of multiple regions with
different scattering length densities.Gommes (2013) Finally, it allows the
Gaussian field to be time-dependent in order to make the model suitable for
inelastic scattering data analysis.Gommes  _et al._ (2021); Gommes (2022)

The structure of the paper is the following. The first section presents the
experimental scattering data that motivate the theoretical development of the
paper. For the purpose of comparison, a state-of-the-art analysis of the
elastic and inelastic scattering data is presented. The Gaussian membrane
model is presented in the following section, and a general expression is
derived for the intermediate scattering function. The most technical parts of
the derivations are presented in the Supporting Information. Finally, the
discussion section is focused on the application of the model to analyze the
scattering data presented in the experimental section.

##  II Experimental Section

Figure 1:  Neutron Spin Echo (NSE) data of unilamellar vesicle (a), together
with the corresponding Small-angle neutron (b) and X-ray (c) scattering
patterns. The surface in (a) is a fit of the Zilman-Granek model, and the red
lines in (b) and (c) result from the joint fitting of the SANS and SAXS data
with a polydispersed three-layer slab model.

The unilamellar vesicles (ULV) are prepared from phospholipids, extracted from
porcine brain tissues as described by Folch et. al.Folch  _et al._ (1951)
Figure 1 displays small-angle neutron (SANS) and X-ray(SAXS) scattering
patterns measured on in-house SAXS instrument KWS-X (a Xeuss 3.0 XL from
Xenocs), at the Forschungszentrum JÃŒlich, Garching and D22 SANS
diffractometer at the Institut Laue-Langevin (ILL),Dubey  _et al._ (2023)
respectively. The SANS measurements were performed at wavelength 11.5 Å and 6
Å and detector distances 17.6 m and 1.4 m, resulting in a $q$ range from
0.0026 to 0.64 Å-1.

To provide a reference for the methodological developments of the paper, the
SAXS and SANS data in Fig. 1 were first analyzed through a state-of-the-art
method. A classical approach for analyzing small-angle scattering patterns of
unilamellar vesicles, consists in assuming a spherical structure with radius
much larger than the membrane thickness. In that case, the scattering can be
approximated by the following separated form factorKiselev  _et al._ (2002);
Pencer  _et al._ (2006)

$I(q)\simeq\left|4\pi
R^{2}\frac{\sin(qR)}{qR}\int_{-\infty}^{+\infty}b(z)\cos(qz)\textrm{d}z\right|^{2}$
(1)

where $b(z)$ is the scattering-length density profile along a coordinate $z$
orthogonal to the membrane. This expression can be further simplified in the
limit where $qR\gg 1$, i.e. when the $q$ range of interest concerns the inner
structure of the membrane rather than the vesicle as a whole. Assuming any
statistical distribution of the vesicle radius $R$, the scattering by the
membrane becomes

$I(q)\simeq\frac{2\pi
A}{q^{2}}\left|\int_{-\infty}^{+\infty}b(z)\cos(qz)\textrm{d}z\right|^{2}$ (2)

where $A=4\pi\langle R^{2}\rangle$ is the average area of the vesicle.

The SANS and SAXS data in Fig 1 were fitted through Eq. (2) to a three-layer
scattering-length density profile, corresponding to two phospholipid molecules
facing each other tail-to-tail. The neutron and x-ray scattering-length
densities of the head and chain segments assumed for the fitting are reported
in Tab. 1. In practice, the logarithm of the SANS and SAXS intensities were
fitted by least-square both independently and jointly, with the lengths of the
chain and of the head $l_{C}$ and $l_{H}$ as fitting parameters. As shown in
Tab. 1, fitting the SANS and SAXS data independently from one another leads to
incompatible values of the parameters, in particular for the head size
$l_{H}$. Reasonable values of the parameters are obtained by minimizing
jointly the error on both the SANS and SAXS. Polydispersity had to be assumed
in the fit, to dampen the oscillations of the form factor that are not present
in the data. This is quantified here as the standard deviation of the chain
length relative to the mean; the value is $\sigma_{C}\simeq 15$ %.

Table 1: Scattering-length densities relevant to neutron and x-ray scattering,
and fitted parameters of the slab model in the case of SANS, SAXS and joint
SANS/SAXS fitting. The errors are standard deviations of the fitted parameters
observed through Monte Carlo simulations.

| $b_{H}$ | $b_{C}$ | $b_{W}$ | $l_{H}$ | $l_{C}$ | $\sigma_{C}$
---|---|---|---|---|---|---
| ($10^{-6}$ Å-2) | ($10^{-6}$ Å-2) | ( $10^{-6}$ Å-2) | (Å) | (Å) | (%)
SANS | 1.87 | -0.07 | 6.37 | 0.1 $\pm$ 0.0 | 19.3 $\pm$ 0.1 | 11.1 $\pm$ 1.1
SAXS | 14.2 | 8.30 | 9.37 | 5.2 $\pm$ 0.0 | 18.7 $\pm$ 0.0 | 14.7 $\pm$ 0.6
SANS & SAXS | id. | id. | id. | 4.5 $\pm$ 0.4 | 17.1 $\pm$ 0.1 | 15.3 $\pm$ 0.3

$b_{H/C/W}$: scattering length densities of the head and chain segments of the
phospholipid, and of the outer water medium; $l_{H/C}$: length of the head and
chain segments; $\sigma_{C}$: standard deviation of the chain length, relative
to the average.

Neutron spin-echo (NSE) spectroscopy is well suited to measure slow
fluctuations of biological membranes such as phospholipid double layers.Nagao
_et al._ (2017); Nagao and Seto (2023) The intermediate scattering function
$I(q,\tau)$, i.e. the Fourier transform of the real-space Van-Hove correlation
function,Squires (2012) is measured on time scales up to several 100 ns on
mesoscopic length scales (about 10-500 Å). Neutron spin-echo measurements were
carried out at the IN15 spectrometer at ILL,Dubey  _et al._ (2023) with three
different wavelengths 12, 10 and 8 Å, with probing time scale up to 335 ns and
$q$ range from 0.046 to 0.122 Å-1.

The NSE data are shown in Fig. 1a, under the form of the intermediate
scattering function $I(q,\tau)$, normalized by the SANS intensity expressed in
this context as $I(q,0)$. To compare the developments of the paper with a
classical approach, the data were fitted with the Zilman-Granek functionZilman
and Granek (1996); Nagao and Seto (2023), namely

$\frac{I(q,\tau)}{I(q,0)}=\exp\left[-Bq^{2}\tau^{2/3}\right]$ (3)

with parameter $B$ as the only fitting parameter for the entire dataset. The
latter parameter can be given a physical interpretation in terms of the
bending rigidity of the membrane but these values are notoriously dependent on
the assumed viscosity of the medium.Nagao and Seto (2023) Assuming the
viscosity of water ($\eta\simeq$ 10-3 Pa.s), the value of $B$ extracted from
the NSE data in Fig. 1a converts to an effective bending modulus of 175 kBT.

##  III The Gaussian membrane model

###  III.1 Definition of the model

A classical approach to model disordered multiphasic structures consists in
defining the various regions that make up a material based on whether the
values taken locally by a given Gaussian random field is larger or smaller
than a predefined threshold. Quiblier (1984); Berk (1987, 1991); Teubner
(1991); Lantuéjoul (2002) In the context of scattering studies, this approach
has notably been used to analyze elastic scattering from porous
materials,Levitz (1998) polymer blends,Chen  _et al._ (1996) gels,Roberts
(1997); Gommes and Roberts (2008) confined liquids,Gommes (2013); Gommes and
Roberts (2018) etc. Recently, this approach has been generalized to analyze
inelastic scattering data as well through the use of time-dependent Gaussian
fields.Gommes  _et al._ (2021)

The approach we propose here to model the structure and dynamics of
fluctuating membranes is based on non-homogeneously clipped Gaussian
fields,Gommes and Pirard (2009); Gommes  _et al._ (2020) as sketched in Fig.
2. In that spirit, the space- and time-dependent structure of the various
regions that make up the membrane (head, chain, as well as outer solution) is
modelled based on the combination of a stochastic Gaussian field
$W(\mathbf{x},t)$ and of a set of deterministic rules used to convert the
continuous values of the field into geometrical sets that describe the
regions. In Fig. 2, two Gaussian fields are shown in grey. The conversion
rules are sketched in the colored flag, which can be thought of as a lookup
table: whether one point $\mathbf{x}$ is assigned to any particular region at
time $t$ depends on the position (here the distance $z$ to a given plane which
controls the average orientation of the membrane) and on the local value of
the random field $W(\mathbf{x},t)$.

Figure 2:  Sketch of the membrane model, with (a) the model’s flag (heads in
red, chains in orange, and outer solution in white) and two Gaussian random
fields having different correlation lengths $l_{Z}$ in direction $z$ (b1 and
c1). In both cases, the correlation length in the orthogonal direction is
$l_{XY}=30$ Å. The membrane structures resulting from these specific fields
and flag are shown in b2 and c2.

A Gaussian field is comprehensively characterized by its space- and time-
dependent correlation function $g_{W}(\mathbf{r},\tau)$, which makes it a
particularly useful concept in the context of scattering studies. The field
correlation function is defined as

$g_{W}(\mathbf{r},\tau)=\langle
W(\mathbf{x},t)W(\mathbf{x}+\mathbf{r},t+\tau)\rangle$ (4)

where the brackets $\langle\rangle$ refer throughout the text to ensemble
averages. In other words, $g_{W}(r,\tau)$ describes the statistical
correlations between the values of $W(\mathbf{x},t)$ at two points at distance
$\mathbf{r}$ from each other, with a time lag $\tau$.

A variety of procedures to construct Gaussian fields with specific space- and
time-correlation functions are discussed for example by Gommes et al.Gommes
_et al._ (2021) The methods developed in the present paper are quite general
and they apply to any of them. For the purpose of illustration, the two fields
shown in the top row of Fig. 2 correspond to a correlation function of the
type

$g_{W}(\mathbf{r},0)=\exp\left[-\frac{r_{x}^{2}+r_{y}^{2}}{l_{XY}^{2}}-\frac{r_{z}^{2}}{l_{Z}^{2}}\right]$
(5)

where $r_{x}$, $r_{y}$, $r_{z}$ are the components of $\mathbf{r}$ parallel or
perpendicular to the membrane; $l_{XY}$ and $l_{Z}$ are correlation lengths in
those specific directions, which are two parameters that comprehensively
describe the Gaussian field. In the cases of Figs. 2b1 and 2c1 the same
correlation length is used in the $XY$ plane ($l_{XY}=30$ Å), and the out-of
plane correlation lengths are $l_{Z}=30$ Å and 100 Å, respectively. The
correlation function in Eq. (5) is written for the value $\tau=0$, and it
characterizes therefore the spatial properties of an instantaneous snapshot of
a possibly time-dependent field. To keep the discussion as general as
possible, we postpone the discussion of specific time dependence to Sec.
III.2.3.

The limits between the layer boundaries are straight lines in the model’s flag
(Fig. 2a). To keep the discussion general, we allow an arbitrary number of
layers $N$ in the model, and we introduce $N+1$ linear threshold functions

$\alpha_{n}(z)=(z-Z_{n})/l_{\alpha}$ (6)

with $n=0,\ldots N$, with parameters $Z_{0}<Z_{1}<\ldots<Z_{N}$ and
$l_{\alpha}$. Any point of space $\mathbf{x}$ is assigned to the $n^{th}$
layer at time $t$ if the value of the Gaussian field satisfies

$\alpha_{n}(z)\leq W(\mathbf{x},t)<\alpha_{n-1}(z)$ (7)

In the case of Fig. 2, with $N=3$, the values are $Z_{0}=-24.5$ Å,
$Z_{1}=-16.5$ Å, $Z_{2}=+16.5$ Å, $Z_{3}=+24.5$ Å, and $l_{\alpha}=10$Å. The
values of $Z_{n}$ control the most probable positions of the interfaces and
the parameter $l_{\alpha}$ controls the slope of the boundary limits in the
flag space, which converts to the amplitude of the interface fluctuations in
real space. Flat interfaces centred exactly at $z=Z_{n}$ are obtained in the
limiting case $l_{\alpha}\to 0$.

In addition to the clipping parameters $l_{\alpha}$ and $Z_{n}$’s, the
structure of the layers is also controlled by the field correlation function.
In the particular case of Fig. 2 and Eq. (5) the correlation length $l_{XY}$
controls the size of the bending-like patterns in the direction parallel to
the membrane. The correlation length in the orthogonal direction $l_{Z}$
controls how the surface fluctuations at the different interfaces are
correlated with each other. In the limit where $l_{Z}\gg|Z_{i}-Z_{j}|$ the
values of the Gaussian field at $z=Z_{i}$ and $z=Z_{j}$ are strongly
correlated ($g_{W}\simeq 1$) so that the fluctuations of the interfaces $i$
and $j$ are almost identical (Fig. 2 right). By contrast, if
$l_{Z}\ll|Z_{i}-Z_{j}|$ the fluctuations of interfaces $i$ and $j$ are
statistically independent of one another (Fig. 2 left). The value of $l_{Z}$
therefore controls the membrane thickness fluctuations.

An important geometrical characteristic of the membrane is the specific
surface area of the individual interfaces. This is conveniently expressed as a
dimensionless roughness factor $a_{A}=a/A$, defined as the ratio of the area
of the distorted interface $a$ to the area of the project flat surface $A$. A
general expression for the roughness factor as a function of dimensionless
numbers $l_{Z}/l_{XY}$ and $l_{\alpha}/l_{XY}$ is calculated in Sec. SI-I of
the Supporting Information. The values are plotted in Fig. 3, together with
some realisations representative of various sets of parameters. In the
particular case of an isotropic field, i.e. for $l_{Z}=l_{XY}$, the roughness
can be calculated through a general expression derived in earlier work (see
Eq. (6) of Ref. Gommes and Pirard, 2009) and expressed as follows

$a_{A}=\frac{2}{\sqrt{\pi}}\frac{l_{\alpha}}{l_{XY}}\left\\{\exp\left[-\left(\frac{l_{XY}}{2l_{\alpha}}\right)^{2}\right]+\sqrt{\pi}\left(\frac{l_{XY}}{2l_{\alpha}}+\frac{l_{\alpha}}{l_{XY}}\right)\textrm{erf}\left[\frac{l_{XY}}{2l_{\alpha}}\right]\right\\}$
(8)

This expression is plotted in Fig. 3 as a solid red line. Globally, the
dependence of $a_{A}$ on the correlation length $l_{Z}$ is strong only for
$l_{Z}<l_{XY}$. This region corresponds to highly distorted and even
disconnected interfaces (e.g. Fig. 3b), which are unrealistic in the context
of membranes. For parameters relevant to membranes, Eq. (8) provides a fair
approximation of the roughness factor.

Figure 3:  Left panel: Dimensionless roughness factor $a_{A}$ of the
interfaces in the membrane model, as a function of the field and clipping
characteristic lengths. The solid red line at $l_{Z}/l_{XY}=1$ is the
analytical result in Eq. (8), which is also plotted at larger values of
$l_{Z}/l_{XY}$ to highlight the weak dependence on $l_{Z}$. Panels a to d
point at specific realizations, also highlighted in the left panel. The red
dot is the value inferred in Sec. IV from the scattering data.

In addition to strictly geometrical aspects, the material characteristic
relevant to scattering is the space- and time-dependent scattering-length
density map $b(\mathbf{x},t)$, which quantifies the total scattering length in
an infinitesimal volume centred on point $\mathbf{x}$ at time $t$. Sivia
(2011); Squires (2012) If we assume that each layer $n$ in the membrane model
has a specific composition-dependent scattering-length density $b_{n}$, the
scattering length density map can be expressed as the following piecewise-
constant functionGommes  _et al._ (2021)

$b(\mathbf{x},t)=\sum_{n=0}^{N+1}b_{n}\mathcal{I}_{n}(\mathbf{x},t)$ (9)

where $\mathcal{I}_{n}(\mathbf{x},t)$ is the indicator function of layer $n$,
which is equal to 1 if point $\mathbf{x}$ is in the $n^{th}$ layer at time $t$
and to 0 otherwise. Values of $b_{n}$ relevant to the head and chain segments
of the phospholipid membrane in Fig. 1, and the outer water solution are given
in Tab. 1.

In the context of the non-homogeneously clipped Gaussian field of Eq. (7), the
indicator function of the $n^{th}$ layer is formally defined as

$\mathcal{I}_{n}(\mathbf{x},t)=H\left[W(\mathbf{x},t)-\alpha_{n}(z)\right]-H\left[W(\mathbf{x},t)-\alpha_{n-1}(z)\right]$
(10)

where $H[]$ is Heaviside’s step function (equal to 1 if the argument is
positive and to 0 otherwise). The indicator functions of the outer regions
($n=0$ and $n=N+1$) are

$\mathcal{I}_{0}(\mathbf{x},t)=H\left[W(\mathbf{x},t)-\alpha_{0}(z)\right]$
(11)

and

$\mathcal{I}_{N+1}(\mathbf{x},t)=1-H\left[W(\mathbf{x},t)-\alpha_{N}(z)\right]$
(12)

corresponding to the water solution on both sides of the membrane.

Note finally that although all illustrations (e.g. in Fig. 2) assume a three-
region model - with head, chain and water - of a single membrane, all the
equations of the paper are general and the number of layers can be arbitrarily
increased. This makes the present approach suitable also for refined models
wherein the head and chain regions are decomposed into subregions with
specific scattering-length densities reflecting possibly different water
concentrations, Kiselev  _et al._ (2008); Lewis-Laurent  _et al._ (2021) or
for multi-layer structures. Lemmich  _et al._ (1996); Krugmann  _et al._
(2020); Maiti  _et al._ (2021)

###  III.2 Scattering properties

In the general context of stochastic models, the intermediate scattering
function at wavevector $\mathbf{q}$ and time $\tau$ is calculated as the
following Fourier transformGommes  _et al._ (2021)

$I(\mathbf{q},\tau)=\int dV_{1}\int dV_{2}\
e^{-i\mathbf{q}\cdot(\mathbf{x}_{1}-\mathbf{x}_{2})}\langle
b(\mathbf{x}_{1},t)b(\mathbf{x}_{2},t+\tau)\rangle$ (13)

where the integrals are on $\mathbf{x}_{1}$ and $\mathbf{x}_{2}$, and the
brackets $\langle\rangle$ stand for an ensemble average. In the particular
case of a structure with a piecewise-constant scattering-length density like
in Eq. (9), this can be written as

$I(\mathbf{q},\tau)=\sum_{n=0}^{N+1}\sum_{m=0}^{N+1}b_{n}b_{m}P_{m,n}(\mathbf{q},\tau)$
(14)

with

$P_{m,n}(\mathbf{q},\tau)=\int\textrm{d}V_{1}\int\textrm{d}V_{2}\quad
e^{-i\mathbf{q}\cdot\left(\mathbf{x}_{1}-\mathbf{x}_{2}\right)}S^{(2)}_{m,n}(\mathbf{x_{1}},\mathbf{x_{2}},\tau)$
(15)

In this expression $S^{(2)}_{m,n}(\mathbf{x_{1}},\mathbf{x_{2}},\tau)$ is a
two-point correlation function, equal to the probability for point
$\mathbf{x}_{1}$ to belong to layer $m$ at time $t$, and for $\mathbf{x}_{2}$
to belong to layer $n$ at later time $t+\tau$. In terms of the indicator
functions, this is defined formally as

$S^{(2)}_{m,n}(\mathbf{x_{1}},\mathbf{x_{2}},\tau)=\langle\mathcal{I}_{m}(\mathbf{x}_{1},t)\mathcal{I}_{n}(\mathbf{x}_{2},t+\tau)\rangle$
(16)

As Gaussian random fields are second-order stationary,Lantuejoul (1991);
Lantuéjoul (2002) the time-dependence of the two-point functions in Eq. (15)
is only through the difference $\tau=t_{2}-t_{1}$. Such simplification does
not apply to the $\mathbf{x}_{1}$ and $\mathbf{x}_{2}$ dependence because the
clipping procedure in Eq. (10) depends explicitly on the position $z$, so that
the spatial stationarity of the Gaussian field is broken by the space-
dependent clipping procedureGommes and Pirard (2009).

In the limit where the time lag $\tau=t_{2}-t_{1}$ is much larger than any
correlation time of the Gaussian field, the fluctuations at
$(\mathbf{x}_{1},t_{1})$ and $(\mathbf{x}_{2},t_{2})$ are statistically
independent. In that limit, the two-point correlation function simplifies to
the following product

$\lim_{\tau\to\infty}S^{(2)}_{m,n}(\mathbf{x}_{1},\mathbf{x}_{2},\tau)=S_{m}^{(1)}(\mathbf{x}_{1})S_{n}^{(1)}(\mathbf{x}_{2})$
(17)

where the one-point function

$S_{n}^{(1)}(\mathbf{x})=\langle\mathcal{I}^{(n)}(\mathbf{x},t)\rangle$ (18)

is the probability for point $\mathbf{x}$ to belong to layer $n$, at any given
time. The limit in Eq. (17) suggests decomposing the two-point correlation
functions as follows

$S_{m,n}^{(2)}(\mathbf{x}_{1},\mathbf{x}_{1},\tau)=S^{(1)}_{m}(\mathbf{x}_{1})S^{(1)}_{m}(\mathbf{x}_{2})+\tilde{S}_{m,n}^{(2)}(\mathbf{x}_{1},\mathbf{x}_{1},\tau)$
(19)

where the first term on the right-hand side is the time-independent
contribution from the average structure, and $\tilde{S}_{m,n}^{(2)}$ is the
contribution from the fluctuations.

The splitting of the two-point correlation function into average and
fluctuating contributions, results in the following decomposition of the
intermediate scattering function

$I(\mathbf{q},\tau)=\bar{I}(\mathbf{q})+\tilde{I}(\mathbf{q},\tau)$ (20)

where the two contributions $\bar{I}(\mathbf{q})$ and
$\tilde{I}(\mathbf{q},\tau)$ are calculated from Eqs. (14) and (15) applied to
the first and second terms on the right-hand side of Eq. (19). We explicit
these two contributions in the following two sections.

####  III.2.1 Scattering contribution of the average structure

The contribution of the average structure to the intermediate scattering
function is calculated by applying Eqs. (14) and (15) to the first term in Eq.
(19). This results in

$\bar{I}(\mathbf{q})=\left|\int e^{-i\mathbf{q}\cdot\mathbf{x}}\
\bar{b}(\mathbf{x})\ \textrm{d}V_{x}\right|^{2}$ (21)

where

$\bar{b}(\mathbf{x})=\sum_{n=0}^{N+1}b_{n}S_{n}^{(1)}(\mathbf{x})$ (22)

is the average scattering-density map. Equation (21) is classical, and it
reduces to Eq. (2) in the case of a large vesicle. It has to be stressed,
however, that in the present context, $\bar{b}(z)$ is not the true
instantaneous scattering-length density profile, but its average value
smoothed over time.

In the case of the membrane model in Eq. (10), the one-point probability
function of layer $n$ is calculated as

$S^{(1)}_{n}(\mathbf{x})=\Lambda_{1}[\alpha_{n}(z)]-\Lambda_{1}[\alpha_{n-1}(z)]$
(23)

where $\Lambda_{1}[\alpha]$ is the probability for a centred Gaussian variable
with unit variance to take values larger than $\alpha$. This can also be
written as

$\Lambda_{1}[\alpha]=\frac{1}{2}\left(1-\textrm{erf}\left[\frac{\alpha}{\sqrt{2}}\right]\right)$
(24)

where erf[] is the error function. The one-point functions of the two outer
regions are

$\displaystyle S^{(1)}_{0}(\mathbf{x},t)$ $\displaystyle=$
$\displaystyle\Lambda_{1}[\alpha_{0}(z)]$ (25) $\displaystyle
S^{(1)}_{N+1}(\mathbf{x},t)$ $\displaystyle=$ $\displaystyle
1-\Lambda_{1}[\alpha_{N}(z)]$ (26)

The one-point probability functions corresponding to Fig. 2 are illustrated in
Fig. 4, together with the average scattering-length density profiles relevant
to x-ray and neutron scattering.

Figure 4:  One-point probabilities $S_{n}^{(1)}(z)$ of the various layers of
the membrane model in Fig. 2 (a1 to a5), and (b) corresponding average
scattering-length density profile relevant to x-ray (red) and neutron (blue)
scattering. The vertical dashed lines are the most probable positions of the
water/head and head/chain interfaces.

The average scattering-length density profiles do not depend on any
characteristic of the Gaussian field, so that the profiles are identical for
the compressible and incompressible cases shown in Figs. 2b2 and 2c2. The
corresponding scattering contribution $\bar{I}(q)$ relevant to x-ray and
neutron are plotted as black lines in Fig. 6c and Fig. 7c, respectively.
Because $\bar{I}(q)$ is the Fourier transform of a smooth function, with no
singularity in any of its derivatives, it decreases asymptotically with $q$
faster than any power law.Lighthill (1958) In particular, with the splitting
of the scattering in Eq. (20) it is the fluctuation contribution that is
responsible for Porod’s $q^{-4}$ scattering.

####  III.2.2 Scattering by the fluctuations

Figure 5: Function $G[\eta,g]$ defined in Eq. (30)

The contribution of the fluctuations to the scattering
$\tilde{I}(\mathbf{q},\tau)$ is specific to the present stochastic model, and
it has no equivalent in the deterministic slab model of Eq. (2). It accounts
for the time-dependent bending-like deformation of the membrane as well as for
its thickness fluctuations. Mathematically, it is calculated by applying Eqs.
(14) and (15) to the second term in Eq. (19). We show in Sec. SI-II of the
Supporting Information that its value per unit area $A$ of the membrane is
calculated as the Fourier transform of the following correlation function

$\displaystyle\tilde{C}_{b}(\mathbf{r},\tau)$ $\displaystyle=$
$\displaystyle\left[\sum_{n=0}^{N}(b_{n}-b_{n+1})^{2}\right]\Gamma_{0}(\mathbf{r},\tau)$
(27) $\displaystyle+$
$\displaystyle\sum_{m>n}(b_{n}-b_{n+1})(b_{m}-b_{m+1})\left[\Gamma_{Z_{m}-Z_{n}}(\mathbf{r},\tau)+\Gamma_{Z_{n}-Z_{m}}(\mathbf{r},\tau)\right]$
(28)

where the function $\Gamma_{L}(\mathbf{r},\tau)$ is given by

$\Gamma_{L}(\mathbf{r},\tau)=l_{\alpha}G\left[\frac{|r_{z}-L|}{l_{\alpha}},g_{W}\left(\mathbf{r},\tau\right)\right]$
(29)

with

$\displaystyle G[\eta,g]$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{\pi}}\left(\exp\left[-\frac{\eta^{2}}{4}\right]-\exp\left[-\frac{\eta^{2}}{4(1-g)}\right]\sqrt{1-g}\right)$
(30) $\displaystyle+$
$\displaystyle\frac{\eta}{2}\left(\textrm{erf}\left[\frac{\eta}{2}\right]-\textrm{erf}\left[\frac{\eta}{2\sqrt{1-g}}\right]\right)$
(31)

The function $G[\eta,g]$ is plotted in Fig. 5 against its two arguments $\eta$
and $g$. Note that in the particular case where $l_{\alpha}$ becomes
vanishingly small, the stochastic model reduces to a deterministic slab model.
In that case, all $\Gamma$’s in Eq. (27) converge to zero, and the only
contribution left to the scattering is from the average structure
$\bar{I}(q)$, as it should.

Figures 6a1 and 6a2 display the correlation functions
$\tilde{C}_{b}(\mathbf{r},0)$ relevant to the two membrane models in Fig. 2,
as a function of the in-plane and out-of-plane components of $\mathbf{r}$,
namely $r_{xy}$ and $r_{z}$. The x-ray contrast is assumed for the figure. The
$r_{xy}$ dependence of the correlation function, corresponding to in-plane
vectors $\mathbf{r}$, are very similar for the two models. It characterises
the wavelike bending of the membrane with a typical size $l_{XY}=30$ Å in the
case of the figure, independently of the out-of-plane correlation length
$l_{Z}$. The $r_{z}$ dependences, however, are distinctly different in Figs.
6a1 and 6a2. Strong oscillations appear in $\tilde{C}(\mathbf{r},0)$ for
$l_{Z}=100$ Å, because in that case the membrane bends while keeping a
constant thickness (Fig. 2 right). By contrast, the membrane thickness
fluctuates in the case with $l_{Z}=30$ Å (Fig. 2 left) which leads to blunter
features in the correlation function.

Figure 6:  Fluctuation contribution to the correlation function
$\tilde{C}_{b}(\mathbf{r},0)$ (top) and scattered intensity
$\tilde{I}(\mathbf{q},0)$ (middle), together with the rotationally-averaged
scattered intensity (bottom), in the case of x-ray contrast. The parameters
are those relevant to the membrane models in Fig. 2, with $l_{XY}=30$ Å and
$l_{Z}=30$ Å (a1, b1) or $l_{Z}=100$ Å (a2, b2). In c, the black line is the
scattering from the average structure. The colored dashed and solid lines are
the fluctuation contributions and the total scattering, for $l_{Z}=30$ Å (red)
and $l_{Z}=100$ Å (blue).

The 3D Fourier transforms $\tilde{I}(\mathbf{q},0)$ exhibit the same
differences in reciprocal space. They are shown in Figs. 6b1 and 6b2 against
in-plane and out-of-plane components of the scattering vector $q_{xy}$ and
$q_{z}$. Because of the cylindrical symmetry of the system, they were
calculated as

$\tilde{I}(q_{xy},q_{z},\tau)=4\pi\int_{0}^{\infty}\textrm{d}r_{z}\int_{0}^{\infty}r_{xy}\textrm{d}r_{xy}\
J_{0}(q_{xy}r_{xy})\cos(q_{z}r_{z})\tilde{C}_{b}(r_{xy},r_{z},\tau)$ (32)

where $J_{0}()$ is the Bessel function of order 0. The rotationally-averaged
intensity, relevant to SANS and SAXS analysis, is then obtained as

$\tilde{I}(q,\tau)=\int_{0}^{\pi/2}\sin(\theta)\textrm{d}\theta\
\tilde{I}(q\sin(\theta),q\cos(\theta),\tau)$ (33)

The so-obtained value of $\tilde{I}(q,0)$ are shown in Fig. 6c.

The membrane thickness fluctuation is found to have two distinctly different
effects on the scattering patterns. First and foremost, it enhances the
scattering at low $q$. This is visible in the different values of the low-$q$
plateaus in the fluctuation contributions $\tilde{I}(q,0)$ in Fig. 6, which is
about 100 times smaller for the larger value of $l_{Z}$. This effect is
physically identical to the well-known enhancement of forward scattering by
compressibilityGlatter (2018). As a second, weaker, effect: thickness
fluctuations add polydispersity to the distance between the hydrophilic heads
on both sides of the membrane, and this dampens the oscillations at
intermediate $q$. By contrast, the high-$q$ scattering is almost blind to
thickness fluctuations, because it is dominated by surface-dependent $q^{-4}$
scattering and the areas are only weakly dependent on $l_{Z}$ (see Fig. 3).

Figure 7:  Fluctuation contribution to the correlation function
$\tilde{C}_{b}(\mathbf{r},0)$ (top) and scattered intensity
$\tilde{l}(\mathbf{q},0)$ (middle), together with the rotationally-averaged
scattered intensity (bottom), in the case of neutron contrast. The parameters
are those relevant to the membrane models in Fig. 2, with $l_{XY}=30$ Å and
$l_{Z}=30$ Å (a1, b1) or $l_{Z}=100$ Å (a2, b2). In c, the black line is the
scattering from the average structure. The colored dashed and solid lines are
the fluctuation contributions and the total scattering, for $l_{Z}=30$ Å (red)
and $l_{Z}=100$ Å (blue).

Despite the strong effect of thickness fluctuations on forward scattering, one
has to consider also that the total scattering at low $q$ is dominated by the
average-structure contribution. This is manifest through Eq. (2), which
reduces to

$\bar{I}(q)\simeq\frac{2\pi
A}{q^{2}}\times\left|\int_{-\infty}^{+\infty}\bar{b}(z)dz\right|^{2}$ (34)

at low $q$, where the second factor is the squared total excess scattering
length of the membrane. The latter $q^{-2}$ scattering hides the
compressibility-related forward scattering, which becomes visible only at
intermediate $q$. This effect is much stronger for neutron-like chain contrast
(Fig. 7) than for x-ray like head contrast (Fig. 6).

####  III.2.3 Inelastic scattering by time-dependent fluctuations

In all mathematical expressions discussed in Sec. III.2.2, the dependence of
the scattering intensity on the Gaussian field is through the field
correlation function $g_{W}$. This is notably the case for the fluctuation
contribution to the scattering in Eqs. (27) and (29), which are valid whether
the Gaussian field is static or time-dependent. In the latter case, the field
correlation function depends on the time lag $\tau$ in line with its
definition in Eq. (4), and so does the intermediate scattering function
$I(q,\tau)$. We consider now four different time-dependent Gaussian random
fields, all of which correspond to the same spatial statistics as the
correlation function $g_{W}(\mathbf{r},0)$ in Eq. (5).

When it comes to adding time dependence to a static Gaussian random field, it
is convenient to think of it as a superposition of a large number of
elementary wave packets $w(\mathbf{x})$ localized in space, namely

$W(\mathbf{x})=\sum_{s}A_{s}w(\mathbf{x}-\mathbf{x}_{s})$ (35)

where the points $\mathbf{x}_{s}$ are uniformly distributed in space,
according to a Poisson point process with density $\theta$, and $A_{s}$ are
random and statistically independent amplitudes. The so-constructed field has
correlation function Lantuéjoul (2002); Gommes  _et al._ (2021)

$g_{W}(\mathbf{r},0)=\theta\langle A^{2}\rangle K(\mathbf{r})$ (36)

where

$K(r)=\int\textrm{d}V_{x}\ w(\mathbf{x})w(\mathbf{x}-\mathbf{r})$ (37)

is the self-convolution of the elementary wave packet.

In the limit of infinitely large density of waves $\theta$, the value of
$W(\mathbf{x})$ is the sum of a large number of independent contributions. It
then results from the central-limit theorem that the local values of the field
are Gaussian distributed. Moreover, as a consequence of Eq. (36), one has to
impose $\theta\langle A^{2}\rangle K(0)=1$ to ensure that the field has
variance equal to one. The specific wave corresponding to the field
correlation function in Eq. (5) is

$w(\mathbf{x})=\exp\left[-2\left(\frac{x^{2}+y^{2}}{l_{XY}^{2}}+\frac{z^{2}}{l_{Z}^{2}}\right)\right]$
(38)

with self-convolution

$K(\mathbf{r})=\frac{\pi^{3/2}l_{XY}^{2}l_{Z}}{8}\exp\left[-\left(\frac{r_{x}^{2}+r_{y}^{2}}{l_{XY}^{2}}+\frac{r_{z}^{2}}{l_{Z}^{2}}\right)\right]$
(39)

A variety of other wave shapes and their corresponding field correlation
functions are given in Tab. 1 of Ref. Gommes  _et al._ , 2021.

In the first type of time-dependent models we consider, the amplitudes $A_{s}$
of the waves are let to oscillate with frequency $\omega_{s}$ and phase
$\varphi_{s}$, namely

$W(\mathbf{x},t)=\sqrt{2}\sum_{s}A_{s}w(\mathbf{x}-\mathbf{x}_{s})\cos(\omega_{s}t-\varphi_{s})$
(40)

where the factor $\sqrt{2}$ keeps the variance of the field equal to one. The
phases are chosen uniformly in $[0,2\pi)$, and the frequencies are distributed
according to any user-specified spectral density $f(\omega)\textrm{d}\omega$.
With such a construction, the space- and time-dependent correlation function
takes the following separable form Gelfand  _et al._ (2010); Gommes  _et al._
(2021)

$g_{W}(\mathbf{r},\tau)=g_{W}(\mathbf{r},0)g^{\prime}(\tau)$ (41)

where the time-dependent factor depends on the spectral density as follows

$g^{\prime}(\tau)=\int_{0}^{\infty}\cos[\omega\tau]f^{\prime}(\omega)\textrm{d}\omega$
(42)

This cosine transform can in principle be inverted to calculate the spectral
density corresponding to any desired correlation function. Not all correlation
functions, however, are realizable as some may lead to non-positive spectral
densities.

For modelling purposes, a natural choice is an exponential correlation
$g^{\prime}(\tau)=\exp(-\tau/\tau_{c})$ with correlation time $\tau_{c}$ as a
parameter. This corresponds to the following spectral density

$f^{\prime}(\omega)=\frac{1}{\pi}\frac{2\tau_{c}}{1+(\omega\tau_{c})^{2}}$
(43)

which is indeed positive and realizable. A realization of this field is shown
in Fig. 8a, together with the time dependent values sampled at three different
points. This specific field is not differentiable with respect to time, which
is a direct consequence of the linearity of $g^{\prime}(\tau)$ for small
$\tau$.Gommes  _et al._ (2021) A smoother dynamics - with a differentiable
field - can be obtained by choosing a time correlation of the type
$g^{\prime}(\tau)=1/\cosh[\tau/\tau_{c}]$. This function is identical to an
exponential for large times, but it is quadratic for small values of $\tau$.
It corresponds to the spectral density

$f^{\prime}(\omega)=\frac{2\tau_{c}\cosh[\pi\omega\tau_{c}/2]}{1+\cosh[\pi\omega\tau_{c}]}$
(44)

A realization of this field is shown in Fig. 8b. Its differentiability is
apparent when looking at the values sampled at specific points in space.

Figure 8: Realizations of time-dependent Gaussian fields with
$l_{XY}=l_{Z}=30$ Å, and various time dependences: (a and b) exponential and
hyperbolic secant correlation with correlation time $\tau_{c}$; (c) diffusive
wave motion with diffusion coefficient $D$; (d) ballistic wave motion with
velocity $c$. In each case, the values of the field $W$ sampled at three
points are plotted against time. Note how the two models on the right are
differentiable with respect to time, and the models on the left are not.

An alternative approach for creating time-dependent Gaussian fields starting
from the wave superposition in Eq. (35), consists in keeping the amplitudes of
the waves $A_{s}$ constant but letting the waves move in space with specific
probabilistic laws, namely

$W(\mathbf{x},t)=\sum_{s}A_{s}w(\mathbf{x}-\mathbf{x}_{s}-\mathbf{j}_{s}(t))$
(45)

where $\mathbf{j}_{s}(t)$ is the displacement jump of wave $s$ at time $t$,
away from its initial position $\mathbf{x}_{s}$. The field correlation
function corresponding to this construction is Gommes  _et al._ (2021)

$g_{W}(\mathbf{r},\tau)=\theta\langle A^{2}\rangle\int\textrm{d}V_{j}\
K(\mathbf{r}-\mathbf{j})f_{\tau}(\mathbf{j})$ (46)

where $f_{\tau}(\mathbf{j})\textrm{d}V_{j}$ is the jump distribution
probability over a duration $\tau$.

In the present context of membrane modelling with the clipping procedure in
Fig. 2, it is only the waves close to $z=0$ that contribute to the structure
fluctuations. Accordingly, we focus on two-dimensional displacements of
elementary waves within the $xy$ plane. We consider two distinctly different
in-plane displacement laws. The first is diffusive, and it corresponds to a
jump probability Berg (1993)

$f_{\tau}(\mathbf{j})=\left(4\pi
D\tau\right)^{-1}\exp\left[-\frac{j_{x}^{2}+j_{y}^{2}}{4D\tau}\right]\delta(j_{z})$
(47)

where $D$ is a diffusion coefficient, which parameterizes the dynamics. In
this equation, the exponential corresponds to a two-dimensional random walk,
the Dirac function $\delta(j_{z})$ restricts the motion to being parallel to
$xy$, and the first factor normalizes the probabilities to one. According to
Eq. (46), this specific displacement law leads to field correlation function

$g_{W}(\mathbf{r},\tau)=\left(1+\frac{4D\tau}{l_{XY}^{2}}\right)^{-1}\exp\left[-\frac{r_{xy}^{2}}{4D\tau+l_{XY}^{2}}\right]\exp\left[-\frac{r_{z}^{2}}{l_{Z}^{2}}\right]$
(48)

where we have used the values of $K(\mathbf{r})$ from Eq. (39). Note that this
time-dependent correlation function reduces to Eq. (5) for $\tau=0$, as it
should. A particular realization of this Gaussian field is shown in Fig. 8c.
This field is also non-differentiable with respect to time because its
correlation function $g_{W}(0,\tau)$ in Eq. (48) is linear in $\tau$ at the
origin.

A smooth time-differentiable field can be generated with moving waves assuming
their ballistic motion with velocity $c$ and random directions parallel to the
$xy$ plane. This corresponds to the following jump probability distribution

$f_{\tau}(\mathbf{j})=\frac{\delta(j_{xy}-c\tau)}{2\pi j_{xy}}\delta(j_{z})$
(49)

where $j_{xy}=\sqrt{j_{x}^{2}+j_{y}^{2}}$ and $j_{z}$ are the in-plane and out
of plane displacements. Here too the delta function restricts the motion to
being parallel to $xy$, and the denominator $2\pi j_{xy}$ normalizes the
probabilities to one. The field correlation function calculated from Eq.(46)
is

$g_{W}(\mathbf{r},\tau)=\exp\left[-\frac{r_{z}^{2}}{l_{Z}^{2}}\right]\exp\left[-\frac{(r_{xy}-c\tau)^{2}}{l_{XY}^{2}}\right]\exp\left[-\frac{2r_{xy}c\tau}{l_{XY}^{2}}\right]I_{0}\left[\frac{2r_{xy}c\tau}{l_{XY}^{2}}\right]$
(50)

where $I_{0}()$ is the modified Bessel function of the first kind. This
expression of $g_{W}(\mathbf{r},\tau)$ also reduces to Eq. (5) for $\tau=0$. A
realization of the ballistic field is shown in Fig. 8d.

Figure 9: Intermediate scattering functions calculated from the parameters of
Fig. 2 ($l_{Z}=30$ Å) and various time dependences: (a and b) exponential and
hyperbolic secant correlation with correlation time $\tau_{c}$; (c) diffusive
wave motion with diffusion coefficient $D$; (d) ballistic wave motion with
velocity $c$. In all cases, the solid red line is the scattering from the
average structure. Note the qualitatively different time dependence of
$I(q,\tau)$ for large $q$ and small $\tau$ in the left and right graphs.

The intermediate scattering functions calculated from Eq. (33) with the same
parameters as in Fig. 2 are plotted in Fig. 9 for the four specific types of
dynamics considered in Fig. 8. In all cases, the intermediate scattering
functions - expressed as $I(q,\tau)/I(q,0)$ \- decrease towards
$\bar{I}(q)/I(q,0)$ for large values of $\tau$, where $\bar{I}(q)$ is the
scattering by the average structure, calculated through Eq. (21).
Mathematically, this results from the fact that the fluctuation contribution
to scattering $\tilde{I}(q,\tau)$ vanishes when the field correlation $g_{W}$
is equal to zero, as is asymptotically the case for $\tau\to\infty$. The two
types of non-differentiable or smooth time-dependences identified in Fig. 8,
lead to qualitatively different shapes of the intermediate scattering function
at high $q$ and low $\tau$. Smooth dynamics lead to quadratic $\tau$
dependence that does not seem to be observed experimentally.Gommes  _et al._
(2021) Other differences between the various dynamics are discussed later,
when using the models to fit the NSE data from Fig. 1.

##  IV Discussion

Figure 10 displays the same SANS and SAXS data as in Fig. 1, fitted here with
the Gaussian membrane model. The fitting parameters are (i) the clipping
characteristics $Z_{n}$’s and $l_{\alpha}$ (see Fig. 2a) and (ii) the field
characteristic lengths $l_{XY}$ and $l_{Z}$. The four $Z$’s are symmetric with
respect to the origin (the membrane is arbitrarily centred on $z=0$), and are
therefore only two degrees of freedom, equivalent to specifying the
thicknesses of the chain and head sections $l_{C}$ and $l_{H}$. The least-
square fit yields the values $Z_{2}=22$ Å, $Z_{3}=29$ Å (with $Z_{0}=-Z_{3}$
and $Z_{1}=-Z_{2}$), $l_{\alpha}=38$ Å, $l_{XY}=80$ Å, and $l_{Z}=233$ Å. A
realization of the model with those sets of parameters is displayed in Fig.
10c. The overall amplitude of the fluctuations (controlled by parameter
$l_{\alpha}$) is found to be of the order of 10 nm, which is globally in line
with independently-determined values for vesicles.Monzel and Sengupta (2016)

Figure 10: Fitting of the SANS (a) and SAXS (b) data of the phospholipid
vesicles (same as in Fig. 1) with the Gaussian membrane model. The dots are
the data, and the solid red line is the fitted model. The solid black line is
the scattering from the average structure, and the dashed red line is the
contribution from the fluctuations. The horizontal dashed black line is the
background. A particular realization of the fitted model is shown in c.

It is interesting to compare the values of $Z_{2}=22$ Å and $Z_{3}-Z_{2}=7$ Å
to the slightly smaller values $l_{C}=17$ Å and $l_{H}=4.5$ Å inferred from
the deterministic slab model (see Tab. 1). To make a meaningful comparison,
however, one has to note that the values of $Z$’s control the volumes of the
layers and not directly their thicknesses. The volume of a layer can indeed be
calculated as the integral of its one-point probability function defined in
Eq. (23), namely

$\frac{V_{n}}{A}=\int_{-\infty}^{+\infty}S^{(1)}_{n}(z)\textrm{d}z$ (51)

which is here expressed per unit of area of the projected membrane $A$,
orthogonal to $z$. Based on Eq. (24) the relation is found to be simply
$V_{n}/A=Z_{n}-Z_{n-1}$. In other words, increasing the amplitude of the
membrane fluctuations through parameter $l_{\alpha}$ does not modify the
membrane volume. However, the average thickness has to be calculated by
dividing the volume by the actual surface area of the distorted membrane,
which can be significantly larger than $A$. The average thickness is therefore
evaluated as

$\langle l_{n}\rangle=\frac{Z_{n}-Z_{n-1}}{a_{A}}$ (52)

where $a_{A}$ is the dimensionless roughness factor calculated in Sec. III.1.
With the values of the fitted parameters, the roughness factor is $a_{A}\simeq
1.35$, shown as a red dot in Fig. 3. With this correction, the values from the
Gaussian membrane fit are $\langle l_{C}\rangle=16.3$ Å and $\langle
l_{H}\rangle=5.2$ Å, in fair agreement with the flat slab model.

In addition to the average thickness, another important characteristic of the
membrane is its compressibility, which can be quantified as the variability
(in time and space) of the membrane thickness. In the case where $l_{\alpha}$
is smaller than $l_{Z}$, as is the case for the structure in Fig. 10, the
position $z_{n}$ of the $n^{th}$ interface is approximately a Gaussian
variable with mean value $Z_{n}$ and standard deviation $l_{\alpha}$.
Moreover, the correlation between the fluctuations of the $n^{th}$ and
$m^{th}$ interfaces is $g_{W}(r_{xy},r_{z}=Z_{m}-Z_{n},\tau)$. The variance of
the thickness fluctuation - for one given position in the $xy$ plane
$r_{xy}=0$, and one given time $\tau=0$ \- is therefore obtained by
considering the difference of two correlated Gaussian variables, namely

$\langle(z_{n}-z_{m})^{2}\rangle=(Z_{n}-Z_{m})^{2}+2l_{\alpha}^{2}\left[1-g_{W}(r_{xy}=0,r_{z}=Z_{n}-Z_{m},\tau=0)\right]$
(53)

where the first and second terms account for the average and fluctuating
contributions. Using the expression in Eq. (5) for the field correlation
function $g_{W}$ then provides the following expression for the standard
deviation of the thickness fluctuation relative to the mean

$\frac{\sigma_{z_{n}-z_{m}}}{Z_{n}-Z_{m}}\simeq\sqrt{2}\
\frac{l_{\alpha}}{l_{Z}}$ (54)

in the limit where $|Z_{n}-Z_{m}|\ll l_{Z}$. Using the values of $l_{\alpha}$
and $l_{Z}$ from the SAXS and SANS fit, one obtains the numerical value
$\sigma_{z_{n}-z_{m}}/(Z_{n}-Z_{m})\simeq 0.23$. This compares resonably with
the 15 % polydispersity of the slab model fitted in Fig. 1, but the
interpretation is different here. In the slab model, one assumes that all the
membranes have permanently slightly different thicknesses. In the Gaussian
model, all the membranes are indeed identical but they fluctuate, also in
thickness. Because SAXS and SANS capture only an instantaneous snapshot of a
fluctuating structure, they cannot discriminate between the two scenarios.

The fluctuations, however, are central to the inelastic scattering underlying
the neutron spin echo (NSE) data in Fig. 1a. In classical approaches for
membrane scattering data analysis, the elastic (SANS, SAXS) and inelastic
(NSE) data are analyzed independently of one another. Elastic scattering data
are typically analyzed through static models, which can be structurally
sophisticated Pabst  _et al._ (2000); Heberle  _et al._ (2013); Chakraborty
_et al._ (2020); Chappa  _et al._ (2021); Lewis-Laurent  _et al._ (2021) and
possibly account also for the diffuse scattering from membrane
deformation.Hamley (2022) In classical approaches, the inelastic scattering
data is analyzed independently of the SANS or SAXS, e.g. through a Zilman-
Granek approach Zilman and Granek (1996, 2002) whereby the intermediate
scattering function is modelled as stretched exponentials at given $q$’s, with
parameters dependent on the elastic moduli and viscosities of the membrane
constituents.Nagao and Seto (2023) The present Gaussian random layer model
offers the possibility of building the inelastic scattering data analysis on
the results of the SAXS or SANS.

Figure 11: One-parameter fits of the Neutron Spin Echo data with the Gaussian
membrane model, for the exponential (a, $\tau_{c}=925$ ns) and diffusive (b,
$D=1.65$ Å2/ns) time dependences. The same data are reproduced for clarity in
c1 ($q=0.046$ Å-1), c2 ($q=0.082$ Å-1), and c3 ($q=0.122$ Å-1), together with
the exponential (red) or diffusive (blue) dynamic models. Panels d1 to d3 are
the same data on logarithmic scales to highlight stretched exponentials; the
solid black line is a standard exponential. The dashed black lines is the
Zilman-Granek model (same fit as in Fig. 1a). The solid black lines shown at
$\tau\simeq 400$ ns in a and b is the asymptotic value
$I(q,\tau\to\infty)/I(q,0)$.

Figure 11 compares the NSE data from Fig. 1c (dots) with the intermediate
scattering function calculated from the Gaussian membrane model. We consider
only the two non-differentiable dynamics, corresponding to Figs. 8a and 8c.
For the calculations, the values of the structural parameters inferred from
the SAXS and SANS are kept constant; the only additional and fitting parameter
is the one describing the time-dependence of the Gaussian field through the
$\tau$ dependence of the correlation function $g_{W}(\mathbf{r},\tau)$. Two
models are considered here, namely: the exponential correlation function (Eqs.
41 and 43) parameterized through the correlation time $\tau_{c}$, and the
diffusive model (Eq. 48) parameterized through diffusion coefficient $D$. The
two models fit equally well the NSE data with either $\tau_{c}=925$ ns or
$D\simeq 1.65$ Å2/ns, as shown in Figs. 11a and 11b, respectively. To compare
the two models, one has to notice that the correlation time of the diffusive
Gaussian field can be estimated as the time needed for a two-dimensional
random walker to travel a distance comparable with the correlation length
$l_{XY}$, i.e. as $l_{XY}^{2}/(4D)\simeq 970$ ns. Expectedly, this value is
comparable with $\tau_{c}$ obtained from the fit of the exponential
correlation. A time-dependent realization of the exponential model is shown in
Fig. 12, over a duration of approximately $2\tau_{c}$. Note how both the
bending and thickness fluctuations are captured in the latter realizations.

Figure 12: Time-dependent realization of the Gaussian random layer model, with
parameters fitted from the SAXS and SANS data (Fig. 10c), and exponential
time-dependence of $g_{W}(\mathbf{r},\tau)$ with correlation time from the NSE
data (Fig. 11).

Figures 11c1 to 11c3 compare the fits of the two Gaussian membrane models with
that of the Zilman-Granek model in the form of Eq. (3), i.e. with a single
fitting parameter $B$. To understand the qualitative differences between the
models, the same data are plotted again in Figs. 11d1 to 11d3 in the form of
$\log[-\log[I(q,\tau)/I(q,0)]]$ versus $\log(\tau)$ in order for stretched
exponentials to appear as straight lines. On such plots, it appears that the
intermediate scattering function of the Gaussian membrane follows the same
stretched-exponential of $\tau$ as the Zilman-Granek model (with exponent
$2/3$) only for large $\tau$ (see e.g. Fig. 11d3). For short times $\tau$ the
dependence is closer to a classical exponential (with exponent 1), which
incidentally seems to describe the NSE data better. The two smooth models
considered in Sec. III.2.3 \- the hyperbolic secant correlation (Fig. 8b) and
the ballistic model (Fig. 8d) - lead to quadratic scaling
$I(q,\tau)\simeq\exp(-\tau^{2})$ for short $\tau$ (not shown), and they are
therefore not suitable for the NSE data analysis. This is in line with earlier
observations with emulsion scattering.Gommes  _et al._ (2021)

In the context of NSE studies, the question of membrane compressibility is
often addressed through head-contrast experiments, which exhibit enhanced
dynamics at intermediate values of $q$.Nagao (2009); Woodka  _et al._ (2012)
Such data are usually analyzed by fitting the intermediate scattering function
with a stretched exponential $I(q,\tau)/I(q,0)=\exp[-(\Gamma\tau)^{2/3}]$ to
extract a $q$-dependent relaxation rate $\Gamma(q)$. When plotting the results
as $\Gamma(q)/q^{3}$ \- to visualize deviations from the $q^{3}$ scaling of
the Zilman-Granek theory - the data often exhibit a peak, which is considered
as a signature of thickness fluctuations. Chakraborty  _et al._ (2020); Nagao
and Seto (2023); Luo  _et al._ (2023) The present Gaussian model offers the
possibility of exploring this data analysis procedure, based on an a
mathematically exact model.

Figure 13: Intermediate scattering functions (a and d) and SANS intensities (c
and f) calculated from the Gaussian model with the same parameter as inferred
from the SANS and NSE fits, only with different values of the correlation
length $l_{Z}$, equal to 1000 Å (a1 and d1), 200 Å (a2 and d2) and 100 Å (a3
and d3). The left column is for chain contrast (a to c) and the right column
is for head contrast (d to f). The corresponding relaxation rates are plotted
in b and e.

Figure 13 plots the intermediate scattering function of the Gaussian membrane
model for both chain and head contrasts, using the same parameters as inferred
from the SANS and NSE fits in Figs. 10 and 11, only with different values of
the field correlation length $l_{Z}$ perpendicular to the membrane. The latter
parameter controls the thickness fluctuations of the membrane through Eq.
(54): the values used in the figure are $l_{Z}=1000$ Å, 200 Å and 100 Å,
corresponding to 5 %, 25 % and 50 % thickness fluctuations, respectively.
Because the intermediate scattering function is well described by a simple
exponential for small $\tau$ (see Fig. 11d1 to 11d3), the relaxation rate
$\Gamma(q)$ was evaluated by fitting the calculated values as
$I(q,\tau)/I(q,0)\simeq\exp[-\Gamma\tau]$ and the values are reported in Fig.
13b and 13e. In the case of chain contrast the relaxation rate follows an
overall scaling close to $q^{2}$, with no visible impact of the thickness
fluctuation. In that respect, the Gaussian membrane model exhibits a typical
diffusion phenomenology with relaxation rate scaling linearly with $\tau$ and
quadratically with $q$. The case of head contrast is distinctly different; the
relaxation rate is found to exhibit a peaked maximum, which is clearly
correlated with thickness fluctuations. The $q$ position of the maximum
approximately corresponds to the minimum of the SANS intensity (see Fig. 13f).

##  V Conclusions

We have developed a mathematical model of a fluctuating membrane, for the
purpose of jointly analyzing small-angle scattering of x-rays and/or neutrons,
as well as neutron spin-echo data, with arbitrary contrasts, within a single
theoretical framework. The model builds on a non-homogeneously clipped time-
dependent Gaussian random field, which modelling procedure provides one with
general analytical expressions for the intermediate scattering function. The
availability analytical results significantly reduces the computational burden
of data analysis; all analyses and simulations presented in this paper were
conducted on a laptop computer.

Our Gaussian membrane model realistically captures the structure and dynamics
of fluctuating membranes with a handful of meaningful parameters. Although the
model is descriptive, exploring the influence of its parameters provides
direct insight into the physics of elastic and inelastic scattering by
membranes. In particular, our analysis enabled us to understand how membrane
compressibility contributes to forward elastic scattering. It also confirmed
through exact mathematical calculations how compressibility significantly
enhances the relaxation rates in neutron spin echo data, in the case of head-
contrast only.

Moreover the synthetic description of membrane structure and dynamics, enables
one to comprehensively analyze large datasets comprising both elastic and
inelastic scattering patterns measured for arbitrary contrasts with one and a
single model. All the structural and dynamic information is then condensed in
a coherent and robust way into a few meaningful parameters. This contrasts
with classical approaches in this context, where elastic and inelastic
scattering patterns are traditionally analyzed independently of one another
with data-specific models. Analyzing jointly different datasets, also offers
the prospect of extracting meaningful information that is not present in the
individual datasets when they are considered separately.

We illustrated our approach with the analysis of SAXS, SANS and NSE data
measured at 37 ∘C on vesicles prepared from porcine brain phospholipids, close
to physiological conditions. The elastic scattering data were analyzed first,
to identify the structural parameters. The latter are equivalent to (i) the
thicknesses of the hydrophobic chains and hydrophilic heads of the membrane
molecules, (ii) the amplitude of the membrane deformations as well as (iii)
their lateral sizes, and (iv) the thickness fluctuation. The latter three
parameters are absent from classical membrane models, but they can be inferred
from the joint SAXS and SANS analysis. The membranes investigated in this work
are found to undergo bending fluctuations with amplitude around 10 nm and
thickness fluctuations of the order of 8 Å.

Building on these values of the structural parameters, the complete NSE
dataset was fitted with a single additional parameter characterizing the
dynamics of the fluctuations. The two specific dynamic models we considered
fit equally well the data. Interestingly, they both predict a $\tau^{2/3}$
scaling of the relaxation rate (similar to the Zilman-Granek model) for
asymptotically large times. At finite times, however, the calculated
relaxation rate is diffusion-like and it scales linearly with $\tau$ and
quadraticaly with $q$. Our data incidentally exhibit this specific type of
scaling.

###### Acknowledgements.

CJG is grateful to the Funds for Scientific Research (F.R.S.-FNRS, Belgium)
for a Research Associate position, and for supporting part of this work
through grant PDR T.0100.22.

## References

  * Watson (2015) H. Watson, Essays Biochem. 59, 43 (2015).
  * Discher and Eisenberg (2002) D. E. Discher and A. Eisenberg, Science 297, 967 (2002).
  * Guinart  _et al._ (2023) A. Guinart, M. Korpidou, D. Doellerer, G. Pacella, M. C. Stuart, I. A. Dinu, G. Portale, C. Palivan, and B. L. Feringa, Proc. Nat. Acad. Sci. USA 120 (2023).
  * Krugmann  _et al._ (2020) B. Krugmann, A. Radulescu, M.-S. Appavou, A. Koutsioubas, L. R. Stingaciu, M. Dulle, S. Förster, and A. M. Stadler, Scientific Reports 10, 16691 (2020).
  * Krugmann  _et al._ (2021) B. Krugmann, A. Koutsioubas, L. Haris, S. Micciulla, D. Lairez, A. Radulescu, S. Förster, and A. M. Stadler, Frontiers in Chemistry 9 (2021), 10.3389/fchem.2021.631277.
  * Büldt  _et al._ (1978) G. Büldt, H. Gally, A. Seelig, J. Seelig, and G. Zaccai, Nature 271, 182 (1978).
  * Pabst  _et al._ (2010) G. Pabst, N. Kucerka, M.-P. Nieh, M. Rheinstädter, and J. Katsaras, Chem. Phys. Lipids 163, 460 (2010).
  * Van Hove (1954) L. Van Hove, Phys. Rev. 95, 249 (1954).
  * Sivia (2011) D. S. Sivia,  _Elementary Scattering Theory_ (Oxford University Press, Oxford, 2011).
  * Squires (2012) G. L. Squires,  _Introduction to the Theory of Thermal Neutron Scattering_ , 3rd ed. (Cambridge University Press, 2012).
  * Glatter (2018) O. Glatter,  _Scattering Methods and their Application in Colloid and Interface Science_ (Elsevier, 2018).
  * Luo  _et al._ (2023) X. Luo, T. Cui, and X. Chu, Front. Phys. 11, 1279007 (2023).
  * Pedersen (1997) J. S. Pedersen, Adv. Coll. Interf. Sci. 70, 171 (1997).
  * Svergun (1999) D. Svergun, Biophys. J. 76, 2879 (1999).
  * Gommes (2018) C. J. Gommes, Microp. Mesop. Mater. 257, 62 (2018).
  * Kiselev  _et al._ (2008) M. Kiselev, E. Zemlyanaya, N. Ryabova, T. Hauss, S. Dante, and D. Lombardo, Chem. Phys. 345, 185 (2008).
  * Kiselev  _et al._ (2002) M. A. Kiselev, P. Lesieur, A. M. Kisselev, D. Lombardo, and V. L. Aksenov, Appl. Phys. A 74, S1654 (2002).
  * Pencer  _et al._ (2006) J. Pencer, S. Krueger, C. P. Adams, and J. Katsaras, J. Appl. Crystallogr. 39, 293 (2006).
  * Chappa  _et al._ (2021) V. Chappa, Y. Smirnova, K. Komorowski, M. Muller, and T. Salditt, J. Appl. Crystallogr. 54, 557 (2021).
  * Monzel and Sengupta (2016) C. Monzel and K. Sengupta, J. Phys. D 49 (2016).
  * Hamley (2022) I. W. Hamley, Soft Matter 18, 711 (2022).
  * McMahon and Gallop (2005) H. T. McMahon and J. L. Gallop, Nature 438, 590 (2005).
  * Zidovska and Sackmann (2006) A. Zidovska and E. Sackmann, Phys. Rev. Letters 96, 048103 (2006).
  * Venable  _et al._ (1993) R. M. Venable, Y. Zhang, B. J. Hardy, and R. W. Pastor, Science 262, 223 (1993).
  * Sodt  _et al._ (2014) A. J. Sodt, M. L. Sandar, K. Gawrisch, R. W. Pastor, and E. Lyman, J. Am. Chem. Soc. 136, 725 (2014).
  * Chakraborty  _et al._ (2020) S. Chakraborty, M. Doktorova, T. R. Molugu, F. A. Heberle, H. L. Scott, B. Dzikovski, M. Nagao, L. Stingaciu, R. F. Standaert, F. N. Barrera, J. Katsaras, G. Khelashvili, M. F. Brown, and R. Ashkar, Proc. Nat. Acad. Sci. USA 117, 21896 (2020).
  * Hayward  _et al._ (2022) D. W. Hayward, S. Jaksch, M. Fomina, P. S. Dubey, H. Frielinghaus, O. Holderer, and M. Monkenbusch, Acta crystallogr. D 78, 1249 (2022).
  * Zilman and Granek (1996) A. G. Zilman and R. Granek, Phys. Rev. Lett. 77, 4788 (1996).
  * Zilman and Granek (2002) A. G. Zilman and R. Granek, Chem. Phys. 284, 195 (2002).
  * Nagao and Seto (2023) M. Nagao and H. Seto, Biophys. Rev. 4 (2023), 10.1063/5.0144544.
  * Granek  _et al._ (2024) R. Granek, I. Hoffmann, E. G. Kelley, M. Nagao, P. M. Vlahovska, and A. Zilman, Eur. Phys. J. E 47 (2024), 10.1140/epje/s10189-023-00400-9.
  * Nagao (2009) M. Nagao, Phys. Rev. E 80, 031606 (2009).
  * Woodka  _et al._ (2012) A. C. Woodka, P. D. Butler, L. Porcar, B. Farago, and M. Nagao, Phys. Rev. Lett. 109, 058102 (2012).
  * Berk (1987) N. Berk, Phys. Rev. Lett. 58, 2718 (1987).
  * Berk (1991) N. F. Berk, Phys. Rev. A 44, 5069 (1991).
  * Teubner (1991) M. Teubner, Europhys. Lett. 14, 403 (1991).
  * Chen  _et al._ (1996) S. H. Chen, D. D. Lee, K. Kimishima, H. Jinnai, and T. Hashimoto, Phys. Rev. E 54, 6526 (1996).
  * Roberts (1997) A. P. Roberts, Phys. Rev. E 55, R1286 (1997).
  * Gommes and Roberts (2008) C. J. Gommes and A. P. Roberts, Phys. Rev. E 77, 041409 (2008).
  * Levitz (1998) P. Levitz, Adv. Coll. Interf. Sci 76, 71 (1998).
  * Prehal  _et al._ (2017) C. Prehal, C. Koczwara, N. Jäckel, A. Schreiber, M. Burian, H. Amenitsch, M. A. Hartmann, V. Presser, and O. Paris, Nature Energy 2, 16215 (2017).
  * Gommes and Pirard (2009) C. J. Gommes and J. P. Pirard, Phys. Rev. E 80, 061401 (2009).
  * Gommes  _et al._ (2020) C. J. Gommes, R. Chattot, and J. Drnec, J. Appl. Crystallogr. 53, 811 (2020).
  * Gommes (2013) C. J. Gommes, J. Appl. Crystallogr. 46, 493 (2013).
  * Gommes  _et al._ (2021) C. J. Gommes, R. Zorn, S. Jaksch, H. Frielinghaus, and O. Holderer, J. Chem. Phys. 155 (2021).
  * Gommes (2022) C. J. Gommes, Gels 8 (2022), 10.3390/gels8040236.
  * Folch  _et al._ (1951) J. Folch, I. Ascoli, M. Lees, J. Meath, and F. LeBaron, J. Bio. Chem. 191, 833 (1951).
  * Dubey  _et al._ (2023) P. S. Dubey, O. Czakkel, H. Frielinghaus, I. Hoffmann, O. Holderer, S. Jaksch, O. Matsarskaia, J. Ollivier, S. Prevost, and D. Saha, “Influence of NSAIDs on the phospholipid mammalian neuronal membrane,”  (2023).
  * Nagao  _et al._ (2017) M. Nagao, E. G. Kelley, R. Ashkar, R. Bradbury, and P. D. Butler, The journal of physical chemistry letters 8, 4679 (2017).
  * Quiblier (1984) J. A. Quiblier, J. Coll. Interf. Sci. 98, 84 (1984).
  * Lantuéjoul (2002) C. Lantuéjoul,  _Geostatistical Simulations_ (Springer, Berlin, 2002).
  * Gommes and Roberts (2018) C. J. Gommes and A. P. Roberts, Phys. Chem. Chem. Phys. 20, 13646 (2018).
  * Lewis-Laurent  _et al._ (2021) A. Lewis-Laurent, M. Doktorova, F. A. Heberle, and D. Marquardt, Biophys. J. 120, 4639 (2021).
  * Lemmich  _et al._ (1996) J. Lemmich, K. Mortensen, J. H. Ipsen, T. Honger, R. Bauer, and O. G. Mouritsen, Phys. Rev. E 53, 5169 (1996).
  * Maiti  _et al._ (2021) S. Maiti, H. Frielinghaus, D. Grässel, M. Dulle, M. Axer, and S. Förster, Sci. Rep. 11 (2021), 10.1038/s41598-021-92995-2.
  * Lantuejoul (1991) C. Lantuejoul, J. Microsc. 161, 387 (1991).
  * Lighthill (1958) M. J. Lighthill,  _An Introduction to Fourier Analysis and Generalised Functions_ , Cambridge Monographs on Mechanics (Cambridge University Press, 1958).
  * Gelfand  _et al._ (2010) A. E. Gelfand, P. Diggle, P. Guttorp, and M. Fuentes,  _Handbook of Spatial Statistics_ (CRC Press, Boca Raton, 2010).
  * Berg (1993) H. Berg,  _Random Walks in Biology_ , 2nd ed. (Princeton University Press, 1993).
  * Pabst  _et al._ (2000) G. Pabst, M. Rappolt, H. Amenitsch, and P. Laggner, Phys. Rev. E 62, 4000 (2000).
  * Heberle  _et al._ (2013) F. A. Heberle, R. S. Petruzielo, J. Pan, P. Drazba, N. Kuerka, R. F. Standaert, G. W. Feigenson, and J. Katsaras, J. Am. Chem. Soc. 135, 6853 (2013).