Multi-scale heat and mass transfer modelling of cell and tissue cryopreservation

Feng Xu, Sangjun Moon, Xiaohui Zhang, Lei Shao, Young Seok Song, Utkan Demirci


Cells and tissues undergo complex physical processes during cryopreservation. Understanding the underlying physical phenomena is critical to improve current cryopreservation methods and to develop new techniques. Here, we describe multi-scale approaches for modelling cell and tissue cryopreservation including heat transfer at macroscale level, crystallization, cell volume change and mass transport across cell membranes at microscale level. These multi-scale approaches allow us to study cell and tissue cryopreservation.

1. Introduction

Cryopreservation aims to preserve cells/tissues without significantly impacting their function (e.g. viability, mechanical properties; Whittingham et al. 1972; Ludwig et al. 1999; Agca 2000; Stachecki & Cohen 2004). Biological activity is first slowed down or even stopped through cooling down to subzero temperatures. This activity is retrieved after warming back to physiological temperature. Cryopreservation has been used for many important applications such as in vitro fertilization (i.e. oocyte (Porcu et al. 1997; Isachenko et al. 2005; Antinori et al. 2007) and sperm (Guthrie & Welch 2005; van den Berg et al. 2007) preservation); stem cell research (Bakken 2006; Demirci & Montesano 2007a; Hunt & Timmons 2007; Kashuba Benson et al. 2008); preservation of organs for transplantation surgery (Ishine et al. 2000); and storage and transportation of tissue engineered products (Nerem 2000). Recent advances in nano- and micro-technologies have allowed new applications in the field of bioengineering (Khademhosseini et al. 2005; Zhu et al. 2005; Demirci 2006; Burg et al. 2007; Cheng et al. 2007a,b; Demirci & Montesano 2007b; Ling et al. 2007), where novel cell cryopreservation techniques have also emerged over the past decade (Lane et al. 1999; Hyttel et al. 2000; Arav et al. 2002; Cho et al. 2003) such as quartz capillary channels (Arav et al. 2002; Risco et al. 2007; He et al. 2008) and cell encapsulation using microdroplets (Demirci & Montesano 2007b; Murua et al. 2009). Figure 1 shows the channel-based and droplet-based methods for cryopreservation. The common feature among these methods is the decreased sample size (e.g. channel volume and droplet diameter; see figure 1b) to achieve higher cooling rates. In vivo tests, e.g. implantation of myoblasts cryopreserved in microcapsules into mice (Murua et al. 2009), have shown that there is little difference between cells cryopreserved using droplet-based methods and non-cryopreserved controls (figure 2).

Figure 1.

Channel-based and droplet-based methods for cryopreservation. (a) In channel-based methods, cells are mixed with media, e.g. CPAs, inside a channel. The whole channel is then immersed in liquid nitrogen for freezing. In droplet-based methods, cell-laden droplets are ejected into nitrogen for freezing (Demirci & Montesano 2007a). (b) A comparison of the devices used in channel-based methods (Arav et al. 2002; He et al. 2008): the traditional straw (TS), the open-pulled straw (OPS), the superfine open-pulled straw (SOPS) and the quartz micro-capillary (QMC). Adapted from He et al. (2008).

Figure 2.

Droplet-based cryopreservation is a promising method as shown from in vivo results. (a) Post-thawed morphology of microencapsulated myoblasts in alginate gel droplet with 10% DMSO. Histological analysis (hematoxylin and eosin staining) 180 days post-implantation shows (c) the myoblasts cryopreserved in microcapsules explanted from the subcutaneous tissue of individuals compared with (b) non-cryopreserved control. Freezing protocol used: 1 h at −20 °C; 23 h at −80 °C; liquid nitrogen (−196 °C) for 44 days. Adapted from Murua et al. (2009). Scale bars: (a) 100 μm, (b,c) 200 μm.

The cellular response to thermal loading at cryogenic temperatures is a complex thermophysical process, i.e. heat transfer process coupled with phase change, moving phase interface (Stefan problem), mass transport (e.g. water) owing to osmotic pressure difference and volume change on freezing (Zhang et al. 2006). These thermophysical events (e.g. ice formation, vitrification) are generally manipulated with chemical adjuvants (e.g. cryoprotective agents (CPAs), antifreeze protein) to improve the cryopreservation outcome. For example, CPAs (e.g. trehalose, dimethyl sulphoxide (DMSO), glycerol, propanediol) have been found to improve the survivability of cryopreserved mammalian cells (Eroglu et al. 2000) and reduce the impact of the freezing process on tissue function (Neidert et al. 2004).

