An Improved Model for Nucleation-Limited Ice Formation in Living Cells during Freezing

Ice formation in living cells is a lethal event during freezing and its characterization is important to the development of optimal protocols for not only cryopreservation but also cryotherapy applications. Although the model for probability of ice formation (PIF) in cells developed by Toner et al. has been widely used to predict nucleation-limited intracellular ice formation (IIF), our data of freezing Hela cells suggest that this model could give misleading prediction of PIF when the maximum PIF in cells during freezing is less than 1 (PIF ranges from 0 to 1). We introduce a new model to overcome this problem by incorporating a critical cell volume to modify the Toner's original model. We further reveal that this critical cell volume is dependent on the mechanisms of ice nucleation in cells during freezing, i.e., surface-catalyzed nucleation (SCN) and volume-catalyzed nucleation (VCN). Taken together, the improved PIF model may be valuable for better understanding of the mechanisms of ice nucleation in cells during freezing and more accurate prediction of PIF for cryopreservation and cryotherapy applications.


Introduction
Cryopreservation and cryotherapy are the two typical biomedical applications of cryogenics. On one hand, cryopreservation, the long-term banking of cells and tissues at cryogenic temperatures, is an enabling technology for the eventual success of cell-based medicine including tissue engineering, regenerative medicine, and assisted reproduction [1][2][3][4][5][6][7][8]. On the other hand, cryotherapy (also called cryosurgery or cryoablation) that utilizes freezing to destroy undesired tissues, is attracting more and more attention as a minimally invasive alternative to radical surgical intervention for treating cancer and many other diseases [8][9][10][11]. However, if the freezing protocol is not well designed, cryotherapy could end up with incomplete tumor destruction and cancer recurrence ensues [11,12]. Therefore, a good understanding of the response of tumor cells to freezing is of importance and significance for optimizing the freezing protocols to improve the treatment outcome of cryotherapy.
When freezing cell suspension, ice forms in the extracellular water first. Consequently, the concentration of extracellular solution increases, which breaks the balance in chemical potential between water inside and outside of cells. As a result, water either osmoses out of the cells through their plasma membrane (and cells dehydrate) or undergoes phase change to form ice in the cells. Freezing injury emerged in this process is generally believed to be related to these two fates of intracellular water: dehydration or solution effect [13] and intracellular ice formation (IIF) [14,15]. Among the two, IIF is considered to be the most deleterious and it damages the cells in both mechanical and non-mechanical fashions [8]. Therefore, it is important to understand IIF during freezing to either minimize it (for cryopreservation) or maximize it (for cryotherapy).
As with all phase change phenomena, ice formation in cells is a two-step process: nucleation of ice embryos and their subsequent growth. Moreover, in the absence of cryoprotective agents (CPAs), once ice embryos are nucleated in a cell, ice will instantly propagate throughout the whole cell [15,16]. This is because the viscosity is low (or diffusion coefficient is high) in physiologic solutions such as that in living cells, which allows instantaneous growth of ice embryos throughout the tiny cellular space. As a result, IIF in the absence of CPAs (e.g., for cryotherapy applications) is nucleation-limited. Therefore, ice nucleation in cells during freezing has been extensively studied in the past decades both experimentally and by modeling [15][16][17][18][19][20][21][22][23][24]. For the latter, the Toner's PIF (Probability of IIF) model [15] has been the most widely used for quantifying IIF that is nucleation-limited [14,25]. This model is based on the catalytic mechanisms including both surface-catalyzed nucleation (SCN) and volumecatalyzed nucleation (VCN) [15,22]. However, this model assumes that apparent ice should form in all cells (i.e., PIF = 1) if they are frozen to low enough temperature, which is not necessary tenable according to several recent experimental observations [14,26,27]. Therefore, it is of significance to further improve the PIF model for better prediction and understanding of ice nucleation in living cells during freezing.
To achieve this goal, we present a modified PIF model in this study by incorporating a new parameter, the critical cell volume (V f ), into the Toner's model. To evaluate the validity and accuracy of the modified PIF model, an extensive set of experiments were conducted to study the freezing response of HeLa cells in terms of both water transport across the cell plasma membrane and IIF. The mechanisms of ice nucleation and their dependence on key parameters such as temperature, cooling rate, and membrane permeability (to water) of HeLa cells were quantitatively investigated. Moreover, we reveal that the critical volume is different for the SCN and VCN mechanisms. This new model will be a valuable tool for understanding and predicting ice nucleation in living cells during freezing for both cryopreservation and cryotherapy applications.

