Originally published as Biophys J. BioFAST on October 8, 2004.
doi:10.1529/biophysj.104.048355
Biophysical Journal 88:647-660 (2005)
© 2005 The Biophysical Society
Kinetics of Intracellular Ice Formation in One-Dimensional Arrays of Interacting Biological Cells
Daniel Irimia * and
Jens O. M. Karlsson 
* Center for Engineering in Medicine and Surgical Services, Massachusetts General Hospital, Shriners Hospital for Children, and Harvard Medical School, Boston, Massachusetts 02129; and Woodruff School of Mechanical Engineering, and Petit Institute for Bioengineering and Bioscience, Georgia Institute of Technology, Atlanta, Georgia 30332
Correspondence: Address reprint requests to Dr. Jens O. M. Karlsson, Georgia Institute of Technology, Woodruff School of Mechanical Engineering, Atlanta, GA 30332-0405. Tel.: 404-385-4157; Fax: 404-385-1397; E-mail: karlsson{at}alum.mit.edu.
 |
ABSTRACT
|
|---|
Although cell-cell interactions are known to significantly affect the kinetics of intracellular ice formation (IIF) during tissue freezing, this effect is not well understood. Progress in elucidating the mechanism and role of intercellular ice propagation in tissue freezing has been hampered in part by limitations in experimental design and data analysis. Thus, using rapid-cooling cryomicroscopy, IIF was measured in adherent cells cultured in micropatterned linear constructs (to control cell-cell interactions and minimize confounding factors). By fitting a Markov chain model to IIF data from micropatterned HepG2 cell pairs, the nondimensional rate of intercellular ice propagation was found to be 
= 10.4 ± 0.1. Using this measurement, a new generator matrix was derived to predict the kinetics of IIF in linear four-cell constructs; cryomicroscopic measurements of IIF state probabilities in micropatterned four-cell arrays conformed with theoretical predictions (p < 0.05), validating the modeling assumptions. Thus, the theoretical model was extended to allow prediction of IIF in larger tissues, using Monte Carlo techniques. Simulations were performed to investigate the effects of tissue size and ice propagation rate, for one-dimensional tissue constructs containing up to 100 cells and nondimensional propagation rates in the range 0.1 
1000.
|
INTRODUCTION
|
|---|
The destruction of biological tissue by freezing has implications for attempts to develop minimally damaging cryopreservation protocols for storage of tissues and organs (Karlsson and Toner, 1996
), or maximally effective procedures for cryoablation of cancerous or other undesirable tissue (Gage and Baust, 1998). At rapid rates of cooling, cell injury is associated with intracellular ice formation (IIF) resulting from supercooling of the cytoplasm (Mazur, 1984; Toner, 1993). Thus, significant effort has been directed toward the development of mathematical models that can predict the probability of IIF (Mazur, 1965; Toner et al., 1990; Pitt et al., 1991; Muldrew and McGann, 1990, 1994), for the purpose of guiding the design of optimal freezing protocols (Toner et al., 1993; Karlsson et al., 1996). Until recently, most attempts to develop theoretical models of IIF have focused on the freezing of cell suspensions. The purpose of this study is to validate our model of IIF for confluent tissue (Irimia and Karlsson, 2002), and to present de novo predictions of the behavior of one-dimensional (1-D) tissue systems during freezing.
In the last decade, it has become apparent that the IIF process in tissue differs significantly from that in isolated cells, both qualitatively and quantitatively. For instance, freezing of intact tissue constructs results in higher rates and probabilities of IIF than is observed during freezing of cell suspensions obtained by disaggregating the tissue (McGrath et al., 1975; Porsche et al., 1991; Yarmush et al., 1992; Armitage and Juss, 1996; Acker et al., 1999; Irimia and Karlsson, 2002). Moreover, early anecdotal evidence, based on cryomicroscopic observation of the spatial pattern of IIF in plant and insect tissue, led to the hypothesis that the increased probability of IIF in tissue was due to the ability of intracellular ice to propagate from cell to cell (Stuckey and Curtis, 1938; Brown, 1980; Berger and Uhrik, 1996). Specifically, Berger and Uhrik (1996) hypothesized that gap junctions play a role in intercellular ice propagation, based on qualitative differences in the IIF pattern when tissue was treated with the gap junction inhibitor heptanol before freezing. Cryomicroscopy experiments after decoupling of gap junction communication using dinitrophenol (Berger and Uhrik, 1996) or low-calcium media (Acker et al., 2001) have also been reported; however, the nonspecific cytotoxicity associated with such treatments limited the ability to ascribe the resulting effects to gap junction inhibition alone. Even heptanol does not act directly on the connexon, but is thought to reduce the fluidity of cholesterol-rich membrane domains in which gap junctions are localized (Bastiaanse et al., 1993). To more definitively demonstrate that the increased propensity for IIF in tissue is due to cell-cell contact and gap junctional intercellular communication, we developed a cell micropatterning technique to directly elucidate the effects of cell-cell interactions, while allowing independent control of confounding factors such as cell-substrate interactions, cell geometry, and time in culture (Irimia and Karlsson, 2002, 2003). By comparing the IIF states in micropatterned one-cell and two-cell constructs during freezing, we were able to demonstrate that cell-cell contact increases the rate of IIF by enabling intercellular ice propagation (Irimia and Karlsson, 2002). Furthermore, we observed a reduction in the average rate of ice propagation after exposure to the highly specific gap junction blocker 18ß-glycyrrhetinic acid (Davidson and Baumgarten, 1988), supporting the hypothesis that active gap junctions are required for some, but not all mechanisms of propagation of intracellular ice (Irimia and Karlsson, 2002).
Despite the evident differences between IIF in isolated and cultured cells, theoretical models of IIF originally developed for cell suspensions (Toner et al., 1990; Karlsson et al., 1994) have been used without modification to describe the kinetics of IIF during freezing of tissue (Karlsson et al., 1993; Tsuruta et al., 1998; de Freitas and Diller, 2004). The ability of cell IIF models to fit experimental data for tissue has been limited in some systems (Tsuruta et al., 1998), but remarkably good in other systems (de Freitas and Diller, 2004). Nonetheless, because previous models of cell freezing do not account for intercellular ice propagation, they must be regarded simply as phenomenological descriptions of the kinetics of IIF when applied to tissue freezing; thus, the model parameters (e.g., nucleation rate coefficients) no longer have physical significance, but are useful only for adjusting the fit between predictions and experimental data. Acker and co-workers have presented a mechanistic model of the thermodynamics of ice propagation through gap junction channels (Acker et al., 2001). Their model predicts that there exists a critical level of cytoplasmic supercooling, below which intercellular ice propagation is thermodynamically possible. However, the predicted magnitude of the critical supercooling was found to be sensitive to the connexon biophysical properties, many of which are not accurately known (Acker et al., 2001). Furthermore, Acker and colleagues analyzed only static (equilibrium) effects, and did not attempt to model the kinetics of the intercellular ice propagation process. Thus, there is still a need for practical models of IIF in tissue, which will afford the ability to predict (and hence, improve) the results of freezing during cryopreservation or cryosurgery.
In a previous study, we have presented the first theoretical model of tissue freezing that explicitly describes the kinetics of ice propagation from cell to cell (Irimia and Karlsson, 2002). Our model assumed that IIF in a given cell within the tissue arose from a combination of two independent stochastic processes, spontaneous IIF (e.g., due to intracellular ice nucleation or membrane rupture) and intercellular ice propagation; the latter mechanism was dependent on the presence of intracellular ice in a neighboring cell, whereas the former was independent of such interactions. A Markov chain model was derived to describe the kinetics of IIF in an ensemble of cell pairs, and predictions were shown to be consistent with our experimental observations of IIF in micropatterned two-cell constructs (Irimia and Karlsson, 2002).
Here, we demonstrate how to scale up our previous model to predict the kinetics of IIF in tissue constructs containing any number of cells, using either Markov chain or Monte Carlo techniques. We present new cryomicroscopy measurements of IIF in micropatterned linear tissue constructs consisting of two or four HepG2 cells each. The rate of intercellular ice propagation was determined by fitting our theoretical model to cryomicroscopic data obtained using two-cell constructs; this represents the first quantitative measurement of the magnitude of this rate. The model was then validated by comparing a priori theoretical predictions with independent experimental IIF data obtained from four-cell constructs. Finally, a parametric analysis was undertaken to investigate the interaction between propagation rate and tissue size in larger constructs. Whereas, for the sake of simplicity, this study is restricted to 1-D tissue constructs, the methods and results described here can readily be extrapolated for analysis of intercellular ice propagation in two- or three-dimensional constructs (Irimia, 2002; Stott et al., 2004).
|
THEORETICAL BACKGROUND
|
|---|
Modeling assumptions
We have previously presented a continuous-time Markov chain model of the kinetics of IIF in an ensemble of cell pairs (Irimia and Karlsson, 2002). Here, we extend our analysis to 1-D cell arrays containing more than two cells. As before, each cell in the tissue construct is assumed to admit only two IIF states (unfrozen and frozen), and the average rate of IIF (J) in an unfrozen cell is assumed to be the sum of the rates associated with two independent stochastic processes: Jp, the average rate of intercellular ice propagation across the interface between a frozen and an unfrozen cell, and Ji, the total rate of IIF from mechanisms independent of the state of neighboring cells (e.g., intracellular ice nucleation); for simplicity, we will refer to the latter class of mechanisms as "spontaneous" or "interaction-independent" IIF. In an unfrozen cell flanked by k frozen cells, the total rate of IIF is
 | (1) |