Cryopreservation involves four steps: CPA loading, freezing, thawing and CPA unloading. Slow freezing (approx. −1 °C min−1) and vitrification (approx. −100 °C min−1, transition from liquid state to glass state without forming crystals) are two currently widely used methods (Karlsson & Toner 2000). Both methods aim to minimize or eliminate cell damage during cryopreservation by minimizing the intracellular ice crystal formation. Slow freezing usually takes advantage of low concentrations of CPAs (1–2 M), offering low chemical toxicity and osmotic shock to cells (Parkening et al. 1976; Karlsson & Toner 1996). Conventional slow freezing methods function in part by the extracellular ice formation, which gradually increases the solute concentration and dehydrates cells. In contrast, vitrification uses high CPA concentrations (4–8 M) coupled with rapid cooling rates (Luyet & Hodapp 1938; Crowe et al. 1998; Arav et al. 2002; Demirci & Montesano 2007a; Gook & Edgar 2007). Vitrification using different tools, including open-pulled straws, electron microscopy grids and cryoloops, all rely on the use of short-term exposure to high levels of CPAs with resulting rapid dehydration of cells prior to freezing (Parkening et al. 1976; Karlsson & Toner 1996; Martino et al. 1996; Lane et al. 1999; Lane & Gardner 2001). The concern with the vitrification approach is the toxic effects and osmotic shock associated with high CPA levels (Arav et al. 2002). The thawing procedure could cause similar problems, such as recrystallization and osmotic shock. Most protocols adopt rapid thawing to prevent intracellular ice formation (IIF) by limiting time for crystallization (Crowe et al. 1998; Demirci & Montesano 2007a).

Cell damage during cryopreservation is generally correlated to biophysical changes including cellular dehydration as well as intracellular and extracellular ice crystal formation (Mazur 1984; He & Bischof 2003). The dehydration events at low cooling rates (less than −1 °C min−1) owing to extracellular ice formation (solution effect) induce elevated intracellular and extracellular water concentration differences that cause both lipid (e.g. thermotropic phase transformations; Caffrey 1987) and protein (e.g. cold denaturation; Privalov 1990) changes at the molecular level. These changes could cause variation of lipid organization and fluidity inside the cells (Crowe et al. 1989), and thus the damage to mammalian cell membranes (Bischof 2006). With increased cooling rate (approx. −1 °C min−1), the transport of water across cell membranes decreases to much lower than even the intracellular cytoplasm supercooling condition. Supercooling, ΔT=(TphsT) with Tphs being the equilibrium phase change temperature, is the driving force of the IIF initiation, i.e. ice crystal nucleation. The decrease in water transport across cell membranes results in IIF, which has been found to highly correlate to cell death in various cell types (e.g. 50% IIF in many cell populations yields 50% survival; Toner 1993). The cell damage is mainly caused by disruption of cell membranes and organelles owing to the volume expansion of intracellular ice crystals. When cells are dehydrated (i.e. Tphs decreases), the driving force for nucleation decreases. At extremely high cooling rate (approx. −100 °C min−1), intracellular liquid remains within the cell and IIF can be avoided. The effects of slow freezing and vitrification on cellular structures are shown in figure 3.

Figure 3.

Morphology of frozen and vitrified pancreatic substitute beads at −90 °C. The beads comprise TC3 cells encapsulated in alginate gel. Beads frozen using a conventional controlled-rate (−1 °C min−1) protocol with 1 M DMSO show considerable ice formation throughout the construct (white spaces a,c). In contrast, beads vitrified with VS55 appear to be ice free (b,d). At higher magnification (c,d), the encapsulated individual cells and clusters of cells appear to be shrunken and compressed within the frozen matrix (c) compared with the more normal appearance of the cells encapsulated in the vitrified matrix (d). Adapted from Song et al. (2005).

Studies have been carried out using numerical modelling to explore the underlying mechanisms of the physical phenomena described above (Rubinsky et al. 1980; Thom et al. 1983; Lin et al. 1990; Zabaras et al. 1991; Pegg 1996; Rabin & Steif 1996, 1998; Kandra & Devireddy 2008; Trelea et al. 2009; Zhmakin 2009). It has been demonstrated that mathematical models can be used to correlate thermophysical phenomena and biological outcomes to optimize experimental methods (Bischof 2006; Song et al. 2009). In this paper, we review methods used to model cryopreservation with a focus on mathematical formulation used during numerical analysis. First, we explain the heat transfer phenomena within the cells and tissues (macroscopic perspective) during cryopreservation processes. Crystallization of CPAs is illustrated and relevant moving boundary problems are covered. Second, we explain the mass transfer, cell dehydration and membrane transport from the microscopic perspective.

