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

doi:10.1098/rsta.2009.0248

Philos Trans A Math Phys Eng Sci. 2010 February 13; 368(1912): 561–583.

doi: 10.1098/rsta.2009.0248

PMCID: PMC3263795

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

Feng Xu,¹ Sangjun Moon,¹ Xiaohui Zhang,¹ Lei Shao,¹ Young Seok Song,²^* and Utkan Demirci^1,³^*

¹Bio-Acoustic-MEMS in Medicine (BAMM) Laboratory, Center for Bioengineering, Department of Medicine, Brigham and Women’s Hospital, Harvard Medical School, Boston, MA, USA

²Polymer System Division, Fiber System Engineering, Dankook University, Yongin-si, Gyeonggi-do, Korea

³Harvard-Massachusetts Institute of Technology Health Sciences and Technology, Cambridge, MA, USA

^*Author for correspondence (Email: ysong/at/dankook.ac.kr; Email: udemirci/at/rics.bwh.harvard.edu).

This article has been cited by other articles in PMC.

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.

Keywords: multi-scale modelling, heat and mass transfer, 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 (more ...)

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 (more ...)

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⁻¹) and vitrification (approx. −100

°C

min⁻¹, 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⁻¹) 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⁻¹), 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⁻¹), 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⁻¹ (more ...)

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

A mathematical equation, expression, or formula.
Object name is rsta20090248-e1.jpg

2.1

where ρ (kg

m⁻³), c (J

kg⁻¹

K⁻¹) and λ (W

m⁻¹

K⁻¹) 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, An external file that holds a picture, illustration, etc.
Object name is rsta20090248-i1.jpg

An external file that holds a picture, illustration, etc.
Object name is rsta20090248-i1.jpg

m⁻³) is the metabolic heat generation and An external file that holds a picture, illustration, etc.
Object name is rsta20090248-i2.jpg

m⁻³) 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

A mathematical equation, expression, or formula.
Object name is rsta20090248-e2.jpg

2.2

Here, 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 ( An external file that holds a picture, illustration, etc.
Object name is rsta20090248-i3.jpg

An external file that holds a picture, illustration, etc.
Object name is rsta20090248-i3.jpg

). With these assumptions, equation (2.2) simply changes into the following form:

A mathematical equation, expression, or formula.
Object name is rsta20090248-e3.jpg

2.3

where α=λ/ρc_p (m²

s⁻¹) is the thermal diffusivity.

By solving equation (2.3), the resulting transient temperature distribution is given as (Zhmakin 2009)

A mathematical equation, expression, or formula.
Object name is rsta20090248-e4.jpg

2.4

where T₀

(K) is the initial temperature and An external file that holds a picture, illustration, etc.
Object name is rsta20090248-i4.jpg

(K) is the final balanced temperature, An external file that holds a picture, illustration, etc.
Object name is rsta20090248-i5.jpg

An external file that holds a picture, illustration, etc.
Object name is rsta20090248-i6.jpg

, J_n is nth kind of the Bessel function and ζ_n is the positive roots of the transcendental equation of ζ_nJ₁(ζ_n)/J₀(ζ_n)=2hr₀/k with h

m⁻²

K⁻¹) 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

A mathematical equation, expression, or formula.
Object name is rsta20090248-e5.jpg

2.5

where q (W

m⁻²) 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^¹ (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)

A mathematical equation, expression, or formula.
Object name is rsta20090248-e6.jpg

2.6

Here,

An external file that holds a picture, illustration, etc.
Object name is rsta20090248-i7.jpg

is defined as the thermal relaxation time, α (m²

s⁻¹) and C_t (m

s⁻¹) 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

A mathematical equation, expression, or formula.
Object name is rsta20090248-e7.jpg

2.7

Various 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

A mathematical equation, expression, or formula.
Object name is rsta20090248-e8.jpg

2.8

This 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:

A mathematical equation, expression, or formula.
Object name is rsta20090248-e9.jpg

2.9

where τ_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)

A mathematical equation, expression, or formula.
Object name is rsta20090248-e10.jpg

2.10

Upon substituting equation (2.1) into this equation, we get

A mathematical equation, expression, or formula.
Object name is rsta20090248-e11.jpg

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^² and τ_q=14

s, τ_T=0.056

s for experiment 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)

A mathematical equation, expression, or formula.
Object name is rsta20090248-e12.jpg

2.12

where κ=L/πδ²vT_mr_f and a₁(χ)=χ^2/3(1−χ). Here, Q (J

mol⁻¹) is the activation energy, T_melting

(K) is the final temperature of the freezing process, R (J

mol⁻¹

K⁻¹) is the gas constant, L (J

kg⁻¹) 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²

s⁻¹) 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):

A mathematical equation, expression, or formula.
Object name is rsta20090248-e13.jpg

2.13

where

An external file that holds a picture, illustration, etc.
Object name is rsta20090248-i8.jpg

arctan

An external file that holds a picture, illustration, etc.
Object name is rsta20090248-i9.jpg

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:

