High-throughput mechanical characterization of giant unilamellar vesicles by real-time deformability cytometry

Abstract

Real-time deformability cytometry (RT-DC) enables high-throughput, contact-free mechanical characterization of soft microscopic objects. Here we apply this technique to giant unilamellar vesicles (GUVs). To interpret vesicle deformation in RT-DC, we present a simulation-based model taking into account the area expansion modulus as the dominant mechanical parameter. Using phase-field simulations over a wide parameter space, we find GUV deformation to depend linearly on GUV area. Based on these results, we derive two complementary fitting strategies for extracting the area expansion modulus K from RT-DC data: a direct model-based fit for single-vesicle characterization and a noise-resistant collective approach that enables robust population-level estimates. Furthermore, we introduce a combined fitting method that integrates both approaches to filter outliers and improve accuracy in heterogeneous or noisy datasets. All methods scale across varying flow rates, channel geometries and buffer viscosities, and produce predictions of K consistent with literature values for different lipid compositions. Compared to traditional techniques such as micropipette aspiration, our approach offers orders of magnitude higher throughput without mechanical contact, making it particularly suitable for GUV population studies. Beyond mechanical phenotyping, this framework opens new avenues for sorting vesicle populations based on membrane mechanics, a capability of growing interest in synthetic biology and soft matter research.

---

1 Introduction

Giant unilamellar vesicles (GUVs) serve as versatile model membrane systems in biophysical research for studying biological processes such as lipid phase separation,³³ membrane fusion⁸ and protein–lipid interactions.^30,37 In addition, GUVs are used in biomedical applications such as targeted drug delivery⁴ or synthetic cell design.³¹ Their versatility stems from their cell-mimicking size (diameter >1 µm), which allows for direct observation of the membrane responses to external factors such as ions, molecules, or hydrodynamic flows under a microscope.⁷ In addition, GUVs share many structural similarities with cellular membranes, while providing precise control over factors such as membrane composition, mechanics and functionality.^3,7,18,37

Mechanical parameters, including deformability, elasticity, and membrane tension are fundamental factors influencing the behavior of GUVs and their various applications. These properties govern how the vesicles undergo dynamic shape changes, maintain structural integrity under external forces, and respond to environmental stresses. At the same time, the mechanical parameters are closely linked to the vesicle composition and are particularly influenced by factors such as lipid types, bilayer thickness and the presence of proteins.^26,29 The mechanical properties of cells, on the other hand, are often attributed to their cytoskeleton rather than the membrane itself. GUVs, in the presence and absence of cytoskeletal elements, provide the unique opportunity to decouple these effects.⁹

However, determining mechanical parameters of GUVs experimentally remains challenging and time consuming. For example, micropipette aspiration experiments,²⁶ response measurements to compression between parallel plates²⁸ or indentation with a conical AFM-tip²⁹ have been used to measure the elastic behaviour of GUVs. A common limitation of these methods is their relatively slow throughput, as they typically allow for the characterization of only one vesicle at a time.

Microfluidic deformability cytometry offers rapid characterization of thousands of GUVs within minutes, potentially enabling ultra-fast screening of large vesicle populations.

One advanced method is real-time deformability cytometry (RT-DC), which provides high-throughput, real-time measurements by subjecting the vesicles to controlled flow forces within a microfluidic channel.²³ As the vesicles traverse the channel, their deformation can be measured and analyzed to potentially quantify the mechanical parameters of GUVs such as elasticity and membrane tension. As an added advantage, the GUVs are not in contact with surfaces during the measurement, which, depending on the type of surface interaction can strongly influence the measurement and make the comparison of different conditions challenging.

Predictions of mechanical properties in RT-DC rely on matching experimentally observed cell deformation with theoretical models. For cells, such theoretical models exist and have been used extensively with RT-DC to study the mechanical properties of cells, allowing rapid phenotyping and mechanical profiling.^{10,16,17,23,27,35,36}

These models consider elastic and viscoelastic materials and are able to link shape deformation with concrete elastic moduli. However, all models consider either volumetric structures (i.e. solid bodies) or elastic shells. Lipid vesicle membranes are neither and hence cannot be evaluated with the current models.

Here, we provide a theoretical method to describe GUV deformation and to measure lipid elastic properties in an ultrafast manner in RT-DC. To this end, we present first numerical simulation results of fluid vesicles in RT-DC channels. To eliminate the influence of experimental noise we introduce a new fitting algorithm which allows real-time extraction of the area expansion modulus of a vesicle sample despite errors in the imaging data. Finally, we illustrate the generality of the approach to deliver accurate estimates when using RT-DC setups of arbitrary flow rates, channel sizes and buffer viscosities.

2 Results

2.1 Deformation of a GUV in a microfluidic channel

Lipid vesicles like GUVs are fluid filled compartments with a bilayer of densely packed lipid molecules at their surface. The lipid molecules (in the liquid disordered state) can freely move in the tangential direction along the surface, which implies a vanishing surface shear modulus. On the other hand, the lipid layer resists bending, areal compression and stretching. It has been shown that the stiffness against bending deformations is negligible in RT-DC measurements due to the large forces exerted.¹⁷ However, the layer's resistance to compression and stretching contributes significantly. The corresponding area expansion modulus K is hence the dominant mechanical parameter and the only determinant of vesicle shapes.

To extract K from RT-DC image data, we simulate GUV deformation in dependence of K. We therefore model GUVs as initially spherical thin, freely-shearable elastic surfaces with radius r immersed in a buffer liquid of viscosity η, within an RT-DC channel with a square-shaped cross-section of side length L, see Fig. 1A. A flow rate Q is imposed to set the GUV in motion, such that it subsequently deforms by pressure and shear forces until it assumes a stationary state. The full numerical method is based on a phase-field method presented in ref. 12, with more details to be found in the methods Section 4.1. The used simulations parameters (Q′ = 0.04 µL s⁻¹, L′ = 20 µm, η′ = 0.015 Pa s) can be scaled to generalize results to arbitrary Q, L, η (see Section 2.2).

