## Abstract

Cryopreservation requires that stored materials be kept at extremely low temperatures and uses cryoprotectants that are toxic to cells at high concentrations. Lyopreservation is a potential alternative where stored materials can remain at room temperatures. That storage process involves desiccating cells filled with special glass-forming sugars. However, current desiccation techniques fail to produce viable cells, and researchers suspect that incomplete vitrification of the cells is to blame. To explore this hypothesis, a cell is modelled as a lipid vesicle to monitor the water content and membrane deformation during desiccation. The vesicle is represented as a moving, bending-resistant, inextensible interface and is tracked by a level set method. The vesicle is placed in a fluid containing a spatially varying sugar concentration field. The glass-forming nature is modelled through a concentration-dependent diffusivity and viscosity. It is found that there are optimal regimes for the values of the osmotic flow parameter and of the concentration dependence of the diffusivity to limit water trapping in the vesicle. Furthermore, it is found that the concentration dependencies of the diffusivity and viscosity can have profound effects on membrane deformations, which may have significant implications for vesicle damage during the desiccation process.

## 1. Introduction

Preservation of biological material can be achieved by reducing the molecular mobility within cells to a specified level. The lack of mobility, of reactants in particular, inhibits the chemical processes that would otherwise deplete cellular metabolites during storage. Cryopreservation achieves this by keeping the internal cellular temperature well below the freezing point of the intracellular fluid, effectively reducing the mobility to zero. However, the cryoprotectants that are necessary to prevent ice nucleation within the cell during the freezing process are toxic at high concentrations, as observed by Fahy *et al*. [1]. This, in addition to the logistical issues with maintaining cryopreservation temperatures, has resulted in recent research into alternative preservation techniques. Lyopreservation is one such alternative and seeks to duplicate a waterless hibernative state used by certain organisms such as brine shrimp and certain nematodes. The resulting state, called anhydrobiosis, is described in detail by Crowe *et al*. [2]. Achieving this state is thought to depend on the presence of a ‘glass-forming’ sugar both in and around the organism. The ‘glass-forming’ property describes how large concentrations of solute cause the viscosity of the surrounding water to reach that of a glass (10^{13} Pa s), described by Aksan *et al*. [3]. Because the mobility within a glass is sufficiently low, these organisms can preserve themselves by raising their sugar concentration through the shedding of all their water, and letting the naturally occurring, glass-forming sugar induce anhydrobiosis.

Because viscosity is thought to be key to achieving anhydrobiosis, studies have been conducted on various glass-forming sugars, including sucrose, fructose and trehalose, including one by Aksan *et al*. [3]. Trehalose is found to have the strongest viscosity response, and thus further studies have investigated the drying of mammalian cells with trehalose inside and out. A spin-drying technique, recently developed by Chakraborty *et al*. [4], is able to produce samples with a water fraction near that required for glass formation. Unfortunately, the mammalian cells do not survive rehydration after the desiccation. Chakraborty *et al*. [4] believe that incomplete glass formation within the cells is the cause; however, this cannot be measured directly due to limitations in experimental observation techniques. Therefore, modelling of the drying process can provide insight in place of experimental measurements, including what aspects of glass-forming sugars are beneficial and what approaches are optimal for proper preservation.

One biological material of interest for lyopreservation is the red blood cell, as reviewed by Scott *et al*. [5]. If blood can be preserved at room temperature with lyopreservation, hospitals would not require costly back-up power systems to keep cryopreserved blood frozen. Understanding the drying dynamics of red blood cells is necessary for obtaining this storage. Because the underlying force driving these dynamics is the flow of water through the cell membrane, and the membranes of lipid vesicles and the membranes of red blood cells exhibit similar osmotic behaviour, described by Lipowsky & Sackmann [6] and Menger & Angelova [7], the red blood cells can be modelled as vesicles. Additionally, both vesicle and red blood cell membranes exhibit a resistance to bending; thus how a vesicle deforms during the drying process can provide insight into how red blood cells might react.

Lipid vesicles have been extensively studied under flow, including the work of Blount *et al*. [8], Cantat & Misbah [9], Schwalbe *et al*. [10], Salac & Miksis [11,12], Veerapaneni *et al*. [13], Zhao *et al*. [14] and Zhao & Shaqfeh [15]. By contrast, much less work has been done on modelling vesicles in a solute–water solution, particularly with a flow of water across the membrane due to the solute. Fettiplace & Haydon [16] empirically measured the dependence of the flow on the solute concentration and solution pressure. Both Chen & Keh [17] and Nardi *et al*. [18] studied vesicle drift from a solute gradient using a spherical vesicle assumption. Bealle *et al*. [19] observed that the loss or gain of water leads to buckling of the membrane, owing to inextensibility, and potential damage. Blount *et al*. [8] asymptotically modelled membrane deformation with a spatially homogeneous solute concentration. While each of these studied behaviours contribute to the drying dynamics in lyopreservation, a study coupling the effects of vesicle deformation, dynamic solute gradients and proper glass formation is key to avoiding cell damage.