2. Heat transfer modelling

(a) Heat transfer during cryopreservation

The lumped model, in which temperature variation across CPAs is negligible, has often been widely used for slow freezing methods. However, this model is not applicable for fast freezing methods like vitrification owing to the non-uniform temperature distribution around CPAs (Han et al. 2008). The frequently used heat equation is given as Embedded Image2.1where ρ (kg m−3), c (J kg−1 K−1) and λ (W m−1 K−1) are the density, specific heat and thermal conductivity of tissue, respectively, t (s) is time, T (K) is the temperature, ΔT is the temperature gradient, Embedded Image (W m−3) is the metabolic heat generation and Embedded Image(W m−3) is the external heat source.

Equation (2.1) has been used to describe cryopreservation processes. For example, in channel-based cryopreservation methods (figure 1a) with no external energy sources, equation (2.1) can be written using the cylindrical coordinate system (Jiao et al. 2006) as Embedded Image2.2Here, several assumptions can be made: (i) negligible latent heat owing to low ice crystallization (Jiao et al. 2006; Han et al. 2008), (ii) ignorable heat transfer in z direction (straw axis) due to the geometry of a very thin and long channel, (iii) temperature-independent (ρ,c and λf(T)) and geometry-independent (ρ,c and λf(r,z)) physical parameters, and (iv) metabolic heat generation is negligible (Embedded Image). With these assumptions, equation (2.2) simply changes into the following form: Embedded Image2.3where α=λ/ρcp (m2 s−1) is the thermal diffusivity.

By solving equation (2.3), the resulting transient temperature distribution is given as (Zhmakin 2009) Embedded Image2.4where T0 (K) is the initial temperature and Embedded Image (K) is the final balanced temperature, Embedded Image, Embedded Image, Jn is nth kind of the Bessel function and ζn is the positive roots of the transcendental equation of ζnJ1(ζn)/J0(ζn)=2hr0/k with h (W m−2 K−1) being the individual heat transfer coefficient.

It should also be noted that the above heat transfer equation, equation (2.1), is based on the classical Fourier law for heating, which assumes that the propagating speed of any temperature disturbance or thermal wave is infinite. The classical Fourier law is given as Embedded Image2.5where q (W m−2) is the heat flux vector representing heat flow per unit time, per unit area of the isothermal surface in the direction of the decreasing temperature and r stands for the position vector.

The Fourier heat transfer equation has wide-ranging applicability. However, it still has limitations due to the Fourier assumption, which is that the temperature disturbance or thermal wave is assumed to propagate at an infinite speed through the medium. This assumption has been shown to be physically unrealizable since any equilibrium state in thermodynamic transition takes time to establish (Liu & Lu 1997). Fourier’s law has been shown to fail during the short duration of an initial transient of heat transfer (e.g. in microscale laser heating of thin metal films and in laser surgery techniques; Peshkov 1944), or when the thermal propagation speed of the thermal wave is not high enough (e.g. in biomaterials; Morse & Feshbach 1953; Cattaneo 1958; Vernotte 1958).

(i) Hyperbolic heat equation

The thermal wave phenomenon was first experimentally observed in liquid helium with a finite speed (Peshkov 1944; Brown et al. 1966), and later in short-pulse laser processing of thin-film engineering structures (Glass et al. 1985; Brorson et al. 1987; Yang 1991; Li 1992; Qiu & Tien 1992, 1993; Banerjee et al. 2005). Similar phenomena have also been experimentally observed in materials with non-homogeneous inner structure (e.g. sand; Kaminski 1990). The non-homogeneous inner structure of biological tissues suggests the existence of unusual heat conduction behaviour, such as temperature oscillation1 (Richardson et al. 1950; Roemer et al. 1985) and wave-like behaviour (Mitra et al. 1995).