Fig. 1 (A) Microfluidic setup of RT-DC and illustration of the physical model. A GUV with an area expansion modulus K is immersed in a surrounding liquid of viscosity η and traverses a channel of side length L at a flow rate Q. (B) Snapshots for RT-DC measurements of GUVs showing a spherical initial state and a bullet-shaped stationary state at the end portion of the channel. (C) Simulation of a single GUV. Exemplary evolution for initially spherical vesicle of radius r = 6.5 µm and area expansion modulus K = 0.4 N m⁻¹, colored by amount of surface stretching (see main text). (D) Stationary deformation and approximate capillary numbers for different r and K colored by amount of surface stretching. Simulation parameters: Q′ = 0.04 µL s⁻¹, L′ = 20 µm, η′ = 0.015 Pa s.

Fig. 1C shows an exemplary evolution of GUV deformation. The initially spherical surface is stretched (and transiently compressed) due to pressure and shear forces until the vesicle reaches a stationary bullet-shape, similar to the shapes observed for RT-DC measurements of GUVs (see Fig. 1B) or cells.¹⁷ To quantify deformation we use the measure.

where A denotes the area and P the perimeter of the two-dimensional image of the GUV shape, as it would be seen by a camera. Note that D = 0 for a circle and D > 0 for any other shape. Note that the GUV relaxes back to its spherical shape once it leaves the channel. Along with the time evolution of D, we illustrate the relative local stretching of the GUV surface, ( − Â₀)/Â₀, with  being the area of a surface element in the current configuration and Â₀ the area of this element at the initial time. Hence, stretch values >0 and <0 correspond to stretching and compression of the surface, respectively, see Fig. 1C.

An illustration of the stationary deformation shows that larger and softer GUVs are more stretched, while smaller GUVs with higher K remain close to spherical (Fig. 1D). The maximum stretching is observed at the lateral front of the vesicle and reaches in our simulations up to 9% for r = 8 µm, K = 0.1 N m⁻¹. This is approximately in the range of stretching tolerated by real lipid bilayers, which is typically limited to a few percents (on the order of 5%^13,22) depending on the lipid type.

For larger stretching, the lipid bilayer may rupture – a process which is not included in our simulations.

Additionally, we computed for each simulation a rough estimate of the capillary number Ca = η × [small gamma, Greek, dot above] × r/K as a measure of the ratio between the external viscous stress and the membrane's resistance to area expansion, see Fig. 1D. Using the simulation parameters Q′ = 0.04 µL s⁻¹, L′ = 20 µm, η′ = 0.015 Pa s gives an approximate shear rate of [small gamma, Greek, dot above] ≈ 35 000 s⁻¹. We find that the computed capillary numbers are fairly close to the deformation values, confirming that GUVs deform in the RT-DC setting primarily in the area-dilatational regime, where the membrane's resistance is dictated mostly by its area expansion modulus K.

Note that small fluctuations and excess area are common in GUVs.² Computing analogously the capillary numbers Ca_κ = η × [small gamma, Greek, dot above] × r³/κ based on the bending rigidity κ, with κ typically being in the order of 10⁻¹⁹ N m for GUVs,²⁵ gives values >10⁵, implying that bending contributions are negligible and possibly available excess area is smoothed out immediately. This is also confirmed by considering the observed deformation energies, which are >10⁶k_BT, much larger than typical energies of fluctuations.

2.2 Dependence on elastic modulus

The goal of this work is to describe the GUV deformation according to its mechanical properties, i.e. the area expansion elastic modulus K, for a given sample of GUVs from RT-DC experiments. The elastic behaviour of lipid bilayers differs strongly depending on the lipid compositions, with values ranging from K ∈ [0.15, 0.26] N m⁻¹ (for DOPC vesicles)^14,25 to K ∈ [1.56, 3.61] N m⁻¹ (for 1 : 1 SM/Chol vesicles).^15,26

Performing simulations for a wide range of elastic moduli and GUV sizes, we find that the deformation in the stationary state depends almost linearly on the GUV size, see Fig. 2A. Linear fits through the simulation data are included to further confirm this observation. Furthermore, an exemplary set of experimentally gathered RT-DC measurements of GUVs is depicted, showing an almost linear relation between D and A for the majority of the data as well. We note that the almost linear relation differs significantly from RT-DC results which modeled the deformation of volumetric elastic bodies¹⁷ or elastic shells^12,17 instead of fluid vesicle surfaces. The dependency between D and K can be well described by a power law D = αK^β, Fig. 2B. The good agreement of fits and simulation data is confirmed by the coefficient of determination R² > 0.99 for all curves. Using these fitted functions now allows for extrapolations outside of the used parameter range, s. t. also vesicles with smaller or larger expected values of K can be treated.

Fig. 2 Simulation data fits. (A) Simulation data (dots) fitted by linear functions (lines) illustrating the linear relation between area and deformation. Exemplary experimental data showing an almost linear relationship in GUV measurements as well. (B) Simulation data (dots) with exponential fits (lines) illustrating the relation between area expansion modulus and deformation. (C) Function fit for the relation between K, area and deformation based on eqn (3) (direct fitting). (D) Fitted function between K and the slope of the linear simulation curves used for the collective fitting method. Parameters: Q′ = 0.04 µL s⁻¹, L′ = 20 µm, η′ = 0.015 Pa s.

Moreover, the quasi-stationary state of the GUVs is determined by a single non-dimensional parameter K × (L′)²/(η′ × Q′).¹⁶ Hence, analogously to the numerical model presented in ref. 17 our computed results (see Fig. 2A) can be transformed to other flow rates, fluid viscosities and channel sizes. More precisely, a change of the parameters from η′, Q′, L′ to the values η, Q, L can be realized by simply rescaling the area expansion modulus K with factor (ηQ(L′)²)/(η′Q′L²) and rescaling the vesicle area with a factor (L′/L)².

2.2.1 Direct fitting. Combining both observed dependencies from Fig. 2(A) and (B), we postulate the following fitting function to describe an explicit relationship between A, D and K:

K̂ = γ × D^−βÂα × Â^−δ, (2)