at time t. In the general case, both the spontaneous IIF rate and the propagation rate may depend on temperature, time, position within the tissue, and a number of state variables (e.g., cell volumes and intracellular cryoprotectant concentrations) which must be described by additional dynamic models (Mazur, 1963; Kleinhans, 1998). For one- and two-dimensional tissue constructs, one can assume that all cells experience the same temperature and extracellular solute concentrations; thus, spatial variations in Ji and Jp can be neglected. Furthermore, by defining a nondimensional time
 | (2) |
the state equations that describe the system dynamics become invariant with respect to the mathematical model of Ji(t), including its coupling with the cell dehydration process (Irimia and Karlsson, 2002). The nondimensional form of the intercellular propagation rate is
 | (3) |
where, as a zeroth order approximation, can be assumed to be approximately constant.
Analytical solution of Markov chain model
With the above assumptions, the freezing of a tissue construct with N cells can be described by a discrete-state, continuous-time, homogeneous Markov chain. The IIF states X of the tissue construct comprise all distinct combinations of states (frozen or unfrozen) of the individual cells in the tissue. For example, all possible states and state transitions for a four-cell array are shown in Fig. 1 in the form of a directed graph. The total number of distinct states in a linear cell array is
 | (4) |
(where the notation int{z} indicates truncation to the nearest integer z).
View larger version (18K):
[in this window]
[in a new window]
|
FIGURE 1 Directed graph showing the possible IIF states (X = 0,...,9) and state transitions in a four-cell linear array. Open symbols represent unfrozen cells, whereas solid symbols represent frozen cells.
|
|
In an ensemble of identical tissue constructs, the elements of the state probability vector P represent the probabilities pX of each IIF state (X = 0,...,M – 1) at a given time, and the time evolution of this vector is given by a Kolmogorov differential equation
 | (5) |