There are physical viewpoints to explain the fundamental wave behaviour in heat conduction (Cattaneo 1958; Vernotte 1958). A heat conduction equation using the concept of a finite heat propagation velocity has been proposed. This equation is a linear extension of the unsteady Fourier equation, equation (2.5), with introduction of an additional parameter, τq (s), to account for the thermal wave behaviour not captured by Fourier’s theory. This equation is given as (Cattaneo 1958; Vernotte 1958) Embedded Image2.6Here, Embedded Image is defined as the thermal relaxation time, α (m2 s−1) and Ct (m s−1) are the thermal diffusivity and the speed of thermal wave in the medium (Mitra et al. 1995; Tzou 1997), respectively. The first-order Taylor expansion of equation (2.6) gives Embedded Image2.7Various physical points of view have been proposed for the thermal relaxation time τq (Tzou 1993), e.g. τq results from the phase lag between the heat flux vector and temperature gradient in a high-rate response, or τq represents a physical constant at which the intrinsic length scales of diffusion behaviour and wave propagation merge together. Most biomaterials have non-homogeneous structure composed of cells, extracellular matrix of superstructures, liquids and solid/soft tissue. This feature results in higher thermal relaxation times compared with engineering materials. For example, τq for biological tissues was measured to be in the range of 10–1000 s at cryogenic temperatures and 1–100 s at room temperature (Vedavarz et al. 1994).

Substituting the heat conduction equation (2.1) into the thermal wave equation, we can get the thermal wave model of heat transfer as Embedded Image2.8This equation is known as a hyperbolic heat equation because there appears a two-double-derivative term (called the wave term mathematically) that modifies the parabolic Fourier heat equation into a hyperbolic partial differential equation (Tang & Araki 1996).

The thermal wave model has been used to explain interesting phenomena (Tzou 1992). However, its validity has been questioned since: (i) it is not built upon the details of energy transport in the material; (ii) although the thermal wave model can capture the macroscale response in time, the wave concept does not capture the microscale response in space (Ozisik & Tzou 1994; Tzou 1997); and (iii) the thermal wave model can lead to unusual solutions (Taitel 1972; Godoy & García-Colín 1997; Körner & Bergmann 1998).

(ii) Dual-phase-lag model

To address the limitations of the thermal wave model, i.e. ignoring the effect of microstructural interactions in the fast transient process of heat transport (e.g. fast cooling during cryopreservation processes), the phase lag time constant for temperature gradient, τT (s), has been introduced to equation (2.6) (Ozisik & Tzou 1994; Tzou 1995, 1997). This gives the heat conduction equation called the dual-phase-lag (DPL) equation: Embedded Image2.9where τq and τT can be interpreted as the periods arising from ‘thermal inertia’ and ‘microstructural interaction’, respectively (Tzou 1995).

The simplest example of the DPL model is its first-order Taylor expansions for both q and T, given as (Xu et al.2008a, 2009) Embedded Image2.10Upon substituting equation (2.1) into this equation, we get Embedded Image2.11

Antaki (2005) pointed out that the DPL model combines the wave features of hyperbolic conduction with a diffusion-like feature not captured by the hyperbolic case. By fitting the experimental data of Mitra et al. (1995) on muscle tissue to the model, it was found that τq=16 s, τT=0.043 s for experiment 12 and τq=14 s, τT=0.056 s for experiment 3.3

(b) Crystallization

Crystallization occurs through ice crystal nucleation and growth in both intracellular and extracellular regions when the temperature of the liquid approaches its crystallization temperature. The crystallization is the phase transition of the first order including energy release owing to the latent heat of fusion. A number of factors such as cooling rate, homogeneity and pressure affect the phase transition from liquid to solid.

The kinematic equation of ice crystallization (χ being the ratio of the total quantity of crystallized ice on cooling to the maximum crystallizable ice) is expressed as (Baudot et al. 2000) Embedded Image2.12where κ=L/πδ2vTmrf and a1(χ)=χ2/3(1−χ). Here, Q (J mol−1) is the activation energy, Tmelting (K) is the final temperature of the freezing process, R (J mol−1 K−1) is the gas constant, L (J kg−1) is the latent heat, δ (m) is the thickness of the transition layer between liquid and solid, rf (m) is the radius of the ice region, and v (m2 s−1) is the kinematic viscosity. Note that the heat transfer equation and crystallization equation are solved in an uncoupled manner. Integration of equation (2.12) leads to the following implicit equation of χ (Baudot et al. 2000): Embedded Image2.13where Embedded Image arctan Embedded Image.