with the dimensionless quantities K̂ = K × L²/(ηQ) and  = A/L². The four dimensionless unknowns α, β, γ, δ are obtained by a least squares fit with the simulation data (i.e. minimizing the squared difference between K̂ from eqn (2) and the scaled, dimensionless values from the simulation). We obtain an optimal fit for α ≈ 0.5, β ≈ 2.16, γ ≈ 3.40 × 10⁻⁴, δ ≈ 1. For arbitrary experimental parameters η, Q, L which may differ from the ones used in our simulations (η′ = 0.015 Pa s, Q′ = 0.04 µL s⁻¹, L′ = 20 µm) and when fixing α = 0.5 and δ = 1 we obtain the final formula

Equation (3) now allows to extract the area expansion modulus K directly for each vesicle imaged in RT-DC. The formula is valid for any flow rate, any medium viscosity, and any channel size, as long as the deformation D and the relative area A/L² are fairly in the range of the here presented simulations (D ≈ ∈ [0, 0.06], A/L² ≈ ∈ [0.12, 0.56]). Fig. 2C shows that the fitting function describes the data quite well with a coefficient of determination of R² ≈ 0.97 and a sum of squared residuals of 0.091.

However, the high-throughput analysis of GUVs in RT-DC is subject to significant data scatter stemming from both inherent physical heterogeneity of the vesicles and various measurement errors. These errors include noise from the imaging system (e.g., low contrast, focus problems, electronic noise) which corrupts the precise determination of the vesicle contour, and additional perturbations (e.g. hydrodynamic interactions of nearby vesicles, flow instabilities, temperature fluctuations) which randomly modulates the effective hydrodynamic stress. Averaging the obtained K-values can cancel out noise. To get a representative, statistically less biased value for the whole data set we will introduce a second fitting method in the following section.

2.2.2 Noise independent collective fitting. Knowing that the theoretical deformation D depends linearly on the vesicle area A (cf. Fig. 2A), we find that the slope of the area-deformation curves is a robust marker for the area expansion modulus K. In fact, postulating the power law K̂ = a× (slope × L²)^−b we find an excellent fit (R² > 0.999) for the dimensionless parameters a ≈ 0.0015 and b ≈ 1.75, see Fig. 2D. Hence, the formularobustly relates the area-deformation slope of data points to their expansion modulus K for any given Q, L, and η.

This insight opens the possibility to averaging out measurement errors directly in the data. The precise idea is to apply a linear regression to the raw experimental data in the area-deformation diagram, and then use the obtained slope to extract the area expansion modulus K via eqn (4). This approach, which we will refer to as collective fitting, promises a more robust identification of K, since noise is averaged out directly on the data set. Note, that this simple way of averaging the data is only enabled, by the observed linearity between theoretical area-deformation curves, which does not apply to RT-DC of volumetric elastic bodies.

We point out, that this idea to first average, then compute K, is obviously different from the approach in the last chapter (compute K for each data point, then average over all K). In fact, due to the nonlinear dependency of K and D (cf. eqn (3)), Jensen's inequality applies, which states that the estimations of both variants must be different.

With eqn (4) we can directly infer K from the entire given experimental data set by simply evaluating the slope of a linear fit through the data. Analogously to eqn (3) and (4) is universal and can be used for arbitrary flow rates, buffer viscosities and channel sizes as long as the obtained slopes and relative vesicle areas (A/L²) remain in the range performed in this study.

2.3 Extraction of dilational elasticity from DOPC vesicles

Next, we compared our theoretical results to experimental data. RT-DC was carried out by flushing electroformed GUVs of 100% DOPC through a standard RT-DC chip with a 20 µm crossection at a flow rate Q = 0.04 µL s⁻¹ and a corresponding buffer viscosity η ≈ 6.3 mPa s, which was estimated using the Buyukurganci-2022 model⁵ for a 0.59% MC-PBS solution (CellCarrier B). GUVs were imaged using high-speed brightfield microscopy and their contour was detected. Considering the imaged area and the perimeter of the contour, the deformation D of the GUVs, as a typical output of RT-DC measurements, was calculated using eqn (1) and plotted against the GUV area, Fig. 3A (left).

Fig. 3 Extraction of K using the combined fitting approach. (A) Comparison of the three different fitting methods for a given data set with snapshots of measured GUVs (here pure DOPC vesicles). Left: Extraction of K using collective fitting (eqn (4)). Middle: Extraction of K using direct fitting (eqn (3)). Computation and visualization of Gaussian mixture model (GMM) with two components. Right: Classification of input data based on the Gaussian mixture model (filtering for outliers). Refitting and noise elimination for the usable data with collective fitting approach (eqn (4)). (B) Extraction of K for different sets of pure DOPC vesicles measured at different flow rates (Q = 0.04 µL s⁻¹, Q = 0.08 µL s⁻¹). (C) Extraction of K for different sets of binary composed SM-Chol vesicles (1 : 1) measured at different flow rates (Q = 0.04 µL s⁻¹, Q = 0.08 µL s⁻¹).

Due to the short passage times of GUVs (∼1 ms), we observe almost constant GUV volume across the channel (volume change <1% in the first 200 µm behind the entry of the channel), see Section 4.4.

Computing the slope of the linear fit through this data enabled us to compute via eqn (4) an area expansion modulus of K = 0.186 N m⁻¹ for the GUV sample. Alternatively, using the direct fitting, we computed K for every data point. The corresponding histogram shows a peak around K = 0.19 N m⁻¹, Fig. 3A (middle). Hence, the extracted K values of both the direct and the collective methods are in perfect agreement with the literature values (K ∈ [0.15, 0.26] N m⁻¹,^14,25).

Since the measured deformation D, is inherently non-negative, the error distribution should be best approximated by a Lognormal distribution. Given the power law relationship between K and D (cf. eqn (3)), the expansion modulus K should also follow a Lognormal distribution. However, the histogram of K, when plotted on a logarithmic scale (or, equivalently, the probability density of ln K), exhibits a clear asymmetry around the peak, deviating from the expected symmetric Gaussian form. This asymmetry suggests the presence of two distinct data populations. Closer inspection of data points corresponding to very low K (high D) revealed those to be unreliable measurements or artifacts (e.g., vesicle doublets or misplaced contour lines). We therefore concluded that the dataset contains two statistically separable groups: the true physical measurements (usable GUVs) and the unreliable outliers. To statistically separate these populations, we model the observed probability density of ln K using a Gaussian Mixture Model (GMM), assuming with mixture weights ω_A and ω_B. The resulting GMM provides an excellent fit to the observed probability density (Fig. 3A middle, orange curve). The expansion modulus extracted from the dominant component, K_A = 0.19 N m⁻¹, shows robust agreement with the previously reported values.