Cell Culture and Sample Preparation
HeLa cells were maintained in culture medium (90% DMEM) supplemented with 10% FBS, 100 U L 21 penicillin G, and 100 U L 21 streptomycin. The cells were cultured at 37uC in 5% CO 2 environment. After reaching ,90% confluence, the culture plate to which the HeLa cells adhered was washed twice with 16 phosphate-buffered saline (PBS). After 5 minutes of trypsinization (in 2 ml trypsin), HeLa cells were centrifuged (100 6 g) for 5 min and re-suspended in 16 PBS for further experiments.

Cryomicroscopy Studies of Transmembrane Water Transport and IIF
As illustrated in Fig. 1, the cryomicroscopy system consists of the Olympus BX53 optical microscope mounted with a temperaturecontrollable FDCS196 cryostage (2196uC to 125uC) assembly. The control accuracy of temperature, cooling/warming rates, and holding time are 60.1uC, 61uC min 21 , and 61 s, respectively (provided by the product manual). The sample chamber was filled with nitrogen vapor during the full course of experiments to minimize condensation on the glass window for visualization. A 506long working distance objective and a Qimaging (Survey, BC, Canada) MicroPublisher 5.0 RTV CCD camera were used for real-time monitoring and recording of all experimental data.
For all cryomicroscopy experiments, a small drop (10 mL) of cell suspension was pipetted onto the center of a coverslip. Another coverslip was then gently applied on top. Since the cells were observed to flow under the coverslip, we inferred that the cells were not compressed by the coverslips presumably because the thickness of the solution film between the two coverslips was greater than the diameter of the cells. HeLa cells were cooled from room temperature to 21uC at 10uC min 21 , followed by holding at the temperature for ,30 s to seed extracellular ice using copper wires with one end winding around the stainless steel cooling tube (as depicted in Fig. 1c). To minimize the initial cell dehydration during the ice-seeding process, the cells were warmed up after iceseeding to 20.5uC at 10uC min 21 . After equilibrating at 20.5uC for 3 min, the cells were cooled to a deep subzero temperature (2 50uC or 2100uC) at various cooling rates. Finally, the cells were heated to 40uC. Although the freezing point of the PBS solution is around 20.5uC, 3 min of equilibrium is not enough for the extracellular ice to melt completely. The sizes of the extracellular crystals are comparable to those of the cells.  Five cooling rates (5,10,15,30, and 60uC min 21 ) were used to quantitatively study the freezing-induced transmembrane water transport of HeLa cells at temperatures ranging from 20.5uC to 2 45uC. For each cooling rate, spherical-shaped cells (n $28) with clear boundary were analyzed to obtain the average volumetric information at various temperatures. Image analysis was performed using gray histogram threshold segmentation with OpenCV (Open Source Computer Vision Library). The total pixel area inside a cell was tallied using the built-in function in OpenCV and then converted to square micrometers. The obtained cell area was then used to calculate the equivalent radius of the cell. Cell volume was further estimated using the equivalent cell radius.
Four cooling rates (45, 60, 75, and 100uC min 21 ) were used to study IIF in HeLa cells at temperatures ranging from 20.5uC to 2 50uC. For each cooling rate, at least 165 cells with clear boundary were analyzed to obtain the IIF information for estimating PIF in the cells. Additional PIF data (60 cells) for freezing at 60uC min 21 from 20.5uC to 2100uC were obtained to validate the effectiveness of the modified PIF models.
The experimental data of transmembrane water transport and IIF were then fitted using mathematical models to extract the model parameters for further prediction and understanding the ice formation process. The fitting was done by pooling all the experimental data of cell volume (or IIF) obtained at the different cooling rates together that are fitted simultaneously using the water transport (or IIF) model, which is therefore called ''pooled fitting'' in this study. Detailed information of the models is given below.

