research papers
The effect of polydispersity, shape fluctuations and curvature on small unilamellar vesicle smallangle Xray scattering curves
^{a}Faculty of Physics, University of Göttingen, FriedrichHundPlatz 1, Göttingen 37077, Germany
^{*}Correspondence email: ysmirno@gwdg.de
Small unilamellar vesicles (20–100 nm diameter) are model systems for strongly curved lipid membranes, in particular for cell q features, whereas thermal fluctuations predominantly affect the scattering at larger q, related to membrane bending rigidity. Finally, it is shown that simulation of fluctuating vesicle ensembles can be used for analysis of experimental SAXS curves.
Routinely, smallangle Xray scattering (SAXS) is employed to study their size and electrondensity profile (EDP). Current SAXS analysis of small unilamellar vesicles (SUVs) often employs a factorization into the (vesicle shape) and the form factor (lipid bilayer electrondensity profile) and invokes additional idealizations: (i) an effective polydispersity distribution of vesicle radii, (ii) a spherical vesicle shape and (iii) an approximate account of membrane asymmetry, a feature particularly relevant for strongly curved membranes. These idealizations do not account for thermal shape fluctuations and also break down for strong salt or proteininduced deformations, as well as vesicle adhesion and fusion, which complicate the analysis of the lipid bilayer structure. Presented here are simulations of SAXS curves of SUVs with experimentally relevant size, shape and EDPs of the curved bilayer, inferred from coarsegrained simulations and elasticity considerations, to quantify the effects of size polydispersity, thermal fluctuations of the SUV shape and membrane asymmetry. It is observed that the factorization approximation of the scattering intensity holds even for small vesicle radii (∼30 nm). However, the simulations show that, for very small vesicles, a curvatureinduced asymmetry arises in the EDP, with sizeable effects on the SAXS curve. It is also demonstrated that thermal fluctuations in shape and the size polydispersity have distinguishable signatures in the SAXS intensity. Polydispersity gives rise to lowKeywords: small unilamellar vesicles; coarsegrained simulations; elastic simulations; smallangle Xray scattering; SAXS curves.
1. Introduction
The shape of fluctuating membranes has received abiding attention in the context of measuring bending rigidity (Gompper & Kroll, 1997) and explaining frequent, but surprising, membrane shapes such as the discoid shape of red blood cells (Canham, 1970; Helfrich, 1973; Seifert et al., 1991; Discher et al., 1994; Safran, 1994; Lim et al., 2002; Li et al., 2005). One of the most widely used experimental techniques to study the structure and shape of lipid bilayers and vesicles is smallangle Xray scattering (SAXS) (Zubay, 1999; Pabst et al., 2000; Kiselev et al., 2002; Brzustowicz & Brunger, 2005; Pencer et al., 2006; Kučerka et al., 2008; Székely et al., 2010; Heberle et al., 2012). It allows information to be obtained starting from the molecular level, such as the location of a particular chemical entity (e.g. carbon double bonds), up to the overall vesicle shape. However, a quantitative analysis of experimental scattering data that includes the effects of (i) size polydispersity, (ii) vesicle shape fluctuations and (iii) the curvaturedependent electrondensity profile (EDP) has remained challenging.
Theoretical methods, such as et al., 2009) and continuum elastic vesicle models (Gompper & Kroll, 1996), can help to interpret and understand experimental data by analysing contributions from different physical factors. For small unilamellar vesicles (SUVs) with radii in the range of 10–50 nm, such methods as coarsegrained MD simulations and the elastic Helfrich (1973, 1986) model can be used efficiently. SUVs are also valuable model systems for studying membrane adhesion and fusion (Komorowski et al., 2018). It is therefore worthwhile to combine simulation techniques and SAXS in order to better understand small unilamellar vesicle shapes and their membrane EDPs. In this work, we use a hierarchical simulation framework that employs different models and simulation techniques on different scales (Müller et al., 2003; Müller et al., 2006): On the scale of the bilayer thickness, we obtain a curvaturedependent radial EDP of a lipid membrane from MD simulations using the coarsegrained MARTINI force field (Marrink et al., 2007). We then use this profile to dress various vesicle shapes generated using the elastic Helfrich model (Seifert, 1997). The resulting threedimensional electrondensity map is used to calculate the scattering intensity via a threedimensional (3D) fast Fourier transform (FFT) and subsequent powder averaging, taking into account many realizations for an ensemble averaging over thermal shape fluctuations and size polydispersity. We compare the simulation results of our numerical model with existing analytical SAXS models (Brzustowicz & Brunger, 2005), and demonstrate its capability to be used for leastsquares fitting of experimental SAXS curves obtained from small unilamellar lipid vesicles in the fluid phase. To this end, we present two examples: (i) SAXS data for a (1:1) mixture of dioleoylphosphatidylcholine and dioleoylphosphatidylethanolamine, formed after extrusion through 30 nm membranes in ultrapure water (MilliQ), and (ii) a (1:1) mixture of dioleoylphosphatidylcholine and dioleoylphosphatidylserine.
(MD) simulations (MarrinkThe simulation and FFT approach allows us to analyse the interplay between size polydispersity, vesicle shape fluctuations and curvaturedependent EDP. Since in the simulation we can control these different phenomena independently, we compare the relative effects of these phenomena on the distinct wavevector regions of the SAXS intensity, and thereby decouple their effect. The goal of our paper is to propose a strategy to analyse scattering data with respect to these three effects. Thereby, we address the limitations of current SUV SAXS analysis paradigms, namely idealization of one or more of the following effects: (i) polydispersity distribution of vesicle radii, (ii) thermal vesicle shape fluctuations and (iii) curvature dependence of the EDP. At the same time we can identify the range of parameters under which common idealizations are justified. In particular, we consider the validity of the factorization approximation of the SAXS intensity, i.e. the common assumption in analytical SAXS models that the scattering function can be described as a product of the bilayer form factor and the of the vesicle shape (Kiselev et al., 2002, 2006). We hope that this contribution will facilitate the analysis of scattering data, which is particularly relevant given the currently increasing complexity in smallangle scattering experiments (Semeraro et al., 2021). Beyond the specific effects considered here, we also want to develop the approach of simulating vesicle ensembles and subsequent 3D FFT in view of the analysis of more complex vesicles and vesicle shape transitions, as required for example in studies of synaptic vesicles. While here we primarily consider SAXS, the approach can equally be used for smallangle neutron scattering studies of vesicles, after minor modification regarding the scattering lengths.
2. Methods
In the following we detail our hierarchical modelling approach, illustrated in Fig. 1. The system has two scales – the larger one, R_{0} ≳ 10 nm, is associated with the vesicle radius and the smaller scale, d_{0} ≃ 4 nm, with the thickness of the membrane. We partition our model according to these two scales. The vesicle shape and its fluctuations are described by the Helfrich Hamiltonian (Helfrich, 1986; Milner & Safran, 1987) which characterizes a vesicle by its average size, R_{0}, whose polydispersity obeys a distribution P(R_{0}) and bending rigidity κ. The membrane structure and the electron density, ρ(d), in turn, characterize the small scale and are investigated by MD simulations of the coarsegrained MARTINI model (Marrink et al., 2007). These simulations provide detailed density profiles across the bilayer, from which we extract the electrondensity contrast that dictates the SAXS intensity. The two scales are coupled via the curvature dependence of the EDP.
First, we describe the study of the curvaturedependent membrane structure by MD simulations, comparing a planar membrane with a very small vesicle of radius ≲10 nm. Subsequently, we explain how we efficiently sample the largescale shape fluctuations by simulating the Helfrich Hamiltonian, using an expansion of the radius in spherical harmonics of the polar and azimuthal angles. Finally, dressing the vesicle shape with a curvaturedependent EDP, we combine these two descriptions to calculate SAXS scattering curves.
2.1. MARTINI model simulation and EDPs
2.1.1. Simulation protocol
Coarsegrained MD simulations of a highly curved vesicle and planar lipid bilayer were performed using the GROMACS simulation package (Abraham et al., 2015) in conjunction with the MARTINI force field (Marrink et al., 2007). A small 1palmitoyl2oleoylsnglycero3phosphocholine (POPC) vesicle was formed by spontaneous aggregation, following the protocol of Risselada et al. (2008, 2014). The vesicle is composed of 1447 in the outer leaflet and 770 in the inner leaflet, embedded in a solvent that contains a total of 97 217 coarsegrained solvent particles. The vesicle radius is about 8 nm, which is twice as large as the membrane thickness, d_{0} = 4 nm. A planar bilayer patch was simulated at full hydration and contained 2048 POPC All systems were equilibrated in the NPT ensemble and simulated for 1 µs in the NVT ensemble at T = 300 K to calculate EDPs. The simulations provide detailed beadtype profiles across the planar or highly curved bilayer, as presented in Fig. 2.
2.1.2. Electrondensity profiles
From these beadtype profiles we extracted the EDPs of the planar membrane and the highly curved vesicle. The numbers of electrons per coarsegrained bead are listed in Table 1 [see also Wanga et al. (2016)]. We used a simplified approach in which the centre of mass of the electron cloud coincides with the centre of mass of the coarsegrained bead.

Fig. 3 summarizes the results from the coarsegrained MARTINI simulations. Fig. 3(a) presents the EDP of the POPC bilayer patch and the small vesicle of radius R_{0} = 8 nm, where we have subtracted the electron density of the water. We fitted these EDPs using two Gaussian functions and a fourthorder polynomial:
where d is the displacement from the bilayer's midplane, and μ_{l} and μ_{r} denote the positions of the left (inner) and right (outer) peaks, respectively. These maxima in the electron density correspond to the headgroup region of the bilayer membrane. σ_{l} and σ_{r} are the and ρ_{l} and ρ_{r} characterize the peak heights. The polynomial fits the electrondensity distribution in the bilayer's hydrophobic interior. a and b are two roots of the fourth polynomial, and c is its value at the midplane, d = 0.
The parameters extracted by a nonlinear leastsquares fit are compiled in Table 2. The fit was made with the Python package SciPy using the method optimize.curve_fit. The fit depends on the initial guess of the parameters and we selected parameters for the fourth polynomial such that the two real roots a and b are close to the vanishing tails of the electron density. The profiles of the planar membrane and highly curved vesicle differ mainly in two aspects: the heights and the positions of the two peaks. In the case of the planar membrane the EDP has two symmetric maxima, whereas for the highly curved vesicle these peaks are shifted closer to each other and differ in height. The higher peak corresponds to the inner monolayer of the small vesicle, where the lipid head groups are tightly packed because of the high curvature, whereas the outer monolayer is characterized by a lower headgroup density due to the smaller curvature (Fig. 2). For the lipid tails the effect is opposite: when referred to the midplane of the bilayer, the area per molecule is lower in the inner monolayer of the vesicle than in the outer, and thus the volume density is also slightly lower in the inner tail region than in the hydrophobic portion of the outer monolayer. Additionally, the high curvature of the vesicle imparts a tension onto the membrane. This results in a thinning of the bilayer that is observable in the inwards shift of the position of the maxima of the electron density.

By virtue of symmetry, a model for the curvature dependence of the EDP of a vesicle can be written in the form of a curvature expansion:
where d denotes the displacement from the centre of the vesicle's midplane, and d_{0} and R_{0} are the bilayer thickness and vesicle radius, respectively. Note that this expansion of the EDP up to second order in curvature is compatible with the smallcurvature expansion that underlies the Helfrich Hamiltonian, i.e. both the largescale description of the vesicle shape and its fluctuations and the response of the smallscale bilayer structure to curvature include effects up to second order. ρ_{m} is the density profile of the planar bilayer, whereas the curvature dependence is described by Δρ_{a}(d) and Δρ_{s}(d) which are odd (antisymmetric) and even (symmetric) functions with respect to their argument, d, respectively. The odd part describes the differences between the inner and outer monolayers in response to the membrane curvature, 1/R_{0}. The even function, Δρ_{s}(d), affects both monolayers in the same way, e.g. the curvatureinduced thinning of the membrane. Since these effects do not depend on the sign of the curvature, they scale like 1/R_{0}^{2} to leading order. If the symmetric contribution, Δρ_{s}(d), chiefly stems from the curvatureinduced thinning of the bilayer membrane, we can approximately relate its functional form to the profile of a planar bilayer by assuming that the profile of a bilayer of thickness d_{0} obeys the affine scaling relation = with the dimensionless = d/d_{0}. Thinning the bilayer by δd_{0}, we obtain for the change of profile
i.e. by symmetry, the thinning of the membrane scales like and the spatial dependence of the symmetric contribution is . Fig. 4 includes a comparison between Δρ_{s}(d), obtained by symmetrizing the difference ρ_{v}(d, R_{0}) − ρ_{m}(d) for the vesicle with R_{0} = 8 nm, and , where the amplitude has been adjusted because δd_{0}(R_{0}) is unknown. The good agreement indicates that (i) Δρ_{s} is, indeed, dominated by the curvatureinduced thinning of the bilayer membrane and (ii) the affine scaling relation is appropriate for the range of thinning induced by the curvature.
Using this model we generated EDPs for vesicles with different radii and these are shown in Fig. 3(c). Note that for vesicle radii R_{0} ≲ 10 nm we observe significant deviations from the EDP of the planar bilayer.
A popular approximation for the EDP of vesicles is the sum of three Gaussian peaks (Brzustowicz & Brunger, 2005):
where the three maxima represent the headgroup regions of the inner and outer monolayers, and the hydrophobic membrane interior, respectively. As in equation (1), the parameters ρ_{k}, μ_{k} and σ_{k} characterize the height, position and widths of the kth peak, respectively. The electron density of this threeGaussian model with symmetric inner and outer maxima is compared in Fig. 3(d) with the curvaturedependent model [equation (2)] for a small vesicle, R_{0} = 30 nm, and a highly curved vesicle, R_{0} = 8 nm. The parameters of the symmetric model [equation (4)] have been adjusted to fit the maximum of the electron density of the inner monolayer and the parameters of the fit are listed in Table 3. This comparison illustrates that our curvaturedependent model and the symmetric threeGaussian model provide a similar representation of the electron density for large radii, R_{0} ≳ 30 nm. Adjusting the parameters of the symmetric threeGaussian model to our curvaturedependent model, the electrondensity maxima coincide (by construction), but the electrondensity contrast is slightly underestimated in the tail region and at the interface between the head groups and the solvent. For smaller radii, R_{0} ≲ 10 nm, however, the two models for the electron density differ significantly because the asymmetry of the profile becomes important.

2.2. Generation of vesicle shapes and SAXS scattering intensity
2.2.1. Size polydispersity and shape fluctuations
In order to obtain the electron density of a vesicle in 3D space, we convolute the position of the bilayer's midplane with the curvaturedependent electrondensity contrast discussed in the previous section. Apart from translations, the position of the bilayer's midplane is dictated by the size of the vesicle and the thermal fluctuations of its shape. In order to consider size polydespersity we assume the vesicle sizes R_{0} are Gaussian distributed, according to
where characterizes the mean vesicle radius of the ensemble and denotes the (dimensional) standard deviation. The latter quantity characterizes the radius polydispersity. The Gaussian distribution is a common approximation for an equilibrated size distribution. Our computational scheme, however, can be straightforwardly generalized to different distributions.
Given the vesicle radius R_{0}, thermal fluctuations will result in deviations from a spherical shape. The freeenergy costs of deviations from a spherical shape are proportional to the bending rigidity, κ, and quantified by the Helfrich Hamiltonian. In the following we assume that the thermally excited deviations from the spherical shape remain small and we parameterize the position of the bilayer's midplane by the distance R(θ, ϕ) from the vesicle's centre of mass. In this spherical coordinate system, we expand R(θ, ϕ) in terms of spherical harmonics,
with , where are associated Legendre polynomials. Thus the shape of a vesicle, i.e. the position of the bilayer's midplane up to translation, is characterized by R_{0} and the set of u_{l,m}. In thermal equilibrium, the Helfrich Hamiltonian asserts that the fluctuation amplitudes, u_{l,m}, are statistically independent and Gaussian distributed with zero mean and variance
where κ is the bending rigidity of the membrane. The case where l = 1 merely corresponds to a translation of the centre of mass and is therefore omitted. We consider modes up to l_{max} = 6. In order to resolve shape fluctuations on the vesicle with a wavelength that is comparable to the membrane thickness d_{0} – the smallest wavelength where the Helfrich Hamiltonian is applicable – the order of spherical harmonics l_{max} should be chosen such that the ratio is of order unity.
The assumption of small deviations from a spherical shape in conjunction with the parameterization via spherical harmonics allows us to generate independent equilibrated vesicle shapes by (i) choosing a radius R_{0} according to equation (5) and (ii) drawing the coefficients u_{l,m} from Gaussian distributions with widths given by equation (7). The configurations thus generated are uncorrelated, i.e. unlike molecular simulations of particlebased models or simulations of dynamically triangulated surfaces we obtain a new membrane shape at each generation step.
2.2.2. Calculation of SAXS scattering intensity
In order to compute the SAXS scattering intensity from an ensemble of vesicles with different sizes and shapes, we combine the generation of configuration snapshots of the location of the membrane's midplane with the curvaturedependent EDPs. The electron density in 3D space is obtained by
where ρ_{v}(d, R_{0}) is the curvaturedependent EDP according to equation (2) and R(θ, ϕ) denotes the local distance of the bilayer's midplane from the vesicle's centre of mass [equation (6)]. This procedure neglects (i) the difference between the local curvature of the bilayer's midplane and the inverse size, 1/R_{0}, of the vesicle, and (ii) the difference between the local normal of the bilayer's midplane and the radial vector to the vesicle's centre of mass. Both approximations are consistent with the assumed small deviation of the fluctuating vesicle from its lowest energy state, a spherical vesicle.
In order to calculate the scattering intensity, we collocalize the electron density, ρ(r), on a regular cubic grid with N × N × N sites. Typically, we use N = 400 or 600 in each Cartesian direction, and the spatial extent of the collocation grid is L = . Thus the collocation grid can resolve spatial distances = . For small vesicles, with = 8, 10 and 20 nm, we use N = 400, resulting in δ/d_{0} = 0.1 for = 20 nm. The same resolution is obtained for the largest vesicle, = 30 nm, for N = 600. There are two scales, and d_{0}, and we choose the dimensionless wavevector to present our results. The scattering intensity is numerically obtained by FFT of the electron density ρ(r) on the regular cubic grid according to
The average runs over all orientations of scattering vectors q – the `powder average' – as well as independent snapshots of vesicle shapes. We use N_{v} = 1000 independent snapshots of the vesicle shape to compute the scattering intensity. The simulations were run on a parallel cluster with IvyBridge Intel E52670 v2 CPUs, 2.5 GHz 2 × 10 cores and 64 GB memory, requiring 2.5 s per vesicle per core. To compare the deviation between two scattering intensities, I(q) and I_{0}(q), we use the meansquared variance:
where N_{Q} denotes the number of q values in the considered interval.
3. Results and discussion
3.1. and form factor
To understand qualitatively the features of the scattering intensity I(q), we factorize I(q) into a which describes the shape of the vesicle and its fluctuations, and a form factor of the membrane, which contains information about the EDP of the bilayer membrane (Kiselev et al., 2002). Such a factorization approximation is justified if the two length scales, vesicle radius R_{0} and bilayer thickness d_{0}, are well separated, i.e. d_{0}/R_{0} << 1. Using the assignment of equation (8), we obtain for the scattering amplitude of a vesicle with shape R(θ, ϕ)
The scattering intensity is obtained by averaging over the orientations of q and the vesicle shapes [equation (9)]. For completeness, we recall the steps that result in the factorization approximation (Kiselev et al., 2002).
(i) In the absence of size polydispersity and thermal fluctuations, the vesicle shape is simply a sphere of radius R(θ, ϕ) = R_{0}, and we obtain
(ii) According to equation (2), the EDP ρ_{v}(d) is a sum of even and odd contributions. In the following, we keep all terms up to second order in the vesicle's curvature, d_{0}/R_{0}:
Inserting this expansion into equation (12), we obtain for the scattering amplitude of a thin spherical vesicle
Using this expression we obtain for the scattering intensity I(q) up to first order in curvature
The leading term is the popular factorization approximation, where I_{TS} is the of an infinitely thin spherical shell and A_{0}^{2} (q) is the powder average of the planar bilayer form factor. Equation (18) indicates under which conditions the factorization approximation is accurate.
The for R_{0} = 30 nm and d_{0} = 4 nm. Even for this vesicle size, one can clearly observe the separation of the two length scales R_{0} and d_{0}. Most notably, the radius of the vesicle can simply be obtained from the wavevector, q_{v} = π/R_{0}, at which the vanishes. The form factor of the planar bilayer in turn, [A_{0}(q)/d_{0}]^{2}, is only a function of qd_{0} and basically remains wavevector independent for q ≲ q_{m} = 2π/d_{0}. Fig. 5 also depicts the result for a spherical vesicle with curvaturedependent EDP. The curvature dependencies, Δρ_{a} and Δρ_{s}, and the firstorder correction term = 2cot(qR_{0}) (d_{0} / R_{0} ) A_{0} (q) B _{1}(q) appear to be negligibly small, even for .
of a spherical shell and the powderaveraged form factor of a planar bilayer are depicted in Fig. 53.2. Models of electrondensity profiles
In this part we will analyse the effect of different EDPs on the scattering intensity. The scattering amplitude of an ideally spherical vesicle with radius R_{0}, whose EDP is given by three Gaussian peaks [equation (4)], takes the form (Brzustowicz & Brunger, 2005)
where δ_{k} = R_{0} + μ_{k} denotes the distance of the headgroup and tail regions from the centre of mass of the vesicle. We consider an asymmetric EDP, where the headgroup peaks differ in their amplitudes, ρ_{1} ≠ ρ_{3}. They are assumed (i) to be symmetrically displaced by a distance μ = −μ_{1} = μ_{3} from the bilayer's midplane, μ_{2} = 0, and (ii) to have the same width, σ = σ_{1} = σ_{3}. Moreover, we assume that (iii) the widths of the Gaussian peaks and the bilayer thickness are small compared with the vesicle's radius, and , respectively. Because of condition (iii) we expand equation (21) up to second order in σ_{k}/R_{0} and μ_{k}/R_{0} and obtain
For giant vesicles (or, equivalently, extremely thin bilayer membranes), μ_{k}/R_{0} → 0 and σ_{k}/R → 0 and equation (22) predicts , i.e. the first minimum of the scattering intensity occurs at qR_{0} = π, in agreement with the factorization approximation. Fig. 6(a) shows the scattering intensity for vesicles of different radii R_{0} and a bilayer thickness of d_{0} = 4 nm for an EDP given by the symmetric threeGaussian model, i.e. ρ_{1} = ρ_{3}. The parameters are compiled in Table 3. We indeed observe that the first minimum is close to q_{v} = π/R_{0} for R_{0} = 30 nm. For smaller values of R_{0} comparable to the bilayer thickness, d_{0} ≃ 2μ, this minimum is shifted to larger wavevectors, i.e. estimating the vesicle radius from the first minimum of the scattering intensity will underestimate R_{0}. This shift of the minimum position, q_{v}R_{0} − π, scales as μ/R_{0}. Fig. 6(b) presents plots for ideally spherical vesicles with radius R_{0} using our curvaturedependent EDP [equation (2)] that has been parameterized from the MARTINI model simulations. The main difference from the symmetric threeGaussian model is the asymmetry of the inner and the outer electrondensity maxima of the head groups, as displayed in Fig. 3(d). Most notably, there is no discernible shift in the location of the first minimum as a function of qR_{0}, i.e. the first minimum remains an accurate estimate of the vesicle's radius even for small vesicles. Apparently, corrections to the factorization approximation and curvatureinduced asymmetry of the EDP cancel out. Using equation (22) which describes the scattering intensity of the asymmetric threeGaussian model up to second order in curvature, we can determine which asymmetry, ρ_{1} ≠ ρ_{3}, results in such a cancellation. Such a cancellation occurs if the coefficient in front of the cosine term of equation (22) vanishes, i.e.
i.e. the necessary asymmetry of the EDP is proportional to the curvature of the vesicle. Clearly, the electron density ρ_{3} of the head groups in the outer monolayer located at r = R_{0} + μ is smaller than ρ_{1} of the inner monolayer at r = R_{0} − μ, in agreement with the curvaturedependent electrondensity model parameterized by results of the MARTINI model simulations.
3.3. Distinction of thermal fluctuations and size polydispersity
Fig. 7(a) presents the scattering intensities obtained for ideally spherical vesicles with fixed radius, R_{0} = 30 nm, ideally spherical vesicles with polydisperse radii using the distribution [equation (5)] with σ_{R} = 0.1, vesicles with thermal shape fluctuations characterized by κ = 10k_{B}T and l_{max} = 6, and vesicles that have both a polydisperse radius R_{0} and thermal shape fluctuations. We observe that both size polydispersity and shape fluctuations smooth out the features of I(q), in particular the first minimum at q_{v}. In experiments, the idealization of an ensemble of monodisperse ideally spherical vesicles breaks down – due to the limitations of experimental vesicle preparation, the vesicle sizes are polydisperse, and the thermal fluctuations that give rise to shape fluctuations are a hallmark of soft matter (Székely et al., 2010; Seifert, 1997). Fig. 7(b) presents the relative effects of size polydispersity and vesicle shape fluctuations on the simulated SAXS curves. To highlight the differences, we show the ratio of the scattering intensity with either only thermal fluctuations or only size polydispersity to a reference curve including both phenomena. This representation corresponds to the case where, if it were possible in an experiment, one or other of these phenomena could be `switched off'. In the absence of polydispersity, i.e. for thermally fluctuating vesicle shapes with constant R_{0}, the ratio exhibits rapidly decaying oscillations that are only sizable in the wavevector range , i.e. switching off polydispersity in the simulation chiefly affects the smallq regime, in accordance with the expectation that vesicle size is the largest scale. In the absence of thermal shape fluctuations, i.e. for ideally spherical vesicles of varying size, we observe that the rapid oscillations are modulated by a sawtooth pattern that extends to much larger q vectors, comparable to the minima of the bilayer's form factor, . Again, this corroborates the expectation that thermal fluctuations also affect the smaller length scales.
As illustrated above, size polydispersity and thermal shape fluctuations give rise to smallq and largeq signatures in the scattering intensity. In order to gauge the ability to extract the dimensionless parameters σ_{R} and κ/k_{B}T that characterize the size distribution and the strength of thermal fluctuations, respectively, from the scattering intensity, we generated scattering data with known values of , /k_{B}T = 10 and = 30 nm, and compared the result with scattering intensities at other parameter doublets (σ_{R}, κ/k_{B}T). This illustrates how sensitively the scattering intensity depends on the two parameters and identifies correlations in their estimates. As before, the scattering intensity is computed from N_{v} = 1000 independent vesicle configurations. Fig. 8 shows a contour plot of the χ^{2} deviation [see equation (10)]. In panel (a) the comparison of the scattering intensities is extended over the lowq wavevector regime , whereas panel (b) displays the comparison for an extended interval of wavevectors, . In both cases, there appear to be no significant correlations between the estimates of σ_{R} and κ/k_{B}T, i.e. in the vicinity of the true values and /k_{B}T, the minor and major axes of the elliptical contour lines of constant χ^{2} are aligned with the σ_{R} and κ/k_{B}T axes in Fig. 8. Extending the wavevector regime to large q significantly increases the accuracy of the estimate of the bending rigidity but does not significantly affect the estimate for the relative variance σ_{R} of the vesicle size distribution. This observation corroborates the discussion of Fig. 7.
3.4. Application to experimental data
Now we use the approach presented above for leastsquares fitting of experimental vesicle SAXS data. First, we demonstrate that simulating an ensemble of vesicles and computing the
by 3D FFT followed by radial averaging can be implemented practically. Second, we test whether the effects of thermal fluctuations and/or static vesicle shape deformations are relevant, in the sense that the achievable experimental data quality suffices to distinguish these effects. Third, we compare the fitting parameters for the bilayer EDP obtained from the new approach with a conventional analytical model.To this end, we compare our method with the vesicle SAXS model of Brzustowicz & Brunger (2005), which is commonly used for SAXS analysis. This model uses a sum of N Gaussians with amplitude ρ_{i} and width σ_{i} to describe ρ_{m}. It assumes spherical symmetry for the vesicle shape and, in addition to the treatment of Kiselev et al. (2002) described in Section 3.1, also accounts for polydispersity of the vesicle sizes by a Gaussian distribution with standard deviation σ_{R}. The latter is integrated analytically, with the integration limits being set to ± ∞ – an approximation which breaks down at high relative polydispersity σ_{R} ≃ R_{0}. Furthermore, the approximation is employed to obtain the final analytical function for leastsquares fitting,
where the wavevectordependent coefficients A_{ij}(q), B_{ij}(q) and C_{ij}(q) are given by
and
For leastsquares fitting of experimental data, the present approach based on simulating fluctuating or deformed vesicles and subsequent Fourier transformation (3D FFT), is used as follows. As in the work of Brzustowicz & Brunger (2005), we use the separatedformfactor (SFF) approximation and model the EDP ρ_{m} in terms of three Gaussians. For all fits we also impose mirror symmetry of the EDP, neglecting curvatureinduced effects, since R_{0} ≫ d_{0}. Making use of the SFF approximation, we can hence factor out F(q) in the fitting function and simulate vesicles with a box profile of constant width . The simulated `δvesicle' hence captures only the vesicle shape, not its EDP, significantly reducing the simulation runs. For each κ, we simulate N = 1000 vesicles for constant (unit) radius to compute and tabulate the vesicle S(Q) in natural units Q = qR_{0}. These tabulated data are scaled in the fitting function with the fitting parameter R_{0}, using interpolation for the sampling points in both q and κ. In the same manner, if applicable, we model a static vesicle shape deformation, 〈u_{l, m}〉 ≠ 0 for l = 2, m = 0, to represent a nonspherical, oblate or prolate average vesicle shape.
Two different lipid vesicle data sets are used for these tests: (a) a (1:1) mixture of dioleoylphosphatidylcholine (DOPC) and dioleoylphosphatidylethanolamine (DOPE), formed after extrusion through 30 nm membranes in ultrapure water (MilliQ), and (b) a (1:1) mixture of DOPC and dioleoylphosphatidylserine (DOPS), which was also formed in ultrapure water but then immersed in 4 mM glucose solution to induce and vesicle deformation. Note that it is well known from phase diagrams of vesicle shapes that a surplus of vesicle surface over volume results in a transition from a spherical to a prolate shape (Seifert et al., 1991; Seifert, 1997). SAXS data were collected on the bending magnet beamline BM29 (BioSAXS) at the European Synchrotron Radiation Facility (ESRF) in Grenoble, France, at photon energy E = 12.5 keV using a multilayer monochromator with ΔE/E ≃ 0.01, and a pixel detector (Pilatus 1M, Dectris) at a sampletodetector distance of 2.867 m to cover a q range of approximately 0.036–4.95 nm^{−1}. The sample suspension was automatically loaded into a vacuummounted quartz capillary of 1.8 mm in diameter for exposure to the beam. For details of the sample preparation, experiment and data correction we refer readers to the article by Komorowski et al. (2018), where the first of the two curves was published (with a standard fitting workflow).
Fig. 9 shows the leastsquares fitting of experimental SAXS curves, by the present approach (red) as well as by the Brzustowicz & Brunger model (orange), for the extruded vesicles of (a) the equimolar DOPC:DOPE mixture and (b) the equimolar DOPC:DOPS mixture in 4 mM glucose solution. In both cases the simulation/FFT fit captures in particular the region around the first formfactor minimum better, resulting in a smaller norm of the residual (reduced χ^{2}). This effect is particularly pronounced for the DOPC:DOPE data. Here κ = 7k_{B}T provides the best fit.
For the DOPC:DOPS system, the static deformation u_{2,0} was varied to account for the expected shape transition due to Indeed, the results indicate a transition to a prolate shape (u_{2,0} > 0) after the vesicle volume decreases as a result of water permeation induced by the osmotic gradient, in line with the theoretical phase diagram (Seifert et al., 1991; Seifert, 1997). Note also that the EDPs obtained from the fits using the current model are more plausible (red versus orange profiles in the inset), i.e. they exhibit in particular a higher headgroup density and a smaller width, more similar in shape to those obtained from multilamellar planar membranes. All fitting parameters are presented in Table 4.
Thus, while we clearly see that polydispersity masks vesicle shape effects to a large extent, a residual benefit of our model is found in terms of χ^{2} (at the cost of one extra fitting parameter), and some residual sensitivity for the effects of interest remains. For preparation methods yielding smaller σ_{R}, such as purification or size fractionation with highperformance columns, and for softer bilayers, the effects would be correspondingly stronger.
4. Conclusions
In this paper we have analysed, by means of computer simulations, the curvature dependence of the electrondensity profile of small unilamellar vesicles, as well as the effects of thermal fluctuations and size polydispersity on smallangle scattering curves. Curvature changes the equilibrium bilayer structure, as captured by coarsegrained molecular dynamics simulations, and results in a thinning of the inner leaflet and a decrease in the headgroup density of the outer leaflet. Hence the EDP becomes asymmetric, even for a membrane with the same lipid composition in the inner and outer leaflets. This curvature effect becomes relevant for radii R_{0} < 30 nm, which can be encountered in experimental preparations of extruded or sonicated lipid vesicles, as well as in biological compartments such as synaptic vesicles. Importantly, the curvatureinduced asymmetry is high enough that it may be observed experimentally in future for very small vesicles. Highly resolved experimental bilayer profiles could thus, in principle, also provide information on the interplay of asymmetric lipid partitioning and curvature. To this end, we have also investigated how the smallangle scattering distribution changes for small R_{0} ≲ 10 nm when the common factorization approximation breaks down. In this regime, the scattering function can no longer be modelled as the product of a powderaveraged bilayer form factor and the transform of the vesicle shape, for example a thin hollow sphere. We find that for a symmetric bilayer the first minimum of the scattering function is positively shifted with respect to the value q_{v} = π/R_{0}, which underestimates the vesicle radius for small vesicles. Interestingly, this shift is eliminated for curvatureinduced bilayer asymmetry. Hence, the curvaturedependent electrondensity model fixes the position of the first minimum of the scattering intensity at q_{v} = π/R_{0}, despite the fact that the factorization approximation is invalid. In summary, the curvature corrections become important for very small vesicles, resulting, for example, in a 20% density difference between inner and outer leaflets for R_{0} ≲ 10 nm. In contrast, for vesicles with R_{0} ≳ 30 nm the curvature has very little effect and, at the same time, the factorization approximation becomes valid.
In addition to curvature, we have studied the effects of size polydispersity and thermal fluctuations, which both result in a smearing out of the scattering curve minima. However, the exact functional forms differ. In fact, polydispersity and thermal fluctuations modify the scattering curve in different q regimes: thermal fluctuations predominantly affect the scattering intensity around the membrane formfactor minima, corresponding to the membrane structure, while the vesicle size polydispersity mainly contributes in the smallq range, reflecting the vesicle size. Hence, our results show that thermal fluctuations are not completely masked by polydispersity, which is almost always unavoidable in experiments. At the same time, this study also explains when the classical vesicle models that assume perfect spherical symmetry work relatively well. However, if further progress is made in purification or size fractionation, more details of the shape fluctuations will become visible and should be accounted for. This also holds for small changes in the average vesicle shape, which are sensitive indicators for changes in vesicle volume and membrane area.
More generally, we have presented an approach which helps to free experimental investigations of vesicles from idealizing assumptions, by directly modelling and simulating vesicle structures for the relevant parameters, and calculating the scattering function based on 3D FFT on a suitable numerical grid. Since this approach is computationally efficient, it is also suitable for the analysis of experimental SAXS data. In future, it could easily be augmented to accommodate either more complicated configurations, including adhering vesicles, or vesicles with
and large membrane proteins. In particular, it could be used to study the structure of synaptic vesicles, including the effects of protein clusters and/or shape transitions. All of the abovementioned systems can be modelled on a numerical grid, while the 3D FFT and radial averaging provide the link to the experimental data.Acknowledgements
We thank Kilian Frank for contributions to the initial work of this project, including code development and many stimulating discussions. Computing time was provided by the NIC Jülich, HLRN Göttingen/Berlin and the GWDG Göttingen. Open access funding enabled and organized by Projekt DEAL.
Funding information
Financial support by CRC 1286 `Quantitative Synaptology', funded by Deutsche Forschungsgemeinschaft, and CRC 803/TP B03 is gratefully acknowledged.
References
Abraham, M. J., Murtolad, T., Schulz, R., Páll, S., Smith, J. C., Hess, B. & Lindahl, E. (2015). SoftwareX, 1–2, 19–25. Google Scholar
Brzustowicz, M. R. & Brunger, A. T. (2005). J. Appl. Cryst. 38, 126–131. Web of Science CrossRef CAS IUCr Journals Google Scholar
Canham, P. B. (1970). J. Theor. Biol. 26, 61–81. CAS PubMed Google Scholar
Discher, D. E., Mohandas, N. & Evans, E. A. (1994). Science, 266, 1032–1035. CAS PubMed Google Scholar
Gompper, G. & Kroll, D. M. (1996). J. Phys. I, 6, 1305–1320. Google Scholar
Gompper, G. & Kroll, D. M. (1997). J. Phys. Condens. Matter, 9, 8795. Google Scholar
Heberle, F. A., Pan, J., Standaert, R. F., Drazba, P., Kučerka, N. & Katsaras, J. (2012). Eur. Biophys. J. 41, 875–890. Web of Science CrossRef CAS PubMed Google Scholar
Helfrich, W. (1973). Z. Naturforsch. Teil C, 28, 693–703. CAS Google Scholar
Helfrich, W. (1986). J. Phys. Fr. 47, 321–329. Google Scholar
Kiselev, M. A., Lesieur, P., Kisselev, A. M., Lombardo, D. & Aksenov, V. L. (2002). Appl. Phys. Mater. Sci. Process. 74, s1654–s1656. Web of Science CrossRef CAS Google Scholar
Kiselev, M. A., Zemlyanaya, E. V., Aswal, V. K. & Neubert, R. H. H. (2006). Eur. Biophys. J. 35, 477–493. PubMed CAS Google Scholar
Komorowski, K., Salditt, A., Xu, Y., Yavuz, H., Brennich, M., Jahn, R. & Salditt, T. (2018). Biophys. J. 114, 1908–1920. CAS PubMed Google Scholar
Kučerka, N., Nagle, J. F., Sachs, J. N., Feller, S. E., Pencer, J., Jackson, A. & Katsaras, J. (2008). Biophys. J. 95, 2356–2367. Web of Science PubMed Google Scholar
Li, J., Dao, M., Lim, C. T. & Suresh, S. (2005). Biophys. J. 88, 3707–3719. PubMed CAS Google Scholar
Lim, H. W. G., Wortis, M. & Mukhopadhyay, R. (2002). Proc. Natl Acad. Sci. USA, 99, 16766–16769. PubMed Google Scholar
Marrink, S. J., de Vries, A. H. & Tieleman, D. P. (2009). Biochim. Biophys. Acta, 1788, 149–168. PubMed CAS Google Scholar
Marrink, S. J., Risselada, J., Yefimov, S., Tieleman, P. & de Vries, A. H. (2007). J. Phys. Chem. B, 111, 7812–7824. PubMed CAS Google Scholar
Milner, S. T. & Safran, S. A. (1987). Phys. Rev. A, 36, 4371–4379. CAS Google Scholar
Müller, K., Katsov, K. & Schick, M. (2003). J. Polym. Sci. B Polym. Phys. 41, 1441–1450. Google Scholar
Müller, M., Katsov, K. & Schick, M. (2006). Phys. Rep. 434, 113–176. Google Scholar
Pabst, G., Rappolt, M., Amenitsch, H. & Laggner, P. (2000). Phys. Rev. E, 62, 4000–4009. Web of Science CrossRef CAS Google Scholar
Pencer, J., Krueger, S., Adams, C. P. & Katsaras, J. (2006). J. Appl. Cryst. 39, 293–303. Web of Science CrossRef CAS IUCr Journals Google Scholar
Risselada, H. J., Bubnis, G. & Grubmüller, H. (2014). Proc. Natl Acad. Sci. USA, 111, 11043–11048. CAS PubMed Google Scholar
Risselada, H. J., Mark, A. E. & Marrink, S. J. (2008). J. Phys. Chem. B, 112, 7438–7447. PubMed Google Scholar
Safran, S. A. (1994). Statistical Thermodynamics of Surfaces, Interfaces and Membranes. Reading: Addison Wesley. Google Scholar
Seifert, U. (1997). Adv. Phys. 46, 13–137. CAS Google Scholar
Seifert, U., Berndl, K. & Lipowsky, R. (1991). Phys. Rev. A, 44, 1182–1202. CAS PubMed Google Scholar
Semeraro, E. F., Marx, L., Frewein, M. P. K. & Pabst, G. (2021). Soft Matter, 17, 222–232. CAS PubMed Google Scholar
Székely, P., Ginsburg, A., BenNun, T. & Raviv, U. (2010). Langmuir, 26, 13110–13129. Web of Science PubMed Google Scholar
Wang, Y., Gkeka, P., Fuchs, J. E. R., Liedl, K. & Cournia, Z. (2016). Biochim. Biophys. Acta, 1858, 2846–2857. CAS PubMed Google Scholar
Zubay, G. (1999). Biochemistry, 4th ed. Dubuque: Brown Publishers. Google Scholar
This is an openaccess article distributed under the terms of the Creative Commons Attribution (CCBY) Licence, which permits unrestricted use, distribution, and reproduction in any medium, provided the original authors and source are cited.