To capture these effects, a lipid vesicle immersed in a sugar–water solution is modelled, where the sugar is glass forming. The adaptable Cartesian-grid, level set, immersed boundary method of Salac & Miksis [11,12] is used to numerically solve for the vesicle shape and fluid flow. The Lagrange multiplier used by Salac & Miksis [11,12] to enforce the area constraint is also used; however, the resulting surface partial differential equation is treated differently by a method described herein. The sugar concentration jump across the membrane is handled by the immersed interface method developed by LeVeque & Li [20]. With this computational machinery in place, the vesicle desiccation is simulated in two dimensions for various drying, membrane and sugar parameters.

## 2. Problem formulation

A two-dimensional lipid vesicle is modelled as a closed interface *Γ*(*t*) immersed in a sugar–water solution. The solution is assumed to be Newtonian with stress , where *p* is the pressure, *μ* is the viscosity and **u** is the solution velocity. The solution is also assumed to be incompressible so that ∇ · **u** = 0. Thus, the solution satisfies
2.1where *ρ* is the density of the solution. The density is spatially varying and given by
2.2where *ρ*_{water} is the density of water, *ρ*_{sugar} is the density of sugar and *c* is the sugar mass concentration (mass per unit volume) in the solution. Following the work of Seifert [21], the membrane force on the solution is derived from the Helfrich free energy , where *κ*(**s**, *t*) is the curvature, *b* is a material bending resistance parameter and *γ*(**s**, *t*) is a quantity that accounts for the membrane tension. This results in the following jump in normal stress:
2.3where ∇_{s} = (*I* − **n ⊗ n**)∇ and Δ_{s} = ∇_{s} · ∇_{s}.

The vesicle membrane velocity is denoted **u*** _{Γ}*. The inextensibility of the vesicle membrane is incorporated by enforcing , where is an arbitrary section of the membrane. This implies that
2.4everywhere along the membrane. There is no tangential slip between the vesicle membrane and the solution: (

**u**−

**u**

*) ·*

_{Γ}*τ*= 0. However, there can be flow through the membrane due either to hydrostatic stress from a pressure jump across the membrane or to osmotic stress from the presence of the sugar. Denoting the volume flux of this flow as

*Q*(

**s**,

*t*), this gives (

**u**−

**u**

*) ·*

_{Γ}**n**=

*Q*and therefore

**u**

*=*

_{Γ}**u**−

*Q*

**n**. The form for the volume flux

*Q*(

**s**,

*t*) is determined empirically by Fettiplace & Haydon [16], discussed by Blount

*et al*. [8], and given as 2.5Here,

*k*

_{o}is the osmotic flow parameter associated with the jump in sugar mass concentration

*c*(

**x**,

*t*) across the membrane and

*k*

_{h}is the hydrostatic flow parameter.