Modeling of Water Transport during Freezing
Transmembrane water transport or cell dehydration during freezing at various cooling rates was modeled using the following equations [14,25,28,29]: where V C is cell volume, A is surface area of the cell, R is the universal gas constant, DH f is molar heat of fusion of water, V b is osmotic inactive volume of the cell, w s is dissociation constant of salt, n s is molar amount of intracellular salt, n w is partial molar volume of water, T R is reference temperature (273.15K), L p is water permeability of the cell plasma membrane, L pg is water permeability of the cell plasma membrane at T R , and E Lp is activation energy for water transport across the cell plasma membrane. This model is based on the assumption that cell membrane is permeable only to water and DH f is independent of temperature. The model parameters including L pg , E Lp , and V b were determined by ''pooled fitting'' the model to our experimental data of cell volume obtained from the cryomicroscopy studies. The L pg and E Lp were then further utilized to predict the cell volume and surface area, which are needed to predict the probability of ice formation (PIF) in cells using the PIF model detailed below.

Modeling of PIF in Cells during Freezing
The probability of ice formation (PIF) in cells accounting for both SCN and VCN is as follows [15,25]: where B is cooling rate in uC min 21 , T 0 is equilibrium freezing temperature of isotonic solutions (20.5uC), A and V are cell surface area and volume (that can be determined from the transmembrane water transport studies), respectively, and I SCN and I VCN are nucleation rates due to SCN and VCN, respectively, which have been widely estimated using the Toner's original model as follows [15]: where the superscript XCN represents the nucleation mechanism for SCN or VCN, N is number of water molecules in contact with the plasma membrane (for SCN) or in the cells (for VCN), T f is equilibrium freezing temperature of cytoplasm, the subscript '0' refers to isotonic condition, V o XCN and k o XCN are kinetic and thermodynamic parameters of IIF, and g is viscosity of the cytoplasm. In order to give a good fitting to our experimental data of IIF, we further modified the Toner's original model for nucleation rates by introducing a critical volume (V f ) into Eq. 6 as follows: The viscosity of the cytoplasm can be estimated using the free volume model as follows [25,30]: where Q s is volume fraction of salts, g w is viscosity of water, g w,o is pre-exponential constant, v Ã w is specific volume of water at 0 K, K 11 /l and K 21 are two free volume parameters for water, and T gw is glass transition temperature of water. The equilibrium freezing temperature T f is defined as the following, where x w is water mole fraction in the cytosol.
The model parameters including V o , k o , and V f for SCN and VCN were determined by ''pooled fitting'' the model to our experimental data of IIF obtained from the cryomicroscopy studies.

Transmembrane Water Transport
Typical images of HeLa cells at room temperature, after ice seeding at 21uC, after seeding ice and warming back to 20.5uC, and after further cooling down to 219.2uC at 30uC min 21 are shown in Fig. 2a, b, c, and d, respectively. When ice was seeded using the copper wire, the differences in chemical potential of water across the cell membrane induced exosomosis of intracellular water and apparent cell dehydration was observable ( Fig. 2b  versus a). After heating back and equilibrating the cells for three minutes at 20.5uC (Fig. 2c), their size and shape restored due to melting of most of the extracellular ice [25]. When further freezing to 219.2uC, extensive dehydration of the cells is clearly visible (Fig. 2d).
The data of normalized cell volume and the corresponding ''pooled fitting'' results at five different cooling rates are shown in Fig. 3a. The initial volume that we used here is the average volume at the equilibrium temperature (20.5uC) after ice-seeding, which was done because cells move and slightly deform during seeding ice in the extracellular solution and it was difficult to identify a specific cell before and after ice-seeding. The average radius of HeLa cells at 20.5uC was measured in this study to be 7.5161.21 mm (n = 166). The parameters in the transmembrane water transport model were determined to be 0.1166 mm min 21 atm 21 and 11.9236 kcal mol 21 for L pg and E Lp , respectively (R 2 = 0.9702).
As shown in Fig. 3a, the curves of normalized cell volume under different cooling rates are very different from each other,  indicating a strong dependence of the water transport process on cooling rate. During the initial phase of the cooling process, cells undergo rapid dehydration. As the temperature continues to decrease, the kinetics of water transport slows down and asymptotically approaches a constant value (Fig. 3b). This may be attributed to the temperature dependence of the plasma membrane permeability to water [31] since L p at 240uC is only 2.4% of its value at 20.5uC. Figure 3b also depicts the starting temperature (pointed by arrows) where the water transport process has entered the equilibrium phase, showing its strong dependence on cooling rate, as well.