A mathematical equation, expression, or formula.
Object name is rsta20090248-e14.jpg

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:

A mathematical equation, expression, or formula.
Object name is rsta20090248-e15.jpg

2.15

where L (J

kg⁻¹) 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

A mathematical equation, expression, or formula.
Object name is rsta20090248-e16.jpg

2.16

where B (K

s⁻¹) is the cooling rate, A

(m²) 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).

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)

A mathematical equation, expression, or formula.
Object name is rsta20090248-e17.jpg

2.17

Here, 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)

A mathematical equation, expression, or formula.
Object name is rsta20090248-e18.jpg

2.18

where Γ₀ is the interface front. The phase function (or the smooth function) is governed by

A mathematical equation, expression, or formula.
Object name is rsta20090248-e19.jpg

2.19

where v is the interface velocity. The unit normal vector to the interface and the curvature of the interface are expressed as n₀=− [nabla]

ξ/|

ξ| and K_interface= [nabla]

n₀, 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 (more ...)

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

A mathematical equation, expression, or formula.
Object name is rsta20090248-e20.jpg

3.1

and

A mathematical equation, expression, or formula.
Object name is rsta20090248-e21.jpg

3.2

where η (Pa

s⁻¹) is the fluid viscosity depending on the CPA concentration, u

s⁻¹) 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:

A mathematical equation, expression, or formula.
Object name is rsta20090248-e22.jpg

3.3

where c (M) is CPA concentration and D (m²

s⁻¹) 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

A mathematical equation, expression, or formula.
Object name is rsta20090248-e23.jpg

3.4

where V (m³) is the cell volume, t (s) is time, L_p is the hydraulic permeability (m³

N⁻¹

s⁻¹) and can be found in He & Bischof (2003) for a variety of cell types, A (m²) is the cell surface area, and π_i (N

m⁻²) and π_o (N

m⁻²) are the intracellular and extracellular osmotic pressures, respectively.

To examine the water transport across the cell membrane, salt concentration ( An external file that holds a picture, illustration, etc.
Object name is rsta20090248-i10.jpg

An external file that holds a picture, illustration, etc.
Object name is rsta20090248-i10.jpg

) and actual salt concentration ( An external file that holds a picture, illustration, etc.
Object name is rsta20090248-i11.jpg

) can be defined as (Collado 2007)

A mathematical equation, expression, or formula.
Object name is rsta20090248-e24.jpg

3.5

and

A mathematical equation, expression, or formula.
Object name is rsta20090248-e25.jpg

3.6

where V_b (m³) denotes the osmotically inactive volume. Consideration of the fraction of the osmotically inactive cell volume ϕ=V_b/V derives the following equation:

A mathematical equation, expression, or formula.
Object name is rsta20090248-e26.jpg

3.7

When taking into account the presence of another semipermeable solute (e.g. CPA), the mean cellular concentration becomes An external file that holds a picture, illustration, etc.
Object name is rsta20090248-i12.jpg

An external file that holds a picture, illustration, etc.
Object name is rsta20090248-i12.jpg

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

A mathematical equation, expression, or formula.
Object name is rsta20090248-e27.jpg

3.8

where α=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

A mathematical equation, expression, or formula.
Object name is rsta20090248-e28.jpg

3.9

where c(r,t)

(M) is the solute concentration. Also, the equation can be rewritten in the derivative form as

A mathematical equation, expression, or formula.
Object name is rsta20090248-e29.jpg

3.10

where

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

A mathematical equation, expression, or formula.
Object name is rsta20090248-e30.jpg

3.11

and the rate of solute or salt concentration variation in a cell is given as

A mathematical equation, expression, or formula.
Object name is rsta20090248-e31.jpg

3.12

The term above in parentheses represents the flux of solutes.

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)

A mathematical equation, expression, or formula.
Object name is rsta20090248-e32.jpg

3.13

A mathematical equation, expression, or formula.
Object name is rsta20090248-e33.jpg

3.14

A mathematical equation, expression, or formula.
Object name is rsta20090248-e34.jpg

3.15

where L_p (m³

N⁻¹

s⁻¹) 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⁻¹

K⁻¹) is the universal gas constant; T (K) is the temperature; c (M) is the concentration of CPAs; A (m²) is the cell surface area; E (Pa) is the cellular elastic modulus; and V_w (m³) 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

A mathematical equation, expression, or formula.
Object name is rsta20090248-e35.jpg

3.16

and

A mathematical equation, expression, or formula.
Object name is rsta20090248-e36.jpg

3.17

where ω (mol

N⁻¹

s⁻¹) 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 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 (more ...)

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’.

¹An unusual oscillation of tissue temperature with heating.

²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.

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

References