The sugar concentration itself undergoes standard advection–diffusion throughout the solution
2.6The presence of sugar affects both the viscosity of the solution and the diffusivity of the sugar. Galmarini *et al*. [22] measured the viscosity of trehalose–water and sucrose–water solutions at various concentrations. The relationship between viscosity and concentration shows an exponential dependence and motivates the following expressions for both the viscosity in (2.1) and the diffusivity in (2.6):
2.7and
2.8where *μ*_{0} and are the viscosity and diffusivity values when the sugar–water solution is near pure water. The expression for the diffusivity is motivated by in the Stokes–Einstein relation. It should be noted that the measurements taken by Galmarini *et al*. [22] involve a different non-dimensional concentration (*m*_{sugar}/*m*_{water}) from that used here. However, the results obtained using their form of concentration for the viscosity and diffusivity are qualitatively identical to those using (2.7) and (2.8). To complete the system for the sugar concentration, there is zero flux of sugar enforced at the vesicle membrane to model the impermeability of mammalian cell membranes to trehalose:
2.9and
2.10The domain is considered to be finite with three solid walls that are impermeable (figure 1). Thus, the solution velocity along the solid walls is zero (**u** = 0), as well as the flux of sugar . The boundary condition at the remaining side is responsible for the evaporation of water
2.11where *M*_{w} and *M*_{s} are the molar weights of water and the sugar, respectively. Note that, as the domain contains less water (*c* → *ρ*), the right-hand side of (2.11) tends to zero. This captures how the rate of water evaporation will slow as there is less water present. The right-hand side of (2.11) is also equal to zero if there is no sugar in the domain (*c* → 0). Without any sugar present, the sugar concentration itself should not rise (the flux of sugar into the domain should be zero). Further discussion over the choice of (2.11) can be found in appendix A.

The system is solved in dimensionless quantities, where the characteristic length is the vesicle radius *R*_{0} = *L/*2*π*, with *L* denoting the length of the vesicle membrane. The diffusion rate of the sugar is used for the characteristic speed , and the density of the sugar is used for the characteristic mass . This results in the dimensionless parameters found in table 1. The dimensionless equations themselves can be found in appendix B.

Because the problem involves many coupled and nonlinear equations, the solution is obtained numerically. The location of the membrane is first determined using a level set method as implemented by Salac & Miksis [11,12]. The level set is updated using the same Lagrangian update and reinitialization schemes, with the addition of velocity extension to embed the membrane velocity **u*** _{Γ}* into the domain. The jump in normal stress is then enforced using the immersed boundary method, as was done by Salac & Miksis [11,12]. Similar to their method, the curvature

*κ*(

**s**,

*t*), tension

*γ*(

**s**,

*t*) and their derivatives must be extended off the membrane. While the extension is similar to the closest-point method described by Ruuth & Merriman [23], the approach here avoids making the approximation that ∇

_{s}= ∇ near the membrane. Because the tension is responsible for enforcing the inextensibility of the membrane, it is important to minimize approximation errors when calculating it. Thus, the surface-gradient operator is discretized directly. The jump in sugar concentration along the interface is handled with the immersed interface method developed by LeVeque & Li [20]. The resulting method is found to be first-order accurate overall in both time and space. The interior sugar mass and the membrane length of the vesicle are conserved to within a relative error that is less than 1%, on average, and at most 3%. Full details of the numerical method are presented by Vogl [24], with the appropriate terms added to the method here to account for the spatially varying density

*ρ*(

**x**,

*t*). A summary of the convergence analysis by Vogl [24] can be found in appendix C.

## 3. Results and discussion

In order to investigate the effects of glass-forming sugars on vesicle morphology and water content during desiccation, results are obtained for a vesicle in a sugar–water solution following (2.1)–(2.3), a sugar concentration governed by (2.6), (2.9)–(2.11), and the interaction between the two governed by (2.5), (2.7) and (2.8). The material parameters are chosen by considering trehalose as the sugar and are found in table 2. These values give the following for the non-dimensional quantities: , and . Because the jump in pressure across the membrane has such a small contribution to the flow of water through the membrane , the value for is taken to be zero.

The remaining parameters *A*_{D}, *A*_{V}, and are varied to observe their effect on the drying process. The choices of initial vesicle shape and sugar concentration are motivated by the equilibrium shapes. Equilibrium vesicle shapes in two dimensions can be computed by minimizing the Helfrich free energy , as was done by Veerapaneni *et al*. [28]. These shapes depend only on the reduced area *v* = 4*πA/L*^{2}, where *A* is the vesicle area and *L* is the vesicle membrane length. Figure 2 shows how the equilibrium shapes start as a circle for *v* = 1.0 and become more and more biconcave as *v* decreases from 1.0. Note that, for reduced areas below 0.3, the equilibrium shape crosses itself and becomes non-physical. In three dimensions, Seifert *et al*. [29] find that the corresponding equilibrium oblate spheroids approach a similar self-intersecting shape for reduced volumes near 0.5, which is equivalent to a reduced area of approximately 0.32 (using from Salac & Miksis [12]). However, for reduced volumes near or below 0.5 in three dimensions, the equilibrium shapes become stomatocytes. The stomatocyte shape results from the balancing of the two principal curvatures in three dimensions, thus it is not an equilibrium shape in the two-dimensional model here. As such, the self-intersecting vesicle is avoided here by keeping the reduced area above 0.3. The initial shape is chosen to be a circle where the reduced area is maximized (*v* = 1.0). The impermeability of the vesicle membrane to the sugar ensures that sugar mass in the vesicle is conserved: (*cA*)_{final} = (*cA*)_{initial} or, equivalently, (*cv*)_{final} = (*cv*)_{initial}. Thus, the initial concentration of the domain is set to *c* = 0.4 so that the final reduced area of the vesicle is *v* = 0.4.

### 3.1. Glass formation

To observe the desiccation progression of the domain and vesicle, the value *c** = 0.6 is chosen to represent a drying front. The time it takes the drying front (*c* = *c**) to travel through the vesicle, so that *c* > *c** everywhere in the vesicle, is denoted the ‘vesicle drying time’. The time it takes the drying front to travel through the domain, so that *c* > *c** everywhere outside of the vesicle, is denoted the ‘domain drying time’. First, the drying number is varied while holding *A*_{D} = 0.0, *A*_{V} = 0.0 and fixed. Figure 3 shows how increasing speeds the drying, leading to lower drying times both for the domain and for the vesicle. However, this effect diminishes for larger , as the drying becomes limited by the diffusion through the domain.

In order to study the effect of the osmotic number, the values *A*_{D} = 0.0, *A*_{V} = 0.0 and are fixed, while is varied. Figure 4 shows how the increased flow of water out of the vesicle lowers the vesicle drying time with increasing . The domain drying time remains essentially flat with increasing as the flow of water out of the vesicle slows the drying of the domain only slightly. Unlike increasing , increasing causes the difference between the domain drying time and vesicle drying time to decrease significantly. This suggests the possibility of a value of where the vesicle dries before the domain.

The value of the glassy viscosity *A*_{V} does not have a significant effect on either drying time for the parameter values explored here. This is because the convective contribution of the solution flow is minimal and has a negligible effect on the sugar distribution, and therefore the drying times do not significantly depend on the solution viscosity. Thus, the next range of parameters explored is *A*_{V} = 0.0, and , with varying glassy diffusivity parameter *A*_{D}. Figure 5 shows that increasing *A*_{D}, which decreases the sugar diffusivity throughout the domain, slows down the drying of both the domain and the vesicle. An interesting observation here is that the vesicle dries faster than the domain for large enough *A*_{D}.

This suggests that the osmotic number and glass diffusion *A*_{D} values determine whether the domain or vesicle dries first. The order of drying can have significant implications in cells undergoing desiccation. Specifically, the values of these parameters might determine whether or not undesirable water content remains in the cell after the desiccation procedure is stopped. Thus, the model is simulated for a variety of and *A*_{D} values to find regimes where the vesicle dries first, limiting the trapping of water. The results are summarized in figure 6, where the solid line divides the (, *A*_{D}) parameter space into two regions: one where the vesicle dries first and one where the domain dries first. Figure 6 also shows how larger values of or *A*_{D} encourage faster drying of the vesicle. Figure 7*a* shows the water content of the vesicle almost matching that of the surrounding domain, at the domain drying time, for and *A*_{D} = 2.5. The large glassy diffusion value means lower overall diffusivity; thus, the drying of the domain is slowed enough to allow the water content in the vesicle to equilibrate. This is in contrast to figure 7*c*, with *A*_{D} lowered to 0.5, where the water content in the vesicle is much higher than that of the domain. A similar effect is observed if *A*_{D} = 0.5 and is decreased from 0.08 in figure 7*d* to 0.04 in figure 7*c*. Whereas the water content of the vesicle in figure 7*d* may not match the surrounding domain as well as figure 7*a*, the phase diagram in figure 6 suggests that a further increase in will result in even faster vesicle drying. Alternatively, the value of *A*_{D} can be increased from 0.5 in figure 7*d* to 2.5 in figure 7*b* to increase vesicle drying.

These results suggest that higher values of and *A*_{D} are beneficial to the desiccation process in that water trapping is less likely. However, the deformation of vesicles during desiccation is also important as high stress on the vesicle membrane is undesirable. Noting the different vesicle shapes in figure 7, it is clear that the osmotic number and glassy diffusivity *A*_{D} have important influence on vesicle deformation. The influence of these material parameters, as well as other parameters, on the vesicle shape during desiccation is explored in the next subsection.

### 3.2. Drying shapes

During drying, the vesicle area decreases from the loss of water. For the vesicle to remain a circle, the membrane length would also have to decrease appropriately (so that *v* = 1.0). However, the inextensibility of the membrane keeps the membrane length constant. This excess membrane length (EML) prevents the vesicle from remaining a circle and causes the membrane to buckle. The EML is then redistributed, owing to the bending resistance, to be more energetically favourable. This is analysed by decomposing the excess membrane into Fourier modes. Given the membrane position **r**(*s*, *t*), the EML is defined as *e*(*s*, *t*) = *|***r**(*s*, *t*) − **x*** _{c}*(

*t*)

*|*−

*R*(

*t*), where

**x**

_{c}is the centre of mass of the vesicle and

*R*(

*t*) is the radius of a circle with the current area of the vesicle. The amplitudes of the Fourier modes of

*e*(

*s*,

*t*) are then computed. As an example, the dominant modes of a vesicle with ,

*A*

_{D}= 0.0,

*A*

_{V}= 0.0 and are shown in figure 8. The figure indicates that the succession of dominant modes of EML is

*k*= 8,

*k*= 4,

*k*= 3 and

*k*= 2. The maximum amplitude of the dominant mode increases from mode to mode due to the continued generation of EML over time. Figure 9 shows that a majority of the area loss occurs while the mode

*k*= 4 is dominant.

The vertical lines in figures 8 and 9 identify a physical shape change between *n*-star shapes, defined as having *n* local maxima in curvature that visually appear as *n* nodules in the shape. Specifically, the lines denote a change to 8-star, 4-star, 3-star and 2-star shapes. In figure 9, the lines for the 3-star and 2-star shapes lie almost on the top of each other because these shape changes occur so late in the desiccation process. They also show that succession of shapes does not immediately occur at the same time as the succession of dominant nodes. The rising node has to first overcome the falling node before the resulting shape change occurs. Figure 10 shows the shapes corresponding to the maximum amplitudes for the dominant modes. Because these maxima occur after the rising node has overtaken the falling node, the shapes are the 8-star, 4-star, 3-star and 2-star shapes. The break in top-to-bottom symmetry of the shapes is due to the drying from one end. Note that similar shapes have been observed by Veerapaneni *et al*. [28] at various energy states.

In order to understand the effect of the osmotic number on the succession of dominant modes and vesicle shapes, the values for are varied while the values for , *A*_{D} = 0.0 and *A*_{V} = 0.0 are held fixed. The results in figure 11 show the amplitudes of the dominant modes *k* = 8, *k* = 4 and *k* = 3. The mode *k* = 2 has been left out as it only becomes dominant after the vesicle is essentially dry. The figure shows that increasing increases the maximum amplitude of the *k* = 8 mode. This is because more EML is generated at the beginning of the desiccation process due to the increased area loss with a higher osmotic flow. This additional EML feeds the mode *k* = 8 before the bending resistance can transition the EML to the mode *k* = 4.

The increasing of values has a more subdued effect on the modes *k* = 4 and *k* = 3. While the increased osmotic flow feeds the mode *k* = 4 in the same way as it feeds the mode *k* = 8, the effect is lessened because the mode *k* = 4 extends late into the desiccation process. Because there is less water inside the vesicle late in the process, there is lower osmotic flow in general, and therefore a weaker effect of the osmotic number . Thus, the maximum of mode *k* = 4 is relatively unaffected by changes in . The same holds for the mode *k* = 3, as this mode becomes dominant even later in the desiccation process than the mode *k* = 4.

The vesicles shapes corresponding to the maximum amplitudes of the modes *k* = 8, *k* = 4 and *k* = 3, as well as the shapes corresponding to the maximum bending energy , are compared to further explore the effect of . The shapes in figure 12*a* show an increasingly defined 8-star shape as is increased. This mirrors the increase in the mode *k* = 8 amplitudes from additional EML. The mode *k* = 4 in figure 12*b* and mode *k* = 3 shapes in figure 12*c* show the effect of greater top-to-bottom asymmetry with increasing . The vesicle translates down as the asymmetry increases and causes the EML to migrate towards the bottom of the vesicle. The maximum bending energy shape in figure 12*d* shows how increasing causes high curvature regions in the top nodules of the vesicle. Thus increasing may help eliminate trapped water in a vesicle; it can also lead to vesicle membrane damage through these more pronounced nodules.

Increasing the drying number has a similar effect to that of increasing . Higher values of increase the initial EML leading to a rise in the mode *k* = 8 amplitude. The result on the mode *k* = 4 and *k* = 3 amplitudes is negligible as with . The reason is the same: these modes occur later in the desiccation process when the drying number has a lessened effect. The 4-star and 3-star shapes show similar increases in top-to-bottom asymmetry and more pronounced top nodules. Thus, increasing the drying number may speed up the drying process, but comes with the same drawback of increased likelihood of membrane damage as .

To understand the effect of the glassy diffusivity *A*_{D} on the vesicle shape, the values for *A*_{D} are varied while , *A*_{V} = 0.0 and *K*_{o} = 0.05 are fixed. Figure 13 shows how the vesicle transitions through the modes significantly earlier in the desiccation process for higher *A*_{D}. This is because the decrease in diffusivity from high *A*_{D} slows the flow of water out of the vesicle. The resulting slower generation of EML allows the bending resistance of the membrane to relax the EML to a lower mode earlier in the drying process. Thus, the maximum amplitude of all the modes is decreased with higher *A*_{D}. Additionally, the vesicle shapes are all smoother with higher *A*_{D}, as shown in figure 14. The maximum bending energy shapes in figure 14*d* for higher *A*_{D} are much less pronounced, as well, owing to the earlier transitioning of modes. Therefore, increasing *A*_{D} not only decreases the possibility of water trapping, but also decreases vesicle deformation during drying.

To understand the effect of the final parameter on the vesicle shape, the values for glassy viscosity *A*_{V} are varied while , *A*_{D} = 0.0 and *K*_{o} = 0.05 are fixed. Figure 15 shows that increasing *A*_{V} slows down the relaxation of all the modes in the desiccation process. This is because the increased viscosity of the surrounding solution resists the force from the bending resistance and inextensibility of the membrane. Thus, the maximum mode *k* = 8 amplitude increases with greater *A*_{V}, because the vesicle is held in the 8-star shape further and further into the desiccation process. Similar to increasing , the mode *k* = 4 does not exhibit an increase in maximum amplitude with increasing *A*_{V}. Again, this is because the maximum amplitude occurs late in the desiccation process where almost all the area loss has occurred. The mode *k* = 3 amplitudes, however, do experience a significant decrease in maxima with increasing *A*_{V}. This is because the vesicle is held in the 4-star shape so long into the desiccation process with larger *A*_{V} that the vesicle transitions almost immediately to the mode *k* = 2 from the mode *k* = 4.

The delaying of the mode relaxation with increased *A*_{V} is also evident in the vesicle shapes. The longer dominance of the mode *k* = 8 results in the stronger defined nodules shown in figure 16*a*. Additionally, the persistence of the 8-star shape leads to different orientation for *A*_{V} = 1.0 and *A*_{V} = 2.0 than seen previously. Whereas the top-to-bottom dried vesicle has, with previous parameters, preferred the 8-star shape with nodules aligned off the *x-* and *y*-axes, here the vesicle aligns the nodules with the axes. This results in different 4-star shapes for the different *A*_{V} values, as shown in figure 16*b*. Note the more pronounced nodules with the higher *A*_{V} 4-star shape than with the *A*_{V} = 0.0 4-star shape. Figure 16*d* shows that these 4-star shapes occur when the bending energy is the highest, and thus is when the likelihood for membrane damage is the highest. As mentioned earlier, the 3-star shape never appears with higher *A*_{V} values because the vesicle transitions almost directly to the 2-star shape, with the maximum mode *k* = 3 amplitude occurring before the shape change. Regardless, the high deformation in the 4-star shape corresponding to high values of the glassy viscosity *A*_{V} can lead to unwanted stress during desiccation.

## 4. Conclusion

A model is developed to explore the behaviour of vesicles undergoing desiccation with a glass-forming sugar present. The vesicle is immersed in a sugar–water solution that is governed by Navier–Stokes equation (2.1), with a jump in normal stress (2.3) from the inextensibility and bending resistance of the vesicle. The vesicle membrane's permeability to the water is driven by osmotic stress from jumps in sugar concentration and from hydrostatic stress from jumps in pressure (2.5); however, the hydrostatic stress is negligible here. The sugar is governed by the advection–diffusion equation (2.6), with its glass-forming properties included through a concentration-dependent viscosity (2.7) and diffusivity (2.8). The drying of the domain occurs through a computational boundary condition (2.11) that approximates the effect of a moving liquid–air interface from evaporation. The system is non-dimensionalized and solved by expanding on the method of Salac & Miksis [11]. This includes a level set method to track the interface, a projection step to enforce the divergence-free conditions, and an immersed interface method to enforce the impermeability of the membrane to sugar. Results are focused on the effects of the concentration-dependent diffusivity, concentration-dependent viscosity, drying rate of the domain and osmotic flow through the membrane. Note that, while a vesicle lacks the cytoskeleton of red blood cells, insight into the desiccation dynamics of red blood cells can be drawn from the qualitative vesicle behaviour observed here. Future work will explore the effect of the cytoskeleton and its elastic contribution to the membrane during drying.

Various desiccation techniques in the lyopreservation field have indicated that incomplete desiccation hinders the successful anhydrobiotic hibernation of cells. In particular, Chakraborty *et al*. [4] hypothesize that trapped water leads to non-viable cells after rehydration in their experiments. Thus, the first series of results here focuses on parameter regimes that minimize the possibility that a vesicle has a vitrified exterior with a liquid interior. It is found that both the glassy diffusivity parameter *A*_{D} and osmotic number play significant roles in determining whether a vesicle vitrifies with, or significantly after, the surrounding solution. If the vesicle membrane has a weak reaction to the sugar, resulting in a low , the flow of water through the membrane lags behind the removal of water from the domain. Thus, the exterior of the vesicle vitrifies while the small flow through the membrane leaves significant water inside. Additionally, if the diffusion of the sugar is heavily dependent on the sugar concentration, the desiccation of the domain is significantly slowed. This gives the flow through the membrane more time to remove water from the vesicle before the exterior is vitrified. Thus, optimal values of and *A*_{D} are suggested to ensure the interior of the vesicle vitrifies along with its exterior.

Trehalose has been identified by Aksan & Toner [30], Aksan *et al*. [3] and Chakraborty *et al*. [4] as a key to the survival of cells during desiccation. Although various hypotheses exist as to the reason, including the protection of the cell membrane against buckling, this has not directly been observed at the cellular level. Thus, the second series of results focuses on the deformation of the vesicle membrane during desiccation. In general, it is found that a vesicle will deform as it loses area to drying, and then relax through 8-star, 6-star, 4-star and 3-star shapes until it reaches the equilibrium bi-concave shape. Increasing the values of the osmotic number or drying number leads to a more pronounced 8-star shape. This is due to increased area loss at the beginning of the desiccation process with high values of or . The greater top-to-bottom asymmetry that results from higher values leads to more pronounced 4-star and 3-star shapes. In particular, the top nodules of the 4-star and 3-star shapes feel the effect of increased or the most, and thus can become strongly defined, leading to areas of high curvature.

The glassy parameters *A*_{D} and *A*_{V} are found to have the most profound effects on vesicle shapes. Just as larger values of the glassy diffusivity parameter *A*_{D} provide more time for the membrane to remove water in the glass formation results, larger values of *A*_{D} also allow the vesicle to relax through the intermediate shapes earlier in the drying process. Thus, the 8-star, 4-star and 3-star shapes all have larger interior area with higher *A*_{D} values and consequently less pronounced overall shapes. The glassy viscosity parameter *A*_{V} has the opposite effect at large values. The high viscosity of the surrounding solution inhibits the relaxation of the vesicle through the intermediate shapes. Thus, the 8-star, 4-star and 3-star shapes all have less interior area with large *A*_{V} values. This results in each nodule having a larger membrane-to-area ratio and therefore being more pronounced. Additionally, for *A*_{V} = 0.0, the vesicles tend to orient the nodules off the *x-* and *y*-axes. However, with higher *A*_{V} values, the vesicles align the nodules with the axes. This leads to more pronounced nodules in the 4-star shape as *A*_{V} increases.

The results therefore support that the presence of trehalose aids in preventing membrane damage during desiccation. With a large value of *A*_{D}, the results suggest that trehalose allows a vesicle undergoing desiccation to relax through the 8-star, 4-star and 3-star shapes more quickly than it would without trehalose. This means the local curvature of the vesicle membrane will remain low, and the vesicle is therefore less likely to incur damage from bending. Similar membrane damage prevention might be expected with a cell membrane, despite the added complexities of the cell itself. The results do, however, run contrary to the general hypothesis that it is the increased viscosity from the trehalose glass that protects the membrane. Here, the increased viscosity prevents the vesicle from relaxing through the intermediate shapes. This means the membrane develops pronounced nodules with high curvature, making the vesicle or cell very susceptible to damage. Thus, it may be that the ideal lyopreservation additive is not necessarily one that causes high viscosity in the surrounding solution, but instead one whose diffusion is highly dependent on the additive concentration.

### 4.1. Suggested experiments

One key result here is that the values of and *A*_{D} have a significant effect on whether the interior of a vesicle vitrifies along with its exterior. This suggests that an ideal lyopreservation additive, sugar or otherwise, will have an appropriate combination of strong concentration-dependent diffusivity and osmotic flow rate. Although a glass-forming sugar, such as trehalose, is likely to have a strongly dependent diffusivity, the osmotic effect of the sugar on the membrane may vary between cell types. Thus, while trehalose may work well in the desiccation of one type of cell, it may work poorly on the desiccation of other cell types whose membranes react less strongly to trehalose. A study comparing the osmotic flow parameters for various sugars in anhydrobiotic cells with those in mammalian cells could investigate whether a weak osmotic response in mammalian cells is responsible for the shortfalls in lyopreservation.

Two other important observations are that vesicles undergo membrane deformation when subjected to a large solute gradient and that solutes with higher viscosity responses amplify these deformations. A study observing vesicles in large concentration gradients of various solutes could verify this. The microchannel set-up of Aksan *et al*. [3] has issues achieving a large gradient in the vicinity of the vesicle, and observing vesicle deformation with the spin-drying set-up of Chakraborty *et al*. [4] would be difficult. However, the dialysis tube set-up of Nardi *et al*. [18] is ideal for observing vesicle deformation during drying. It could be used to observe vesicles undergoing large concentration gradients of fructose (low-viscosity response) and trehalose (high-viscosity response) to verify the effect of viscosity on deformation. It should be noted that no deformation of the vesicles was reported in the Nardi *et al*. [18] experiment. However, this is likely to be due to the relatively small concentration gradient induced in that experiment. Future work with the model presented here may include simulating a vesicle in a small concentration gradient to reproduce their zero-deformation drift.

## Funding statement

C.J.V., M.J.M. and S.H.D. acknowledge partial financial support from NIH grant no. 1R01GM086886, NSF RTG grant no. DMS-0636574. M.J.M. also acknowledges partial support from NSF grant no. DMS-1312935. D.S. acknowledges funding from the State University of New York at Buffalo.

## Appendix A. Evaporation boundary condition

One approach for an evaporation boundary condition is to consider the entire microchannel of Aksan *et al*. [3] as the domain, with the vesicle at one end and a moving liquid–air interface at the other (figure 17). However, the focus of this problem is the vesicle dynamics where a 20 µm diameter vesicle occupies only one-tenth of the microchannel length (200 µm). Thus, instead of considering the full microchannel with a moving liquid–air interface, the domain is set to the 60 µm of microchannel containing the vesicle. The effect of the moving liquid–air interface is incorporated with a computational boundary condition at the fixed evaporation boundary. The solution at the evaporation boundary is considered to be far enough away from the vesicle to remain essentially static (**u** = 0). The moving liquid–air interface causes sugar to enter the domain from the left, and thus the computational boundary condition for the sugar concentration *c*(**x**, *t*) takes the form of

As water is removed, the evaporation process will slow, and thus it is expected that *F*(*c*) → 0 as *c* → *ρ*. Additionally, if the microchannel starts with no sugar, there is no flux of sugar into the domain, and thus it is expected that *F*(*c*) → 0 as *c* → 0. To find a suitable form for *F*(*c*), the sugar dynamics at the liquid–air interface are considered.

If the liquid–air interface is assumed to be planar and moving with a velocity −*h′*(*t*)**n**, where **n** points towards the air, the impermeability of the interface to sugar yields the following condition:
A1The speed *h′*(*t*) is determined from the conservation of water mass during the evaporation process. The mass flux of water (MW) across the liquid–air interface is MW = *ρ*_{water}*h′*(*t*), which is obtained by taking the time derivative of the total mass of water evaporated (*ρ*_{water}*h*(*t*)). As described by Asano [31], the mass flux of water is also proportional to the jump in water vapour pressure *v* across the interface: MW = *k _{p}*(

*v*

^{liquid}−

*v*

^{air}) where

*k*is a mass transfer coefficient. The air in the microchannel is assumed to be dry enough to consider

_{p}*v*

^{air}= 0. Raoult's law, described by Asano [31], states that the vapour pressure of a water–sugar mix follows

*v*

^{liquid}=

*x*

_{water}

*v*

_{0}, where

*x*

_{water}is the molar fraction of water and

*v*

_{0}is the vapour pressure of pure water. Denoting

*m*

_{s}and

*M*

_{s}as the molar concentration and mass of sugar and

*m*

_{w}and

*M*

_{w}as the molar concentration and mass of water, the mass flux of water is written in terms of the mass concentration of sugar (

*c*), Using this expression for MW—that MW also equals

*ρ*

_{water}

*h′*(

*t*) from the conservation of water mass, and the sugar impermeability condition (A 1)—the following is the boundary condition at the liquid–air interface: A2Using the right-hand side of (A 2), the form

*F*(

*c*) =

*λ*[

*ρ*−

*c*]/[

*ρ*– (1 −

*M*

_{w}/

*M*

_{s})

*c*]

*c*has the desired behaviour for both

*c*→

*ρ*and

*c*→ 0. Thus the computational boundary condition is (A 2), but with the material quantity

*k*

_{p}v_{0}/

*ρ*

_{water}replaced by a general model parameter

*λ*to account for the fact that (2.11) is enforced at the computational boundary instead of at the liquid–air interface. While future work may involve incorporating the liquid–air interface directly, this computational boundary condition is sufficient to observe the qualitative vesicle dynamics during desiccation.

## Appendix B. Non-dimensionalization

Denote the following non-dimensional variables: **x*** = **x/***R*_{0}, *t** = *t*/[*R*_{0}/*U*], **u*** = **u**/*U*, *μ**(**x**, *t*) = *μ*(**x**, *t*)/*μ*_{0}, , *p*^{*} = *p/*[*μ*_{0}*U*/*R*_{0}], *γ** = *γ*/[*μ*_{0}*U*], *κ** = *κ*/[1/*R*_{0}], and *σ** = *σ*/[*μ*_{0}*U*/*R*_{0}]. These variables are substituted into the solution equations (2.1) and (2.3) to give
B1
B2
B3where and . In order to non-dimensionalize the sugar equations (2.6) and (2.11), the following quantities are also defined: , and *Q** = *Q*/*U*. This gives
B4
B5
B6where , and .

## Appendix C. Convergence analysis

The following is an excerpt of the convergence analysis conducted by Vogl [24]. Convergence of the numerical scheme is verified with the following parameter values: , , *K*_{o} = 0.05, , *A*_{D} = 0.0, *A*_{V} = 0.0 and . The computational domain is [−2.5, 2.5] × [−2.5, 2.5] with an initial shape of , where . The perturbation off a circle is required because a circle is rotationally symmetric, and thus the orientation of the initial buckling is not well defined. Without a perturbation off a circle, the numerical error and grid orientation create the perturbation and define the orientation. While this is acceptable for simulation runs where the grid resolution remains constant, the convergence analysis requires this perturbation to remain constant across multiple grid resolutions. Thus, the perturbation is chosen so that the shape remains symmetric across the *x-* and *y*-axes and so that the highest mode number is large enough to avoid influencing the final result.

The results are then compared across the spatial and temporal resolutions (Δ*x*, Δ*t*) = (5/2^{6}, 2 × 10^{−3}), (5/2^{7}, 1 × 10^{−3}) and (5/2^{8}, 5 × 10^{−4}) . Because the problem is nonlinear and highly coupled, a numerical solution with high resolution (Δ*x*, Δ*t*) = (5/2^{9}, 1 × 10^{−4}) is used in place of an exact analytic solution. Table 3 shows both the *L*_{2} and *L _{∞}* errors for the sugar concentration

*c*(

*x*,

*t*) at

*t*= 0.5. Note that both error norms are indicative of second-order convergence, with some fluctuations due to the adaptive meshing.

A similar process is used to analyse the convergence of the level set *ϕ*(**x**, *t*). Because the membrane location is the quantity of interest (*ϕ* = 0), the error is computed on a neighbourhood of the membrane (|*ϕ*| < 0.5). The *L*_{2} and *L*_{∞} errors are again compared at *t* = 0.5 in table 4. Unlike the sugar concentration errors, the level set errors show convergence closer to first order. This is despite all the numerical stencils involved being of the appropriate accuracy to produce a second-order scheme overall. However, the numerical delta function in the immersed boundary method, which has a height of 1/Δ*x*, degrades the overall order of the scheme. Thus, the scheme shows first-order convergence, again with some fluctuations from the adapting meshing.

- Received June 17, 2014.
- Accepted July 28, 2014.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.