There are three main approaches that macroscopically model the solidification of liquids: uncoupled method, Stefan approach (sharp interface method) and zone model. As explained through equations (2.1)–(2.5) and (2.12), the uncoupled method is based on the assumption that the latent heat generation of liquid during cooling process is negligible. Consequently, the energy equation is decoupled from the kinetics of the liquid (Ren et al. 1994). Basically, the Stefan problem has characteristics of a moving boundary between liquid and solid with a sharp interface (Friedman 1968). In the classical Stefan problem addressing the melting of ice, the contribution of diffusion and convection to solidification can be neglected. In this case, crystallization is determined only by the temperature evolution (Friedman 1968). Assuming that the thermophysical parameters are constant, the following simple equation is used in both liquid and solid domains: Embedded Image2.14 The interface between liquid phase and solid phase is determined by the continuity of the heat flux, including the latent heat, at the interface as boundary conditions: Embedded Image2.15where L (J kg−1) is the latent heat during crystallization and melting, nΓ is the outer normal vector to the solid domain and vΓ is the normal component of the interface velocity. In Stefan’s approach, the entire domain is sharply divided into solid and liquid sub-domains. However, it was reported that the sharp interface assumption is not always reasonable since smooth change of the phase is more accurate to describe the underlying physics, i.e. smooth formation of ice crystals (Benard & Advani 1995). In §2c, we will explain how to deal with the moving boundary. As an alternative, the zone model uses a degree of crystallization χ (equation (2.12)). It is found that one can model the crystallization process as a propagating zone that can be either very large or very narrow (sharp interface), depending on the processing conditions and material parameters (Benard & Advani 1995). Furthermore, the zone models are based on the hypothesis that the total heat content can be calculated by an enthalpy function (Le Bot & Delaunay 2008).

There is another parameter named the probability of intracellular ice formation (PIF), where when PIF is larger than the threshold that the entire intracellular volume will freeze at once. A model has also been proposed for PIF, using heterogeneous nucleation theory (Toner 1993) as Embedded Image2.16where B (K s−1) is the cooling rate, A (m2) is the cellular surface area, Tseed (K) is the ice seeding temperature, ΔT is the supercooling defined as (TphsT), Ω is a heterogeneous term of the nucleation related to the ability of water molecules to cross the interface and join the ice phase and κ is the thermodynamic term of the nucleation related to the Gibbs free energy for the formation of a critical nucleus of ice (Toner 1993; Devireddy et al. 2002).

(c) Moving boundary problems for cryopreservation

Once the formulation of solidification is determined, we have to select a numerical method to calculate the moving boundary between liquid and solid phases during cryopreservation. The methods for the moving boundary problem are classified mainly into two groups (Crank 2005): front tracking method and front capturing method. First, the front tracking method that possesses the Lagrangian feature uses a set of marker points indicating the front interface. In this method, there is a discontinuity of solution across the interface, which is modelled by a differential equation with lower dimensionality. This method also requires interface updating, followed by interface smoothing. One of the promising methods using this Lagrangian approach is the arbitrary Lagrangian Eulerian method, which moves meshes in the domain. Second, the front capturing method is based on the Eulerian concept (Braess & Wriggers 2000). Since the front capturing method uses a fixed mesh, it is more suitable for a system with a complex, three-dimensional geometry. Furthermore, the front capturing method commonly encompasses three popular sub-methods: phase field (volume of fluid (VOF)), level set, and marker and cell method. For instance, in the VOF method, the phase function f(x) is given as (Zhmakin 2009) Embedded Image2.17Here, the entire domain consists of a solid sub-domain Ωs and a liquid sub-domain Ωl. The interface is determined by a function with a value in the range of 0 through 1. On the other hand, the level set method uses the smooth function ξ(x, t) instead of f(x), defined as (Sussman et al. 1994) Embedded Image2.18where Γ0 is the interface front. The phase function (or the smooth function) is governed by Embedded Image2.19where v is the interface velocity. The unit normal vector to the interface and the curvature of the interface are expressed as n0=−∇ξ/|∇ξ| and Kinterface=∇⋅n0, respectively.

(d) Case study

A case study is performed here to investigate the effect of the non-Fourier effect on tissue freezing (figure 4a). At t=0, the tissue that is initially of temperature T=37 °C (e.g. maintained in an incubator) is frozen by suddenly immersing into liquid nitrogen. The skin is modelled as a layered structure (i.e. epidermis, dermis and fat layers) mimicking the real skin. The problem is solved by using three models, i.e. Fourier model (equation (2.1)), thermal wave model (equation (2.8)) and DPL model (equation (2.11)). The relevant parameters used for heat transfer analysis can be found in Xu et al. (2008b).

Figure 4.

