Slow Sedimentation and Deformability of Charged Lipid Vesicles

The study of vesicles in suspension is important to understand the complicated dynamics exhibited by cells in in vivo and in vitro. We developed a computer simulation based on the boundary-integral method to model the three dimensional gravity-driven sedimentation of charged vesicles towards a flat surface. The membrane mechanical behavior was modeled using the Helfrich Hamiltonian and near incompressibility of the membrane was enforced via a model which accounts for the thermal fluctuations of the membrane. The simulations were verified and compared to experimental data obtained using suspended vesicles labelled with a fluorescent probe, which allows visualization using fluorescence microscopy and confers the membrane with a negative surface charge. The electrostatic interaction between the vesicle and the surface was modeled using the linear Derjaguin approximation for a low ionic concentration solution. The sedimentation rate as a function of the distance of the vesicle to the surface was determined both experimentally and from the computer simulations. The gap between the vesicle and the surface, as well as the shape of the vesicle at equilibrium were also studied. It was determined that inclusion of the electrostatic interaction is fundamental to accurately predict the sedimentation rate as the vesicle approaches the surface and the size of the gap at equilibrium, we also observed that the presence of charge in the membrane increases its rigidity.


Introduction
Lipid vesicles in the same size range as cells are widely used to investigate the physical properties of cell membranes, providing a simplified model system to investigate the role that the membrane plays in the mechanical deformation of cells. It is important to consider lipid membrane mechanics in order to understand many relevant processes involving cell deformation. For example, red blood cells (RBCs) undergo large deformations in the microcirculation [1] which appears to optimize gas exchange between blood and the surrounding tissues [2,3]. Given that unilamellar lipid vesicles share some features with RBCs, such as the ability to generate equilibrium biconcave shapes and dynamics under shear flow such as tank-threading and tumbling, the study of its behavior under flow conditions can provide useful information in the understanding of blood rheology [4][5][6].
Given the high stretching elastic modulus of lipid membranes, these are generally assumed to be incompressible (constant area per lipid molecule) [7][8][9]. As a consequence, the dynamics of vesicles under flow are sensitive to the value of the reduced volume, t~6 ffiffiffi p p V =S 3=2 (V is the volume of the vesicle and S its surface area) [10], which quantifies the amount of area available for the deformation of the vesicle. The maximum value t~1 corresponds to a sphere and values of tv1 correspond to shapes with an initial excess of area (as compared to a sphere).
In order to enforce local incompressibility a local isotropic tension, in the form of a Lagrange multiplier, can be added to the force density acting on the membrane [11][12][13][14][15][16]. The value of the Lagrange multiplier can be calculated approximately by considering a nearly incompressible elastic membrane and calculating the tension due to deformation, or by enforcing that the resulting flow field has zero divergence on the surface of the membrane. Either method neglects the effect of thermal fluctuations which is generally justified due to the relatively high stresses the hydrodynamics impose on the deforming vesicles.
Lipid membrane deformability influences the movement of vesicles close to a surface [17]. For example, it was shown that deflated vesicles under shear flow at a critical shear rate are lifted from the surface by means of a force of viscous origin. After lifting from the surface, the vesicles acquire the shape of prolate ellipsoids as they move in the direction of the flow. A direct numerical simulation (DNS) using a spectral boundary integral equation method was used to quantitatively determine the lift force acting on deflated vesicles near a wall [16]. It was demonstrated that lift velocity strongly depends on the vesicle's reduced volume and the viscosity ratio.
Despite the fact that the dynamic behavior of vesicles under external flow fields has been widely studied experimentally [18][19][20], computationally [13,16,21], and theoretically [12,22], sedimentation of lipid vesicles had received little attention. Recently, Boedec and coworkers studied the sedimentation of vesicles using computational [15] and theoretical [23] methods. In both cases the contribution of thermal fluctuations is neglected and sedimentation is studied in an infinite space. Experimental observations of pear-like shapes and microtether extrusion of sedimenting vesicles have also been reported [24].
Buoyancy-driven motion of drops has been widely studied in the fluid mechanics literature both experimentally and numerically. For example, the buoyancy-driven interaction of viscous drops was studied experimentally and numerically using the boundary integral method [25]. The effect of surfactants on buoyancydriven motion and interaction of viscous deformable drops was also studied using the boundary-integral method [26][27][28][29]. In all these works, the computational method has proven to be precise and in good agreement with both theoretical predictions and experimental data.
In this work we focus on the gravity induced sedimentation and deformation of initially quasi-spherical vesicles as they approach a flat horizontal no-slip surface. For this purpose we develop a computational algorithm based on the boundary integral method [30]. Our algorithm is mostly based on the method described by Zinchenko and coworkers [31,32] which was developed to study the interaction of deformable drops. More recently a similar algorithm has been used to study the dynamics of a vesicle in shear flow [10] and in a wall-bound shear flow [16].
The assumption of membrane incompressibility is based on the fact that direct area expansion is expensive from an energy point of view. Studies carried out by micropipette aspiration have shown that vesicles can undergo strain at low tensions. For example 1stearoyl-2-oleoyl-phosphatidylcholine (SOPC) vesicles have shown strains of up to 0.06 for applied tensions of 0.5 mN/m [33], this strain is accounted for by considering the smoothing of suboptical thermal fluctuations of the membrane in the low tension regime [34]. Therefore, even though changes in the lipid surface density are not expected to be significant, small membrane strain is likely as a result of smoothing of thermal undulations, which are accessible under low tension conditions. In other words, the thermal undulations act as an area reservoir which allows the membrane to deform under moderate forces.
In our model, membrane mechanics are dictated by the traditional Helfrich Hamiltonian [35] and we take into account the smoothing of thermal undulations by relating area strain to tension by the constitutive function proposed by Evans and Rawicz [34]. This equation considers area strain to be caused by the superposition of the smoothing of thermal undulations plus a direct expansion of the area per molecule. The latter requires much higher energy and is not available at low tensions, such as those exerted by gravity in our experiments.
Many experiments with lipid vesicles are performed in glucose/ sucrose solutions on glass substrates, using fluorescent dyes to label the membrane. These studies do not control for ionic conditions, and usually ignore the influence of electrostatic interactions between different components in their experimental setups. These electrostatic interactions may affect the behavior of a suspended vesicle, for instance in determining the sedimentation rate close to the wall and the gap between the vesicle and the glass surface at equilibrium. In our work electrostatic interaction between the membrane and the planar surface are taken into account by using the Derjaguin, Landay, Verwey, Overbeek (DLVO) theory of colloidal stability [36][37][38], where the linear Derjaguin solution for the Poisson-Boltzmann equation of the system was used to find the interaction force between the vesicle and the substrate.
By using single plane imaging microscopy (SPIM) we study the sedimentation dynamics of vesicles, this technique allows to follow a vesicle as a function of time with a better time resolution than a scanning technique such as confocal microscopy. We study the shapes and the size of the gap between the vesicle and the surface at equilibrium using confocal microscopy.
The objective of this work is to develop a robust computational model, which is able to accurately predict the behavior of charged vesicles sedimenting towards a charged flat surface. The model should contain the basic physics in order to be able to predict variables such as sedimentation rate as the vesicle approaches the surface and the gap between the vesicle and the glass surface at equilibrium.
Giant unilamellar vesicles were prepared via electroformation as described in the literature [39]. First, platinum wire electrodes inserted into Teflon wells are painted with POPC lipids dissolved in chloroform at 2 mg/mL (2-4 mL) and labeled with 1 mol% NBD-PE. The electrodes are then allowed to dry in vacuum overnight. 1 mL of sucrose solution (1 Molar) is then added to the wells and the electrodes are connected to an AC signal generator producing a sinusoidal 2 V peak to peak signal at 10 Hz for 1.5 hrs. 100 mL of solution containing vesicles is taken from the well, and the vesicles are then resuspended in 300 mL of a glucose solution with the same osmolarity as the sucrose solution (1 Molar). This sucrose-glucose system results in the precipitation of the vesicles due to density differences between the glucose and sucrose solutions.
The sedimentation of vesicles was studied through a homebuilt SPIM (UMR5088 Universite Paul Sabatier, Toulouse III and CNRS) [40]. First the glucose solution is placed into a quartz cuvette, followed by careful injection of the sucrose solution containing the vesicles. Once the sample is prepared, precipitating vesicles are localized with a 20X objective in air and a series of images are acquired through a CCD camera with a 15 seconds interval between each image.

