## Abstract

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 2007*a*; 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*. 2007*a*,*b*; Demirci & Montesano 2007*b*; 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 2007*b*; 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 1*b*) 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).

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 2007*a*; 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 2007*a*).

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*=(*T*_{phs}−*T*) with *T*_{phs} 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. *T*_{phs} 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.

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
2.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, (W m^{−3}) is the metabolic heat generation and (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 1*a*) with no external energy sources, equation (2.1) can be written using the cylindrical coordinate system (Jiao *et al*. 2006) as
2.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 (). With these assumptions, equation (2.2) simply changes into the following form:
2.3where *α*=*λ*/*ρc _{p}* (m

^{2}s

^{−1}) is the thermal diffusivity.

By solving equation (2.3), the resulting transient temperature distribution is given as (Zhmakin 2009)
2.4where *T*_{0} (K) is the initial temperature and (K) is the final balanced temperature, , , *J*_{n} is *n*th kind of the Bessel function and *ζ*_{n} is the positive roots of the transcendental equation of *ζ*_{n}*J*_{1}(*ζ*_{n})/*J*_{0}(*ζ*_{n})=2*hr*_{0}/*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
2.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 oscillation^{1} (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)
2.6Here, is defined as the thermal relaxation time, *α* (m^{2} s^{−1}) and *C*_{t} (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
2.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 2.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:
2.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.*2008*a*, 2009)
2.10Upon substituting equation (2.1) into this equation, we get
2.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 1^{2} 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)
2.12where *κ*=*L*/*πδ*^{2}*vT*_{m}*r*_{f} and *a*_{1}(*χ*)=*χ*^{2/3}(1−*χ*). Here, *Q* (J mol^{−1}) is the activation energy, *T*_{melting} (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, *r*_{f} (m) is the radius of the ice region, and *v* (m^{2} 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):
2.13where arctan .

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:
2.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:
2.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 §2*c*, 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
2.16where *B* (K s^{−1}) is the cooling rate, *A* (m^{2}) is the cellular surface area, *T*_{seed} (K) is the ice seeding temperature, Δ*T* is the supercooling defined as (*T*_{phs}−*T*), Ω 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)
2.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)
2.18where Γ_{0} is the interface front. The phase function (or the smooth function) is governed by
2.19where **v** is the interface velocity. The unit normal vector to the interface and the curvature of the interface are expressed as **n**_{0}=−∇*ξ*/|∇*ξ*| and *K*_{interface}=∇⋅**n**_{0}, respectively.

### (d) Case study

A case study is performed here to investigate the effect of the non-Fourier effect on tissue freezing (figure 4*a*). 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*. (2008*b*).

The comparison between results predicted by the three heat transfer models is shown in figure 4*b* for temperature spatial distribution along skin depth and in figure 4*c* 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 4*c*). 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 4*b*. 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
3.1and
3.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:
3.3where *c* (M) is CPA concentration and *D* (m^{2} 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
3.4where *V* (m^{3}) is the cell volume, *t* (s) is time, *L*_{p} is the hydraulic permeability (m^{3} N^{−1} s^{−1}) and can be found in He & Bischof (2003) for a variety of cell types, *A* (m^{2}) 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 () and actual salt concentration () can be defined as (Collado 2007)
3.5and
3.6where *V*_{b} (m^{3}) denotes the osmotically inactive volume. Consideration of the fraction of the osmotically inactive cell volume *ϕ*=*V*_{b}/*V* derives the following equation:
3.7When taking into account the presence of another semipermeable solute (e.g. CPA), the mean cellular concentration becomes with *N*_{s} 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
3.8where *α*=*S*_{m}/*V* (*S*_{m} 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
3.9where *c*(**r**,*t*) (M) is the solute concentration. Also, the equation can be rewritten in the derivative form as
3.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
3.11and the rate of solute or salt concentration variation in a cell is given as
3.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 (*J*_{w}) is given as (Cui *et al*. 2002; Li 2006; Song *et al*. 2009)
3.13
3.14
3.15where *L*_{p} (m^{3} 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* (m^{2}) is the cell surface area; *E* (Pa) is the cellular elastic modulus; and *V*_{w} (m^{3}) 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 (*J*_{c}) is expressed as
3.16and
3.17where *ω* (mol N^{−1} s^{−1}) and *c*_{up} (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 5*a*; 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 5*b*,*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 5*d*,*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 5*d*,*e*).

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.

## Acknowledgements

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.

## Footnotes

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.

- © 2010 The Royal Society