Modelling of tissue cryopreservation. Different models (i.e. Fourier model, thermal wave model and DPL model) are used to indicate the non-Fourier effect on heat transfer process during tissue cryopreservation. (a) Schematic of skin tissue immersed in liquid nitrogen for freezing. (b) Temperature distribution within skin at t=30 s. (c) Temperature variation with time at dermis–fat interface at r=1.6 mm. (b,c) Solid line, τq=0 s and τT=0 s; dashed line, τq=10 s and τT=0 s; and dotted line, τq=10 s and τT=10 s.

The comparison between results predicted by the three heat transfer models is shown in figure 4b for temperature spatial distribution along skin depth and in figure 4c for temperature temporal variation at the dermis–fat interface. The results demonstrate that tissue temperatures during the freezing process predicted by the different models deviate substantially. With the thermal wave model, the temperature inside the tissue was undisturbed during the initial stage of freezing before jumping instantaneously (t<15 s in figure 4c). This may be viewed as the wave front emerging from the finite propagation of the thermal wave, i.e. existence of a relaxation time τq that is the phase lag in establishing the heat flux and the associated conduction through a medium. Unlike the thermal wave model, no wave behaviour is observed in results from the DPL model, but a non-Fourier diffusion-like behaviour exists. This is due to the second thermal relaxation time τT which accounts for the diffusion of heat ahead of sharp wave fronts that would be induced by τq. This relaxation time constant is the phase lag in establishing the temperature gradient across the medium during which conduction occurs through its small-scale structures. The existence of τq weakens the thermal wave and thereby destroys the sharp wave front. It is noticed that a sudden temperature step at both boundaries is associated with the DPL model, as shown in figure 4b. In summary, these results demonstrate that non-Fourier features could play an important role in the tissue cryopreservation process.

3. Mass transport modelling of cryopreservation

(a) Mass transfer at macroscale

To load and unload CPAs, macroscopic flow along with CPAs can be used. For instance, a microfluidic device can be designed and fabricated to load cells with CPAs along a microfluidic channel, where the concentration of CPAs progressively varies along the channel by diffusion (Song et al. 2009). In this case, the macroscopic flow at steady state is governed by the Navier–Stokes equations as Embedded Image3.1and Embedded Image3.2where η (Pa s−1) is the fluid viscosity depending on the CPA concentration, u (m s−1) is the velocity vector and P (Pa) is the fluid pressure. Assuming that CPAs interact only with water molecules, the CPA transport is modelled by using the convection and diffusion equation at steady state: Embedded Image3.3where c (M) is CPA concentration and D (m2 s−1) is the diffusion coefficient of CPAs.

(b) Cell dehydration

As cells are cooled down, the formed ice rejects solute resulting in unfrozen fractions with higher solute concentrations. The reduced amount of water in the intracellular region causes an increase in the intracellular concentration of solutes (Katkov 2002), and thus an osmotic pressure difference across the cell membrane. This pressure difference exposes the unfrozen fraction setting up a driving force for exosmosis of water from the cell, or dehydration. The dehydration can induce serious damage to cells (Fathallah et al. 1995). Mazur’s model is widely used for describing cell dehydration (Mazur 1963), which is given as Embedded Image3.4where V (m3) is the cell volume, t (s) is time, Lp is the hydraulic permeability (m3 N−1 s−1) and can be found in He & Bischof (2003) for a variety of cell types, A (m2) is the cell surface area, and πi (N m−2) and πo (N m−2) are the intracellular and extracellular osmotic pressures, respectively.

To examine the water transport across the cell membrane, salt concentration (Embedded Image) and actual salt concentration (Embedded Image) can be defined as (Collado 2007) Embedded Image3.5and Embedded Image3.6where Vb (m3) denotes the osmotically inactive volume. Consideration of the fraction of the osmotically inactive cell volume ϕ=Vb/V derives the following equation: Embedded Image3.7When taking into account the presence of another semipermeable solute (e.g. CPA), the mean cellular concentration becomes Embedded Image with Ns being moles of this solute. By applying Henry’s law of absorption with coefficient k (Batycky et al. 1997) to describe the absorption of semipermeable solute to internal organelle membranes, we can get Embedded Image3.8where α=Sm/V (Sm being the specific surface area of organelle membranes) and K means the partition coefficient into organelles.