While both methods (direct and collective fitting) give equally good results for the analyzed data set, the presence of large data scatter and unreliable measurements (outliers) may corrupt the extracted K values for smaller sample sets. We therefore propose a combined method using the statistically unbiased collective fitting after largely eliminating outliers with the direct fitting approach, as described in the next section.

2.3.1 Combined fitting gives robust estimates of area expansion modulus. As described before, the presence of wrongly measured outliers is observed and expected, even after data postprocessing (see Section 4.5). The corresponding variations in deformation might strongly influence the predictions, especially for smaller vesicle populations.

To improve results in such cases, we propose a combined fitting strategy: first, remaining outliers are filtered out of the data set. Using the Gaussian Mixture Model we classify data points by calculating how likely the sample belongs to each of the mixture components (posterior probability of each class given the point) and assigning them to the component with the highest probability. Hence, the given data set is divided into two clearly distinct clusters. The cluster with the lower elasticity (inferred from higher deformations) is marked as outliers and eliminated from the further analysis. In a second step, applying the collective fitting method to the remaining data points gives a noise independent prediction of the area expansion modulus of the vesicles, see Fig. 3A (right). The obtained modulus K = 0.18 N m⁻¹ is again comparable to literature values.¹⁴

The approximately equal results for all three fitting strategies (direct, collective, and combined) indicate robustness of all these approaches for the large data set used. However, given the inherent variability present for smaller datasets or stiffer GUVs, the extracted K-values show moderate differences between these models. The combined fitting method is anticipated to be the most reliable, as it effectively filters outliers and exhibits greater insensitivity to noise. Therefore, all subsequent data analysis and reported parameters will be based on the combined fitting strategy.

2.3.2 Extraction for different lipid compositions. We illustrate the power of collective fitting with supplementary filtering step for different lipid compositions in the following. It is known that the lipid composition influences K, however literature values vary. RT-DC would be the first high-throughput method to measure K of vesicle populations.

First we analyzed different data sets of 100% DOPC vesicles for the two flow rates Q = 0.04 µL s⁻¹ and Q = 0.08 µL s⁻¹. The usage of different flow rates also influences the viscosity of the surrounding fluid due to shear thinning. To account for this, we computed the effective viscosity again via the model from ref. 5 and obtained η ≈ 6.3 mPa s for Q = 0.04 µL s⁻¹ and η ≈ 4.7 mPa s for Q = 0.08 µL s⁻¹.

Following the combined fitting procedure from Section 2.3.1, we computed a linear fit through the filtered data and computed its slope for two independently measured data sets. Using eqn (4) we obtained elastic moduli between (K = 0.13 N m⁻¹ and K = 0.23 N m⁻¹), see Fig. 3B. All the predictions are fairly in the expected range derived from micropipette aspiration (K ∈ [0.15, 0.26] N m⁻¹,^14,25).

In a second set of experiments we studied SM-Chol vesicles (1 : 1 mixture). These vesicles are known to be much stiffer then DOPC vesicles with an expected value for K in the range K ≈ ∈ [1.5, 3.6] N m⁻¹ (ref. 15 and 26). Again we extracted the area expansion modulus for multiple data sets at the two flow rates Q = 0.04 µL s⁻¹ and Q = 0.08 µL s⁻¹, see Fig. 3C. The extracted values for K are in the range of K ≈ ∈ [2.08, 3.18] N m⁻¹, which agrees perfectly with the expected reference range.

Moreover, the variations of extracted K-values between the different data sets are for both sets of experiments comparable to those between different flow rates, highlighting the universality of the extraction method to give comparable results across different flow rates and buffer viscosities.

In contrast to the controlled membrane tension environment used in aspiration experiments, GUVs in RT-DC are not initially subjected to a consistent prestress. As a result, small excess area, common in GUVs,² must be smoothed out first, leading to an initial deformation response that does not reflect the true elastic stretching of the lipid bilayer. Consequently, the apparent stiffness deduced from RT-DC experiments tends to underestimate the intrinsic area expansion modulus K of the membrane. However, the observation that varying flow rates result in differing deformations (D) but yield similar elastic moduli (K) provides a strong justification that the influence of membrane fluctuations and excess area is secondary. If significant membrane excess area were present, it would stretch out already at very low flow rates, introducing a dominant, flow-rate-independent component to the deformation response. This would in turn lead to calculated modulus K increasing almost linearly with the flow rate (cf. eqn (3) and (4)), which we do clearly not observe.

2.4 Analyzing mixed sets of differently composed vesicles

The previous examples illustrate that we can extract the area expansion modulus K for a given set of identical GUVs using the collective fitting with a supplementary filtering step. In the following, we will apply the combined fitting method to classify and identify mixed sets of different GUVs.

We consider a data set consisting of 1000 randomly selected samples from a given DOPC set as well as 1000 samples from a set of binary composed SM-Chol-vesicles measured at the same flow rate Q = 0.08 µL s⁻¹ with areas in the range A ∈ [60, 225] µm² and deformations D ≤ 0.06. Fig. 4 (top left) shows the merged data set. As can be seen from the plot, the data is spread across the whole domain, while forming two distinguishable clusters as well as an amount of outliers. This suggests a mixture of (at least) two different types of GUVs within the data. The composition of the data set is assumed to be unknown in this example. Instead we will distinguish the given data by the actual mechanical behaviour. This is again done by applying the direct fitting method using function (3), giving an estimate for the area expansion modulus for any vesicle in the given data set. Visualizing all estimates in a histogram (see Fig. 4 (top right)) illustrates that the data (indeed) consists of at least two different types of vesicles. Note that the predictions are again plotted on a logarithmic scale to visualize variations properly. Following the ideas from Section 2.3.1 while assuming three classes (type A, type B, outliers O) within our data, we use a mixture of three Gaussian distributions to classify the distribution of the data, i.e.