Confocal Microscopy and Image Analysis
Vesicles at equilibrium were studied using confocal microscopy using an Olympus FV1000 (Tokio, Japan). Calcein was added to the glucose solution at a 1 mM concentration in order to improve contrast. The calcein provides the external solution with counterions at an approximately 2 mM concentration.
In order to measure the deformation of the vesicles, the confocal stacks were analyzed using an in-house algorithm by the following method: in each image plane, a threshold was applied, and the boundaries of the vesicle were detected using the function cvFindContours of the OpenCV library [41,42]. The cvFitEllipse2 function from the same library was then used on the detected contours to retrieve the dimensions and positions of the vesicles found in the given plane. A KMean algorithm was applied on the vesicle centers to track each vesicle from plane to plane throughout the stack. After this automated detection, the contour of the vesicle along the z axis is defined by its borders in each image plane. The Nelder-Mead simplex algorithm [42] is then used to fit the z contour of the vesicle with a function given by: The center of the polar coordinates in the z axis corresponds to the plane where the vesicle width is maximum. Best fit is obtained by finding the best values of parameters a and b and minimizing the least-square error. From this fit, area strain and curvature are determined.
To measure the equilibrium distance we waited for 15 minutes after putting the sample to start taking images of the vesicles. After waiting for this period of time no significant movement of the vesicles was observed. The whole experiment typically lasted for an hour. We performed these experiments for unlabeled POPC vesicles and for NBD-PE labeled POPC vesicles in solutions with no added salt and in solutions with 3 mM NaCl concentration.