With the help of the Reynolds transport theorem also known as the Leibniz–Reynolds transport theorem (Collado 2007) and considering the conservation of solutes within the moving control volume V (t), the total moles of salt N(t) can be calculated as Embedded Image3.9where c(r,t) (M) is the solute concentration. Also, the equation can be rewritten in the derivative form as Embedded Image3.10where ∂V is the volume boundary moving with velocity u (or surface of the control volume) and v is the mass average velocity of convection across the volume boundary. In the case of a spherical cell geometry Embedded Image3.11and the rate of solute or salt concentration variation in a cell is given as Embedded Image3.12The term above in parentheses represents the flux of solutes.

(c) Mass transport across cell membrane

Models have been proposed to describe the mass transport across cell membranes by diffusion during freezing (Kleinhans 1998; Woods et al. 1999; Cui et al. 2002; Ateshian et al. 2006; Li 2006), such as Fick’s equation (Rubinsky & Pegg 1988; Rubinsky 1989; Bischof & Rubinsky 1993) and Mazur’s equation (Batycky et al. 1997; Zhmakin 2009). These models can commonly be classified into three subtypes, i.e. one-parameter, two-parameter and three-parameter models, depending on whether solute permeability, water permeability and/or the interaction between the solute and the solvent are considered (Kleinhans 1998; Woods et al. 1999; Ateshian et al. 2006). In these models, it is assumed that a permeable CPA diffuses in/out of cells and water channels interact with low-molecular-weight CPA channels. The three-parameter model, also known as the Kedem–Katchalsky model, adopts the reflection coefficient to couple the water and solute transport across the cell membrane. As a result, the model can help to understand the change in fluid fluxes, volume changes, and CPA molarity of cells due to osmotic pressure. Assuming that cells are submerged in the CPAs, the water flux across the cell membrane (Jw) is given as (Cui et al. 2002; Li 2006; Song et al. 2009) Embedded Image3.13 Embedded Image3.14 Embedded Image3.15where Lp (m3 N−1 s−1) indicates the membrane hydraulic conductivity; σs is the membrane reflection coefficient of CPAs (related to how the cell membrane can reflect solute particles from passing through); R (J mol−1 K−1) is the universal gas constant; T (K) is the temperature; c (M) is the concentration of CPAs; A (m2) is the cell surface area; E (Pa) is the cellular elastic modulus; and Vw (m3) is the water volume of cells. In addition, superscripts i and e in the above equations denote the intracellular and extracellular regions, respectively. Next, the CPA flux (Jc) is expressed as Embedded Image3.16and Embedded Image3.17where ω (mol N−1 s−1) and cup (M) are the membrane permeability of CPAs and the upstream concentration of CPAs, respectively. On the other hand, since these microscopic equations are coupled with the macroscopic transport phenomena (equations (3.1)–(3.3)), we have to take into account all of them simultaneously (Friedman 1968).

It should be noted that the current approaches as reviewed in this paper have limitations. They cover macroscale (greater than 1 mm, tissue level) and microscale (approx. 1 μm to 1 mm, cellular level), where the continuum assumption is still valid. However, when the scale goes down to nanoscale (less than 100 nm, molecular level, e.g. protein, lipid and DNA), the continuum assumption is not applicable. In this range, on which we do not focus in this paper, mathematical models follow molecular dynamics (Bromfield et al. 2009; Wang et al. 2009).

(d) Case study

During cryopreservation, cells may undergo some damage due to osmotic shock, dehydration and extracellular/intracellular ice crystallization (Mazur et al. 1992; Hyttel et al. 2000; Ogonuki et al. 2006). Therefore, it is of great importance to minimize the cell damage through the cryopreservation processes. In particular, osmotic shock before freezing and after thawing (i.e. in CPA loading and unloading processes) acts as a primary factor that reduces cell viability during cryopreservation (Song et al. 2009). In this case study, the mass transport in a microfluidic channel is investigated by carrying out macroscopic analyses with equations (3.1)–(3.3). The model, equations (3.4)–(3.13), was used to simulate the CPA loading and unloading processes in a microfluidic device (figure 5a; Song et al. 2009). Cells are infused into the middle channel with simultaneous injection of CPAs during CPA loading process or phosphate buffered saline (PBS) buffer during CPA unloading process into the two side channels. The Kedem and Katchalsky model (equations (3.13)–(3.17)) is used here to characterize osmotic shock at a microscale, i.e. water and CPA transport across the cell membrane, figure 5b,c. Cells experience the variation in the CPA and water concentration (e.g. water goes into cells and CPA comes out of cells in the CPA unloading step) in a progressive manner while flowing through the channel due to the diffusion process and laminar flow (Reynolds number Re∼3; figure 5d,e). This smooth variation of intracellular water concentration prevents cells from abrupt change in water flux and thus from high levels of osmotic shock. The model simulation results are consistent with the experimental results for both CPA loading and unloading processes (figure 5d,e).