Fig. 4 Extraction of area expansion modulus K for a mixture of 1000 DOPC- and 1000 SM-Chol-vesicles. Top (left): Given (merged) data set. Top (right) – Step 1: Extraction of K using direct fitting (eqn (3)). Computation and visualization of Gaussian mixture model with three components. Bottom (left) – Step 2: Classification of input data based on the Gaussian mixture model. Refitting and noise elimination for the two classes with collective fitting approach (eqn (4)). Bottom (right) – Step 3: Comparison of data classifications and K-extractions with given ground truth (DOPC-vesicles are colored in red, SM-Chol-vesicles are colored in blue), containing possible unreliably results.

For the computed histogram we get a Gaussian mixture distribution with mean values K_O ≈ 0.08 N m⁻¹, K_A ≈ 0.24 N m⁻¹ and K_B ≈ 1.34 N m⁻¹ and corresponding standard deviations of σ_A ≈ 0.58, σ_B ≈ 1.17 and σ_O ≈ 0.55 respectively. Also the mixture weights are ω_A ≈ 0.50, ω_B ≈ 0.46 and ω_O ≈ 0.04. The ratio of the dominant weights (ω_A and ω_B) accurately indicates the real half-and-half mixing ratio of the two populations.

In a second step, we use these results to classify the given data to any of the three classes based on the probability distribution and apply again collective fitting to eliminate noise within the data and to get more precise K-predictions for the data set.

Fig. 4 (bottom left) shows the data set after classification and filtering for outliers. Also we applied collective fitting for both vesicle types. For type A we get an estimate of K_typeA ≈ 0.16 N m⁻¹. For the vesicles of type B we get an estimate of K_typeB ≈ 2.68 N m⁻¹.

Finally, since the vesicle compositions of the given data are known, we can compare in a third step the estimates of our fitting methods with the ground truth. Fig. 4 (bottom right) shows the merged data set with 1000 DOPC-vesicles (colored in red) and 1000 SM-Chol-vesicles (colored in blue) and their respective reference values for K.²⁶

Firstly, we observe that the direct fitting with function (K) in step one gives clearly distinguishable estimates for K. This suggests that function (3), despite its sensitivity to noise, can be used to indicate the presence of different types of GUVs in a data set and even help to distinguish them by their estimates of the area expansion modulus.

Despite possible misclassification of some of the GUVs due to their similar measured behavior, we get reasonable predictions for both DOPC (K_typeA ≈ 0.16 N m⁻¹ vs. K_ref ≈ ∈ [0.15, 0.26] N m⁻¹) and SM-Chol vesicles (K_typeB ≈ 2.68 N m⁻¹ vs. K_ref ≈ ∈ [1.5, 3.6] N m⁻¹) when using the combined fitting.

Note, that the extracted K-values from the initial direct fitting method deviate from the final results of the combined fitting (K_A ≈ 0.24 N m⁻¹ vs. K_typeA ≈ 0.16 N m⁻¹ and K_B ≈ 1.34 N m⁻¹ vs. K_typeB ≈ 2.68 N m⁻¹). This discrepancy is a direct consequence of the different statistical treatments: the GMM-derived component means reflect a soft, posterior-weighted expectation that accounts for contributions of all data points according to their classification probabilities (hence points in regions of overlap influence multiple components), whereas the combined method uses a hard assignment (each vesicle is assigned to the most likely class) and then extracts K for each class.

This example illustrates that the combined fitting approach can also be applied to detect the presence of different types of GUVs within a data set and extract the corresponding elastic moduli.

3 Conclusion

In this study, we applied real-time deformability cytometry to characterize giant unilamellar vesicles. To quantify GUV mechanical properties, we introduced a numerical model to describe GUV deformation during this contact-free high-throughput microfluidic process. Including the dilational elasticity of the vesicle surface as the dominant mechanical element, we find that vesicles assume a stationary bullet shape, akin to experimentally observed shapes of cells. Interestingly, we find that GUV deformation depends linearly on the imaged area of the GUV, with a slope that depends on the area expansion modulus K. Combined with experimental measurements of various vesicle types, we provided evidence that K is the primary factor governing vesicles deformation, while bending contributions are negligible and effects of possibly available excess area are small.

Based on simulation results for a large number of parameter combinations, we developed fitting functions to extract K. We showed that our results enable reasonable predictions of K for arbitrary flow rates, buffer viscosities and channel sizes, even compensating noise. The direct fitting approach (Section 2.2.1) enables K-predictions for single vesicles from imaging their deformation and size. While this approach is more sensitive to noisy image data, it can be used to filter for outliers in wide spread data (see Section 2.3.2) or classify different types of vesicles in a mixed data set, as shown in Section 2.4, as long as the vesicles are distinguishable by their mechanical behavior.

To obtain statistically less biased estimates, we introduced the collective fitting method (see Section 2.2.2) based on the slope of area-deformation curves. This method yields the averaged area expansion modulus for the whole data set. It works well for noisy image data as it combines information of a whole collective of vesicle deformations to eliminate noise effects in the measured vesicle contours.

The presence of strong variations within the experimental data still influences the results of the collective fitting method. To eliminate the influence of these outliers and ensure the applicability of our method, we introduced the combined fitting approach (see Section 2.3.1). It combines the advantages of both fitting methods leading to more reliable and robust predictions, as shown in Section 2.3.2. Also, we showed the applicability of this method for mixed sets of differently composed vesicles (see Section 2.4).

All fitting approaches are presented here such that they scale for arbitrary experimental settings, i.e. arbitrary flow rates, channel sizes and buffer viscosities. Compared to micropipette aspiration, RT-DC tends to slightly underestimate the area expansion modulus, probably due to the absence of a controlled membrane tension environment. Nevertheless, RT-DC reproduces the mechanical properties of vesicles reasonably well and provides significantly higher throughput in a contact-free manner. Therefore, we believe that it is particularly suited for GUV population analysis. RT-DC provides additional information, like GUV numbers and size distributions which are useful for downstream experiments. For synthetic cell research and beyond, it will become exciting to sort GUVs dependent on their mechanical properties which would lie within the capabilities of the existing RT-DC method.^20,21