where Q is the generator matrix containing the Poisson process rates for all state transitions. In practice, only M – 1 elements of P need to be included in Eq. 5, whereas the remaining state probability can be determined from the constraint 
pX = 1. Linear state equations such as Eq. 5 can be solved analytically using standard methods (Bay, 1999).
The solutions P(
) for a system of cell pairs (N = 2) were previously derived (Irimia and Karlsson, 2002). In this system, the probability p1 of the singlet state (the state in which only one of the two cells is frozen) was found to be sensitive to variations in over three orders of magnitude. Thus, the value of the nondimensional propagation rate can be estimated by fitting the analytical solution
 | (6) |
to experimental data obtained from two-cell constructs. Another important result from our previous analysis, generalized here to an N-cell tissue construct, is that the probability of the state X = 0, representing a tissue with all cells unfrozen, is given by
 | (7) |
Thus, for purposes of comparing theoretical predictions of P() to experimental measurements of the state probabilities as a function of temperature, pX(T), it is possible to estimate the value of from experimental data without knowing Ji:
 | (8) |
In this article, we present experimental measurements of the probabilities of IIF states for a four-cell linear array. The generator matrix describing the state transitions for states X = 0,...,8 of such a construct, as defined in Fig. 1, is
 | (9) |
The analytical solution to Eq. 5 for the four-cell system is presented in the Appendix.
Monte Carlo simulation
It is evident from Eq. 4 that the computational complexity associated with the Markov chain model is exponential in the tissue size N, and thus becomes prohibitive for larger tissue constructs. Therefore, we used a lattice Monte Carlo technique to simulate IIF in tissues with N > 4, as described below.
The tissue is described by a regular lattice in which each lattice site represents a cell, and intercellular ice propagation is allowed only between nearest neighbors. If the tissue is in a given IIF state X at time t, the probability that an unfrozen cell at lattice site j will freeze within a subsequent time interval 
t is described by a Poisson process, with a rate given by Eq. 1:
 | (10) |
where kj is the number of frozen nearest neighbors of cell j. Nondimensionalizing the above equation using Eqs. 2 and 3, and making the assumptions that and kj are constant during the time interval t, one obtains a homogeneous Poisson process:
 | (11) |
where 
(t + t) – (t). To enforce the validity of the assumption that kj is constant during , we required that the probability of an IIF event occurring anywhere in the tissue during this time interval be bounded by some small probability 
, resulting in the constraint
 | (12) |
where the sum is taken over all unfrozen cells in the tissue. For time intervals satisfying Eq. 12, the probability of multiple mechanisms of IIF occurring in the same cell j is much smaller than , and thus the probabilities of IIF by spontaneous and propagative mechanisms are independent, and are given by
 | (13a) |
 | (13b) |
respectively. Similarly, the overall probability of IIF in cell j can be approximated by
 | (14) |
For the Monte Carlo simulations in this study, time was discretized into uniform steps of size
 | (15) |
