The effect of polydispersity, shape fluctuations and curvature on small unilamellar vesicle small-angle X-ray scattering curves

It is a challenge to distinguish the effect of shape fluctuations and size polydispersity on experimental small-angle X-ray scattering curves of small unilamellar vesicles. Here it is shown that both effects have distinguishable spectral patterns, and an efficient simulation tool is presented for simulating and analysing experimental data. The importance of curvature-induced electron-density profile asymmetry for estimating the vesicle size from SAXS scattering curves is also demonstrated.


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 small-angle X-ray 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 curvature-dependent electron-density profile (EDP) has remained challenging.
Theoretical methods, such as molecular dynamics (MD) simulations (Marrink 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 coarse-grained MD simulations and the elastic Helfrich (1973Helfrich ( , 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 curvature-dependent radial EDP of a lipid membrane from MD simulations using the coarse-grained 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 three-dimensional electron-density map is used to calculate the scattering intensity via a three-dimensional (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 least-squares 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 (Milli-Q), and (ii) a (1:1) mixture of dioleoylphosphatidylcholine and dioleoylphosphatidylserine.
The simulation and FFT approach allows us to analyse the interplay between size polydispersity, vesicle shape fluctuations and curvature-dependent 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 structure factor of the vesicle shape (Kiselev et al., 2002(Kiselev et al., , 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 small-angle neutron scattering studies of vesicles, after minor modification regarding the scattering lengths.

Methods
In the following we detail our hierarchical modelling approach, illustrated in Fig. 1. The system has two characteristic length 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 coarse-grained MARTINI model (Marrink et al., 2007). These simulations provide detailed density profiles across the bilayer, from which we extract the electron-density 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 curvature-dependent 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 curvature-dependent EDP, we combine these two descriptions to calculate SAXS scattering curves.
2.1. MARTINI model simulation and EDPs 2.1.1. Simulation protocol. Coarse-grained 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 1-palmitoyl-2-oleoylsn-glycero-3-phosphocholine (POPC) vesicle was formed by spontaneous aggregation, following the protocol of Risselada et al. (2008Risselada et al. ( , 2014. The vesicle is composed of 1447 lipids in the outer leaflet and 770 lipids in the inner leaflet, embedded in a solvent that contains a total of 97 217 coarse-grained 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 lipids. All systems were equilibrated in the NPT ensemble and simulated for 1 ms in the NVT ensemble at T = 300 K to calculate EDPs. The simulations provide detailed bead-type profiles across the planar or highly curved bilayer, as presented in Fig. 2. 2.1.2. Electron-density profiles. From these bead-type profiles we extracted the EDPs of the planar membrane and the highly curved vesicle. The numbers of electrons per coarse-grained 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 coarse-grained bead. Fig. 3 summarizes the results from the coarse-grained 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 fourth-order 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 head-group region of the bilayer membrane. l and r are the peak widths, and l and r characterize the peak heights. The polynomial fits the electron-density 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 least-squares 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

Figure 2
Plots of the number density of lipid heads (NC 3 and PO 4 ), tails (C1A-C4A and C1B-C5B) and water (W) versus displacement from the bilayer's midplane for a planar membrane (solid lines) and highly curved vesicle (dotted lines). The schematic of a curved membrane at the top illustrates how lipid packing responds to curvature: lipid heads (smaller circles) in the inner leaflet are tightly packed, whereas lipid heads in the outer leaflet have a larger hydration shell (larger circles). 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 small-curvature expansion that underlies the Helfrich Hamiltonian, i.e. both the large-scale description of the vesicle shape and its fluctuations and the response of the small-scale 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 (anti-symmetric) 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 curvature-induced thinning of the membrane. Since these effects do not depend on the sign of the curvature, they scale like 1=R 2 0 to leading order. If the symmetric contribution, Á s (d), chiefly stems from the curvature-induced thinning of the bilayer membrane, we can approximately relate its functional form to the profile of a planar bilayer by assuming that  (a) Electron-density contrast of a planar bilayer patch (orange dashed line) and an 8 nm vesicle (blue solid line) of POPC lipids obtained by MD simulations, and fitted curves using equation (1) (dashed-dotted lines). The difference between vesicle and planar membrane electron densities is shown by the green dotted line. (b) Anti-symmetric, d 0 /R 0 Á a (d) (red dashed line), and symmetric, (d 0 /R 0 ) 2 Á s (d) (violet solid line), responses of the electron density to curvature (green dotted line). (c) Electron-density model of a vesicle of arbitrary radius R 0 > 8 nm, as given by equation (2). (d) Electrondensity contrast of vesicles with R 0 = 8 nm and R 0 = 30 nm obtained with the curvature-dependent model, equation (2) (blue and red dashed lines, respectively), and the curvature-independent three-Gaussian model, equation (4) (blue and red solid lines, respectively). the profile of a bilayer of thickness d 0 obeys the affine scaling relation m ðdÞ = ðd dÞ with the dimensionlessd d = 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 Àd 0 ' 1=R 2 0 and the spatial dependence of the symmetric contribution is Á s ðdÞ ' dðd m =ddÞ. 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 dðd m =ddÞ, 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 curvature-induced 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 head-group 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 k th peak, respectively. The electron density of this three-Gaussian model with symmetric inner and outer maxima is compared in Fig. 3(d) with the curvature-dependent 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 curvature-dependent model and the symmetric three-Gaussian model provide a similar representation of the electron density for large radii, R 0 > $ 30 nm. Adjusting the parameters of the symmetric three-Gaussian model to our curvature-dependent 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.

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 curvature-dependent electron-density 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 R characterizes the mean vesicle radius of the ensemble and R R 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 free-energy 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, Electron-density contrast for the symmetric function, Á s , and the function d(d m /dd), plotted versus the rescaled parameterd d ¼ d=d 0 , with d 0 = 4 nm for the vesicle and 4.4 nm for the flat membrane. The discretization of the EDP used in the code is shown as points. Table 3 Parameters of the symmetric three-Gaussian model [equation (4)], shown in Fig. 3(d). 1 = 3 , = À 1 = 3 , 1 = 3 and 2 = 0.
with Y l;m ð; Þ ¼ P m l ðcos Þ exp ðimÞ, where P m l ðcos Þ 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 2R 0 =l max d 0 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 particle-based 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 curvature-dependent EDPs. The electron density in 3D space is obtained by where v (d, R 0 ) is the curvature-dependent 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 = 8R. Thus the collocation grid can resolve spatial distances =d 0 = 8R=ðNd 0 Þ. For small vesicles, with R = 8, 10 and 20 nm, we use N = 400, resulting in /d 0 = 0.1 for R = 20 nm. The same resolution is obtained for the largest vesicle, R = 30 nm, for N = 600. There are two characteristic length scales, R and d 0 , and we choose the dimensionless wavevector qR 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 hÁ Á Ái 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 Ivy-Bridge Intel E5-2670 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 mean-squared variance: where N Q denotes the number of q values in the considered interval.

Structure factor and form factor
To understand qualitatively the features of the scattering intensity I(q), we factorize I(q) into a structure factor, 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 AðqÞ ¼ BðqÞ Using this expression we obtain for the scattering intensity I(q) up to first order in curvature dd Á a ðdÞ sinðqdÞ: The leading term is the popular factorization approximation, where I TS is the structure factor of an infinitely thin spherical shell and A 2 0 ðqÞ is the powder average of the planar bilayer form factor. Equation (18) indicates under which conditions the factorization approximation is accurate.
The structure factor of a spherical shell and the powderaveraged form factor of a planar bilayer are depicted in Fig. 5 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 structure factor 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 first-order correction term I 1 ðqÞ=I TS ðqÞ = 2 cotðqR 0 Þ ðd 0 =R 0 Þ A 0 ðqÞB 1 ðqÞ appear to be negligibly small, even for d 0 =R 0 ¼ 0:13.

Models of electron-density 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 head-group and tail regions from the centre of mass of the vesicle. We consider an asymmetric EDP, where the head-group peaks differ in their amplitudes, 1 6 ¼ 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, k = k ( 1 and k =R 0 ( 1, Scattering intensity of an infinitely thin spherical shell of R 0 = 30 nm (green dotted line), a planar membrane bilayer, R 0 = 1, of thickness d 0 = 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ðqÞ $ sin 2 ðqR 0 Þ, 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 three-Gaussian 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 curvature-dependent EDP [equation (2)] that has been parameterized from the MARTINI model simulations. The main difference from the symmetric three-Gaussian model is the asymmetry of the inner and the outer electron-density 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 curvature-induced asymmetry of the EDP cancel out. Using equation (22) which describes the scattering intensity of the asymmetric three-Gaussian model up to second order in curvature, we can determine which asymmetry, 1 6 ¼ 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 curvature-dependent electron-density model parameterized by results of the MARTINI model simulations. 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 IðqRÞ=I both ðqRÞ of the scattering intensity with either only thermal fluctuations or only size polydispersity to a reference curve including both Scattering intensity I as a function of the dimensionless wavevector qR 0 obtained (a) from the symmetric three-Gaussian EDP with 1 = 3 [cf. equation (4) and Table 3] and (b) with the curvature-dependent EDP, equation (2), parameterized from the MARTINI model simulation (cf. Fig. 3).

Distinction of thermal fluctuations and size polydispersity
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 < $ qR < $ 40, i.e. switching off polydispersity in the simulation chiefly affects the small-q regime, in accordance with the expectation that vesicle size is the largest characteristic length 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, qR ' 64 ' 2R=d 0 . 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 small-q and large-q 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 Ã R ¼ 0:1, Ã /k B T = 10 and R = 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     Isocontour plots of 2 of simulated intensities for different polydispersities R and bending rigidities /(k B T), calculated using equation (10) with (a) 0 qR < 12:2 and (b) 0 qR < 110. the 2 deviation [see equation (10)]. In panel (a) the comparison of the scattering intensities is extended over the low-q wavevector regime 0 qR < 12:2, whereas panel (b) displays the comparison for an extended interval of wavevectors, qR < 110. 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 Ã R 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.

Application to experimental data
Now we use the approach presented above for least-squares fitting of experimental vesicle SAXS data. First, we demonstrate that simulating an ensemble of vesicles and computing the structure factor 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 AE 1 -an approximation which breaks down at high relative polydispersity R ' R 0 . Furthermore, the approximation 2 R q cosðqR 0 Þ ( R sinðqR 0 Þ is employed to obtain the final analytical function for least-squares fitting, where the wavevector-dependent coefficients A ij (q), B ij (q) and C ij (q) are given by and For least-squares 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 separated-form-factor (SFF) approximation and model the EDP m in terms of three Gaussians. For all fits we also impose mirror symmetry of the EDP, neglecting curvature-induced 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 d sim ( d 0 . 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 structure factor 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, hu l, m i 6 ¼ 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 ultra-pure water (Milli-Q), 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 osmotic pressure 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 sample-to-detector 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 vacuum-mounted 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 least-squares 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 form-factor 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 research papers osmotic pressure. 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 head-group 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 high-performance liquid chromatography columns, and for softer bilayers, the effects would be correspondingly stronger.

Conclusions
In this paper we have analysed, by means of computer simulations, the curvature dependence of the electron-density profile of small unilamellar vesicles, as well as the effects of thermal fluctuations and size polydispersity on small-angle scattering curves. Curvature changes the equilibrium bilayer structure, as captured by coarse-grained molecular dynamics simulations, and results in a thinning of the inner leaflet and a decrease in the head-group 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 small-angle 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 powder-averaged 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 curvature-induced bilayer asymmetry. Hence, the curvature-dependent electron-density model fixes the position of the first minimum of the scattering  (a) SAXS data for 30 nm extruded DOPC:DOPE (1:1) vesicles in ultrapure water (black circles), with least-squares fits (solid lines) to the analytical model according to Brzustowicz & Brunger (2005) (orange) and to simulated and tabulated scattering curves computed for fluctuating vesicles according to the present method (red). The inset shows the respective EDPs, in both cases modelled by three Gaussians (head group, tail region, head group). (b) SAXS data for DOPC:DOPS (1:1) vesicles suspended in 4 mM glucose solution (black circles) to induce a vesicle shape deformation by osmotic pressure. Again, the fit to the Brzustowicz & Brunger model (orange) is compared with the present approach (red), and the corresponding EDPs are shown in the inset. The fitting parameters are summarized in Table 4. Table 4 Fitting parameters obtained by least-squares fitting of the SAXS data for 30 nm extruded DOPC:DOPE (1:1) vesicles and for 100 nm extruded DOPC:DOPS (1:1) vesicles + 4 mM glucose (Fig. 9).
The structural bilayer parameters are h = h1 = h2 and h = h1 = h2 for a symmetric bilayer profile. The amplitude of the Gaussian representing the chain region is set to c = À1, i.e. the electron-density difference from water is fitted (shown in Fig. 9, insets) and normalized to the electron-density difference between water and the bilayer centre plane. u 2,0 is the static deformation. 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 form-factor minima, corresponding to the membrane structure, while the vesicle size polydispersity mainly contributes in the small-q 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 lipids 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 above-mentioned systems can be modelled on a numerical grid, while the 3D FFT and radial averaging provide the link to the experimental data.