Modeling Ice Crystal Formation of Water in Biological System

1. Introduction

Ice crystal formation of water plays a very important role in broadest scientific and engineering fields such as cloud evolution in atmosphere and brine pockets in marine ice. One of the most attractive attentions of this issue is about the mechanisms of cell injury induced by freezing in cryobiology [1]. During the past few decades, a series of cryogenic techniques have been developed for clinical applications in cryosurgery and cryopreservation. The former is aimed to destroy the target diseased tissues (such as tumor) while minimize pain and bleeding of the tissue as much as possible. It has been demonstrated to be a safe, minimally invasive and cost-effective treatment approach. In contrast to the cryosurgery, the task of cryopreservation is however to preserve cell, tissues or organs at a super-low temperature condition with the addition of protective agents, and ensure that they can be normally used after several re-warming stages. Clearly, no matter which process it involves, the formation of the ice crystal is ineluctable and decisive to the final fate of those cells and tissues subject to treatment.

Cryoinjury induced by ice formation is usually considered as a main killer of cells during cryosurgery or cryopreservation. A two-factor hypothesis [2], “intracellular ice formation” (IIF) at rapid cooling rates and “solution effects” at slow cooling rates, proposed by Mazur, has been serving well to explain the mechanism of freeze injury. The biological materials would subject to damage due to extracellular and intracellular ice formation during freezing. For the cryopreservation at super-cooling temperature, the occurrence of IIF is the single most important factor determining whether or not cells will survive cryopreservation [3]. So far, the ice crystal formation problems have been partially addressed in material and atmospheric sciences. However, similar researches are hard to find in cryobiology area due to complexity of the water structure in a biological system. This chapter is dedicated to summarize the typical physical and mathematical model for ice nucleation and growth in cryobiology area. In addition, the heat and solute transport, and interaction between the cell and ice during freezing are also discussed.

2. Ice nucleation

Freezing often denotes a process of new phase defined by ice creation via nucleation and subsequent crystal growth from supercooled water. This process undergoes three stages including the chemical diffusion stage for water molecule transportation, the surface kinetic stage for hydrogen bonds formation, and the final heat conduction stage for heat releasing and diffusing. It implies that the ice nucleation—whose formation determines the growth of ice crystal to a large extent—turns out to be the initial and crucial step during crystallization. In terms of whether there exist foreign bodies or not, the nucleation is usually classified as homogeneous or heterogeneous cases. By definition, homogeneous ice nucleation requests a severe condition that the liquid water is absolutely dust-free, or there exist foreign particles indeed, but whose sizes are small enough to be negligible. Since this rigid demand is hard to be satisfied, the nucleation is in fact regarded as heterogeneous in most cases and may be influenced by many external factors, such as the size of foreign particles, surface roughness and the effective size of the nucleators at different supercoolings. There exist a great many barriers to build up an integrated theory on nucleation, because of the dissimilar properties might exhibit between bulk fluid and the micro/nanoscale molecules. Up to now, the two major theoretical approaches to characterize the formation of ice can be classified as thermodynamic approximation and molecular simulations [4].

2.1. Thermodynamic model of ice nucleation

The ice nucleation does not happen necessarily if only temperature is below the freezing point. This process is determined by the free energy of system, which implies that the ice germ growth or disappearance depends on the free energy decrease or increase. For constant rate cooling protocols, the temperature of system is given by $$T\left(t\right)=T_0-Bt$$, where $$B$$ is the cooling rate. Considering a cluster of ice denoted by $$\beta$$ with $$i$$ molecules from its mother water phase $$\alpha$$, the driving force of nucleation is derived from the Gibbs free energy difference$$\Delta G_i$$ [5]

$$\Delta G_i=i{\nu }^{\beta }\Delta G_t+{(36\pi )}^{1/3}{i}^{2/3}{\sigma }^{\alpha \beta }{({\nu }^{\beta })}^{2/3}$$E1

