Fluctuation spectroscopy of giant unilamellar vesicles using confocal and phase contrast microscopy

A widely used method to measure the bending rigidity of bilayer membranes is fluctuation spectroscopy, which analyses the thermally-driven membrane undulations of giant unilamellar vesicles recorded with either phase-contrast or confocal microscopy. Here, we analyze the fluctuations of the same vesicle using both techniques and obtain consistent values for the bending modulus. We discuss the factors that may lead to discrepancies.

Bending rigidity of cellular membranes plays a key role in membrane remodeling. Knowledge of its value is needed to quantify processes that involve curvature changes such as budding (as in endo-and exocytosis), tubuation and fusion. Various experimental methods have been devised to measure bending rigidity [1], e.g. micropipette aspiration [2][3][4], electrodeformation [5][6][7][8], optical tweezers [9,10], and scattering based techniques [11,12]. One of the most popular methods is the fluctuation spectroscopy, pioneered by Brochard and Lenon [13], due to its ease of implementation [7,[14][15][16]. In essence, a time series of vesicle contours in the focal plane (the equator of the quasi-spherical vesicle) is recorded. The quasi-circular contour is decomposed in Fourier modes. The fluctuating amplitudes have variance dependent on the membrane bending rigidity and tension. Imaging is most commonly done by phase contrast microscopy [1,7,15,[17][18][19][20][21][22] but other methods such as confocal [23][24][25] and light sheet microscopy [26] have also been employed. The increased variety of imaging methods raises the question if they all yield the same results.
Recently, Rautu et al. [25] pointed out that in phase contrast imaging, projections of out-of-focus fluctuations may contribute to the contour statistics leading to systematic overestimation of the bending rigidity value when compared to other methods such as micropipette aspiration and X-ray scattering. Their analysis claimed to resolve the discrepancies and confocal microscopy was suggested as a better imaging technique due to the precise control over the focal depth. However, comparing bending rigidity numbers obtained by different imaging techniques is only meaningful if the same system is probed. It is known that many factors such as sugars (and gravity), salt, buffers, solution asymmetry, concentration of fluorescent lipids, preparation method or type of bilayer configuration (stacked or free-floating), influence the measured mechanical properties of bilayer membranes [1,18,[27][28][29][30][31]. For example, even measurements with the same method can give wide range of values, e.g., the bending rigidity of a DOPC bilayer measured with flickering spectroscopy has been reported form 15 k B T [24] to 27 k B T [32], where k B T is the thermal energy; see Table 1 in the Supporting Information (SI) for a comprehensive list.
In order to compare imaging with phase contrast and confocal microscopy, we measure the bending rigidity of the same giant vesicle with both techniques. We highlight some important issues to be considered to ensure reliable measurements. We also show that results obtained with both methods are consistent.

Equilibrium fluctuations of a quasi-spherical vesicle
First, we summarize the theoretical basis of the fluctuations analysis and the experimental methods (details are provided in SI section S5). We also correct published expressions for the relaxation frequency and cross-spectral density of the shape fluctuations.
The contour in the equatorial plane of a quasispherical vesicle is decomposed in Fourier modes, r(φ, t) = R 0 1 + qmax q=−qmax u q (t) exp(iqφ) , where R 0 = (3V /4π) 1/3 is the radius of an equivalent sphere with the volume V of the GUV. In practice, q max is the maximum number of experimentally resolved modes. The statistical analysis of the fluctuating amplitudes u q yields the values of membrane bending rigidity κ and the tension σ since |u q | 2 ∼ k B T /κ q 3 +σq , whereσ = σR 2 0 /κ. More precisely, the statistics of the two-dimensional circular modes, u q , is derived from the three-dimensional shape modes, f lm , which describe the nearly-spherical shape in terms of spherical harmonics [14,33], Here, l max is an upper cutoff, in the order of the ratio of the GUV radius and bilayer thickness. The contour in the focal plane corresponds to the equator of the quasispherical vesicle, θ = π/2, i.e., r(φ, t) = R( π 2 , φ, t), which leads to the following expression for the mean squared amplitudes (1) where n lq = (2l + 1)(l − q)!/4π(l + q)! and P lq are the associated Legendre polynomials. The short-wavelength shape fluctuations are dominated by the bending rigidity, while the long wavelengths are controlled by tension; the crossover occurs around mode √σ .
To validate our methodology, we have simulated the thermal shape fluctuations of a GUV, see also SI section S5. We have generated a sequence of three-dimensional shapes (and their corresponding equatorial contours) using the evolution equations [33,34] where ζ lm is the thermal noise driving the membrane undulations, η in and η ex are the viscosity of the solution inside and outside the vesicle. Note that the relaxation time given by Eq. 2 in Rautu et al. [25] has incorrect dependence on the viscosities of the enclosed and suspending solutions. The simulated contours were analyzed by our code and the extracted bending rigidity and tension were compared to the input values to confirm accuracy of the contour detection, Fourier decomposition and data fitting algorithms. The time evolution of the modes also enables us to access information provided by the time correlations u q (0)u * q (t) = |u q | 2 exp(−t/τ q ) . If q 1, the correlation time tends to that of a planar membrane τ −1 q = κ q 3 +σq /2R 3 0 (η ex + η in ). If the cross-spectral density |u q (0)||u q (t)| is utilized, the correct time dependence in the exponential includes a factor of 2 and Eq. 3 in Zhao et al. [35] needs to be corrected (see SI section S5, Eqn. 39). Figure 1 shows a typical fluctuations spectrum, given by Eq.(1), fitted to the experimental data for the same vesicle imaged with confocal and phase contrast microscopy using a 40x objective with 0.6 numerical aperture (NA), pinhole size of 1 Airy unit (AU) and polarization correction (see below and SI section S3). The contour was detected with sub-pixel resolution [7]. The experimental data was fitted with Eq.(1) with Levenberg-Marquardt algorithm and yielded bending rigidity κ = (23.9 ± 1.6) k B T and tension σ = 5.1 ± 1.4 × 10 −9 N/m and κ = (22.3 ± 2.1)k B T and σ = 3.1 ± 1.2 × 10 −9 N/m from the confocal and phase contrast microscopy data. Fluctuation spectrum of the same DOPC vesicle (shown in the inset) obtained with confocal and phase contrast microscopy with a 40x/NA 0.6 objective, pinhole size of 1 AU and polarization correction. The dye concentration is 0.2 mol%. Scale bar is 15 µm. The vertical lines denote the cutoff resolution for the modes: optical resolution (solid line), phase contrast (blue dashed) and confocal (red dashed line). The crossover mode √σ is 7.
Bending rigidity obtained from confocal and phase-contrast microcopy: effect of resolution and vesicle size Giant unilamellar vesicles (GUV) were electroformed from DOPC (99.8 mol% dioleoylphosphatidylcholine and 0.2 mol% Texas-Red 1,2-hexadecanoyl-wn-glycero-3-phosphoethanolamine, TR-DHPE) in 20 mM sucrose and subsequently diluted in 22 mM glucose, see SI section S2 for details. Low sugar concentration was used in order to minimize the effects of gravity [29] and effect of sugars [27], but still allow the vesicles to settle to the chamber bottom for easier recording. Low dye content minimizes effects of fluorophores [30]. By imaging a population of 18 vesicles with both methods, the bending rigidity obtained are 22.5±2.0 k B T with confocal and 23.3±1.6 k B T phase contrast microscopy; each vesicle was analyzed with both imaging techniques as in Fig. 1 and then the results were averaged over the population. Figure 2 shows the box and whisker plot for more detailed statistics. Based on F statistics and ANOVA (analysis of variance) test , the probability obtained is p = 0.48 for null hypothesis testing. This indicates no significant difference between the two imaging techniques.
Since only modes with wavenumber q > σR 2 0 /κ are sensitive to the bending rigidity, it is desirable to have more resolved modes, i.e., modes with amplitude greater than optical resolution limit ≈ 250 nm [36]. The average mean fluctuation amplitude scales with the size of the vesicle u q ∼ R 0 k B T /κ, hence larger GUVs admit more spatially resolvable fluctuation modes as shown in Fig. 3. However, even for the same vesicle we find that Imaging with phase contrast and confocal microscopy for objectives of the same numerical aperture (NA) give consistent results. Box and whisker plot comparison for a DOPC vesicle population where each vesicle was analyzed with phase contrast and confocal imaging with 40x objectives with NA 0.6 and NA 1.3. Pinhole size is 1 AU with polarization correction. The dye concentration is 0.2 mol%.
the number of resolved modes is higher for phase contrast than for confocal imaging. Indeed, Fig. 1 shows that the noise level is higher for confocal microscopy, and on average phase contrast imaging resolves 8-10 modes more than confocal one does. The poorer mode resolution with confocal microscopy is likely due to poor contour recognition. The reasons for this are discussed in the next section.

FIG. 3.
Larger vesicles allow resolving more fluctuation modes thus yielding more reliable determination of the bending rigidity. Data are collected on DOPC vesicles with different sizes. The dye concentration is 0.2 mol%. Regression lines are added to guide the eyes. Imaging was done with 40x objectives with different numerical aperture (NA), pinhole size of 1 AU and polarization correction.
We found that the vesicle population needs to have broad range of radii to avoid the size bias in the bending rigidity values we discovered for confocal microscopy. In this case with 40x/NA 0.6 (air) objective the mean Pearson correlation and standard deviation coefficient is 0.65±0.21 (see SI section S4 for the histograms generated with bootstrapping resampling technique). This indicates the bending rigidity obtained from confocal microscopy can be systematically underestimated if the vesicle population contains similar sized vesicles. For example, analysis of vesicles with radius around 10 µm yields κ smaller by roughly 6 k b T than the values obtained for vesicles with radius of 30 µm. The size dependence is insignificant for phase contrast microscopy with a mean correlation coefficient of 0.28±0.18 with 40x/NA 0.6 (air). Analyzing the same vesicle population with 40x/NA 1.3 objective in phase contrast and confocal imaging yields 21.0±2.0 k B T and 21.7±2.0 k B T respectively. Higher numerical aperture in phase contrast leads to negligible correlation coefficient of -0.07±0.34 between bending rigidity and vesicle size and decrease in the correlation coefficient to 0.43±0.14 for confocal imaging with 40x/NA 1.3 objective. The bias originates from out of plane fluorescence which worsens the contour detection. This issue is investigated in the next section.

Out-of-focus fluorescence affects contour detection quality in confocal microscopy
The vesicle contour is detected from radial intensity line profiles, see SI section S2. In confocal cross sections, weak fluorescence from the vesicle membrane located above and below the focal plane may result in signal projected in the interior of the vesicle image which is higher compared to the surrounding background. The resulting asymmetry in the intensity line profile (Fig.4a) leads to an artificial contour displacement, i.e., poor contour detection (note that such an asymmetry is absent in images acquired with phase contrast microscopy). This asymmetry creates a systematic error shifting the vesicle contour by 0.53 µm. The error is larger than the pixel resolution of the system, 0.252 µm, hence the higher modes are averaged out. Smaller vesicles or larger pinholes lead to higher signal inside the vesicle (see inset in Fig. 4b) corresponding to greater asymmetry which increases the error from contour fitting and introduces dependence of the bending rigidity on vesicle size. For imaging with higher numerical aperture objectives (e.g. NA 1.3), the asymmetry in the intensity line profiles is suppressed and contour detection is correct. Note that phase contrast images do not suffer from the asymmetry-induced error irrespective of the objective NA.
We investigated the impact of out-of-focus fluorescence on the fluctuations statistics by varying the pinhole size for confocal imaging on the same vesicle. The standard pinhole size in confocal microscopy is defaulted to 1 Airy unit (AU) (full width at half maximum FWHM=1.6 µm) for 40x/NA 0.6 objective. We analyzed the same vesicles with different optical sectioning at 0.3 AU (FWHM=0.9 µm), and 2 AU (FWHM=2.9 µm). The mean bending rigidity did not show significant differences based on ANOVA testing and post hoc Dunnett test, however the error increases with the pinhole size. The sensitivity to the vesicle size also becomes more pronounced with

Dye related artifacts: vesicle tubulation and polarization
Confocal imaging relies on fluorophores added to the membrane, and some studies have used up to 10 mol% dye [23,24]. To probe the effect of fluorophore on κ, we changed the dye concentration from 0.2 mol% to 2 mol% TR-DHPE. The bending rigidity of this population of vesicles showed non significant difference with κ=20.09±2.49 k B T with one ANOVA testing. However, it was observed that over 2-3 min of recording, around 50 % of the vesicles developed inward structures such as buds or visible tubes as shown in Fig. 5. Vesicles with such defects displayed significantly higher bending rigidity, 25.01 ± 2.11 k B T .
TR-DHPE belongs to a family of polarity-sensitive fluorescent probes. As a result, the signal intensities are different at the pole and equator of the vesicle (see SI section S3). This may lead to errors in the contour detection in these regions. The polarization effect was corrected by using circular rotation plates to have even intensities across the equatorial vesicle plane. The analysis of the same vesicle with and without the polarization correction showed a 3 k B T lower bending rigidity without any correction with 40x/NA 0.6 (air) objective. This softening effect became insignificant with 40x/NA 1.3 (oil) objective (SI section S3). This is likely due to loss of signal at low intensity regions where the higher mode fluctuations intensities are averaged out with background noise due to out-of-focus fluorescence.

Effect of nearby vesicles on fluctuation spectra
The equilibrium shape fluctuations of an isolated GUV are driven by Gaussian thermal noise. Defects such as buds, nanotubes, invaginations or docked LUVs modify the vesicle fluctuations [7] and their effect can be detected in the statistics at each point on the vesicle contour profile using the ensemble-averaged probability density function (PDF) as shown in Fig. 6a In addition to defects attached the membrane, we also found that hydrodynamic flows from nearby vesicles can amplify vesicle fluctuations.
We characterized the Gaussianity of the fluctuations using the fourth PDF moment, Kurtosis, K. For a Gaussian distribution, K = 3. In Fig. 6 we demonstrate how thermal fluctuations maybe modified (see supplementary video). As shown, majority of the contours is characterized by normal distribution. However near other flickering structures, the fluctuation map density is modified. The non-Gaussian enhanced fluctuations are observed with leptokurtic nature (K > 3). This observation serves as a caution to filter out vesicles with sub-optical structures affecting the fluctuations.
For all the experiments above, we find a deviation of ≈ 2 k B T if the experiments are repeated for the same vesicle conditions on a different day. This points out the sensitive nature of the experimental procedure that can slightly vary from day to day.

Conclusions
We compare the bending rigidity of bilayer membranes determined from flickering spectroscopy of GUVs imaged with confocal and phase contrast microscopy. Examining the same vesicle with both imaging techniques shows no significant differences in the bending modulus obtained from the two methods, in contrast to the overestimation reported by Rautu et al [25] when phase contrast microscopy is used. Our analysis indicates that membrane defects such as buds and tubes induced by long exposure to laser in confocal microscopy can significantly stiffen the membrane. Furthermore, we find that errors in contour detection that could impact data interpretation can arise from fluorescence signal "pollution" and dye polarization. The bending rigidity we obtain (∼ 22k B T for DOPC) is in line with the values obtained with other techniques such as micropipette aspiration, X -ray scattering, electrodeformation and neutron spin echo [3,9,11]. Exploring the effect of various parameters, we find that optimal imaging conditions for bending rigidity measurements from confocal imaging include high magnification objective, high numerical aperture, circular polarization correction, minimum dye concentration, small pinhole size, and broad vesicle size distribution.
In conclusion, we demonstrate that phase contrast and confocal microscopy produce the same results if precautions are taken to minimize effects of the dye and improve contour detection. Our study suggests that the many published results obtained by phase contrast microscopy are likely to be unaffected by the projections of out-offocus fluctuations onto the imaging plane in contrast to the claim by Rautu et al [25]. Since dye related artifacts such as laser-induced defects can compromise the data, it is advantageous to use phase contrast imaging as it does not require dyes.

ACKNOWLEDGEMENTS
This research was funded in part by NSF-CMMI awards 1748049 and 1740011. PV acknowledges support from the Alexander von Humboldt Foundation and HAF acknowledges financial support from Prof. Reinhard Lipowsky for visits to the Max Planck Institute of Colloids and Interfaces,

S1. BENDING RIGIDITY VALUE OF DOPC BILAYERS
The bending rigidity values of bilayer membranes made of the same lipid can vary across studies due to different conditions, e.g., sugars, salt, buffers, dye concentration, as well as the preparation method [37]. Table I illustrates the wide range of reported values of the bending rigidity values DOPC bilayers. Refer to Table II for the bending rigidity values obtained in this study for different microscopy setting.

Vesicle preparation
Giant unilamellar vesicles (GUVs) were prepared using the classical electroformation method [51] from DOPC and the fluorescent lipid Texas Red 1,2-hexadecanoyl-sn-glycero-3-phosphoethanolamine (TR-DHPE). The composition of the GUVs explored are 99.8 % DOPC 0.2 % TR-DHPE and 98 % DOPC 2 % TR-DHPE (mole fractions). Stock solutions of DOPC and TR-DHPE at 10 mg/ml and 1 mg/ml in chloroform were diluted to a final concentration of 4 mM for varying proportions. A small volume, 10 µl, of the solution was spread on the conductive surface of two glass slides coated with indium tin oxide (ITO) (Delta Technologies). The glass slides were then stored under a vacuum for 1-2 hours to remove traces of organic solvent. Afterwards, a 2 mm Teflon spacer was sandwiched between the glass slides and the chamber was gently filled with 20 mM sucrose solution. The slides (conductive side facing inward) were connected to an AC signal generator Agilent 33220A (Agilent Technology GmbH, Germany). An AC field of voltage 1.5 V and frequency 10 Hz applied for 2 hours at room temperature, resulting in 10-50 µm sized vesicles. The harvested vesicles were diluted 10 times in 22 mM glucose solution to obtain fluctuating vesicles.

Microscopy and video recording
The equatorial fluctuations for both phase contrast and confocal mode were recorded with Leica TCS SP8 scanning confocal microscope using a HCX PL APO 40x/ Numerical Aperture (NA) 0.6 Ph2 (air) objective and a HC PL APO 40x/ NA 1.3 (oil) objective. The pinhole size during the experiment was fixed to 1 AU (Airy units) unless stated otherwise. Table 1 compiles the pixel size and focal depth for different experimental conditions. The scanning speed was fixed to 1 kHz in bidirectional mode and the polarizer plates were rotated (100%) to remove the polarization effect of the fluorescent dye unless stated otherwise. The dye was excited with a 561 nm laser (diode-pumped solid-state laser) with 1.61% (laser intensity) HyD3 detector (hybrid) and the gain was fixed to 23%. Phase contrast imaging was recorded with PCO CS dimax (PCO AG, Kelheim, Germany)) mounted on confocal microscope. 1500-2000 images were recorded at 3.83 frames per second (fps) with confocal and 60 fps with phase contrast imaging.
In this section, we list different focal depths and pixel sizes for different microscopy and numerical aperture settings for 40x objective. Focal depth or FWHM (full width half maximum) of phase contrast imaging was determined using the standard formula d = λ N A 2 The wavelength of transmission light was assumed to be 550 nm. Sub-pixel contour recognition The intensity profile in the radial direction for N wedges were determined from three different interpolation schemes (Gaussian, parabolic and linear weighting of neighbouring pixel) for sub-pixel contour recognition. This was done to check if different interpolation schemes affects the bending rigidity values due to uncertainty introduced at higher wave-numbers for experimental vesicle contour fluctuations. The mean bending rigidity obtained was similar for all the three schemes for the same vesicle. Figure (S7

S3. POLARIZATION EFFECTS
We analyzed the same vesicle with and without polarization effects. The polarization effects were corrected using circular plates that were rotated 100%. Figure (S8) illustrates the effect of dye polarization for vesicles imaged with different numerical apertures. Using one Anova test, we find a significant difference of 3 k B T for the 40x/0.6 NA case. The difference tends to be negligible for 40x/1.3 NA case.

FIG. 8. Polarization effects. (a)
Confocal images of the same vesicle with and without polarization effects for 40x/0.6 NA case. The polarization effects were removed using circular plates that were rotated 100%. (b, c) Comparison between the same vesicles for different numerical apertures. Using one Anova test, we find a significant difference of 3 kBT for 40x/0.6 NA case. The difference tends to be negligible for 40x/1.3 NA case. Pinhole size is 1 AU.

and red AU 2).
Details about the various statistical techniques can be found in Ref. [52]. Here we explain the bootstrapping sampling technique. A more rigorous reference is the textbook [53]. In practice, the finite amount of data or length of experiment limits the accuracy to infer data confidently. Bootstrapping is an inference method about the population from a given sample. In bootstrap-resamples, the population is in fact the sample and this quantity is known. This allows to measure the quality of inference of the 'true' sample from a re-sampled data. For example, let's consider the average mass of the human population world wide. It is difficult to measure the mass of every individual globally, therefore, a small sample is measured. Let's assume the sample size of N people. From that sample size, only one mean can be measured. In order to have a reasonable estimate about the population statistics, we need to have variability of the mean that we computed. The simplest bootsampling statistics can be considered by taking the original data N individuals and resampling to create a new sample of the same size N (e.g. we might 'resample' 10 times from [60, 61,62,63,64,65,66,67] kg and get [61,64,63,63,60,60,62,65] kg). This process is repeated a large number of times, 100 to 10000, to create a histogram that be applied to any estimator testing. Bootstrap resampling was carried out using MATLAB's bootstrp ().
In the case of our experiments, the finite amount of data or length of experiment limits the accuracy to infer data confidently. The bootstrap resampling requires choosing random replacement from a given data set and examining each sample the same way. This way a particular data point from the original set can reappear randomly multiple times in a particular bootstrap sample. The element size of the bootstrap sample is the same as the element size of the original data. This technique allows to obtain uncertainty of the quantity one estimates.
Bootstrap resampling algorithm for estimating standard error [53]: 1. Obtain N independent bootstrap samples X * 1 , X * 2 , X * 3 , ...X * N , each consisting of n data values drawn with a replacement from x where x = [x 1 , x 2 , x 3 ...x n ]. Note for estimating a standard error, the number N will ordinarily be larger than 30 to satisfy the Central Limit Theorem. Computations allow to use a large number N such as 10 3 to 10 4 . 2. Determine the bootstrap replication for every bootstrap resample: where s() is a statistical function like sample mean. For example, if s(x) is the sample meanx then S(X * ) is the mean of bootstrap data set. 3. Compute the standard error SE by utilizing the standard deviation of N replications In our case we determine the SE of mean Pearson correlation using bootsampling statistics.

Mathematical Model
The total energy of the system is given by the Helfrich model [54] as Eq. (5) where κ is the bending rigidity, c 1 and c 2 are the local radii curvatures, A is the total surface area, V is the interior volume of the vesicle, σ is the surface tension, and p is the pressure difference across the membrane.
For a quasi-spherical vesicle in equilibrium, the shape can be decomposed into spherical harmonics (Y lm ) such that the position of the surface is given by where the characteristic radius R 0 is given by V = 4 3 πR 3 0 . The spherical harmonics are defined as Y lm = n lm P lm (cos θ)e imφ , n lm = (2l + 1)(l − m)! 4π(l + m)!
P lm (cos θ) are the associated Legendre polynomials. As l = 1 account for translational modes, for the sake of this paper, f lm (t) will be restricted to f 1m (t) = 0 for l = 1. Furthermore, volume conservation requires that [54] Assuming there is no external fluid flow, the harmonic coefficients (f lm ) for l > 1 are described by the following stochastic differential equation [54] ∂ where and E l = (l + 2)(l − 1) l(l + 1) +σ .
λ = η in /η ex is the ratio of viscosities of the solutions inside and outside the vesicle. The dimensionless tension is σ = σR 2 0 /κ. ζ lm (t) is a stochastic term accounting for thermal noise; the corresponding time correlation is given as The δ functions are the traditional Kronecker and Dirac delta functions. From Eq. 11, the variance of ζ lm (t) is given by Numerical Method At this point, it is convenient to decompose f lm and ζ lm into real and imaginary components such that f lm (t) = X lm (t) + i Y lm (t) and ζ lm (t) = a lm (t) + i b lm (t). As a lm and b lm are independent of each other then Eq. (9) can then be rewritten as and similarly for Y lm . As Eq. (14) is a simple Langevin equation, the exact time update [? ] is given as such that ∆t is the time step size and n is a sample value from the normal distribution N (0, 1). In order to properly resolve the dynamics of the higher order coefficient, a sufficiently small time step must be chosen so that ∆t << τ lmax . Yet as each harmonic coefficient is independent of each other, Eq. 15 can be evaluated for all X lm and Y lm simultaneously. Given all the harmonic coefficients (f lm ), the cross-section at the equator, R(θ = π/2), can easily be computed using Eq. (6). Here we demonstrate an example of a numerically simulated vesicle with predefined bending rigidity and membrane tension. Figure (S10)b shows a time sequence of equatorial vesicle contours with bending rigidity of κ = 10 −19 J and membrane tension of σ = 10 −9 N/m. The size of the vesicle is R 0 = 10 −5 m. By implementing our image detection technique and fitting algorithm from Gracia et al. [7], we are able to reproduce the bending rigidity and membrane tension respectively as κ = (1.00 ± 0.01) × 10 −19 J and σ = (1.1 ± 0.2) × 10 −9 N/m with the Helfrich spectrum given in Figure (S10). Notably our image detection is able to resolve more than 45 shape fluctuation modes.

Fluctuations statistics: derivations of the basic results
Here we summarize the main results for the dynamics of a quasi-spherical vesicle [33,34,[55][56][57][58][59][60]. Mean Squared Magnitude of the Fourier Modes: The dynamics of the spherical harmonics modes is governed by the following Langevin equation, where the relaxation time τ l is given by Eq.(10) and the noise is The analytic solution to Eq. (16) is given by and therefore The ensemble average of |f lm | 2 of Eq. (19) is then Using Eq. (17), Eq. (20) simplifies to At long times, t >> τ l , Eq. (22) simplifies to Recallσ = σR 2 0 /κ. Since the dynamics of the different spherical harmonics modes are completely decoupled, we can more generally say Next, we consider the contour of the GUV at the equator as a function of the spherical harmonic coefficients: The Fourier coefficient for the q-th mode is then given by as all the other terms integrate to zero. In the above equation, we have inserted the definition of the spherical harmonic, Y(θ, φ) = n lm P lm (cos θ)e imφ (see Eq. (7)), which shows that the dependence on φ cancels out. The mean squared amplitude of u q is then given by f * l q f lq n lq n l q P lq (0)P * l q (0).
Time Correlation for Fourier Modes: Time correlations present another useful metric to analyze the membrane fluctuations. As the different spherical harmonics modes are independent, the average time correlations, can be simplified to Using (19), (30) can be rewritten as Since the first term in (32) has both the smallest decay rate (τ −1 q ) and largest mean-squared amplitude, the time correlation can be approximated to leading order as If we consider limit of undulations with short wavelengths (shorter than the vesicle radius), q 1, then the leading order decay rate can be approximated as which is the decay rate derived using planar fluctuations. However, we suggest using the exact decay rate from the spherical harmonics as it is both more accurate and valid for all Fourier modes. When comparing the time correlations in Fig.11, the exact decay rate, from the full spherical harmonics (SpH), is immediately more accurate than if the planar membrane (PM) decay rate is used. To get the accuracy even better, the higher order terms in Eq. (32) must be included. If all of the terms are included then the time correlation is directly on top of the curve from produced by the numerical simulation. However, as it is not feasible to include all the terms for real membranes, it is of interest to know how many terms are enough to sufficiently reproduce the numerical simulations. As shown in Figure 5, the time correlation produced by including the first two terms in Eq. (32 lies almost directly on top of the true solution. Including more terms would improve the accuracy further, but it is not likely to be significant due to experimental error. Cross-Spectral Density: Similar to time correlations, the Cross-Spectral Density (CSD) is given by |u q (0)||u q (t)| − |u q | 2 (0) .
For the sake of clarity of explanation, in this section we will use the leading order approximation of u q , u q (t) ≈ f qq (t)n qq P q cos π 2 .
From (37), it is clear that u q (t) =ζ q (t) for large values of t. Furthermore, it is worth noting that all Fourier modes, except q = 0, have both real and an imaginary component, u q = A q + iB q , and that these two components are independent of each other. Likewise, the thermal noise can be decomposed into independent real and imaginary components:ζ q =ζ Aq + iζ Bq . The real component of Eq. (37) can then be written as A q (t) = A q (0)e −t/τq +ζ Aq (t) (38) and a similar expression for B q . Therefore, it can be shown that |u q (0)||u q (t)| = |u q (0)| A 2 q (t) + B 2 q (t) 1/2 = |u q (0)| |ζ q | 2 (t) + 2 ζ Aq (t)A q (0) +ζ Bq (t)B q (0) e −t/τq + |u q (0)| 2 e −2t/τq 1/2 If we assume that t >> t q , then we can perform the following expansion |u q (0)||u q (t)| = |u q | 2 (0) + ζ Aq (t)A q (0) + ζ Bq (t)B q (0) e −t/τq The second term in (40) is first order with the thermal noise and thus averages to zero. Therefore, to leading order, the CSD is given as |u q (t)||u q (t)| − |u q (0)| 2 = 1 4 |u q (0)| 2 e −2t/τq + O e −3t/τq (41) and thus the slowest decaying mode is O e −2t/τq . This contradicts H. Zhou et al. [35] who give it as O e −t/τq . This factor of two is a consequence that each Fourier coefficient has both a real and imaginary component that are completely independent of each other. Finally, users are recommended to use time correlations over CSD. CSD requires the same amount of work and contains the same higher order error as the time correlation method. Yet, CSD has an additional layer of truncation error introduced in the expansion in Eq. (40).