Intracellular Ice Formation
Two types of apparent intracellular ice formation during freezing were observed in the present study with HeLa cells (Fig. 4). The darkening IIF occurred with a clear change in gray level in the cells, resulting in sudden darkening of the cells (Fig. 4bd), indicating the instantaneous propagation of the intracellular ice once ice embryos were nucleated. In contrast, the twitching IIF only exhibited a sudden trembling movement that was observable under microscope with minimal change in gray level (Fig. 4e-h). As illustrated in Fig. 4h, bubble-like intracellular materials (red arrows) were observable to seep out of cells after the twitching movement, possibly due to IIF induced damage to the cell plasma membrane and the subsequent leakage of cytoplasm out of the cells. This phenomenon was also observed in cells undergoing darkening IIF.
To verify the validity of the modified PIF model, the parameters for transmembrane water transport across the cell membrane were utilized to predict the cell volume under freezing at four different cooling rates (Fig. 5a), which was used for simultaneous fitting to the experimental PIF data with the modified PIF model and the results are shown in Fig. 5b. The critical volume together with the kinetic and thermodynamic parameters for both SCN and VCN were determined and are given in Table 1. As shown in Fig. 5b, it appears that darkening dominated twitching at slower cooling rates while at higher cooling rate (100uC min 21 ), twitching dominated.
According to the experimentally obtained maximum PIFs at cooling rates of 45, 60, 75 and 100uC min 21 (Fig. 5c), the overall PIF and twitching PIF appears to increase along with the cooling rate. In contrast, the darkening PIF exhibits an inverted parabolic shape, where the probability reaches its maximum value of ,40% between 60 and 75uC min 21 for HeLa cells.
To show the importance of introducing the critical volume in the modified PIF model, the experimental PIF data obtained under different cooling rates (45, 60, 75 and 100uC min 21 ) were fitted using both the original Toner's and the modified PIF models. Fig. 6a illustrates a comparison of the two model predictions at temperature from 0uC to 2100uC, where the PIF data was observed to keep constant at the lower temperatures (2 50uC to 2100uC). While both models were capable of delivering accurate and logical predictions at the beginning of the cooling process with high fitting goodness (R 2 $0.9963), the fitting curves plotted by the original PIF model rapidly increase to 100% in the deep subzero region. In contrast, the modified PIF model gives more accurate prediction of PIF for HeLa cells over the entire temperature range. The data from additional experiments for freezing at 60uC/min from 20.5uC to 2100uC shown in Fig. 6b confirm that no IIF would occur below 250uC and validate the effectiveness of the modified PIF model for accurate prediction of IIF over a wide range of temperature. Note that a few of ice nucleus may still form while not grow to a detectable dimension because of the extremely lower diffusion coefficient. The concentrated residual intracellular solution will be vitrified at the lower temperatures. Since the PIF model we used is phenomenological, we omit such undetectable ice nuleation in our statistics.
The ice nucleation parameters obtained from the curve fitting to data in Fig. 6a using the original and modified PIF model are shown in Table 1. All SCN related parameters generally appears Table 1. Ice nucleation parameters estimated using both original and modified PIF model: O refers to the original PIF model, M refers to the modified version, and N/A represents not applicable. to be in good agreement for the two models. However, a large difference is noticeable in terms of the model parameters for VCN.