where $$\Delta G_t$$ is negative and denotes the Gibbs free-energy difference between $$\alpha$$ and $$\beta$$ phase per volume, according to classical homogeneous nucleation theory; $${\nu }^{\beta }$$is the molecular volume of phase $$\beta$$, and $${\sigma }^{\alpha \beta }$$ is the interfacial free energy per unit area. According to thermodynamics, Gibbs free-energy difference between $$\alpha$$ and $$\beta$$ could be given as

$$\Delta G_t=\frac{\Delta H_f\Delta T}{T_f}\frac{2T}{T_f+T}$$E2

where $$\Delta H_f$$ is the heat of fusion, $$T_f$$ is the equilibrium freezing temperature, and $$\Delta T=T_f-T$$. The free energy of form of a cluster of ice is initially positive and goes through a maximum as the cluster size is increased. Maximizing with respect to $$i$$ and assuming closely packed water molecules, the following relationships for the critical cluster could be derived:

$$r_c=-\frac{2{\sigma }^{\alpha \beta }}{\Delta G_t}\text{and}\Delta G_{i_c}\text{=}\frac{\text{16}\pi }{\text{3}}\frac{{\left({\sigma }^{\alpha \beta }\right)}^3}{{\left(\Delta G_t\right)}^2}$$E3

where $$r_c$$ is the radius of the critical cluster. For heterogeneous nucleation theory $$\Delta G_{i_c}$$ is respectively modified as

$$\Delta G_{i_c\text{,}h}\text{=}f\left(\eta ,x\right)\Delta G_{i_c}$$E4

where $$x$$ is defined as $$x\text{=}{R}_s/r_p$$, $$R_s$$ is the radius of the foreign particle, and $$r_p$$ is the radius of curvature of the critical nuclei. $$\eta$$ is approximately regarded as $$\cos \theta$$ where $$\theta$$ is the contact angle between the nucleation phase and the substrate. The factor f for convex surface of particle could read as

$$f\left(\eta ,x\right)=\frac{1}{2}\left\{1+{\left(\frac{1-\eta x}{g}\right)}^3+x^3\left[2-3\left(\frac{x-\eta }{g}\right)+{\left(\frac{x-\eta }{g}\right)}^3\right]+3\eta x^2\left(\frac{x-\eta }{g}-1\right)\right\}$$E5

where $$g={\left(1+x^2-2\eta x\right)}^{1/2}$$