In summary, we successfully demonstrated the application of RT-DC for characterizing GUVs. Beyond precisely measuring the area expansion modulus (K), this study establishes the capability to simultaneously quantifying vesicle size, membrane viscosity, and mechanical heterogeneity, thus paving the way for rapid screening and sorting of GUVs by RT-DC.

4 Methods

4.1 Numerical simulations

The presented simulations describe a single initially spherical vesicle with radius r and area expansion modulus K placed in the middle of a square channel of length L = 20 µm. Approximating the vesicle shapes as rotationally symmetric and using the concept of the equivalent channel radius,¹⁶ the simulations are reduced to a 2D axisymmetric setting with a cylindrical channel of radius R = 10.94 µm. This approximation by axisymmetric shapes was analyzed in the context of RT-DC for volumetric elastic bodies in ref. 36. Based on their results we expect accurate results for small vesicles, but a slight underestimation of up to 10% in the area expansion modulus of the largest vesicles used in this study (area ≈200 µm). The computational domain is prescribed to move with the vesicle barycenter to well resolve the flows and forces around the vesicle at all times. In this case, as shown in ref. 17 a channel length of 40 µm is sufficient to provide accurate vesicle deformation. Accordingly, the computational domain is chosen Ω = [0, 40] µm × [0, 10.94] µm.

The axisymmetric model is based on a phase-field description of elastic surfaces in fluid flow developed in ref. 12. The resulting system is discretized using a semi-implicit Euler scheme in time and a finite-element-scheme in space (see ref. 12, Section 4). The timestep size is chosen as 0.05 µs. The grid size is approximately 0.5 µm in the fluid phases, while we use an adaptive refinement strategy to ensure a sufficient smooth resolution of the used phase-field function (by at least five grid points at any time around the interface) to describe the vesicle membrane.

The laminar flow with magnitude 0.21 ms⁻¹ is induced through appropriate Dirichlet boundary conditions for the velocity field v. This ensures the typical flow rate Q′ = 0.04 µL s⁻¹. No slip boundary conditions are imposed on the channel wall as well. In the simulations the densities and viscosities within the vesicles as well as in the surrounding fluid are set equally as ρ₀ = ρ₁ = 1000 kg m⁻³ and η′ = η₀ = η₁ = 0.015 Pa s. Note, that stationary results are independent of the viscosity of the vesicles's interior η₁ since there is (almost) no internal flow in the stationary state.

The problems are discretized in the finite element toolbox AMDiS.^24,32,34 The resulting linear system of equations is solved using the direct solver UMFPACK.⁶

4.2 Fitting procedure

Throughout this work we used multiple function fits to describe different relationships within the simulation and experimental data. Each fitting problem is formulated as an unconstrained least-squares optimization problem and solved based on the Levenberg–Marquardt algorithm (LMA,¹⁹), which is a standard solver in many software packages. In case of the one-dimensional function fits, e.g. the relation between K and the slope of the simulation curves, we directly used the computed parameters. For the two-dimensional fit based on eqn (2) we used a two-stepped fitting procedure to reduce the complexity of the problem while maintaining the accuracy of the solution. It turned out, that fixing the parameters α = 0.5 and δ = 1 gives as good results as solving for four unknowns, which led to the more elegant function (3).

The effective viscosities necessary for all our calculations were estimated using the experimentally established model from ref. 5 assuming CellCarrier B (Zellmechanik Dresden) as a buffer medium. The slight dilution of CellCarrier B used in our experimental setup (Section 4.4) tends to marginally influence fluid viscosity, but has only a negligible effect on the predicted K. Local variations in shear thinning were shown to have negligible effect on the mechanical characterization by RT-DC at the considered flow rates.³⁶

4.3 GUV preparation

Giant unilamellar vesicles (GUVs) composed of 100% 18 : 1 (Δ 9-Cis) (PCDOPC), or 50% 18 : 0 SM (d18 : 1/18 : 0)(SM) and 50% Chol (Avanti Polar Lipids) were prepared via electroformation as described in ref. 1, using a Vesicle Prep Pro (Nanion Technologies GmbH). In brief, 40 µL of the lipid mixture dissolved at 5 mM in chloroform were spread on the conductive side of an indium tin oxide (ITO) coated glass slide (Visiontek Systems Ltd) and dried under vacuum for 15 min. Using vacuum grease (Carl Roth), a rubber ring was placed on top of the lipid coated slide. The ring was filled with 275 µL of a 320 mM sucrose solution (Sigma-Aldrich) and a second ITO slide was used to create a closed chamber with the conductive side facing downward. Next, the electrodes of the Vesicle Prep Pro were connected to the conductive side of the ITO slides and the Standard program was started, generating an AC field of 3 V at 5 Hz for 2 h. The temperature was set to 70 °C for GUVs containing SM and to 37 °C for GUVs containing DOPC.

4.4 RT-DC measurements

For experimental determination of vesicle deformation, an AcCellerator (Zellmechanik Dresden) combined with a high-speed CMOS camera (MC1362, Microtron) on an inverted microscope (AxioObserver, Carl Zeiss AG) with an xy-stage and a 40×/0.65 objective was used. 150 µL of 320 mM Glucose solution was added to 90 µL s⁻¹ of electroformed vesicles in 320 mM Sucrose solution in an Eppendorf tube. After 10 min vesicles were collected from the bottom of the tube and mixed with CellCarrier B (Zellmechanik Dresden) in a 1 : 9 volumetric ratio. 1 mL glass syringes filled with CellCarrier B or GUVs in CellCarrier B were mounted into neMESYS syringe pumps (Cetoni) and connected to a rectangular 20-µm-wide channel (Flic20, Zellmechanik Dresden) polydimethylsiloxane chip, via PTFE tubing (S181012, BOLA) and a PTFE plunger (SETonic). The chip was mounted onto the microscope stage and measurements were conducted using Shape-In 2 (Zellmechanik Dresden) at a mean temperature of 23.7 °C, flow rates of 0.04, 0.06, 0.08 or 0.12 µL s⁻¹ and a sheath flow three times higher than the sample flow.