where n is the maximum number of nearest neighbors for the lattice (n = 2 for linear cell arrays). At each time step, a random number rj was drawn from a uniform distribution on the interval (0,1), for every unfrozen cell j. An IIF event was recorded for cell j at this time step if 
calculated using Eq. 14. Alternatively, to distinguish between the two the mechanisms of IIF, two random numbers were drawn for each unfrozen cell, and compared against the probabilities of spontaneous IIF and propagation, respectively, calculated using Eq. 13.
|
MATERIALS AND METHODS
|
|---|
Cell micropatterning
Experiments were performed on tissue constructs consisting of either two or four cells micropatterned in a linear arrangement on a glass substrate, using photolithography and surface modification techniques previously described (Irimia and Karlsson, 2002, 2003). Briefly, glass coverslips (Fisher Scientific, Pittsburgh, PA) were cleaned and coated with positive photoresist 1818 (Shipley, Marlborough, MA), which was patterned by exposure to ultraviolet light through a chrome-on-glass mask. After removal of the exposed photoresist, an aqueous solution of 5 mM PEG-disilane (Shearwater Polymers, Huntsville, AL) was applied to the substrates, which were then dried at 75°C for 5 min before removal of the unexposed photoresist. This procedure resulted in substrates with rectangular regions of unmodified glass, surrounded by a thin layer of PEG; the dimensions of these glass islands were 80 x 30 µm or 160 x 30 µm for two-cell and four-cell constructs, respectively.
The human hepatoma cell line HepG2 (American Type Culture Collection, Manassas, VA) was cultured at 37°C under a humidified 5% CO2 atmosphere, in minimum essential medium (Gibco BRL Life Technologies, Rockville, MD) supplemented with 10% v/v fetal bovine serum (Sigma-Aldrich, St. Louis, MO), 2.2 g/l sodium bicarbonate (Sigma), 1 mM sodium pyruvate (Sigma), 100 µg/ml streptomycin (Boehringer Mannheim, Indianapolis, IN), and 100 U/ml penicillin (Boehringer Mannheim); and subcultivated by disaggregation in a solution of 0.2% w/v trypsin (Gibco), 0.2% w/v glucose (Sigma), and 0.5 mM ethylenediamine-tetraacetic acid (EDTA; Sigma) in isotonic Ca2+ and Mg2+ free Dulbecco's phosphate buffered solution (PBS; Gibco).
For all experiments, HepG2 cells were harvested, and washed by centrifugation for 2.5 min at 200 x g and resuspension in versene solution (5 mM EDTA in Ca2+ and Mg2+ free PBS), followed by a second centrifugation (2.5 min at 200 x g) and resuspension in Ca2+ and Mg2+ free PBS. A 2-ml volume of cell suspension was seeded at a density of 1 x 105 cells/ml onto the patterned substrates in a 35-mm petri dish. The cells were incubated for 20 min at 37°C to allow cells to adhere to the exposed regions of the glass surface, after which nonadherent cells were removed by flushing the coverslip with culture media. After 24 h of culture, the substrates were washed with a solution of PBS, and cell nuclei were stained by exposure to a 2-µM solution of the fluorescent dye SYTO13 (Molecular Probes, Eugene, OR) for 10–15 min at 37°C, to aid in identification of the position and number of cells in the cell arrays. A representative four-cell construct, under bright-field and epifluorescence illumination, is shown in Fig. 2.
View larger version (94K):
[in this window]
[in a new window]
|
FIGURE 2 Representative micrographs of a micropatterned tissue construct containing four HepG2 cells in a linear array: (A) bright field; (B) epifluorescence illumination (SYTO 13 dye). Scale bars are 20 µm.
|
|
Cryomicroscopy
The cryomicroscopy system and techniques are described in detail elsewhere (Irimia and Karlsson, 2002). Briefly, an upright light microscope (E-600; Nikon, Tokyo, Japan) was fitted with a freezing stage (BCS 196; Linkam Scientific Instruments, Tadworth, UK) and an analog video camera (DVC-0A; Digital Video Camera Company, Austin, TX). The sample was cooled by heat conduction to a silver block in contact with a nitrogen gas stream, which was precooled using a coiled copper tubing heat exchanger immersed in liquid nitrogen. The stage temperature was regulated by a feedback control system using an electrical resistance heater and a platinum-resistance temperature detector (Pt100, 1/10 DIN Class B), which were both imbedded in the silver block. The temperature sensor was calibrated using the melting point of ice crystallized from a sample of purified water (Elix water purification system; Millipore, Bedford, MA). The system is able to control temperatures to a precision of ±0.1°C, and the spatial uniformity of the temperature across the microscope field of view is better than ±0.1°C.
Samples were prepared by dispensing a 2-µl droplet of PBS containing 2 µM SYTO13 onto an 18-mm circular glass coverslip (Fisher), then inverting the patterned substrate and placing it on the droplet, such that the cultured cells were sandwiched between the two coverslips. The samples were immediately used for cryomicroscopy experiments. The samples were initially cooled to –1.8°C, and held isothermally at this temperature during seeding of ice in the extracellular liquid. Tissue constructs comprising either two or four cells were located using epifluorescence illumination to identify SYTO13-stained nuclei, and the remainder of the experiment was recorded on video tape under bright-field transillumination. Samples were cooled to –60°C at a controlled rate of 130°C/min to induce IIF in the cell arrays. The IIF temperature for each cell was determined during frame-by-frame playback of the video recording by noting the temperature at which the characteristic darkening of the cytoplasm (resulting from scattering of the transilluminating light by intracellular ice crystals) could be observed. These IIF temperature data could be used to determine the state probabilities, the overall probability of IIF in the cell arrays, and the probabilities of IIF at each location within the arrays, as a function of temperature.
Statistical analysis
A nonlinear curve fit to experimental data from two-cell constructs was performed using the Levenberg-Marquardt algorithm. Agreement between experimental data from four-cell constructs and independent theoretical predictions was evaluated using the Kolmogorov-Smirnov test. Distributions of IIF temperatures for various cell positions were compared using analysis of variance and Student's paired t-test. Effects were considered significant at the 95% confidence level (p < 0.05).
|
RESULTS AND DISCUSSION
|
|---|
Measurement of rates of spontaneous intracellular ice formation and propagation in two-cell constructs
In a previous study, we measured the kinetics of the IIF state probabilities in a population of 164 micropatterned two-cell constructs during rapid freezing (Irimia and Karlsson, 2002). We now present additional cryomicroscopic observations, representing a total of 246 micropatterned HepG2 pairs, and use these data to estimate the magnitude of the rates of IIF by interaction-independent mechanisms (Ji) and intercellular ice propagation (Jp). The rate of interaction-independent IIF was estimated from experimental measurements of the fraction of unfrozen cell pairs (p0), using Eqs. 2 and 8. Fig. 3 A illustrates the relationship between temperature and nondimensional time described by Eq. 8. When comparing or fitting model predictions (in which kinetics are described in terms of the nondimensional time variable ) to experimental data (which describe probabilities as a function of temperature), Fig. 3 A was used to convert between units of temperature and nondimensional time. Moreover, because the temperature history T(t) used in our freezing experiments is known (T = 271.4 K – 2.17 K s–1 x t), it is also possible to use Fig. 3 A to convert between time t and nondimensional time , as shown. The use of dimensionless variables in the model formulation has several benefits. For example, it allowed the theoretical description of intercellular ice propagation to be developed without making any assumptions about the physical mechanisms underlying the interaction-independent IIF process, or the functional form of Ji(t). In addition, it makes possible the extrapolation of our model predictions to cell types other than HepG2, or to different freezing procedures that may alter the kinetics of interaction-independent IIF; whereas the absolute timing of IIF events depends on the exact form of Ji(t) for a given cell type and cooling protocol, the relative sequence and spatial distribution of IIF events depends only on the magnitude of the nondimensional propagation rate (), and will not be affected by changes in Ji. Physically, one can interpret the value of the variable as the average number of interaction-independent IIF events per cell, subject to the hypothetical condition that cells be allowed to freeze multiple times.
Because represents the integral of the rate of interaction-independent IIF (see Eq. 2), it is possible to determine the value of Ji by taking the derivative of the function (t). Thus, the data in Fig. 3 A were numerically differentiated using the central difference method, followed by a 25-point boxcar filter to suppress noise. The resulting estimate of the temperature-dependence of Ji is shown in Fig. 3 B. It is evident that the rate of interaction-independent IIF increases with decreasing temperature, and that its magnitude is on the order 0.1–0.3 s–1 over the temperature range observed. We have previously established that the magnitude and temperature dependence of Ji are the same in two-cell HepG2 constructs as in single-cell constructs of adherent HepG2 (Irimia and Karlsson, 2002). Thus, we believe that the rate of interaction-independent IIF for a given cell type is not affected by the number of cells in the tissue construct, and that Ji(T) will be the same if tissue culture conditions and freezing procedures are reproducible. Therefore, it should be noted that the data in Fig. 3 are applicable only to HepG2 cells cultured 24 h on glass substrates, and frozen at a rate of 130°C/min from an initial temperature of –1.8°C.
In our earlier study, we demonstrated that in an ensemble of two-cell constructs, observations of the transient singlet state (i.e., constructs with one frozen and one unfrozen cell) could be used to estimate the magnitude of the nondimensional intercellular ice propagation rate (Irimia and Karlsson, 2002). At that time, we determined an approximate value of the propagation rate for HepG2 cells ( 
5) by estimating the peak magnitude of the probability of the singlet state. Here, we present a more accurate estimate of for this system, by performing a least-squares curve fit of Eq. 6 to new cryomicroscopy data. Thus, IIF was observed in an ensemble of 246 micropatterned cell pairs; the fraction of cell pairs found in the singlet state is shown as a function of temperature in Fig. 4, together with the best fit of Eq. 6 to the experimental data (where was converted to temperature using Fig. 3 A); the nonlinear regression yielded a best fit value = 10.4 ± 0.1 (R2 = 0.77). Thus, in HepG2 cell pairs, the rate of intercellular ice propagation is an order of magnitude larger than the rate of interaction-independent IIF. Considering the magnitude of Ji as estimated in Fig. 3 B, intercellular ice propagation across the interface between a frozen and an unfrozen HepG2 cell occurs with an approximate frequency of 1–3 Hz per interface.
View larger version (17K):
[in this window]
[in a new window]
|
FIGURE 4 Probability of the singlet state (the state in which one cell is frozen and the other is unfrozen), measured in an ensemble of micropatterned constructs of HepG2 cell pairs during freezing at a rate of 130°C/min (circles), and predicted using the best fit of Eq. 6 (solid line).
|
|
Ice propagation has been shown to be mediated in part by intercellular gap junctions (Berger and Uhrik, 1996; Irimia and Karlsson, 2002). However, HepG2 is a transformed hepatoma cell line that is known to have impaired gap junctional intercellular communication due to low expression levels and aberrant localization of connexins (Yang et al., 2003). Therefore, it is likely that the nondimensional propagation rate in primary hepatocytes and other cell species that express higher numbers of functional gap junctions is significantly >10.4.
Model validation using four-cell construct
The above data for two-cell constructs were used to make predictions about the kinetics of IIF in independent experiments on four-cell constructs, to validate our model of intercellular ice propagation. Equations 5 and 9 were solved analytically subject to the initial condition p0 = 1 (see Appendix), and the predicted probabilities of IIF states X = 0,...,8 (defined in Fig. 1) for = 10.4 are shown in Fig. 5 (the probability of the fully frozen state is not shown, because p9 is not an independent variable). To compare these predictions to cryomicroscopy data, experimentally measured IIF temperatures had to be converted to dimensionless times. As shown in Eq. 2, the conversion to dimensionless time depends only on the interaction-independent rate of IIF (Ji). Based on our previous finding that the magnitude of Ji(t) is indistinguishable in micropatterned one- and two-cell constructs (Irimia and Karlsson, 2002), we assumed that this rate would also be the same in four-cell constructs frozen using the same temperature profile. Thus, we converted IIF temperatures from these four-cell experiments to dimensionless times using the correlation (T) shown in Fig. 3 A, which was constructed using only data from two-cell constructs. This allowed us to make a priori predictions of the IIF kinetics in four-cell arrays, without using any of our experimental data from this system to estimate model parameters. Thus, experimental data from cryomicroscopic observations of 71 four-cell constructs are compared against independent theoretical predictions in Fig. 5. The Kolmogorov-Smirnov test was used to evaluate model accuracy: experimental measurements for all of the IIF states (X = 0,...,8) shown in Fig. 5 were found to conform to the theoretically predicted distribution for the corresponding state probabilities (p < 0.05). The agreement between theoretical predictions and experimental data for the unfrozen state probability p0() validates our assumption that Ji is approximately the same in four-cell arrays as in single-cell and two-cell constructs. Similarly, the agreement between theory and experiment for the intermediate, partially frozen states (X = 1,...,8) lends support to the validity of our model of intercellular ice propagation, and to the accuracy of our measurement of the nondimensional propagation rate .
The experimental observations are consistent with the theoretical result that among the intermediate states, states X = 3 and X = 7 should be the most prevalent, whereas states X = 4 and X = 6 should occur with the lowest peak frequency. The latter are the only two states that can only form via two consecutive spontaneous IIF events in the same tissue construct, which occurs only rarely (because the propagation rate is an order of magnitude larger than the rate of spontaneous IIF). In contrast, among all states with either two or three frozen cells, X = 3 and X = 7 are the only states with a single propagation "front" (interface between an unfrozen and a frozen cell), respectively; thus these states will be consumed more slowly than others, in which multiple sites are available for intercellular ice propagation. Note that if there were no intercellular ice propagation ( = 0), then the relative magnitudes of the state probabilities would be drastically different from the experimentally observed relationships, as shown in Fig. 5 (dashed lines). The fact that our model captures both the quantitative and qualitative features of the experimental data strongly suggests that IIF in tissue constructs is due to the mechanisms we have proposed (Irimia and Karlsson, 2002).
Spatial distribution of intracellular ice in linear tissue constructs
Because cells at the termini of a linear tissue construct have fewer cell-cell interactions than cells in the interior of the construct, we hypothesized that the probability and rate of IIF for edge cells would be lower than for cells in the center of the tissue construct. Indeed, cryomicroscopic observation of the freezing of four-cell constructs revealed that the average IIF temperature in the two distal cells (–20.6 ± 0.5°C) was somewhat lower than the average temperature of IIF in the central cells (–20.0 ± 0.5°C), although the difference between the means was not statistically significant. To further analyze the spatial distribution of IIF, we compared the kinetics of IIF at different absolute positions within the linear cell array. The four cells were labeled as shown in the schematic in Fig. 6, where the orientation of each tissue construct was specified such that either cell A or cell B corresponded to the location of the first IIF event observed. It was possible to uniquely identify the array orientation in all but five cases, in which the first observable state was X = 7; these cases were not included in the analysis. The cumulative incidence of IIF in each of the four cells in the array is shown as a function of temperature in Fig. 6. When the data are analyzed in this way, the cell position does have a statistically significant effect on the mean IIF temperature, as determined by analysis of variance (p < 0.05). Furthermore, analysis of the difference between adjacent cells within each tissue construct revealed that the average paired difference between IIF temperatures in cells A and B was not significantly different from zero, whereas IIF did, on the average, occur at a higher temperature in cell C than in cell D for a given construct (p < 0.05). These experimental observations agreed with theoretical predictions obtained using Eqs. A17–A20 (with dimensionless time converted to temperature using Fig. 3 A), as shown in Fig. 6 and confirmed using a Kolmogorov-Smirnov test (p < 0.05). Thus, a gradient in the probability of IIF is expected within the half of the array in which the first spontaneous IIF event did not occur (i.e., cells C and D), presumably because the directed propagation of intracellular ice from the other array half (cells A and B) dominates any contribution of spontaneous ice formation. Note that if the orientation of the linear arrays is not considered, this spatial dependence will manifest itself as a decreased probability of IIF in both edge cells relative to the probability of IIF in the central cells, due to symmetry.
View larger version (19K):
[in this window]
[in a new window]
|
FIGURE 6 Measured cumulative probability of IIF in cell A (solid squares), cell B (open squares), cell C (open circles), and cell D (solid circles) in micropatterned linear arrays of four HepG2 cells during freezing at a rate of 130°C/min. Also shown are theoretical predictions of the cumulative incidence of IIF in cells A (dotted line), B (dash-dotted line), C (dashed line), and D (solid line). The inset schematic identifies the position of each cell within the array; the orientation of the array was defined such that the first observed IIF event corresponded to either location A or B.
|
|
To explore the phenomenon of nonuniform probability of IIF in larger tissue constructs, we used Monte Carlo techniques to simulate ice propagation in linear arrays containing as many as 100 cells. Whereas the accuracy of the predictions depends on the magnitude of the step size used in the simulation, the effect of the value of used in Eq. 15 was determined in a four-cell linear array with = 10.4. We characterized the resulting predictions by calculating the median nondimensional time of IIF (50), i.e., the instant at which intracellular ice is present in 50% of all cells in an ensemble of tissue constructs with identical properties. For our four-cell system, the theoretical value of the median IIF time is 50 = 0.2820, as determined using the analytical solution for the kinetics of IIF, which is presented in the Appendix. Table 1 shows the step size and the corresponding median times of IIF in an ensemble of 2.5 x 105 identical constructs (i.e., a total of 106 cells), as predicted using various values of . As expected, the accuracy of the predictions improved with decreasing , due to the reduced probability of multiple IIF events within a single time step; note that although nominally represents the probability of IIF within the interval , reasonably accurate results can be obtained even by setting > 1, because Eq. 15 represents a conservative lower bound on the required step size. For = 0.1, the median IIF time could be predicted with an accuracy and precision both <0.5%; thus, this value was used for all further Monte Carlo simulations. Due to the trade-off between accuracy and computational cost, simulations with higher precision were not practical. To overcome this problem, it is possible to use variable-size time steps calculated using Eq. 12; a significant improvement in performance would be expected for systems with high propagation rates. Using an alternative approach, we have recently adapted the stochastic simulation algorithm developed by Gillespie (1976) for prediction of chemical reaction kinetics, to make our Monte Carlo simulations of IIF more accurate and computationally efficient (e.g., Stott et al., 2004). However, because this algorithm assumes that the reaction rate coefficients are constant, a conventional simulation technique like the one presented here may be required for simulation of more complex phenomena such as IIF with temperature-dependent ice propagation rates.
Monte Carlo simulations were run for ensembles of tissue constructs with sizes N = 1, 3, 10, 30, and 100, for values of the nondimensional propagation rate ranging over four orders of magnitude. Each ensemble comprised 1 x 105–2 x 105 cells, corresponding to 2 x 103–3 x 104 constructs, depending on tissue size. The resulting probabilities of IIF at the center and edge of the cell array were compared at the median time of IIF for the entire ensemble. If the kinetics of IIF were spatially uniform, then the probability of IIF at the median time of IIF should be 50% throughout the tissue, including both the central and distal cells, by definition (Fig. 7, dotted line). However, as shown in Fig. 7, the local probability of IIF in the center of the tissue is greater than the cumulative incidence of IIF in the entire tissue, and the local probability of IIF at the tissue edges is smaller than the overall probability of IIF in the cell array. Thus, paradoxically, gradients in the probability of IIF in a tissue construct comprising interacting cells are expected to occur even if the underlying rates of spontaneous IIF and propagation (Ji and Jp) are spatially uniform. This phenomenon may confound analysis of the effects of temperature- and concentration-gradients on the probability of IIF in tissue (Diller, 1992), or the interpretation of the cause of spatially nonuniform tissue injury (Merchant et al., 1996; Rabin et al., 1997; Rupp et al., 2002). It is interesting to note that as the nondimensional propagation rate increases, the disparity in the incidence of IIF between central and distal cells initially increases to a maximum value, and subsequently decreases. For low propagation rates, IIF is dominated by spontaneous freezing events, the probability of which has been assumed to be spatially uniform, and therefore no difference in the response of central and distal cells is expected in the limit 
0. Conversely, in the limit 
, intercellular ice propagation will be instantaneous; thus, because all cells in the tissue construct will freeze simultaneously upon the first spontaneous IIF event, the probability of IIF should be the same at each location within the tissue. The value of that results in the maximum gradient in the incidence of IIF depends on the number of cells in the tissue construct. The larger the construct is, the larger the range of for which nonuniformity in IIF continues to increase with increasing , and the larger the maximum discrepancy between central and distal cells. Therefore, we expect this effect to be highly significant in macroscopic tissues and tissue-engineered devices, in which the number of cells is orders of magnitude larger than in the cell arrays simulated in this study.
View larger version (19K):
[in this window]
[in a new window]
|
FIGURE 7 Predicted probability of IIF in the center (open symbols) and edge cells (solid symbols) at the median IIF time (50), for various values of the dimensionless propagation rate and tissue size: N = 3 (circles); N = 10 (squares); N = 30 (triangles); N = 100 (diamonds). The reference line (dotted line) represents the probability that would be observed if IIF were spatially uniform.
|
|
Parametric analysis: effect of propagation rate and tissue size
Monte Carlo simulations were used to investigate the effect of tissue size and propagation rate on the overall kinetics of IIF in the tissue construct. Whereas the median time of IIF (50) in an ensemble of cell arrays is inversely proportional to the time-averaged total rate of IIF in the cell population, this metric was used to quantify the kinetics of ice formation. To separate effects due to the number of cells in a construct from edge effects such as those demonstrated in Figs. 6 and 7, we simulated freezing of 1-D tissue constructs with both open and periodic boundary conditions. In the latter case, the tissue topology is that of a closed ring, resulting in spatially uniform probability of IIF for any tissue size or propagation rate . In Fig. 8, the median nondimensional time of IIF as a function of construct size is shown for = 0.1, 1, 10, 100, and 1000. Each ensemble consisted of 103-105 constructs (5 x 104-2 x 105 cells). As expected, the median IIF time decreased with increasing propagation rate, approaching the theoretical limit 50 = ln(2)/N as . In this limiting case, the total rate of IIF is proportional to the number of cells in the construct, because a spontaneous IIF event in any one of the cells will instantly cause freezing of all N cells in the array. On the other hand, in the limit 0, the expected median IIF time is 50 = ln(2) = 0.69, independent of the number of cells in the tissue; as seen in Fig. 8, results for < 1 approached this limiting case. For intermediate finite values of , the median IIF time was inversely proportional to N for small construct sizes, but independent of N for larger arrays. The value of N that demarcated the regimes of size-dependent and size-independent IIF kinetics increased with increasing propagation rate. In the size-independent regime (large N), the overall kinetics of IIF in the tissue were rate limited by the intercellular ice-propagation process. Conversely, in the size-dependent regime (small N), the kinetics of IIF were equivalent to those that would result from instantaneous propagation of intercellular ice. A study by Zhang et al. (2003) has claimed that intercellular ice propagation at a rate = 10.4 can be considered "instantaneous" in a four-cell construct. However, it is evident from Fig. 8 that if N = 4, the dimensionless propagation rate must be on the order 103 to be considered approximately instantaneous (i.e., yielding kinetics comparable to the limiting case ). As we have recently demonstrated, this discrepancy is due to errors in the theoretical analysis by Zhang and colleagues (Karlsson, 2004).
View larger version (16K):
[in this window]
[in a new window]
|
FIGURE 8 Predicted median nondimensional time of IIF (50) in one-dimensional cell arrays with open (solid lines) or periodic boundary conditions (dashed lines), as a function of construct size and dimensionless propagation rate. Also shown are theoretical predictions for the limiting case of instantaneous ice propagation (dotted line).
|
|
The boundary conditions on the one-dimensional cell arrays had only a marginal effect on the overall kinetics of IIF (see Fig. 8). The median IIF time was nominally smaller for periodic boundary conditions than for open boundary conditions, but the difference was never >15%. The kinetics of IIF are expected to be slower in open-ended linear arrays, because intercellular ice propagation must stop when the propagation front reaches the tissue edge. Mirroring our observations in Fig. 7, edge effects (i.e., observed differences between the values of 50 for open and periodic boundary conditions) initially increased with increasing propagation rates, then decreased with increasing propagation rates after reaching a maximum magnitude at an intermediate value of , which was dependent on tissue size.
We also investigated the relative importance of the two mechanisms of IIF (propagative versus spontaneous IIF) in linear tissue constructs. As shown in Fig. 9, the percentage of cells frozen as a result of intercellular ice propagation generally increased with increasing and with increasing construct size. Because at least one cell per construct must freeze spontaneously (e.g., via nucleation) before intercellular ice propagation is initiated, the fraction of cells frozen via propagative mechanisms must approach (N – 1)/N in the limit ; thus, for large , the cumulative probability of propagative IIF is primarily a function of tissue size. On the other hand, for N > 10, the fraction of cells frozen due to intercellular ice propagation was relatively insensitive to tissue size at all values of the nondimensional propagation rate. It is noteworthy that in this size-independent regime, over 70% of cells were predicted to freeze as a result of ice propagation if > 10 (the putative physiological range). This suggests that, for most practical applications of tissue freezing, intercellular ice propagation is the dominant mechanism of IIF.
View larger version (24K):
[in this window]
[in a new window]
|
FIGURE 9 Predicted percentage of IIF events resulting from intercellular ice propagation, as a function of the dimensionless propagation rate and the number of cells in the tissue construct. Constructs were linear with open boundary conditions.
|
|
Acker and McGann have recently proposed the hypothesis that propagative IIF is innocuous, and that cell damage is only associated with spontaneous IIF (Acker and McGann, 2002, 2003). On the other hand, our results above suggest that if IIF occurs in a confluent tissue with > 10, then the probability that a freezing event was spontaneous is significantly less likely than the probability that it was due to intercellular ice propagation. Taken together with Mazur's two-factor model of cryoinjury (Mazur et al., 1972), a corollary of these two assertions is that the expected drop in tissue survival when increasing the cooling rate beyond its optimal value should be no more than 30% if 10, and significantly less for larger . This conclusion follows from the fact that the probability of cell injury due to non-IIF mechanisms (so-called "solution effects") decreases with increasing cooling rates, and our result above that the probability of spontaneous IIF in macroscale tissue always remains <30% for > 10. Unfortunately, however, results from the literature do not support the above corollary: in contrast, the difference in viability of tissues frozen at optimal and supraoptimal rates of cooling is routinely found to be >50% (e.g., Song et al., 1995; Pasch et al. 1999; Pegg, 2002); for sufficiently large cooling rates, tissue destruction is complete (Borel Rinkes et al., 1992). Thus, either the rate of intercellular ice propagation in other tissue types is significantly less than in cultured HepG2 constructs (yielding a higher probability of spontaneous IIF at rapid cooling rates), or IIF is associated with cell injury whether it occurs spontaneously or by propagation from neighboring cells. We have already concluded that the rate of intercellular ice propagation presently measured in HepG2 constructs is likely to be lower, not higher, than values expected for other tissues types. In light of this conclusion, the claim by Acker and McGann that propagative IIF is nonlethal appears to be inconsistent with the known patterns of cryoinjury in rapidly frozen tissue. This discrepancy may be due in part to the fact that in the studies by Acker and McGann (2002, 2003) samples were never frozen below –40°C before thawing and viability assessment, whereas in practical cryopreservation applications, the final freezing temperature is typically –196°C. It is no surprise then, that the correlation between IIF and cell injury in frozen tissues (Karlsson et al., 1993; Bischof et al., 1997; Zieger et al., 1997) appears to be just as strong as for frozen suspensions of isolated cells (Mazur, 1984; Toner, 1993).
|
CONCLUSIONS
|
|---|
This study strongly supports the hypothesis that the known correlation between IIF and cell-cell contact in tissue is due to the ability of intracellular ice to propagate from a cell to its neighbors, by demonstrating that the model we previously developed to describe the kinetics of intercellular ice propagation in two-cell constructs yields accurate a priori predictions of the dynamic behavior of a four-cell system. The success of this theoretical model also validates its main assumptions, namely that the rate of spontaneous IIF in each cell is not affected by cell-cell interactions, and that the nondimensional propagation rate can be taken to be approximately constant. The latter conclusion implies that in HepG2 cells frozen rapidly, the temperature dependence of the rates of intercellular ice propagation and spontaneous IIF are equivalent; the two rates differ only by the multiplicative factor . Considering the fact that the rate of spontaneous IIF (Ji) can be described using a nucleation model (Karlsson, 2004), the absolute rate of intercellular ice propagation (Jp) presumably has a strong dependence on cytoplasmic supercooling. Nonetheless, to fully understand the mechanisms of ice propagation, and to develop models that can be used to optimize tissue freezing protocols, further work is required to accurately characterize the dependence of the intercellular ice propagation rate on temperature and cell water content.
By extrapolating model predictions to larger tissue systems, we were able to demonstrate that intercellular ice propagation is the dominant mode of IIF in rapidly frozen tissues, even in those with modest propagation rates, such as HepG2 cultures. In the context of previously published data on viability of rapidly frozen tissues, this observation reinforces the well-established paradigm that IIF causes irreversible cell damage, whether it occurs spontaneously or by propagation. We have also demonstrated that intercellular ice propagation can result in gradients in the probability of IIF (and, hence, in cell viability) during freezing of tissues, even in the absence of any nonuniformities in tissue properties or in temperature or concentration fields. In general, the temporal and spatial distributions of IIF events are determined by a complex interrelationship between tissue size and nondimensional propagation rate, although there also exist regimes in which the response of the tissue is relatively insensitive to either the number of cells or to the magnitude of the dimensionless propagation rate.
|
APPENDIX
|
|---|
Equation 5 was solved using the four-cell generator matrix given by Eq. 9, subject to the initial condition p0(0) = 1 and pX(0) = 0 (X = 1,...,8), and assuming a constant nondimensional propagation rate . If 
1/2, 1, 2, and 3, then the analytical solution for 
0 is:
 | (A1) |
 | (A2) |
 | (A3) |
 | (A4) |
 | (A5) |
 | (A6) |
 | (A7) |
 | (A8) |
 | (A9) |