The homogeneous nucleation rate $$J$$ considering both diffusion barrier$$\Delta G_{i\text{'}}$$and nucleation barrier$$\Delta G_{i_c}$$could be given by

$$J=\frac{nkT}{h}\exp \left(-\frac{\Delta G_{i\text{'}}}{kT}\right)\exp \left(-\frac{\Delta G_{i_c}}{kT}\right)$$E6

where $$h$$ is Planck constant, $$k$$ the Boltzmann constant, and $$n$$ denotes molecules per unit volume of mother phase. Thus the heterogeneous nucleation rate $$J_{het}$$ could be obtained by modifying $$\Delta G_{i_c}$$ as $$\Delta G_{i_c\text{,}h}$$ in the expression of $$J$$.

2.2. Diffusion limited model of ice growth inside cell

Karlssonet al. [6] has developed a physiochemical theory of ice growth inside biological cells by coupling the water transport models, theory of ice nucleation, and theory of crystal growth. Considering a given cell with the control volume of $$V_{cv}\left(t\right)$$, which is composed of water-NaCl-glycerol solution. Evolution of the number of ice nuclei denoted by $$N\left(t\right)$$ in a given cell is a stochastic process given as

$$N\left(t\right)=\mathrm{int}\left[{\displaystyle {\int }_0^{t}J\left(t\right)V_{cv}\left(t\right)}dt\right]$$E7

where int denote the nearest integer truncated from the ensemble average. When the crystallizing phase is of different composition than the bulk solution, such as the case with ice crystallization from aqueous solutions, the rate of crystal growth is believed to be limited by the diffusion water molecules to the ice-solution interface. Considering a spherically symmetric ice crystal under isothermal conditions, its radius at time t could be given as [6]

$$r\left(t,t_i\right)={\left[{\displaystyle {\int }_{T\left(t_i\right)}^{T\left(t\right)}{\phi }^2\overline{D}\left(-\frac{dT}{B}\right)}dt\right]}^{1/2}$$E8

where $$\overline{D}$$ is the effective diffusivity of water and is calculated according to the viscosity of the intracellular solution. $$\phi$$ is the nondimensional crystal growth parameter and is derived from the equation obtained by curve fitting in advance based on the nondimensional supersaturation [6]. If the crystals within a given cell are small enough such that neighboring crystals and their depletion regions do not overlap, the total crystallized volume could be given as

$$V_{ice}\left(t\right)={\displaystyle \sum _{i\text{=}1}^{N\left(t\right)}\frac{4\pi }{3}{r}^3\left(t,t_i\right)}$$E9

Thus the crystallized volume fraction considering the impingement appearance could be given as

$$X_{ice}\left(t\right)=1\text{-}\exp \left(\frac{V_{ice}}{V_{cv}}\right)$$E10

According the water-transport model developed by Mazur [7] to cells containing a ternary water-NaCl-glycerol solution, one could derive the time evolution of cell volume

$$\frac{d{V}_{cv}}{dt}=-\frac{L_{p}A{R}_{g}T}{{\nu }_w}\left[\frac{\Delta H_f}{R_g}\left(\frac{1}{T_0}-\frac{1}{T}\right)-\ln \left(\frac{V_{cv}-\left(n_s{v}_s\text{+}{n}_g{v}_g\right)}{V_{cv}-\left(n_s{v}_s\text{+}{n}_g{v}_g\right)\text{+}{v}_w\left(n_s{\nu }_s\text{+}{n}_g\right)}\right)\right]$$E11

where $$L_p$$ is the water permeability; $$A$$ is the cell membrane surface area; $$R_g$$ gas constant; $$v_w$$, $$v_s$$ and $$v_g$$ denote specific volumes of water, NaCl and the glycerol, respectively; $$n_s$$ and $$n_g$$ are the number of moles of NaCl and glycerol in cell. Coupling the Eq. (10) and Eq. (11), one could calculate the crystallized volume fraction during freezing. Recently more sophisticated models [8-10] of intracellular ice formation and its growth have been developed to obtain more close to the real cell environments.

2.3. Molecular simulation of ice nucleation

Molecular dynamics (MD) is a powerful computer simulation technique to track time evolution of a system of interacting particles, such as atoms, molecules and coarse-grained particles. It has been widely used in the study of the structure, dynamics, and thermodynamics of water phase change and their complexes addressing a variety of issues including ice nucleation. The basic idea of the MD simulations is to generate the atomic trajectories of a system of finite particles by numerical integration, of a set of Newton’s equations of motion for all particles in the system with certain given boundary conditions, an initial set of positions and velocities. In addition, the potential energy of the particles system needs to be specified in order to describe the interaction between atoms or molecules. The detailed potential energy expression and its parameter are derived from both experimental work and high-level quantum mechanical calculations.

The dynamics of the ice–water interface determines the ice nucleation and growth such that considerable MD methods have been developed. Karim and Haymet [11, 12] have used the TIP4P water model to simulate the contact of the basal plane of ice with water and investigate dynamics of the ice–water interface. Their results show that the interface was found to be stable at 240 K. In addition, there appears a gradual change in the orientation and mobility of the molecules in the interfacial region. A very extensive MD simulation by Matsumoto et al. [13] predicted the freezing of pure water at different supercooling. They found that the occurrence of the ice nucleation was impressed extensively by the number of spontaneously developing long-lived hydrogen bonds at the same location. Handel et al. [14] had performed direct calculations of the ice Ih-water interfacial free energy for the TIP4P model by extending the cleaving method to molecular systems, which agrees with the experiments results.

MD simulations also could be used to investigate the effects of salt on ice nucleation, which is more close to biological environments. Smith et al. [15] had studied the potential of mean force of single ions as a function of the distance from the ice-water interface, and the free energy of transfer of $${\text{Na}}^{\text{+}}$$ and $${\text{Cl}}^{\text{-}}$$ ions from bulk water into the ice and to the interface. Their calculations suggested that the interface may have a potential that is attractive for $${\text{Cl}}^{\text{-}}$$ and repulsive for $${\text{Na}}^{\text{+}}$$. The investigation from Vrbka and Jungwirth [16], and Carignano et al. [17] indicated that the growth rate of ice was much lower in solutions with salt than in pure water.

The mechanism of ice growth inhibition by antifreeze proteins (AFP) is still poorly understood, owing to the difficulty in elucidating the molecular-scale structure of an ice–water interface, and in experimentally observing the structure and dynamics of AFP at an ice–water interface. MD is one of most popular methods that currently have the potential to clarify the mechanism of ice growth inhibition by AFP at the molecular level [18-21]. Wathen et al. [19] developed a MD which combines molecular representation and detailed energetics calculations of molecular mechanics techniques with the less-sophisticated probabilistic approach to study systems containing millions of water molecules undergoing billions of interactions. Such model could enable the 3-D shape and surface properties of the molecules to directly affect crystal formation. They have applied this technique to study the inhibitory effects of antifreeze proteins (AFPs) (from both fish and insect) on ice-crystal formation, including the replication of ice-etching patterns, ice-growth inhibition, and specific AFP-induced ice morphologies. Their work suggests that the degree of AFP activity results more from AFP ice-binding orientation than from AFP ice-binding strength. The simulation results were consistent with experimental observations very well. The effects of cryoprotectants such as ethylene glycol and glycerol solutions on ice nucleation have also been investigated [22, 23].

The electrostatic field was recently proposed to improve the output of a cryopreservation process [24-26]. As the molecular dynamics simulated [27], a DC field with magnitude of $$E=5.0\times {10}^9\text{V}/\text{m}$$ could induce crystallization of supercooled water. Recently, experiments [24-26, 28] have demonstrated that the external electric field can effectively affect the water nucleation and ice growth. A lower electrostatic field strength in the range of $$E={10}^3\sim {10}^5\text{V}/\text{m}$$, which is easily obtained in experiments, can induce increase in the nucleation temperature of water [25] and the asymmetric growth of ice [24].

3. Ice crystal growth using phase field method

To better understand the detailed events in cell damage due to extracellular and intracellular ice formation during freezing, it is rather important to model and predict the micro-scale ice crystal formation and growth behavior, and thus characterize the effects of several core factors, such as undercooling, anisotropy, thermal noise, and external electric field etc. Among various ways ever developed, the phase field method is rather flexible in characterizing the crystal formation and has been extensively adopted to model solidification process in material science [29]. Recently some authors [30, 31] begun to pay attention to possible application of this method in cryopreservation. He and Liu [32] extended deeply this theoretical strategy to characterize the ice crystal growth under electrostatic field and describe multiple ice nucleuses’ competitive growth behavior.

A dimensionless phase field model [33] is adopted to characterize the ice growth. A non-conserved phase field $$\varphi$$ ($$\varphi =0$$ is corresponding to ice and $$\varphi =1$$ to water) describes the transition from water to ice, which is driven by undercooling. A conserved dimensionless temperature field $$u$$ is used to characterize the energy variation due to temperature diffusion and phase transition. It reads as

$${\varphi }_t=m\nabla \cdot \left(B\left(\theta \right)\nabla \varphi \right)+m{\epsilon }^{-2}\varphi (1-\varphi )\left[\varphi -0.5+30\epsilon aSu\varphi (1-\varphi )\right],$$E12

$$u_t={\nabla }^{2}u-S^{-1}30{\varphi }^2{\left(1-\varphi \right)}^2{\varphi }_t+\psi \left(x,y,t\right),$$E13

where

$$B\left(\theta \right)=\left(\begin{array}{l}{\gamma }^2\left(\theta \right)-{\gamma }^{\prime }\left(\theta \right)\gamma \left(\theta \right)\\ {\gamma }^{\prime }\left(\theta \right)\gamma \left(\theta \right){\gamma }^2\left(\theta \right)\end{array}\right)$$E14

In the above equations, the term $$\psi$$ represents the thermal noise, which indicates the fluctuation of temperature, and $$u=\left(T-T_M\right)/\Delta T$$. $$\gamma \left(\theta \right)=1+\varsigma \cos n\left(\theta -{\theta }_0\right)$$ denotes dimensionless anisotropic surface tension and $$\theta =\arctan \left({\varphi }_y/{\varphi }_x\right)$$, where $${\theta }_0$$ is the angle between reference direction and $$x$$ axis. Four dimensionless parameters are defined as follow

$$m=\frac{\nu \sigma T_m}{\kappa L},\epsilon =\frac{\delta }{w},a=\frac{\sqrt{2}w}{12d},S=\frac{c_p\Delta T}{L},$$E15

where $$d=c_p\sigma T_m/\rho L^2$$ [33].

3.1. Ice crystal growth from pure water

The physical properties of pure water [30] used in the simulation are listed as: $$\kappa =1.31\times {10}^{-3}c{m}^2/s$$, $$L=333.5J/c{m}^3$$, $$T_M=273.15K$$, $$c_p=4.212J/c{m}^3\cdot K$$ and $$\sigma =7.654\times {10}^{-6}J/c{m}^2$$. The length scale $$w$$ is considered as $$2.7\times {10}^{-4}cm$$. Thus the phase field model parameters for all the simulations, according to Ref. [33], can be fixed with $$\epsilon =0.005$$, $$m=0.035$$, $$a=400$$ and $$S=0.5$$, which is corresponding to the cooling temperature $$233.56K$$.

An explicit time-differencing scheme is adopted to solve $$\varphi$$-equation, and implicit Crank-Nicolson algorithm is chosen for solving $$u$$-equation. All the simulations are performed in a two dimensional domain with $$1500\times 1500$$ grid points ($$\Delta x={10}^{-2}$$ and $$\Delta t={10}^{-4}$$) under the Neumann boundary conditions. The initial condition of the phase field for single ice nucleus is set for $$\varphi =\frac{1}{2}\left[\tanh \left(\frac{r-r_c}{2\sqrt{2}\epsilon }\right)+1\right]$$. The initialization of the dimensionless temperature is $$u=-\varphi$$ for all the simulations, which is corresponding to initial temperature $$T_M$$ in solid and $$T_M-\Delta T$$ in liquid.

Figure 1.

Micro-scale dendritic ice crystal patterns with (a) 4-fold structures without thermal noise, (b) 6-fold structures with thermal noise, (c) two seeds and (d) four seeds competitive growth. The parameters $$\gamma =0.04$$ and $$n=4$$ in (a), (d), and $$\gamma =0.02$$ and $$n=6$$ in (b), (c). The angle $${\theta }_0=0$$ for (a-d). [32]

Simulations indicate that the influence of surface tension anisotropy on the shape of ice is more evident compared with that on the growth rate. The variation of the anisotropy strength does not significantly change the velocity at the dendritic tip, which is strongly affected by the local cooling strength. Four-fold and six-fold dendritic ice crystals are often observed in experiments. Fig. 1(a) shows a four-fold dendritic ice crystal evolving from an initial circle shaped seed [32]. The thermal noise, which often induces the temperature fluctuation on the interface, is considered as the main reason of sidebranching growth [34]. A six-fold dendritic ice with sidebranching is shown in Fig. 1(b). It is necessary to consider the influence of the other crystal seeds on the ice growth. Fig. 1(c), (d) represent two seeds and four seeds competitive growth, which results in the crease of growth for both collided dendrites and preferred growth in directions of unhindered grains.

3.2. The effects of external electric on Ice crystal growth

The growth rate depends on the dynamical state of the ice-water interface, such as the surface tension and the temperature at interface. In addition, the effects of the external electric field on the ice formation have been demonstrated in both theoretical and experimental research. The influence of low electrostatic field strength on water nucleation and ice growth may be attributed to dipole polarization of the water molecules by the electrostatic field. Under the electrostatic field, the phase transition must overcome excess energy $$\Delta G^E=E\Delta \mu$$, where $$\Delta \mu ={\mu }_w-{\mu }_i$$ is the difference of electrical dipole moments between water and ice. However, for the low electrostatic field strength ($$E<{10}^5$$V/m), $$\Delta G^E$$ is much smaller than $$L$$ and can be negligible. In other words, the water molecules with dipole moments may follow the Boltzmann distribution function $$p=A\exp \left({\mu }_{w}E\cos \phi /kT\right)$$ [25], and get in the most stable state with the maximum energy along the direction of the electrostatic field ($$\phi =0$$). Although for the low electrostatic field magnitude, the value variation of $$p$$ with respect to $$\phi$$ is very small in bulk phase, the effects of electrostatic field may become largely strong at the interface and affect dynamical state of interface [24]. The polar water molecules (or clusters) may be torqued, rearranged, and forced to join the ice lattice via a special orientation and position under action of the electrostatic field. Thus the rate of attachment (or detachment) of water molecules on the ice interface is affected by external electrostatic field. In order to characterize such influence of the electrostatic field, the surface tension should be modified as $$\exp \left(\lambda \cos \left(\theta -\phi \right)\right)\sigma$$ during the ice growth, where $$\lambda$$ is proportional to $${\mu }_{w}E/kT$$.

From the above dissection, the effects of electrostatic field with low magnitude can induce change of the surface tension with the form $${\sigma }_E=\exp \left(\lambda \cos \left(\theta -\phi \right)\right)\sigma$$. Thus in phase field model, the dimensionless anisotropic surface tension should be modified as $${\gamma }_E\left(\theta \right)=\exp \left(\lambda \cos \left(\theta -\phi \right)\right)\gamma \left(\theta \right)$$. The value of $$\lambda$$ is proportional to $${\mu }_{w}E/kT$$, which characterizes the strength of electric field. In our simulations, the value of $$\lambda$$ is chosen in the range 0.02~0.2 and $$\phi =0$$, which indicates that the orientation of electric field is along positive $$x$$ axis.

Figure 2.

The influence of electric field on ice crystal growth at (a) $$\gamma =0.04$$, $$n=4$$ and $$\lambda =0.1$$, at (b) $$\gamma =0.02$$, $$n=6$$ and $$\lambda =0.1$$, at (c), (d) $$\gamma =0.02$$, $$n=6$$ and $$\lambda =0.2$$. The angle $${\theta }_0=0$$ for (a-c) and $${\theta }_0=\pi /9$$ for (d). [32]

Under the electrostatic field, the growth of dendritic ice crystals represents asymmetry behavior. Fig 2 shows that the main branches parallel to the electrostatic field grow faster than the other branches [32]. One can also find that the growth of sidebranching is strengthened along the direction of electric field and weakened along the inverse direction. The above influence becomes clearer by increasing the value of $$\lambda$$ (comparing Fig. 2(b) with (c)), which indicates that a more strong magnitude of electrostatic field would strengthen the asymmetric growth of ice. Fig. 2(d) shows how the reference direction affects the ice growth under electric field. The above results concerning the effects of electric field are in good accordance with the qualitative analysis in Ref. [25] and experimental observation in Ref. [24].

It is noteworthy that the biological media is saline solution including various ions, such as $${\text{Na}}^{\text{+}}$$ and $${\text{Cl}}^{-}$$. Under the electrostatic field, these positive and negative ions have completely different movement and gather near the corresponding electrodes, which induce non-uniform distribution of ions. In addition, the heterogeneous structure of ions, such as their different size, also leads to different ions’ behaviors. These behaviors induced by the electrostatic field will affect the dynamical behavior of ice interface and strengthen the asymmetric growth of ice [4]. Thus, the effects of the ions under the electrostatic field are also indicated via modifying the surface tension. The effects of ions coupling with the electrostatic field on the dynamical behavior of ice interface should be incorporated to a general model. The present investigation focuses on the effects of the electrostatic field on the ice interface but is not to completely characterize the effects of ions, which will be addressed in further work.

4. The interaction between ice growths with cell

The interaction between the cell and ice during freezing is very important to investigate the cell cryoinjury during freezing. However, it is still a rather difficult task to track the evolution of the solidification front and characterize the interaction between the cell membrane and ice in detail. Recently, the progress [35, 36] gave some promising research directions to solve this issue.

Mao et al. [35] developed a sharp interface method to simulate the response of a biological cell during freezing. The cell is modeled as an aqueous salt solution surrounded by a semi-permeable membrane. The concentration and temperature fields both inside and outside a single cell are computed taking into account heat transfer, mass diffusion, membrane transport, and evolution of the solidification front. The external ice front is computed for both stable and unstable growth modes. It is shown that for the particular geometry chosen in this study, the instabilities on the front and the diffusional transport have only modest effects on the cell response. For the cooling conditions, solute and cell property parameters used, the low Peclet regime applies. The computational results are therefore validated against the conventional membrane-limited transport (Mazur) model. Good agreement of the simulation results with the Mazur model are obtained for a wide range of cooling rates and membrane permeabilities. A spatially non-isothermal situation is also considered and shown to yield significant differences in the cell response in comparison to the isothermal case.

Another important process is from Chang et al. [36]. They developed a modified level set method for modeling the interaction of biological cells with ice front during a directional solidification. The simulation result has shown that different types of cell interaction with the ice interface can lead to significant differences in the final volume of individual cells. The cell is easily exposed to high concentrations of electrolytes in the interdendritic regions due to cell entrapment. Cell pushing can be used to control the dehydration of the cell more easily, which avoids exposing the cell to high concentrations that develop in the interdendritic spaces. Cell engulfment protects the cell from high concentration regions and leaves the cell vulnerable to intracellular ice formation. However, it could be used to control the dehydration of the cell and provide limited interaction with high electrolyte concentrations through combination of cell pushing and cell engulfment. The numerical results also indicate that the interactions between cells are very important for cryopreservation conditions, where suspensions have higher volume fractions of cells. Level set method could deal with the interaction of clusters of cells with ice by assuming that the clusters act like a single particle with a radius that represents the overall size of the cell cluster. Increasing the radius of the representative particle would increase the probability of engulfment and pushing by the interface. Although the intermolecular forces between cells and the fluid flow around various arrays of single cells and cell clusters have not been included in the current model, these phenomena could be easily included in the simulations in order to study their effects on the partitioning of cells by the ice interface.

5. Summary

The present chapter has presented an overview on the theoretical modeling of the microscale phase change of biological system subject to freezing with special focus on the ice crystal formation and growth. In addition, the heat and solute transport, and interaction between the cell and ice during freezing have also been discussed. Although considerable investigation focus on ice crystal formation and its growth, there still exists some outstanding questions [37-40] including ice structure and ice phase diagram. Here we discuss some important issues and challenges about ice formation in biological system in brief. Firstly, the effects of cell membrane on ice nucleation and growth inside a cell with restricted space during freezing have not been investigated deeply, which is in fact very important to clarify the puzzle about the ice propagation between cells. In addition, the improved understanding on the effects of cell–cell interactions suggested that the mechanism of ice propagation should be mediated to some extent, but not at all, by the gap junction, if uncontrolled factors are eliminated. Secondly, a detailed three dimensional model needs to be developed to characterize the complicated freezing process that coupling micro heat and mass transfer, and ice nucleation and growth, which involves the moving water-ice interface and deformed cell membrane. Thirdly, the advanced experimental measures need to be designed to determine the parameter in the physical model and confirm the model validity. All the efforts will warrant a bright future for better understanding and practice of cryobiology.

Acknowledgments

This work is supported by the NSFC under Grant No. 51006114 & 81071255, the Specialized Research Fund for the Doctoral Program of Higher Education, and Research Fund from Tsinghua University under Grant 523003001.