Vesicles were measured at the end portion of the squared channel in their stationary state with sufficient distance to the channel exit (>30 µm).

Moreover, due to the short passage times of GUVs (∼1 ms), we observed almost constant GUV volume across the channel (volume change <1% in the first 200 µm behind the entry of the channel), see Fig. 5.

Fig. 5 Volume (normalized) over channel position for RT-DC measurements of GUVs (flowing from right to left) showing almost constant GUV volume across the channel.

4.5 RT-DC data post processing

It is known that in RT-DC pixelation effects lead to a systematic deformation offset, such that even perfectly circular objects appear slightly deformed when their contours are pictured on a pixel grid. Herold¹¹ quantified this effect for the camera resolution used in our experimental setup (340 nm per pixel) and showed that the deviation decreases with increasing object size but remains non-negligible, particularly for small objects. To avoid systematic overestimation of deformation, we apply the published correction function (ref. 11, eqn (1)) to all our data sets prior to further analysis. This ensures that our presented results in this work are not strongly biased by pixelation artifacts. Finally, we restricted the experimental data to the range of our simulation data, i.e. we just considered data points for A ∈ [60, 225] µm² and D ≤ 0.1.

Moreover, the contour detection of giant unilamellar vesicles in RT-DC is susceptible to segmentation errors, often stemming from low image contrast or poorly defined vesicle edges. Such inaccuracies yield irregular or non-convex contours, leading to unreliable measurements of vesicle area (A) and deformation (D).

To robustly eliminate these artifacts, we employed a filtering strategy based on the ratio between convexified and non-convexified vesicle areas. This ratio, also known as RT-DC ‘porosity’, should ideally be close to unity (≈1.0). Deviations significantly greater than 1.0 indicate geometric flaws, such as indentations from overlapping vesicles (doublets) or severely jagged contours.

We chose a robust threshold and excluded all data points where the porosity was greater than 1.03. This criterion effectively removed the majority of missegmented contours and ambiguous outliers from the dataset. While additional parameters, such as aspect ratio and brightness metrics, were also evaluated as quality indicators, the convexity criterion proved to be the most effective single filter for isolating geometrically compromised events.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data that support the findings of this study are openly available on figshare at https://doi.org/10.6084/m9.figshare.29480798.

Acknowledgements