Agca Y. 2000. Cryopreservation of oocyte and ovarian tissue. ILAR J 41, 207–220.220. [PubMed]
Antaki P. J. 2005. New interpretation of non-Fourier heat conduction in processed meat. J. Heat. Transf. 127, 189–193.193 (doi:10.1115/1.1844540)
Antinori M., Licata E., Dani G., Cerusico F., Versaci C., Antinori S. 2007. Cryotop vitrification of human oocytes results in high survival rate and healthy deliveries. Reprod. Biomed. Online 14, 72–79.79. [PubMed]
Arav A., Yavin S., Zeron Y., Natan D., Dekel I., Gacitua H. 2002. New trends in gamete’s cryopreservation. Mol. Cell. Endocrinol 187, 77–81.81. [PubMed]
Ateshian G. A., Likhitpanichkul M., Hung C. T. 2006. A mixture theory analysis for passive transport in osmotic loading of cells. J. Biomech 39, 464–475.475 (doi:10.1016/j.jbiomech.2004.12.013) [PMC free article] [PubMed]
Bakken A. M. 2006. Cryopreserving human peripheral blood progenitor cells. Curr. Stem Cell Res. Ther. 1, 47–54.54 (doi:10.2174/157488806775269179) [PubMed]
Banerjee A., Ogale A. A., Das C., Mitra K., Subramanian C. 2005. Temperature distribution in different materials due to short pulse laser irradiation. Heat Transf. Eng. 26, 41–49.49 (doi:10.1080/01457630591003754)
Batycky R. P., Hammerstedt R., Edwards D. A. 1997. Osmotically driven intracellular transport phenomena. Phil. Trans. R. Soc. Lond. A 355, 2459–2488.2488 (doi:10.1098/rsta.1997.0143)
Baudot A., Alger L., Boutron P. 2000. Glass-forming tendency in the system water–dimethyl sulfoxide. Cryobiology 40, 151–158.158 (doi:10.1006/cryo.2000.2234) [PubMed]
Benard A., Advani S. G. 1995. Energy equation and the crystallization kinetics of semi-crystalline polymers:regimes of coupling. Int. J. Heat Mass Transf 38, 819–832.832 (doi:10.1016/0017-9310(94)00205-A)
Bischof J. C. 2006. Micro and nanoscale phenomenon in bioheat transfer. Heat Mass Transf./Waerme- und Stoffuebertragung 42, 955–966.966.
Bischof J. C., Rubinsky B. 1993. Microscale heat and mass transfer of vascular and intracellular freezing. J. Heat Transf 115, 1029–1035.1035 (doi:10.1115/1.2911357)
Braess H., Wriggers P. 2000. Arbitrary Lagrangian Eulerian finite element analysis of free surface flow. Comput. Methods Appl. Mech. Eng 190, 95–109.109 (doi:10.1016/S0045-7825(99)00416-8)
Brazhnikov A. M., Karpychev V. A., Luikova A. V. 1975. One engineering method of calculating heat conduction processes. Inzh. Fiz. Zh 28, 677–680.680.
Bromfield J. J., Coticchio G., Hutt K., Sciajno R., Borini A., Albertini D. F. 2009. Meiotic spindle dynamics in human oocytes following slow-cooling cryopreservation. Hum. Reprod. 24, 2114–2123.2123 (doi:10.1093/humrep/dep182) [PubMed]
Brorson S. D., Fujimoto J. G., Ippen E. P. 1987. Femtosecond electronic heat transfer dynamics in thin gold film. Phys. Rev. Lett 59, 1962–1965.1965 (doi:10.1103/PhysRevLett.59.1962) [PubMed]
Brown J. B., Chung D. Y., Matthews P. W. 1966. Heat pulses at low temperatures. Phys. Lett. 21, 241–243.243 (doi:10.1016/0031-9163(66)90794-3)
Burg T. P., Godin M., Knudsen S. M., Shen W., Carlson G., Foster J. S., Babcock K., Manalis S. R. 2007. Weighing of biomolecules, single cells and single nanoparticles in fluid. Nature 446, 1066–1069.1069 (doi:10.1038/nature05741) [PubMed]
Caffrey M. 1987. The combined and separate effects of low temperature and freezing on membrane lipid mesomorphic phase behavior:relevance to cryobiology. Biochim. Biophys. Acta 896(doi:10.1016/0005-2736(87)90365-8) [PubMed]
Cattaneo C. 1958. A form of heat conduction equation which eliminates the paradox of instantaneous propagation. Comp. Rend 247, 431–433.433.
Cheng X., Liu Y. S., Irimia D., Demirci U., Yang L., Zamir L., Rodriguez W. R., Toner M., Bashir R. 2007a. Cell detection and counting through cell lysate impedance spectroscopy in microfluidic devices. Lab Chip 7, 746–755.755 (doi:10.1039/b705082h) [PubMed]
Cheng X., Irimia D., Dixon M., Sekine K., Demirci U., Zamir L., Tompkins R. G., Rodriguez W., Toner M. 2007b. A microfluidic device for practical label-free CD4(+) T cell counting of HIV-infected subjects. Lab Chip 7, 170–178.178 (doi:10.1039/b612966h) [PubMed]
Chester M. 1963. Second sound in solid. Phys. Rev 131, 2013–2015.2015 (doi:10.1103/PhysRev.131.2013)
Cho B. S., Schuster T. G., Zhu X., Chang D., Smith G. D., Takayama S. 2003. Passively driven integrated microfluidic system for separation of motile sperm. Anal. Chem. 75, 1671–1675.1675 (doi:10.1021/ac020579e) [PubMed]
Chua K. J., Chou S. K., Ho J. C. 2007. An analytical study on the thermal effects of cryosurgery on selective cell destruction. J. Biomech. 40, 100–116.116 (doi:10.1016/j.jbiomech.2005.11.005) [PubMed]
Collado F. J. 2007. Reynolds transport theorem for a two-phase flow. Appl. Phys. Lett. 90, 024101.(doi:10.1063/1.2430675)
Crank J.Free and moving boundary problems. New York, NY: Oxford University Press; 2005.
Crowe J. H., Hoekstra F. A., Crowe L. M., Anchordoguy T. J., Drobnis E. 1989. Lipid phase transitions measured in intact cells with Fourier transform infrared spectroscopy. Cryobiology 26, 76–84.84 (doi:10.1016/0011-2240(89)90035-7) [PubMed]
Crowe J. H., Carpenter J. F., Crowe J. F. 1998. The role of vitrification in anhydrobiosis. Annu. Rev. Physiol. 6, 73–103.103 (doi:10.1146/annurev.physiol.60.1.73) [PubMed]
Cui Z. F., Dykhuizen R. C., Nerem R. M., Sembanis A. 2002. Modeling of cryopreservation of engineered tissues with one-dimensional geometry. Biotechnol. Prog. 18, 354–361.361 (doi:10.1021/bp0101886) [PubMed]
Demirci U. 2006. Acoustic picoliter droplets for emerging applications in semiconductor industry and biotechnology. J. Microelectromech. Syst 15, 957–966.966 (doi:10.1109/JMEMS.2006.878879)
Demirci U., Montesano G. 2007a. Cell encapsulating droplet vitrification. Lab Chip 7, 1428–1433.1433 (doi:10.1039/b705809h) [PubMed]
Demirci U., Montesano G. 2007b. Single cell epitaxy by acoustic picolitre droplets. Lab Chip 7, 1139–1145.1145 (doi:10.1039/b704965j) [PubMed]
Devireddy R. V., Smith D. J., Bischof J. C. 2002. Effect of microscale mass transport and phase change on numerical prediction of freezing in biological tissue. J. Heat Transf 124, 365–374.374 (doi:10.1115/1.1445134)
Eroglu A., Russo M. J., Bieganski R., Fowler A., Cheley S., Bayley H., Toner M. 2000. Intracellular trehalose improves the survival of cryopreserved mammalian cells. Nat. Biotechnol 18, 163–167.167 (doi:10.1038/72608) [PubMed]
Fathallah H., Coezy E., de Neef R. S., Hardy-Dessources M. D., Giraud F. 1995. Inhibition of deoxygenation-induced membrane protein dephosphorylation and cell dehydration by phorbol esters and okadaic acid in sickle cells. Blood 86, 1999–2007.2007. [PubMed]
Friedman A. 1968. The Stefan problem in several space variables. Trans. Am. Math. Soc. 133, 51–87.87 (doi:10.2307/1994932)
Glass D. E., Ozisik M. N., Vick B. 1985. Hyperbolic heat conduction with surface radiation. Int. J. Heat Mass Transf 28, 1823–1830.1830 (doi:10.1016/0017-9310(85)90204-2)
Godoy S., García-Colín L. S. 1997. Nonvalidity of the telegrapher’s diffusion equation in two and three dimensions for crystalline solids. Phys. Rev. E 55, 2127–2131.2131 (doi:10.1103/PhysRevE.55.2127)
Gook D. A., Edgar D. H. 2007. Human oocyte cryopreservation. Hum. Reprod. Update 13, 591–605.605 (doi:10.1093/humupd/dmm028) [PubMed]
Guthrie H. D., Welch G. R. 2005. Impact of storage prior to cryopreservation on plasma membrane function and fertility of boar sperm. Theriogenology 63, 396–410.410 (doi:10.1016/j.theriogenology.2004.09.020) [PubMed]
Han X., Ma H. B., Jiao A. J., Critser J. K. 2008. Investigations on the heat transport capability of a cryogenic oscillating heat pipe and its application in achieving ultra-fast cooling rates for cell vitrification cryopreservation. Cryobiology 56, 195–203.203 (doi:10.1016/j.cryobiol.2008.02.006) [PMC free article] [PubMed]
He X., Bischof J. C. 2003. Quantification of temperature and injury response in thermal therapy and cryosurgery. Crit. Rev. Biomed. Eng. 31, 355–422.422 (doi:10.1615/CritRevBiomedEng.v31.i56.10) [PubMed]
He X., Park E. Y., Fowler A., Yarmush M. L., Toner M. 2008. Vitrification by ultra-fast cooling at a low concentration of cryoprotectants in a quartz micro-capillary:a study using murine embryonic stem cells. Cryobiology 56, 223–232.232 (doi:10.1016/j.cryobiol.2008.03.005) [PMC free article] [PubMed]
Hunt C. J., Timmons P. M. 2007. Cryopreservation of human embryonic stem cell lines. Methods Mol. Biol. 368, 261–270.270. [PubMed]
Hyttel P., Vajta G., Callesen H. 2000. Vitrification of bovine oocytes with the open pulled straw method:ultrastructural consequences. Mol. Reprod. Dev. 56, 80–88.88 (doi:10.1002/(SICI)1098-2795(200005)56:1<80::AID-MRD10>3.0.CO;2-U) [PubMed]
Isachenko V., Montag M., Isachenko E., Zaeva V., Krivokharchenko I., Shafei R., van der Ven H. 2005. Aseptic technology of vitrification of human pronuclear oocytes using open-pulled straws. Hum. Reprod. 20, 492–496.496 (doi:10.1093/humrep/deh605) [PubMed]
Ishine N. B., Rubinsky C., Lee Y. 2000. Transplantation of mammalian livers following freezing:vascular damage and functional recovery. Cryobiology 40, 84–89.89 (doi:10.1006/cryo.1999.2225) [PubMed]
Jiao A. J., Han X., Ma H. B., Critser J. K. 2006. Numerical investigations of transient heat transfer characteristics and vitrification tendencies in ultra-fast cell cooling processes. Cryobiology 52, 386–392.392 (doi:10.1016/j.cryobiol.2006.01.009) [PubMed]
Kaminski W. 1990. Hyperbolic heat conduction equation for materials with a nonhomogeneous inner structure. J. Heat Transf. 112, 555–560.560 (doi:10.1115/1.2910422)
Kandra D., Devireddy R. V. 2008. Numerical simulation of local temperature distortions during ice nucleation of cells in suspension. Int. J. Heat Mass Transf 51, 5655–5661.5661 (doi:10.1016/j.ijheatmasstransfer.2008.04.026) [PMC free article] [PubMed]
Karlsson J., Toner M. 1996. Long-term storage of tissues by cryopreservation:critical issues. Biomaterials 2, 243–256.256 (doi:10.1016/0142-9612(96)85562-1) [PubMed]
Karlsson J. O. M., Toner M., editors. (eds) 2000. Cryopreservation San Diego, CA: Academic Press.
Kashuba Benson C. M., Benson J. D., Critser J. K. 2008. An improved cryopreservation method for a mouse embryonic stem cell line. Cryobiology 56, 120–130.130 (doi:10.1016/j.cryobiol.2007.12.002) [PMC free article] [PubMed]
Katkov I. I. 2002. The point of maximum cell water volume excursion in case of presence of an impermeable solute. Cryobiology 44, 193–203.203 (doi:10.1016/S0011-2240(02)00029-9) [PubMed]
Khademhosseini A., Yeh J., Eng G., Karp J., Kaji H., Borenstein J., Farokhzad O. C., Langer R. 2005. Cell docking inside microwells within reversibly sealed microfluidic channels for fabricating multiphenotype cell arrays. Lab Chip 5, 1380–1386.1386 (doi:10.1039/b508096g) [PubMed]
Kim M. H., Sasaki K., Imai R. 2009. Cold shock domain protein 3 regulates freezing tolerance in Arabidopsis thaliana. J. Biol. Chem. 284, 23.(doi:10.1074/jbc.M109.025791) [PubMed]
Kim W. S., Khunajakr N., Dunn N. W. 1998. Effect of cold shock on protein synthesis and on cryotolerance of cells frozen for long periods in Lactococcus lactis. Cryobiology 37, 86–91.91 (doi:10.1006/cryo.1998.2104) [PubMed]
Kleinhans F. W. 1998. Membrane permeability modeling:Kedem–Katchalsky vs. a two-parameter formalism. Cryobiology 37, 271–289.289 (doi:10.1006/cryo.1998.2135) [PubMed]
Körner C., Bergmann H. W. 1998. Physical defects of the hyperbolic heat conduction equation. Appl. Phys. A:Mater. Sci. Process. 67, 397–401.401 (doi:10.1007/s003390050792)
Lane M., Gardner D. K. 2001. Vitrification of mouse oocytes using a nylon loop. Mol. Reprod. Dev. 58, 342–347.347 (doi:10.1002/1098-2795(200103)58:3<342::AID-MRD13>3.0.CO;2-X) [PubMed]
Lane M., Bavister B. D., Lyons E. A., Forest K. T. 1999. Containerless vitrification of mammalian oocytes and embryos. Nat. Biotechnol. 17, 1234–1236.1236 (doi:10.1038/70795) [PubMed]
Le Bot C., Delaunay D. 2008. Rapid solidification of indium:modeling subcooling. Mater. Charact 59, 519–527.527 (doi:10.1016/j.matchar.2007.03.010)
Li J. D.Reduction of core loss in silicon steel sheets by laser processing. Beijing, China: Tsinghua University; 1992.
Li L. Y. 2006. Numerical simulation of mass transfer during the osmotic dehydration of biological tissues. Comput. Mater. Sci. 35, 75–83.83 (doi:10.1016/j.commatsci.2005.03.006)
Lin S., Gao D. Y., Yu X. C. 1990. Thermal stress induced by water solidification in a cylinder tube. J. Heat Transf 112, 1079–1082.1082 (doi:10.1115/1.2910482)
Ling Y., Rubin J., Deng Y., Huang C., Demirci U., Karp J. M., Khademhosseini A. 2007. A cell-laden microfluidic hydrogel. Lab Chip 7, 756–762.762 (doi:10.1039/b615486g) [PubMed]
Liu J., Lu W. Q. 1997. Dual reciprocity boundary element method for solving thermal wave model of bioheat transfer. Space Med. Med. Eng. 10, 391–395.395. [PubMed]
Ludwig M., Al-Hasani S., Felderbaum R., Diedrich K. 1999. New aspects of cryopreservation of oocytes and embryos in assisted reproduction and future perspectives. Hum. Reprod. 14(Suppl.1), 162–185. [PubMed]
Luyet B. J., Hodapp A. 1938. Revival of frog’s spermatozoa vitrified in liquid air. Proc. Soc. Exp. Biol. 39, 433–434.434.
Martino A., Songsasen N., Leibo S. P. 1996. Development into blastocysts of bovine oocytes cryopreserved by ultra-rapid cooling. Biol. Reprod. 54, 1059–1069.1069 (doi:10.1095/biolreprod54.5.1059) [PubMed]
Mazur P. 1963. Kinetics of water loss from cells at subzero temperatures and the likelihood of intracellular freezing. J. Gen. Physiol 47, 347–369.369 (doi:10.1085/jgp.47.2.347) [PMC free article] [PubMed]
Mazur P. 1984. Freezing of living cells:mechanisms and implications. Am. J. Physiol. 247, 3(Pt 1), C125.–C142. [PubMed]
Mazur P., Cole K. W., Hall J. W., Schreuders P. D., Mahowald A. P. 1992. Cryobiological preservation of Drosophila embryos. Science 258, 1932–1935.1935 (doi:10.1126/science.1470915) [PubMed]
Mitra K., Kumar S., Vedavarz A., Moallemi M. K. 1995. Experimental evidence of hyperbolic heat conduction in processed meat. J. Heat Transf. 117, 568–573.573 (doi:10.1115/1.2822615)
Morse P. M., Feshbach H.Methods of theoretical physics. New York, NY: McGraw-Hill; 1953.
Murua A., Orive G., Hernandez R. M., Pedraz J. L. 2009. Cryopreservation based on freezing protocols for the long-term storage of microencapsulated myoblasts. Biomaterials 30, 3495–3501.3501 (doi:10.1016/j.biomaterials.2009.03.005) [PubMed]
Neidert M. R., Devireddy R. V., Tranquillo R. T., Bischof J. C. 2004. Cryopreservation of collagen-based tissue equivalents. II. Improved freezing in the presence of cryoprotective agents. Tissue Eng. 10, 23–32.32 (doi:10.1089/107632704322791664) [PubMed]
Nerem R. M. 2000. Tissue engineering:confronting the transplantation crisis. Proc. Inst. Mech. Eng 214, 95–99.99 (doi:10.1243/0954411001535273) [PubMed]
Ogonuki N., et al. 2006. Spermatozoa and spermatids retrieved from frozen reproductive organs or frozen whole bodies of male mice can produce normal offspring. Proc. Natl Acad. Sci. USA 103, 13 098–13 103 (doi:10.1073/pnas.0605755103) [PubMed]
Ozisik M. N., Tzou D. Y. 1994. On the wave theory in heat conduction. J. Heat Transf. 116, 526–535.535 (doi:10.1115/1.2910903)
Parkening T., Tsunoda Y., Chang M. 1976. Effects of various low temperatures, cryoprotective agents and cooling rates on the survival, fertilizability and development of frozen-thawed mouse eggs. J. Exp. Zool 197, 369–374.374 (doi:10.1002/jez.1401970310) [PubMed]
Pegg D. E. 1996. Problems in the cryopreservation of tissues and organs. Cryobiology 33, 658–659.659.
Peshkov V. 1944. Second sound in helium II. J. Phys 8, 381.
Porcu E., Fabrri R., Seracchioli R., Chiotti P., Magrini O., Flamigni C. 1997. Birth of a healthy female after intracytoplasmic sperm injection of cryopreserved human oocytes. Fertil. Steril. 68, 724–730.730 (doi:10.1016/S0015-0282(97)00268-9) [PubMed]
Privalov P. L. 1990. Cold denaturation of proteins. Crit. Rev. Biochem. Mol. Biol. 25, 281–305.305 (doi:10.3109/10409239009090612) [PubMed]
Qiu T. Q., Tien C. L. 1992. Short-pulse laser heating on metals. Int. J. Heat Mass Transf. 35, 719–726.726 (doi:10.1016/0017-9310(92)90131-B)
Qiu T. Q., Tien C. L. 1993. Heat transfer mechanisms during short-pulse laser heating of metals. J. Heat Transf. 115, 835–841.841 (doi:10.1115/1.2911377)
Quintanilla R., Racke R. 2006. A note on stability in dual-phase-lag heat conduction. Int. J. Heat Mass Transf. 49, 1209–1213.1213 (doi:10.1016/j.ijheatmasstransfer.2005.10.016)
Rabin Y., Steif P. S. 1996. Analysis of thermal stresses around cryosurgical probe. Cryobiology 33, 276–290.290 (doi:10.1006/cryo.1996.0028) [PubMed]
Rabin Y., Steif P. S. 1998. Thermal stresses in a freezing sphere and its application to cryobiology. ASME J. Appl. Mech 65, 328–333.333.
Ren H. S., Wei Y., Hua T. C., Zhang J. 1994. Theoretical prediction of vitrification and devitrification tendencies for cryoprotective solution. Cryobiology 31, 47–56.56 (doi:10.1006/cryo.1994.1006)
Richardson A. W., Imig C. G., Feucht B. L., Hines H. M. 1950. Relationship between deep tissue temperature and blood flow during electromagnetic irradiation. Arch. Phys. Med. Rehabil 31, 19–25.25. [PubMed]
Risco R., Elmoazzen H., Doughty M., He X., Toner M. 2007. Thermal performance of quartz capillaries for vitrification. Cryobiology 55, 222–229.229 (doi:10.1016/j.cryobiol.2007.08.006) [PubMed]
Roemer R. B., Oleson J. R., Cetas T. C. 1985. Oscillatory temperature response to constant power applied to canine muscle. Am. J. Physiol. 249, R153–R158. [PubMed]
Roetzel W., Putra N., Das S. K. 2003. Experiment and analysis for non-Fourier conduction in materials with non-homogeneous inner structure. Int. J. Thermal Sci. 42, 541–552.552 (doi:10.1016/S1290-0729(03)00020-6)
Romero-Mendez R., Franco W., Aguilar G. 2007. Laser-assisted cryosurgery of prostate:numerical study. Phys. Med. Biol. 52, 463–478.478 (doi:10.1088/0031-9155/52/2/011) [PubMed]
Rubinsky B. 1989. The energy equation for freezing of biological tissue. J. Heat Transf. 111, 988–996 (doi:10.1115/1.3250815)
Rubinsky B. 2000. Cryosurgery. Annu. Rev. Biomed. Eng 2, 157–187.187 (doi:10.1146/annurev.bioeng.2.1.157) [PubMed]
Rubinsky B., Pegg D. E. 1988. A mathematical model for the freezing process in biological tissue. Proc. R. Soc. Lond. B 234, 343–358.358 (doi:10.1098/rspb.1988.0053) [PubMed]
Rubinsky B., Cravalho E. G., Mikic B. 1980. Thermal stress in frozen organs. Cryobiology 17, 66–73.73 (doi:10.1016/0011-2240(80)90009-7) [PubMed]
Song Y. C., Chen Z. Z., Mukherjee N., Lightfoot F. G., Taylor M. J., Brockbank K. G., Sambanis A. 2005. Vitrification of tissue engineered pancreatic substitute. Transplant. Proc. 37, 253–255.255 (doi:10.1016/j.transproceed.2004.11.027) [PubMed]
Song Y. S., Moon S., Hulli L., Hasan S. K., Kayaalp E., Demirci U. 2009. Microfluidics for cryopreservation. Lab Chip 9, 1874–1881.1881 (doi:10.1039/b823062e) [PMC free article] [PubMed]
Spencer J. M. 2004. Cryosurgery for skin cancer. Dermatol. Surg. 30, 1269.(doi:10.1111/j.1524-4725.2004.30391.x) [PubMed]
Stachecki J. J., Cohen J. 2004. An overview of oocyte cryopreservation. Reprod. Biomed. Online 9, 152–163.163. [PubMed]
Sussman M., Smereka P., Osher S. 1994. A level set approach for computing solutions to incompressible two-phase flow. J. Comput. Phys. 114, 146–159.159 (doi:10.1006/jcph.1994.1155)
Taitel Y. 1972. On the parabolic, hyperbolic and discrete formulation of the heat conduction equation. Int. J. Heat Mass Transf. 15, 369–371.371 (doi:10.1016/0017-9310(72)90085-3)
Tang D. W., Araki N. 1996. The wave characteristics of thermal conduction in metallic films irradiated by ultra-short laser pulses. J. Phys. D Appl. Phys 29, 2527–2533.2533 (doi:10.1088/0022-3727/29/10/001)
Thom F., Matthes G., Richter E., Hackensellner H. A. 1983. The ductility of mammalian tissue in dependence on the deformation temperature. Cryo-Lett. 4, 341–348.348.
Toner M. 1993. Nucleation of ice crystals inside biological cells In Advances in low-temperature biology (ed Steponkus P., editor. ) London, UK: JAI.
Trelea I. C., Passot S., Marin M., Fonseca F. 2009. Model for heat and mass transfer in freeze-drying of pellets. J. Biomech. Eng. 131, 074501.(doi:10.1115/1.3142975) [PubMed]
Tzou D. Y. 1992. Thermal shock phenomena under high-rate response in solids. In Annual review of heat transfer (ed. Tien C. L., editor. ), pp. 111–185.185 Washington, DC:Hemisphere.
Tzou D. Y. 1993. An engineering assessment to the relaxation time in thermal wave propagation. Int. J. Heat Mass Transf. 36, 1845–1851.1851 (doi:10.1016/S0017-9310(05)80171-1)
Tzou D. Y. 1995. A unified field approach for heat conduction from macro- to micro-scales. J. Heat Transf. 117, 8–16.16 (doi:10.1115/1.2822329)
Tzou D. Y.Macro- to micro-scale heat transfer:the lagging behavior. Washington, DC: Taylor and Francis; 1997.
van den Berg H., Repping S., van der Veen F. 2007. Parental desire and acceptability of spermatogonial stem cell cryopreservation in boys with cancer. Hum. Reprod. 22, 594–597.597 (doi:10.1093/humrep/del375) [PubMed]
Vedavarz A., Kumar S., Moallemi M. K. 1994. Significance of non-Fourier heat waves in conduction. J. Heat Transf. 116, 221–224.224 (doi:10.1115/1.2910859)
Vernotte P. 1958. Les paradoxes de la theorie continue de l’equation de la chaleur. Comp. Rend 246, 3154–3155.3155.
Wang B., Tchessalov S., Cicerone M. T., Warne N. W., Pikal M. J. 2009. Impact of sucrose level on storage stability of proteins in freeze-dried solids:II. Correlation of aggregation rate with protein structure and molecular mobility. J. Pharm. Sci 98, 3145–3166.3166 (doi:10.1002/jps.21622) [PubMed]
Whittingham D. G., Leibo S. P., Mazur P. 1972. Survival of mouse embryos frozen to -196° and -269°C. Science 178, 411–414.414 (doi:10.1126/science.178.4059.411) [PubMed]
Woods E. J., Liu J., Gilmore J. A., Reid T. J., Gao D. Y., Critser J. K. 1999. Determination of human platelet membrane permeability coefficients using the Kedem–Katchalsky formalism:estimates from two- vs. three-parameter fits. Cryobiology 38, 200–208.208 (doi:10.1006/cryo.1998.2146) [PubMed]
Xu F., Seffen K. A., Lu T. J. 2008a. Non-Fourier analysis of skin biothermomechanics. Int. J. Heat Mass Transf. 51, 2237–2259.2259 (doi:10.1016/j.ijheatmasstransfer.2007.10.024)
Xu F., Wen T., Seffen K. A., Lu T. J. 2008b. Biothermomechanics of skin tissue. J. Mech. Phys. Solids 56, 1852–1884.1884 (doi:10.1016/j.jmps.2007.11.011)
Xu F., Lu T. J., Seffen K. A., Ng E. Y. K. 2009. Bioheat transfer of skin tissue. Appl. Mech. Rev. 62, 050801.(doi:10.1115/1.3124646)
Yang H. Q. 1991. Non-Fourier effect on heat conduction during welding. Int. J. Heat Mass Transf. 34, 2921–2924.2924 (doi:10.1016/0017-9310(91)90252-A)
Yarmush M. L., Toner M., Dunn J. C., Rotem A., Hubel A., Tompkins R. G. 1992. Hepatic tissue engineering. Development of critical technologies. Ann. NY Acad. Sci. 665, 238–252.252 (doi:10.1111/j.1749-6632.1992.tb42588.x) [PubMed]
Zabaras N., Ruan Y., Richmond O. 1991. On the calculation of deformations and stresses during axially symmetric solidification. ASME J. Appl. Mech. 58, 865–871.871.
Zhang A., Xu L. X., Sandison G. A., Zhang J. 2003. A microscale model for prediction of breast cancer cell damage during cryosurgery. Cryobiology 47, 143–154.154 (doi:10.1016/j.cryobiol.2003.08.002) [PubMed]
Zhang A., Xu L. X., Sandison G. A., Cheng S. 2006. Morphological study of endothelial cells during freezing. Phys. Med. Biol. 51, 6047–6060.6060 (doi:10.1088/0031-9155/51/23/007) [PubMed]
Zhmakin A. I.Fundamentals of cryobiology:physical phenomena and mathematical models. Berlin, Germany: Springer; 2009.
Zhu X., et al. 2005. Fabrication of reconfigurable protein matrices by cracking. Nat. Mater. 4, 403–406.406 (doi:10.1038/nmat1365) [PubMed]

Articles from Philosophical transactions. Series A, Mathematical, physical, and engineering sciences are provided here courtesy of
The Royal Society