Transmembrane Water Transport
Since cell volume and area are required for predicting PIF (Eqs. [4][5][6][7], it is important to study cell dehydration as a result of transmembrane water transport during freezing for investigating IIF. Therefore, the dehydration responses of HeLa cells at five different cooling rates were experimentally studied and the experimental data were fitted using the transmembrane water transport model (Eqs. 1-2) to determine L pg , E Lp , and V b for further predicting the cell volume and area during freezing at any cooling rates. Although the two parameters (0.1166 mm min 21 atm 21 for L pg and 11.9236 kcal mol 21 for E Lp ) obtained for HeLa cells are not the same as that reported in a previous study [32], the differences (37.23% for L pg and 8.97% for E Lp ) are not huge either. In addition, our data are considered to be more accurate and suitable for the subsequent prediction of PIF because they were obtained by ''pooled fitting'' to experimental data obtained with five different cooling rates while in the previous study, the measurements were done under only one cooling rate.
As shown in Fig. 3a, there are one and two data points at temperatures between 23 and 27uC for the cooling rates of 30 and 60uC min 21 , respectively, for which the model could not give a good fit. This unusual increase in cell volume during the initial stage of freezing was probably an artifact due to the unclear boundary of cells for determining cell volume using imaging processing with OpenCV, which was a result of freezing induced concentration of extracellular solution that were pushed together around the cells by extracellular ice formed at these two high cooling rates (data not shown). However, no significant effect on the fitting results of model parameters (L pg , E Lp , and V b ) were observed, probably because the impact of these few data points were minimized by the many other data points as a result of the use of the ''pooled fitting'' methods to extract the two model parameters for transmembrane water transport.

Intracellular Ice Formation
Both the twitching and darkening phenomena associated with IIF were observed in the study. As previously reported, the twitching phenomenon was observed to occur as a sudden tremble with minimal change in gray level in cells whereas the darkening event was found to exhibit as either ice formation propagated from one side of the plasma membrane to the other or a sudden opacity flashing throughout the cell [15,16]. The number of twitching occurrence was found to increase along with the cooling rate (Fig. 5c). Faster cooling rate may cause the presence of a higher supercooled degree of intracellular solution, which might result in the formation of fine ice crystals that are transparent to light and therefore, not easily observable under microscope [16]. Consequently, only a brief tremble of the cell is noticeable. The darkening phenomenon of IIF, on the other hand, exhibits an inverted parabolic trend where the probability reaches a maximum value of ,40% between 60 and 75uC min 21 for HeLa cells. This phenomenon suggests that the darkening IIF may be associated with large intracellular ice crystals that can block light from going through them, resulting in a dark appearance in the area occupied by them. Twitching IIF is cell type dependent because the twitching IIF in HeLa cells was found to occur from ,220uC to 245uC which covers the homogeneous nucleation temperature (,240uC) while the twitching IIF of granulocytes [16] was found to occur 6.3-7.4uC above the homogeneous nucleation temperature. Therefore, our finding does not support the hypothesis that twitching IIF is induced only by the VCN mechanism [15,16]. Nevertheless, the observations that darkening IIF initiates from one side of the plasma membrane (Fig. 4) and the incidence of darkening cells is mainly between 210uC and 235uC (Fig. 5c) suggest darkening IIF may be mainly associated with SCN.
The maximum PIF due to both SCN and VCN in HeLa cells were found to be less than 0.5 during freezing at cooling rates below 75uC min 21 . This may contribute to the insufficient destruction of cancer cells in the tumor peripheral where the cooling rate is usually less than 50uC min 21 during cryoablation using a single cryoprobe [33]. Moreover, as shown in Fig. 6, the original model reported by Toner et al. [15] for PIF was unable to predict PIF in cells if the maximum PIF is less than 1. We overcome this problem by modifying the original model using a critical volume (V f ) below which the nucleation rate due to either SCN or VCN is zero (Eq. 7). Interestingly, the V f SCN is constantly higher than V f VCN (Fig. 5a and Table 1) and a change in slope is clearly noticeable in the experimental data of PIF (summation of both twitching and darkening) where the predicted PIF due to SCN only reaches maximum. It could be inferred from these observations that the VCN mechanism might not be a major contributor to PIF only after the SCN mechanism stops. In other words, intracellular ice during freezing is first nucleated primarily by SCN induced by extracellular ice through the cell plasma membrane when the cell volume is above V f SCN . When the cell dehydrates with a volume less than V f SCN , nucleation of intracellular water due to SCN diminishes and may be catalyzed in the whole cellular space heterogeneously by the internal organelles and/or homogeneously by the stochastic thermodynamics when the cell volume is above V f VCN . This assumption is in good agreement with the observation [15] that VCN is mostly found in the absence of external ice. Here, our interpretation may explain why both SCN and VCN mechanisms could occur in the presence of extracellular ice albeit they may occur sequentially rather than in parallel.

Summary and Conclusions
In this study, the transmembrane water transport and IIF in HeLa cells were thoroughly investigated both experimentally and by modeling. Our modified PIF model is shown to be efficient in predicting PIF over a wide range of temperature regardless of the value of the maximum PIF while the original PIF model might give a misleading prediction of PIF if the maximum PIF is less than 1. The critical volume that we introduced in the modified model is found to increase with the increase of cooling rate from 45 to 100uC min 21 . From the data of the critical volume and the change in slope of the PIF data, we infer that the IIF due to VCN mainly occurs after IIF as a result of SCN stops. The new PIF model for nucleation-limited IIF together with the concept of critical volume introduced in this study may be valuable for understanding the mechanism and kinetics of ice formation in living cells and optimizing the protocols of both cryopreservation and cryotherapy applications.

Supporting Information
Video S1 Video for twitching IIF.

(AVI)
File S1 Raw data of the normalized cell volume.