SA acknowledges support from DFG (grant numbers 328170591 and 417223351). The authors gratefully acknowledge computing time on the high-performance computer at the NHR Center of TU Dresden. This center is jointly supported by the Federal Ministry of Education and Research and the state governments participating in the NHR (https://www.nhr-verein.de/unsere-partner). Also the authors acknowledge computing time on the compute cluster of the Faculty of Mathematics and Computer Science of Technische Universität Bergakademie Freiberg, operated by the computing center (URZ) and funded by the Deutsche Forschungsgemeinschaft (DFG) under DFG grant number 397252409. K. G. was funded by the Deutsche Forschungsgemeinschaft (DFG) – SFB-1638/1 – 511488495 – P15 and the ERC starting grant ENSYNC (101076997). K. G. further acknowledges funding from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy via the Excellence Cluster 3D Matter Made to Order (EXC-2082/1 – 390761711).

References

1. T. Abele, T. Messer, K. Jahnke, M. Hippler, M. Bastmeyer, M. Wegener and K. Göpfrich, Two-photon 3d laser printing inside synthetic cells, Adv. Mater., 2022, 34(6), 2106709.
2. M. Aleksanyan, H. A. Faizi, M.-A. Kirmpaki, P. M. Vlahovska, K. A. Riske and R. Dimova, Assessing membrane material properties from the response of giant unilamellar vesicles to electric fields, Adv. Phys.: X, 2023, 8(1), 2125342 PMID: 36211231.
3. B. Ali Doosti, D. Fjällborg, K. Kustanovich, A. Jesorka, A.-S. Cans and T. Lobovkina, Generation of interconnected vesicles in a liposomal cell model, Sci. Rep., 2020, 10(1), 14040.
4. T. M. Allen and P. R. Cullis, Liposomal drug delivery systems: from concept to clinical applications, Adv. Drug Delivery Rev., 2013, 65(1), 36–48 . Advanced Drug Delivery: Perspectives and Prospects.
5. B. Büyükurganci, S. K. Basu, M. Neuner, J. Guck, A. Wierschem and F. Reichel, Shear rheology of methyl cellulose based solutions for cell mechanical measurements at high shear rates, Soft Matter, 2023, 19, 1739–1748.
6. T. A. Davis, Algorithm 832: Umfpack v4.3-an unsymmetric-pattern multifrontal method, ACM Trans. Math. Softw., 2004, 30(2), 196–199.
7. R. Dimova, Giant vesicles and their use in assays for assessing membrane phase state, curvature, mechanics, and electrical properties, Ann. Rev. Biophys., 2019, 48, 93–119.
8. R. Dimova, K. A. Riske and M. Damijan, Electrodeformation, electroporation, and electrofusion of giant unilamellar vesicles, Handbook of Electroporation ( 2016), pp. 235–252.
9. A. Fink, S. Fazliev, T. Abele, J. P. Spatz, K. Göpfrich and E. A. Cavalcanti-Adam, Membrane localization of actin filaments stabilizes giant unilamellar vesicles against external deforming forces, Eur. J. Cell Biol., 2024, 103(2), 151428.
10. B. Fregin, F. Czerwinski, D. Biedenweg, S. Girardo, S. Gross, K. Aurich and O. Otto, Dynamic real-time deformability cytometry – time-resolved mechanical single cell analysis at 100 cells/s, Biophys. J., 2020, 118, 605a.
11. C. Herold, Mapping of deformation to apparent young's modulus in real-time deformability cytometry, arXiv, 2017, preprint, arXiv:1704.00572 DOI:10.48550/arXiv.1704.00572.
12. M. Kloppe and S. Aland, A phase-field model of elastic and viscoelastic surfaces in fluids, Comput. Methods Appl. Mech. Eng., 2024, 428, 117090.
13. R. Lipowsky, Remodeling of membrane shape and topology by curvature elasticity and membrane tension, Adv. Biol., 2022, 6(1), 2101020.
14. L. Lu, W. J. Doak, J. W. Schertzer and P. R. Chiarot, Membrane mechanical properties of synthetic asymmetric phospholipid vesicles, Soft Matter, 2016, 12(36), 7521–7528.
15. T. J. McIntosh, S. A. Simon, D. Needham and C. H. Huang, Structure and cohesive properties of sphingomyelin/cholesterol bilayers, Biochemistry, 1992, 31(7), 2012–2020 PMID: 1536844.
16. A. Mietke, O. Otto, S. Girardo, P. Rosendahl, A. Taubenberger, S. Golfier, E. Ulbricht, S. Aland, J. Guck and E. Fischer-Friedrich, Extracting cell stiffness from real-time deformability cytometry: theory and experiment, Biophys. J., 2015, 109(10), 2023–2036.
17. M. Mokbel, D. Mokbel, A. Mietke, N. Träber, S. Girardo, O. Otto, J. Guck and S. Aland, Numerical simulation of real-time deformability cytometry to extract cell mechanical properties, ACS Biomater. Sci. Eng., 2017, 3(11), 2962–2973 PMID: 33418716.
18. L.-R. Montes, A. Alonso, F. M. Goni and L. A. Bagatolli, Giant unilamellar vesicles electroformed from native membranes and organic lipid mixtures under physiological conditions, Biophys. J., 2007, 93(10), 3548–3554.
19. J. J. Moré, The Levenberg-Marquardt algorithm: implementation and theory, in Numerical Analysis, ed. G. A. Watson, Springer Berlin, Heidelberg, 1978, pp. 105–116.
20. A. A. Nawaz, D. Soteriou, C. K. Xu, R. Goswami, M. Herbig, J. Guck and S. Girardo, Image-based cell sorting using focused travelling surface acoustic waves, Lab Chip, 2023, 23(2), 372–387.
21. A. A. Nawaz, M. Urbanska, M. Herbig, M. Nötzel, M. Kräter, P. Rosendahl, C. Herold, N. Toepfner, M. Kubánková and R. Goswami, et al., Intelligent image-based deformation-assisted cell sorting with molecular specificity, Nat. Methods, 2020, 17(6), 595–599.
22. V. N. Ngassam, W.-C. Su, D. L. Gettel, Y. Deng, Z. Yang, N. Wang-Tomic, V. P. Sharma, S. Purushothaman and A. N. Parikh, Recurrent dynamics of rupture transitions of giant lipid vesicles at solid surfaces, Biophys. J., 2021, 120(4), 586–597.
23. O. Otto, P. Rosendahl, A. Mietke, S. Golfier, C. Herold, D. Klaue, S. Girardo, S. Pagliara, A. Ekpenyong, A. Jacobi, M. Wobus, N. Töpfner, U. F. Keyser, J. Mansfeld, E. Fischer-Friedrich and J. Guck, Real-time deformability cytometry: on-the-fly cell mechanical phenotyping, Nat. Methods, 2015, 12, 199–202.
24. S. Praetorius, The adaptive multi-dimensional simulation toolbox (AMDiS), a discretization module on top of the Dune framework, https://gitlab.com/amdis/amdis.
25. W. Rawicz, K. Olbrich, T. McIntosh, D. Needham and E. Evans, Effect of chain length and unsaturation on elasticity of lipid bilayers, Biophys. J., 2000, 79(1), 328–339.
26. W. Rawicz, B. Smith, T. McIntosh, S. Simon and E. Evans, Elasticity, strength, and water permeability of bilayers that contain raft microdomain-forming lipids, Biophys. J., 2008, 94(12), 4725–4736.
27. P. Rosendahl, K. Plak, A. Jacobi, M. Kraeter, N. Toepfner, O. Otto, C. Herold, M. Winzi, M. Herbig and Y. Ge, et al., Real-time fluorescence and deformability cytometry, Nat. Methods, 2018, 15(5), 355–358.
28. E. Schäfer, T.-T. Kliesch and A. Janshoff, Mechanical properties of giant liposomes compressed between two parallel plates: impact of artificial actin shells, Langmuir, 2013, 29(33), 10463–10474 PMID: 23869855.
29. E. Schäfer, M. Vache, T.-T. Kliesch and A. Janshoff, Mechanical response of adherent giant liposomes to indentation with a conical afm-tip, Soft Matter, 2015, 11, 4487–4495.
30. E. Sezgin and P. Schwille, Model membrane platforms to study protein-membrane interactions, Mol. Membr. Biol., 2012, 29(5), 144–154.
31. L. Van de Cauter, L. van Buren, G. H. Koenderink and K. A. Ganzinger, Exploring giant unilamellar vesicle production for artificial cells – current challenges and future directions, Small Methods, 2023, 7(12), 2300416.
32. S. Vey and A. Voigt, AMDiS: adaptive multidimensional simulations, Comput. Visualization Sci., 2007, 10, 57–67.
33. O. Wesołowska, K. Michalak, J. Maniewska and A. Hendrich, Giant unilamellar vesicles-a perfect tool to visualize phase separation and lipid rafts in model systems, Acta Biochim. Pol., 2009, 56(1), 33–39.
34. T. Witkowski, S. Ling, S. Praetorius and A. Voigt, Software concepts and numerical algorithms for a scalable adaptive parallel finite element method, Adv. Comput. Math., 2015, 41, 1145–1177.
35. L. D. Wittwer, F. Reichel and S. Aland, Numerical simulation of deformability cytometry: transport of a biological cell through a microfluidic channel, in Modeling of Mass Transport Processes in Biological Media, ed. S. Becker, A. V. Kuznetsov, F. de Monte, G. Pontrelli, and D. Zhao, Academic Press, 2022, ch. 3, pp. 33–56.
36. L. D. Wittwer, F. Reichel, P. Müller, J. Guck and S. Aland, A new hyperelastic lookup table for RT-DC, Soft Matter, 2023, 19, 2064–2073.
37. N. H. Wubshet and A. P. Liu, Methods to mechanically perturb and characterize guv-based minimal cell models, Comput. Struct. Biotechnol. J., 2023, 21, 550–562.

---