The cumulative probability of IIF in the four-cell array is given by:
 | (A10) |
To characterize the overall rate of IIF, we computed the median nondimensional time of IIF (50) by finding the time point for which pIIF(50) = 0.5.
The formulation above does not distinguish between the two possible orientations of the asymmetric IIF states. However, to investigate the effect of cell position relative to the initial site of IIF, it was necessary to define three new states X = 4', 7', and 8', corresponding to the mirror images of states X = 4, 7, and 8, respectively, as defined in Fig. 1. Thus, the following modifications were required to the state equations defined by Eqs. 5 and 9:
 | (A11) |
 | (A12) |
 | (A13) |
Furthermore, the following equations describe the evolution of the mirror image states:
 | (A14) |
 | (A15) |
 | (A16) |
After solving the corresponding system of differential equations, the probabilities of IIF at cell positions A, B, C, and D (as defined in Fig. 6) were determined from the state probabilities as follows:
 | (A17) |
 | (A18) |
 | (A19) |
 | (A20) |
|
ACKNOWLEDGEMENTS
|
|---|
This work was supported in part by the National Science Foundation (Faculty Early Career Development Award BES-0242377), and by the Georgia Tech-Emory Center for the Engineering of Living Tissues, a National Science Foundation Engineering Research Center (EEC-9731643).
Submitted on June 28, 2004;
accepted for publication September 27, 2004.
|
REFERENCES
|
|---|
TOP
ABSTRACT
INTRODUCTION