Fluid Mechanics Model
The system considered in this paper is a giant lipid vesicle (diameter&20mm) immersed into a semi-infinite external fluid. The vesicle is settling due to a gravitational force towards a no-slip, infinite planar surface. A schematic of the system is shown in Figure 1. Fluid 1, which is internal to the vesicle, and fluid 2, in which the vesicle is suspended are both assumed to be Newtonian. The fluids are separated by a phospholipid bilayer membrane which introduces a jump in the stress field between both fluids. Gravitational acceleration g is assumed to act along the z-axis, normal to the planar surface. R 0 , is the radius of the quasispherical undeformed vesicle. r i (i~1,2) is the density of each fluid. The viscosities of fluids 1 and 2 are lm and m, respectively, where l is the viscosity contrast between the internal and external fluids.
All equations are presented in dimensionless form unless otherwise noted. Distances are scaled using R 0 . The characteristic time and traction are, respectively, where Dr~r 1 {r 2 and k b is the bending rigidity modulus of the membrane. The analysis presented in this work is based upon the creeping motion approximation in which the inertial terms in the equation of motion are neglected entirely for vanishing small Reynolds numbers. In this limit the Stokes equations can be used to model the flow of fluids 1 and 2. During sedimentation the membrane of the vesicle is a moving interface, whose position is unknown a priori. The fluid-structure interaction problem, the Stokes equations coupled to the membrane mechanics equations, can be solved using the boundary integral formalism [31,43] which yields in non-dimensional form.
where u is the velocity at each node over the membrane, k k~(l{1)=(lz1) and g 0~D rgR 4 0 =k b [44] is analogous to the Bond number. Experiments performed on pure POPC membranes have determined that k b ranges from 10 k B T to 40 k B T [44,45] where k B is the Boltzmann constant and T is absolute temperature. We have chosen the value of 40 k B T as the characteristic value to scale the experimental data and as a parameter in the computer simulations. The traction jump Df:(s (1) {s (2) ) :n n, wheren n is the unit vector normal to the surface, is determined from the configuration of the membrane. This term will be further discussed in the following section.
The Green's function kernels G and T for a wall bounded semiinfinite space were derived by Blake [46] and can be viewed as the fundamental solution to Stokes equations plus additional terms to account for the presence of the planar wall.

Membrane Mechanics
The traction jump Df introduced in Eq. 3 must include all external forces acting on the vesicle. For the problem under consideration Df~Df mem zDf grav zDf elec , where the terms on the right refer to the contribution of the membrane elasticity, the gravitational pull and the electrostatic interactions, respectively. Vesicle deformation at equilibrium results from the balance of all these terms.
There has been extensive theoretical research on lipid membrane deformation that involves the minimization of the free energy functional first proposed by Helfrich [7]. In general, it has been observed that the shape of a vesicle in suspension is limited by its resistance to bending, which is governed by the bending rigidity modulus of the membrane. The shape of the vesicle is determined by minimization of the shape energy which may be written as where C 1 and C 2 are the principal curvatures of the membrane. Both integrals are over the surface area of the membrane. The last term in Eq. 4 takes into account the constraint of constant area [35]. It is customary to assume membrane incompressibility based on the high energy cost of surface area dilation of the lipid membrane as compared to other modes of deformation. Most authors rely on approximations based on the calculation of membrane tensions that are proportional to membrane strain [10,14] which limit deformability of the membrane. Incompressibility can be imposed exactly by explicitly calculating the tension corresponding to the Lagrange multiplier by enforcing a divergence free velocity field on the membrane surface. This has been done for vesicles undergoing small deformations [23,47,48] and for vesicles in a wall-bound shear flow [16]. An algorithm for the simulation of three dimensional vesicle dynamics was developed by Boedec and coworkers [15].
Micropipette aspiration experiments have shown that vesicle strain can occur at low tensions [34]. It has been hypothesized that this strain in the low-tension regime results from the smoothing of thermal membrane fluctuations occurring in the suboptical range, which act as an area reservoir. Area dilation as a function of tension is thus given by [34] a~k where c is a coefficient that depends on the topology of the vesicle, S~sR 2 0 =k b is the effective reduced tension [44] and K ext~KA R 2 0 =k b is the dimensionless membrane's modulus for direct area expansion. The variable A is the dimensionless surface area of the membrane. Given the relatively large values of the membrane dilation modulus K A compared to k b , K ext is a large number of order 10 8 .
At low tensions, stretching of the membrane is ruled by the thermal energy (logarithmic term in Eq. 5) which is characterized by the bending coefficient, while at high tensions membrane strain involves surface dilation and is dominated by K ext (linear term in Eq. 5). While membrane strain due to direct surface dilation is difficult to achieve under normal flow conditions due to the relatively high magnitude of K ext , smoothing of thermal undulations are easily accessible at the low tensions that arise during gravity induced sedimentation.
In the present work we use Eq. 5 to calculate the tension on the membrane for a given area strain. Although we are modeling a dynamic process, the sedimentation is slow and will be regarded as a quasi-equilibrium process for which a uniform membrane tension can be assumed. We have performed simulations which have included the calculation of local tensions to enforce local incompressibility and have found the tangential stresses due to tension gradients to be at least two orders of magnitude smaller than the other forces acting on the membrane at any given time for the conditions being considered.
The contribution of the membrane elasticity to the force density can be obtained through the functional derivative of Eq. 4 as given by [35]. In non-dimensional form the Helfrich force can be expressed as Here, H~(C 1 zC 2 )=2 is the mean curvature, K G~C1 C 2 is the Gaussian curvature and + 2 s is the Laplace-Beltrami operator over the surface of the membrane. Since the contribution of the curvature term to the total energy will depend on the bending modulus, membrane deformation will depend on factors that influence this parameter such as lipid composition [45].
The gravitational pull is incorporated into the calculations through [30] Df where z is the vertical position of each differential element of the membrane with respect to a reference plane, in this case the glass surface, which has been made dimensionless by scaling with respect to R 0 . The calculation of the electrostatic contribution Df elec is explained in the following section.

Electrostatic Interactions
Experiments performed in the absence of salt result in long screening lengths. For this reason the electrostatic interaction of the vesicle with the surface will play an important role in determining the sedimentation dynamics of the vesicle and the equilibrium state. Microscope glass slides were used as substrates in the experiments. The glass, unless treated using methods such as silanization [49], will present a surface charge, where the main mechanism by which the glass surface acquires charge while in contact with water is by dissociation of the silanol groups [50]. The glass surface charge density is negative and has a value around s glass~{ 0:2 mC/m 2 for glass in water at pH 7.5 [51,52]. Vesicles are labeled with 1 mol% NBD-PE fluorescent probe, which has a negative charge and gives the membrane a negative surface charge density [53]. Hence, the electrostatic interaction between the vesicle and the glass will be repulsive. We can obtain an expression for the interaction force between the membrane and the glass surface from the linearized Poisson-Boltzmann (P-B) equation [36,37] Y is the electric potential and k {1 is the Debye length: k~2 where e is the fundamental charge, C is the ion concentration and e w is the dielectric constant of water. Equation 8 is valid for potentials less than 40 mV (about twice the thermal potential, k b T=e) [37]. In our first set of experiments no salts are added to the solutions, resulting in a Debye length of about 300 nm due to the presence of hydroniums and residual salts. When calcein is added to the solution the Debye length is lowered to 214 nm.The interaction between the vesicle and the glass surface can be modeled using the linear Derjaguin approximation which expresses that the interaction force between a small plane segment of the membrane and the plane surface of the glass is given by [36] f where h is the distance between the two surfaces and the potentials, for the glass and the vesicle are given, respectively, by: and Y ves~s ves R 0 (1zkR 0 )e w : The surface charge density of the vesicle can be estimated to be around s ves *{2:9 mC/m 2 and its size ranges from 1 mm to 50 mm in radius. Thus, the potential generated by the vesicle is high enough that the linearized P-B Eq. 8 cannot accurately model its behavior. However, for distances greater than the Debye length the behavior of the solutions of the linear and non-linear forms are very similar, provided the surface potential, Y ves , is replaced by a renormalized surface potential given by [54,55] which is independent of s ves and has a value of around 100 mV in our experimental setup. Substituting Eq. 13 into 10 for the vesicle potential, the expression to calculate the electrostatic repulsion becomes Combining Eqs. 11 and 14 and expressing in terms of nondimensional parameters, the contribution of the electrostatic interaction yields: Where we have introduced the dimensionless parameters k ie~R0 k and Note that in Eq. 15 the distanceh h~h=R 0 is expressed in dimensionless form dividing the distance by the characteristic length R 0 .
A total of six dimensionless parameters have been introduced, two from the fluid dynamics model,k k and g 0 , two from the electrostatic model, k ie and y ie , and two from the membrane mechanics model, K ext and the reduced effective tension S. The latter depends on the instantaneous configuration of the membrane and is recalculated at each time step.
When both the external and internal fluids have the same viscosity,k k~0 which simplifies Eq. 3 by canceling the second integral in the right hand side. Given that we found a very small difference between the viscosities of the sucrose and glucose solutions, we used this simplifying assumption in the solution of the model.
The remaining four dimensionless parameters are calculated in order to match the experimental conditions. The radius R 0 of the vesicle together with all other physical parameters determine the values of k ie , y ie , K ext and g 0 . A range of vesicle radii was chosen to match experimental condition.

Numerical Method
The general boundary integral method has been widely used to model the dynamic behavior of drops, capsules, vesicles and cells suspended in general flows. We refer to [31,43,56] for the general description of the method. However, in the following paragraphs we discuss particular considerations for the simulation of gravityinduced sedimentation.
To obtain the computational domain we use the method of uniform triangulation to discretize the initially spherical membrane into a uniform triangulated mesh. As it has been described by other authors [31] we begin with a regular icosahedron inscribed into the sphere. Each face of the icosahedron is divided into four smaller triangles by dividing each edge at its midpoint, the new vertices are projected radially onto the sphere and the process is applied recursively as many times as necessary to obtain the desired refinement. In the current work we do three refinements to obtain 642 vertices and 1280 triangular elements. An example of the computational mesh is shown in Fig. 2.
Integration of the right hand side terms in Eq. 3 are approximated by a simple surface trapezoidal rule that requires the integrands only at the triangle vertices [31]. This rule can be written, for any function g(x), in the following form, by reassigning the contribution of each triangular element to the vertex: where the summation in (18) is over all flat triangular areas DS with vertex x i . One of the greatest challenges in the numerical algorithm is the calculation of the curvatures and the Laplace-Beltrami operator in Eq. 6 over the discrete triangulated surface. A common algorithm used to calculate the mean and Gaussian curvature is to fit via least-squares a quadratic surface to each node and its neighbors and calculate the curvature from this function. From our tests we have determined that this method introduces sufficient error in the calculation of the curvature to be problematic during the calculation of the Laplace-Beltrami operator. For this reason we have decided to use the discrete differential-geometry operators as presented in [57].
The mean curvature estimate is derived from a discretization of the Laplace-Beltrami operator applied to the 1-ring neighborhood. Given a patch of triangles surrounding point x i as shown in Fig. 3, the estimates for the Gaussian curvature, K i and mean curvature H i , at x i , given by Meyer et al. [57] are: Voronoi area can be calculated by: Finally the Laplace-Beltrami operator of the mean curvature H can be estimated as [57] A similar method was recently used by Boedec and coworkers [15] for the calculation of the curvatures and Laplace-Beltrami operator.
The numerical algorithm advances in a relatively simple fashion and can be described as follows: from the shape of the vesicle the traction Df on the membrane is computed, the boundary-integral Eq. 3 can then be solved for the velocity at the nodes, finally nodes are advected via integration in time using a second order Runge-Kutta method. The algorithm continues until an equilibrium state is achieved.
An important difficulty in the boundary-integral calculations is efficient mesh control. Namely, if the nodes are simply advected with the flow, an initially regular mesh of triangles covering the surface becomes highly irregular after a short simulation time, thus invalidating the calculation. Mesh degeneration is especially problematic in gravity-induced motion given the large displacements of the vesicles. A mesh stabilization method can be developed by, instead of advecting the membrane nodes with the interfacial velocity u from Eq. 3, using the normal velocity (u :n n)n n plus an artificial tangential velocity field which can be constructed to maintain certain uniformity of the mesh. In the  current work we construct the tangential velocity field using the passive mesh stabilization algorithm first introduced by Zinchenko and coworkers [31]. The algorithm is based on the global minimization of the rate of change of the distances between neighboring nodes. To achieve this one must solve an optimization problem over the whole surface of the membrane. We found that by using such algorithm long simulations of gravity-induced sedimentation of vesicles were possible with reasonable time steps.

Results and Discussion
In this section the results from the computer simulations are presented and compared to the experimental data. First, the sedimentation dynamics are studied, then the equilibrium state of the vesicle is analyzed.

Sedimentation of Vesicles Towards a Flat Surface
The sedimentation of vesicles in isotonic conditions was studied using SPIM. The sedimentation process is driven by the density difference between the sucrose and glucose solutions (&35 kg/ m 3 ). The best six image series were analyzed. The criterion for this selection was to ensure that the vesicle did not present significant lateral displacement during sedimentation. Fig. 4 shows sample images from one of the series that were analyzed. The sedimentation rate from the images is calculated by a simple two point backwards finite difference scheme for the position of the centroid between successive images. All experimental results have a time lapse between images of 15 seconds.
Due to the optical configuration of the SPIM, an interference pattern results in an image artifact where some regions of the membrane appear to be brighter. This could be interpreted as these regions having a higher concentration of fluorescente probe. Nevertheless, given the preparation technique and the experimental configuration we do not expect this to be the case. We consider the assumption of uniformly distributed probe to be more accurate.
The figure also shows the sedimentation sequence as predicted by our computer simulations under similar conditions. Computationally, we calculate the sedimentation rate with a method similar to the one used experimentally, but the time between data points is much smaller. The simulations are initialized with the vesicle at a certain distance from the surface (between 6-8 radii). The position of the centroid of the vesicle is calculated at each time step and the sedimentation rate is calculated.
Sedimentation rates calculated experimentally and computationally are shown in Fig. 5. In the figure the distance to the surface is set to non-dimensional units by dividing by R 0 and velocity is also given in non-dimensional form. Far away from the surface, the simulations match the experimental results accurately. The sedimentation velocity decreases as the vesicle approaches the flat surface and the rate of change is well predicted by the simulations. At distances larger than 0:2 the electrostatic interactions are negligible and the whole process is dominated by the hydrodynamics.
As the vesicle approaches the wall, it decelerates. For uncharged vesicles (dark solid line in the figure) the velocity decreases exponentially approaching zero. This behavior is a result of the slow draining of the fluid that separates the vesicle and the surface, which corresponds to a lubrication regime [29]. When the electrostatic interaction is taken into account (light solid line in the figure) the sedimentation velocity decreases rapidly to zero further away from the surface avoiding this regime entirely. Close to the wall the electrostatic interactions become increasingly important causing a rapid deceleration, as the gravity-induced force is balanced by the electrostatic interactions.
Closer to the surface there appears to be a larger discrepancy between the experimental and computational results. Experimental data lies between the simulations with and without electrostatic interactions. We attribute this to the method used in determining the distance between the vesicle and the surface in the SPIM experiments, in which the gap between the vesicle and the surface is determined by calculating the distance between the centroids of the vesicle and its reflection (see Fig. 4) divided by two, minus the radius of the vesicle. This method does not take into account the deformation of the lower part of the vesicle, which induces an underestimation of the size of the gap at equilibrium.
One could also look for an explanation to this discrepancy in the assumptions made when developing our model. The most important simplification affecting the sedimentation rate is the assumption of having a uniform tension on the membrane. Given the symmetry of the problem, during the sedimentation process the membrane will become immobile and flow inside of the vesicle will be restricted. With a uniform tension, the membrane is mobile and there is the possibility of flow inside of the vesicle, as it happens in sedimenting droplets. This phenomenon will directly affect sedimentation rate. We have performed simulations, and compared the velocity of particles with mobile and immobilized interfaces under the same conditions and have found differences of less than 0.1 percent under these conditions. This difference is well below our experimental uncertainty. Also, this effect would be more important far away from the wall where the sedimentation rate is higher, opposite to our observations. Finally, in the following section we show that our model accurately predicts the equilibrium position of the vesicle as measured by the more precise confocal microscopy.
Once the forces acting on the vesicle balance out, the vesicle reaches equilibrium and the sedimentation rate approaches zero. At this point the vesicle acquires a shape determined by the equilibrium of the electrostatic and gravitational forces with the membrane resistance to deform.

Gap between the Vesicle and the Glass Surface at Equilibrium
The gap between the vesicle and the surface at equilibrium was further studied using confocal microscopy images. Confocal images and the corresponding frame from the computer simulation for similar conditions are shown in Fig. 6. The small gap between the vesicle and the surface is seen for both the confocal image and the computer simulation.
The gap at equilibrium measured experimentally and in the computer simulations are shown in Fig. 7 as a function of g 0 . The Debye length for NBD-PE labeled vesicles was calculated to be 214 nm when the ionic concentration of the calcein is accounted for. Computed simulations using that value for the Debye length show good agreement with the experimental results. The confocal microscopy combined with our image analysis algorithm provides a more accurate way of measuring the gap than the SPIM images. The distance to the surface decreases as g 0 increases due to the increase in the effective weight of the membrane as compared to the electrostatic force.
Experimentally, the addition of a small amount of salt results in the screening of the electrostatic forces and a decrease of the Debye length. By adding 3 mM of NaCl the Debye length is decreased to 5.5 nm. This condition can be used as a control, to verify that our observations are due to the electrostatic interactions. The experimental data and computer simulations under these conditions are also shown in Fig. 7. Overall, under these conditions, we observe a decrease of the equilibrium distance between the vesicle and the surface. The computer simulations predict a smaller gap distance than those measured experimentally, especially for smaller vesicles. Before commenting on the possible explanations for this behavior we note in passing that for larger vesicles the gap distance seems to converge for both Debye lengths being considered. This indicates that for heavier vesicles, which can reach closer to the wall, sedimentation is limited by the viscous draining of the fluid between the vesicle and the wall rather than the electrostatic interactions.
The gap size is determined through image analysis using a custom ImageJ [58] algorithm. The process consists in plotting the intensity profile along a line that crosses the water gap between the vesicle and the glass surface. From the plot, the gap size is defined as the distance between points with an intensity corresponding to 70 percent of the highest measured value. Several distances measured for the same vesicle (between 3 and 5) are then averaged. The error bars shown in figures 7 and 8 are obtained from a combination of the deviation from the mean of the measured values with a confidence interval of 95%, together with the bias error due to the resolution of the microscope configuration used (around 0.29 microns).
Our scaling normalizes the weight of all vesicles. It can be observed in Fig. 5 that sedimentation rate is independent of the size of the vesicle in non-dimensional terms. Physically, sedimentation rate does depend on the size of the vesicle, as heavier vesicles sediment faster. From the characteristic length R 0 and the characteristic time defined in Eq. 2 it can be seen that physical velocity varies with R 2 0 . That is, vesicles with radius of order one are sedimenting a hundred times slower that those of radius 10! Due to the stability of the vesicles, experiments are limited to one hour. With this in mind, a better explanation to the observed behavior is that smaller vesicles have not yet reached equilibrium  and are sedimenting very slowly when the measurement is performed. Computationally, we are allowing the vesicles to reach equilibrium defined as the state when displacements are below the numerical resolution of our algorithm, thus the measured gap is smaller.
Overall, the computational model is able to predict in an accurate fashion the size of the gap at equilibrium for a wide variety of conditions. The agreement in the results also indicates an appropriate selection of the values of some important parameters that were not measured in our experiments but relied on calculated values and theoretical predictions, such as the glass surface charge density and the counterions concentration in the solutions.
So far we have presented results corresponding to charged vesicles labeled with NBD-PE in solution. We have considered solutions with no added salt where electrostatic interactions are important (relatively large Debye length) and with salt added to screen the effect of electrostatic interactions. The addition of a small amount of salt (3 mM) generates a slightly hypertonic condition which causes the vesicle to deflate. In order to investigate the effect of this on the measured gap we perform additional experiments with unlabeled vesicles (no charge). The gap distance for charged vesicles in a solution with added salt compared to unlabeled vesicles is shown in Fig. 8. It can be seen that there is no difference in the behavior of both groups and that the same behavior discussed above is observed in unlabeled vesicles. In the following section we investigate the deformation of both charged and uncharged vesicles.
The thorough understanding of the role played by the electrostatic interactions on the sedimentation process, as well as on determining its position at equilibrium might be important in understanding the formation and stability of certain colloidal systems. The results presented in this section will play an important role in this understanding.

Vesicle Deformation at Equilibrium
When the vesicle reaches equilibrium, the sedimentation process stops and the vesicle adopts an equilibrium shape determined by the balance of the forces acting on it. The shape of the vesicle at equilibrium was studied through confocal microscopy and image analysis as explained in the Materials and Methods section. The analysis was applied to images similar to that shown in Fig. 6 (A). Contrast was improved by staining the glucose solution with low concentrations of calcein (1 mM). At these low concentrations, it should not influence the osmotic conditions of the solution. It has been shown [34] that membrane deformation is caused by the superposition of two modes: smoothing of thermal undulations and direct stretching of the membrane. At low external forces, such as those in the present experiments, it is expected that strain be due to the former mode only. Tension on the membrane at equilibrium is proportional to the value of g 0 . By plotting the area strain as a function of g 0 in a semi-log plot (Fig. 9) it is evident that the logarithmic term in Eq. 5 is the dominating term in defining the shape of the membrane at equilibrium as expected.
In Fig. 6 experimental images of deformed vesicles at equilibrium were shown together with the computational prediction at similar conditions. It is interesting to note that those relatively large deformations are achieve with very small area strains (less than 0.01 in all cases), which can be explained due to the smoothing of thermal undulations [33]. Hence, it becomes apparent that considerations of thermal undulations is important in the study of vesicle sedimentation.
The tension-strain relationship (Eq. 5) suggests that strain at a certain tension depends on the membrane bending rigidity. It has been suggested previously [59][60][61] that electrostatic interactions between neighboring NBD-PE molecules should provide a contribution to the bending rigidity. In order to verify this observation we compare the strain of charged and uncharged vesicles.
As mentioned above, salt produces a slightly hypertonic condition which causes the vesicle to deflate, this will certainly affect the deformability of the particles. For this reason in the following we report the deformation of both labeled and unlabeled vesicles suspended in the same solution with no salt added. In the previous section it was shown that not labeling the vesicles had the same effect on the measured gap as adding a small amount of salt to screen the electrostatic interactions. Fig. 9 shows the measured area strain for both NBD-PE labeled and unlabeled vesicles. The simulation results for vesicles with a bending rigidity modulus k b~4 0k B T is also shown in the figure.
For the majority of the unlabeled vesicles the measured strain is in line with the prediction of the computational model. For a smaller fraction of the data points, mostly lighter vesicles, a much smaller strain is measured. All data points correspond to the same sample, hence tonicity effects should be ruled out. The most probable cause for this observation is, as explained above, that this vesicles have not yet reached an equilibrium state and have not yet fully deform to balance the gravitational force.  Charged vesicles undergo a smaller strain at equilibrium. Solid triangular data points in Figs. 7 and 9 correspond to the same vesicles. Based on the comparison between the experimental data and the computer simulations, we consider that these vesicles have reached equilibrium, and electrostatic forces are balanced with the gravitational force. These results suggest that the presence of charged fluorescent probes does have the effect of increasing the bending rigidity as has been reported previously [59].
Measurement of the mechanical properties of vesicles represents one of the possibilities presented by the current computational model. By comparing good experimental data with the computer simulations, one might be able to determine the mechanical properties of membranes which could complement techniques such as micropipette aspiration.

Conclusions
We have reported the implementation of a computational algorithm for the simulation of lipid vesicles in suspension. Experiments using SPIM and confocal microscopy were performed to verify the performance of the algorithm. The algorithm is used to simulate the sedimentation of vesicles due to gravity, taking into consideration the electrostatic interactions between the membrane and a glass surface towards which it is sedimenting. We have shown that our simulations are able to predict, with reasonable accuracy, the rate at which the vesicles sediment and the velocity variations as the vesicle approach the glass surface. The algorithm also predicts the equilibrium gap between the vesicle and the surface at equilibrium. It was shown that the consideration of the electrostatic contribution to surface interactions is essential in order to accurately predict the sedimentation rate, especially at close range from the surface, and the fluid gap between the vesicle and the surface at equilibrium. We have shown that by modeling the mechanical behavior of the membrane with Eq. 5, which superposes the deformation due to the smoothing of thermal undulations and direct membrane stretch, area strains are reasonably small and in line with those observed experimentally. The model used in our simulations also takes into account the intrinsic mechanical properties of the membrane through the membrane bending modulus, which can be modulated in cells through membrane composition. For example, by adding 30% cholesterol the bending rigidity of the membrane can be doubled [62]. It is also suggested that charged fluorescent probes on the membrane have the effect of increasing its bending rigidity.
Our vesicle model has considered the basic physical contributions: Stokes flow, bending rigidity, membrane deformation through smoothing of thermal undulations but limited by the high energy cost of direct membrane deformation and electrostatic interactions between the vesicle and the glass surface. In its current form the algorithm can be used to study the effect that the different physical parameters discussed have on the sedimentation rate, distance to the surface at equilibrium and membrane strain. It can also be extended in order to study important phenomena such as colloidal stability, and biofilm formation.
Finally, regarding future work, we consider our current model to be robust and accurate enough to allow us to investigate the roll that the physics described and studied in the current paper play in the slow motion of a suspended vesicle. Nevertheless, the dynamics and motion of vesicles under shear or general flow pose increasing challenges as the calculation of local tensions becomes mandatory. The method proposed by Boedec et al. is a very good starting point, but the inclusion of viscosity contrasts different than one and other physics, such as those investigated in the current paper, could prove challenging and to the best of our knowledge remain unsolved.