Figure 5.

Modelling of CPA uploading and unloading processes using a microfluidic device for cell cryopreservation. (a) The microfluidic device has a three-input channel with dimensions of 100 μm in height, 100 μm in width and 1.5 m in length. Bright field image is the microfluidic channel inlet during flow. (b) Dimensionless water fluxes and CPA fluxes across cell membrane during the CPA loading step and (c) CPA unloading step. Here, V0 is the initial water volume in the cell, C0 is the final CPA concentration, t0=V0/(A0RTC0Lp) is the characteristic time, Jw0=V0/(A0t0) is the initial water flux and Jc0= C0/(A0t0) is the initial CPA flux. Normalized CPA concentration variation along the microfluidic channel in the (d) loading step and (e) unloading step agrees well with experimental data. (b,c) Solid line, water flux; dotted line, CPA flux. (d,e) Filled circle, experiment run 1; open circle, experiment run 2; inverted triangle, experiment run 3; solid line, numerical. Adapted from Song et al. (2009).

From this case study, we can understand how water and CPA fluxes are correlated with osmotic shock that cells experience during CPA loading and unloading processes. Such a modelling approach gives insight into the macroscopic transport behaviour of CPAs in a microfluidic channel. This can help to optimize the current cryopreservation protocols and to develop new protocols by minimizing osmotic shock during CPA loading and unloading.

4. Summary

We present mathematical formulation and methods regarding heat and mass transfer processes during cryopreservation. The models cover the multiple spatial and temporal scales from cell scale to tissue scale. In the heat transfer section, modelling of heat transfer in tissues and ice crystallization in cells during cryopreservation were elucidated. In the mass transfer section, physiological phenomena related to cells such as cell dehydration and cell membrane transport were outlined. The imbalance of osmotic pressure across cell membranes that serves as a driving force was formulated mathematically. Since these macroscopic analyses for cryopreservation are linked to the microscopic results, it is necessary to investigate the entire cryopreservation process from multiple perspectives, i.e. through combining microscale and macroscale analyses.

The multi-scale modelling approach as reviewed in this paper can also be applied to other fields such as cryosurgery, where malignant or unwanted cells/tissues are destroyed by localized freezing using a fine surgical probe (Rubinsky 2000). Cryosurgery is now used routinely for treatment of cancer and other diseases (Rubinsky 2000; Spencer 2004). The need to design and predict the treatment outcome necessitates mathematical models. For example, both macroscale (Chua et al. 2007; Romero-Mendez et al. 2007) and microscale (Zhang et al. 2003) models have been developed for cancer cryosurgery to achieve maximal cell destruction within the tumour while preserving neighbouring healthy tissue.

The ultimate goal of the mathematical modelling is to improve cryopreservation outcome by supplying reliable prediction of local temperature and crystallization in cells/tissues during cryopreservation. This necessitates accurate information on parameters used in the models such as thermal properties (e.g. thermal conductivity and specific heat), mass transfer properties (e.g. mass diffusivity) and geometry (e.g. cell shape and tissue microarchitecture). These parameters are often temperature dependent. More quantitative measurements of these model parameters are thus needed. Also, the connection between thermophysical parameters and the ultimate biological outcome (e.g. cell viability and tissue function) is still not well understood. Furthermore, most of the mathematical models handle the freezing process and few consider the storage process during cryopreservation, where the tolerance of cells to cryopreservation is induced (e.g. cold shock protein; Kim et al. 1998, 2009). An enhanced understanding of the mechanisms involved in this cell tolerance induction will help further development and application of cryopreservation techniques.


F.X., S.M., X.Z., L.S. and U.D. would like to acknowledge the support from NIH R21 (EB007707) and NIH R01 (AI081534). Y.S.S. would like to thank the research fund of Dankook University in 2009.


  • One contribution of 9 to a Theme Issue ‘Multi-scale biothermal and biomechanical behaviours of biological materials’.

  • 1 An unusual oscillation of tissue temperature with heating.

  • 2 The experiments were designed to show that heat waves take a finite time to reach a particular point inside the sample. In the experiment, two identical meat samples at different initial temperatures were brought into contact with each other.

  • 3 The experiments were designed to show wave superposition: one thin sample is sandwiched by two thicker ones.


View Abstract