A Nonparametric Method for Estimating Rates of Intracellular Ice Formation Due to Independent Competing Mechanisms

Karlsson, Jens O. M.

doi:10.1115/1.4053109

Abstract

The probability of intracellular ice formation (IIF) has conventionally been analyzed by counting the cumulative number of IIF events observed in a cell population, and normalizing to the total cell count to estimate the cumulative IIF probability. However, this method is invalid when attempting to distinguish among multiple, independent IIF mechanisms, because of confounding effects due to competition for a finite pool of unfrozen cells. Therefore, an alternative approach for analyzing IIF data is proposed, based on treating IIF as a marked point process, in which the points represent IIF events and the marks represent different mechanisms of IIF. Using the new method, it is possible to quantify the kinetics associated with any IIF mechanism for which corresponding events can be marked (i.e., experimentally distinguished from competing IIF mechanisms). The proposed approach is nonparametric, making possible characterization of IIF mechanisms that have not yet been fully elucidated. The new analytical approach was compared to the conventional method of IIF analysis using data from a simulated experiment, demonstrating that the new method yielded superior estimates of the cumulative distribution function of IIF times when two competing mechanisms of IIF were active. The proposed algorithm was also applied to cryomicroscopic IIF observations in adherent endothelial cells, yielding rate estimates for two concurrent IIF processes. Furthermore, a proof is presented to demonstrate that when the proposed data analysis algorithm is applied to IIF data from a single mechanism of IIF, the results are equivalent to those obtained by the conventional method of analysis.

Introduction

The rational design of optimal freezing (and warming) procedures for cryopreservation (or cryoablation) of cells or tissue requires prediction of the probability of cell damage associated with intracellular ice formation (IIF). In particular, because IIF is known to be a stochastic process, it is important to model the kinetics of IIF, and to measure the rate constants associated with these kinetics. To ensure that IIF models have predictive power under conditions that differ from those of experiments used to estimate model parameters, a mechanistic approach to modeling is required—i.e., mathematical models should be based on the physicochemical mechanisms that are responsible for the IIF phenomenon. However, various pathways of IIF exist, each with a distinct physicochemical mechanism; thus, in a given system, theoretical models must typically be elucidated (and kinetic parameters estimated) for multiple IIF pathways. A prerequisite for both model development and parameter estimation is the ability to rigorously analyze experimental observations of IIF to obtain valid quantitative estimates of the kinetics associated with each IIF mechanism active in the system.

The existence of multiple mechanisms for ice formation in cells complicates the estimation of the rate associated with any one IIF pathway; all previous attempts to deal with this problem in the cryobiology literature have been limited to parametric regression techniques. For example, Pitt et al. [1] proposed the presence of two distinct IIF pathways (corresponding to a time-dependent and a time-independent mechanism), and quantified the probability distributions associated with each mechanism by fitting a parametric mathematical model to the experimentally observed distribution of IIF times (and temperatures). Similarly, the kinetics of IIF associated with surface-catalyzed nucleation (SCN) and volume-catalyzed nucleation (VCN) have been quantified by fitting parametric models of the two nucleation rates to experimentally measured IIF probability distributions [2,3]. In multicellular tissue, we have previously quantified the relative magnitude of the rates of IIF associated with two hypothesized mechanisms (intercellular ice propagation versus spontaneous IIF) by approximating the ratio of the two rates as a constant parameter, the value of which was estimated by fitting experimental data with a mathematical model that contained this nondimensional parameter [4]. However, when there is insufficient understanding about a putative mechanism of IIF to justify development of a mechanistic mathematical model, parametric regression techniques cannot be used to extract information about the kinetics of the proposed IIF mechanism. For example, Scheiwe and Körber [5] hypothesized that when observing cell freezing using conventional video cryomicroscopy techniques, events manifesting as cell darkening and cell twitching represent two distinct pathways of IIF; however, they did not propose mathematical models corresponding to these two categories of IIF.

Previous attempts at applying nonparametric approaches to the analysis of darkening and twitching IIF (e.g., comparison of the corresponding mean IIF temperatures, or comparison of the total incidence of IIF events by each mechanism) have been faulty, inasmuch as the confounding effects of competition between the two mechanisms have been neglected [5]. Because the two IIF mechanisms act on the same pool of unfrozen cells, a naïve analysis gives the appearance of correlation even between independent mechanisms, e.g., when the fraction of cells exhibiting twitching-type IIF increases, the fraction of cells exhibiting darkening-type IIF necessarily decreases, even if the two manifestations of IIF are caused by mechanistically independent processes. Past attempts at analysis have fallen prey to this paradox, by failing to account for the reduction of the pool of unfrozen cells due to competing IIF mechanisms [5,6].

The lack of appropriate nonparametric techniques for analysis of systems with multiple, competing IIF mechanisms is an obstacle to the understanding of IIF mechanisms in freezing of tissue and tissue engineered constructs. High-speed video cryomicroscopy techniques pioneered by our group have led to the discovery of several previously unknown mechanisms of IIF in tissue constructs [7,8]; however, these new IIF mechanisms are not amenable to parametric analysis, because the processes have yet to be fully elucidated and mathematically modeled. To make possible rigorous analysis of the measured kinetics associated with each of these putative IIF mechanisms, a new nonparametric method for estimating the rates of IIF caused by multiple, competing mechanisms is proposed here. This work summarizes the conventional approach for quantifying IIF kinetics, describes the proposed new technique for data analysis, and demonstrates that the two methods yield equivalent results when there is only one mechanism of IIF (or equivalently, when no attempt is made to distinguish among multiple mechanisms). A simulated experiment is used to demonstrate the advantages of the new algorithm over the conventional method of analyzing IIF kinetics when multiple mechanisms of IIF compete for the same pool of unfrozen cells. As a further illustration of the new analytical technique, our method is applied to experimental data from high-speed video microscopy studies of adherent bovine endothelial cells [7], demonstrating the estimation of kinetic data for independent competing mechanisms of IIF.

Conventional Method for Analysis of Intracellular Ice Formation Kinetics

Experimental measurements of IIF kinetics yield a running tally of the number of cells that have experienced an intracellular freezing event during the course of observation. Typically, these data are reported as a plot of the cumulative probability of IIF (P_IIF) as a function of time elapsed (t) or temperature (either may be used as the random variable when the cooling rate is nonzero, and conversion between time and temperature is straightforward if the experimental temperature profile is known); the conventional method of data analysis uses the following relationship to compute the cumulative IIF probability:

P_{IIF} (t) = \frac{N - N_{u} (t)}{N}

(1)

where N_u(t) represents the number of unfrozen cells remaining at time t; N, the initial value of N_u(t). The probability distribution obtained in Eq. can be interpreted as a nonhomogeneous Poisson process, characterized by a time-dependent rate of IIF (J_IIF):

P_{IIF} (t) = 1 - exp (- \int_{0}^{t} J_{IIF} (τ) d τ)

(2)

If multiple, independent mechanisms of IIF are active, then J_IIF represents the composite arrival rate of IIF events:

J_{IIF} (t) = \sum_{μ} J_{IIF}^{μ} (t)

(3)

where each individual mechanism μ is characterized by a rate

J_{IIF}^{μ} (t)

⁠. Furthermore, for purposes of model development, one can define the hypothetical cumulative probability distribution corresponding to a single mechanism μ (representing the probability of IIF that would result if all other mechanisms were inhibited):

P_{IIF}^{μ} (t) = 1 - exp (- \int_{0}^{t} J_{IIF}^{μ} (τ) d τ)

(4)

Combining Eqs. –, one obtains the expected result for the composite probability distribution resulting from multiple independent mechanisms:

P_{IIF} (t) = 1 - \prod_{μ} (1 - P_{IIF}^{μ})

(5)

For example, when there are only two IIF mechanisms under consideration (e.g., SCN and VCN), Eq. (5) takes the familiar form $P_{IIF} = P_{IIF}^{SCN} + (1 - P_{IIF}^{SCN}) P_{IIF}^{VCN}$ [2].

It is clear that if a parametric model of $J_{IIF}^{μ} (t)$ is available for each mechanism of IIF, then the combination of Eqs. (2) and (3) can be fit to experimental data that have been transformed using Eq. (1). The resulting set of best-fit parameters can then be used with Eq. (4) to analyze or compare the IIF time (or temperature) distributions for the various mechanisms μ. However, as noted above, a parametric model cannot always be postulated a priori for every active IIF mechanism. In this case, as is evident from Eq. (5), it is in general not possible to solve for the individual probability distributions $P_{IIF}^{μ} (t)$ from experimental data that have been reduced to a cumulative probability distribution P_IIF(t) using the conventional method described by Eq. (1). At best, one is sometimes able to directly measure the IIF kinetics associated with individual mechanisms by physically or chemically suppressing all competing mechanisms (e.g., using fast cooling rates, high solute concentrations, or restricted temperature ranges) [1,2]. It is possible to use Eqs. (1) and (5) to solve for at most one unknown probability distribution, if the IIF kinetics of all other mechanisms are known [3].

In this work, we concern ourselves with cases in which experimental conditions make possible the identification of the mechanism μ associated with each individual IIF event (e.g., categorizing IIF events by appearance [5–8]). In such cases, the experimental data can be represented by a marked point process, in which the points are individual IIF events and the marks represent the mechanism of the corresponding IIF event. The cumulative number of IIF events due to a given mechanism μ is given by

N_{IIF}^{μ} (t) = \sum_{k} u (t - t_{k}^{μ})

(6)

where u represents the Heaviside step function;

t_{k}^{μ}

⁠, the time of the kth IIF event associated with the mechanism μ; the sum is taken over all values of k for which

t_{k}^{μ} \leq t

⁠. In the cryobiology literature, experimental data of this nature have sometimes been presented using a transformation analogous to Eq. , as follows [5,6]:

P_{μ} (t) = \frac{N_{IIF}^{μ} (t)}{N}

(7)

However, it is important to realize that the cumulative probability distribution P_μ(t) as defined by Eq. (7) is only equivalent to $P_{IIF}^{μ} (t)$ as defined by Eq. (4) in the trivial case when μ is the only active mechanism of IIF. In the presence of multiple IIF mechanisms, however, a rigorous estimation of $P_{IIF}^{μ} (t)$ requires careful accounting not only of the incidence of IIF events of type μ, but also of the changes in N_u(t) due to depletion of the pool of unfrozen cells via competing mechanism of IIF. Although Eq. (4) can be used to estimate the characteristic rate $J_{IIF}^{μ} (t)$ if the probability distribution $P_{IIF}^{μ} (t)$ is known, neither parametric nor nonparametric approaches can be used to estimate $J_{IIF}^{μ} (t)$ if Eq. (7) has been used for experimental data reduction. Thus, a new approach to analysis of experimental data comprising distinct manifestations of IIF is proposed below.

New Method for Analysis of Intracellular Ice Formation Kinetics

The rate

J_{IIF}^{μ} (t)

can be determined directly (i.e., without the need for a parametric model) by estimating the time-derivative (i.e., slope) of the following quantity:

n_{μ} (t) \equiv \int_{0}^{t} J_{IIF}^{μ} (τ) d τ

(8)

Similarly, if n_μ(t) is known, then the distribution of IIF times (or temperatures) associated with the mechanism μ follows from Eq. :

P_{IIF}^{μ} (t) = 1 - exp (- n_{μ} (t))

(9)

The quantity n_μ is the integrated intensity of the Poisson process associated with IIF events of type μ; it can be interpreted as the population mean of the cumulative number of μ-events per cell, in a hypothetical experiment in which frozen cells are replaced by unfrozen cells immediately after each IIF event (of any mechanism). The goal of the analysis below will be to determine an estimate of n_μ(t) from experimental IIF data, which then makes possible evaluation of $J_{IIF}^{μ} (t)$ and $P_{IIF}^{μ} (t)$ in a rigorous manner.

In a population of cells, IIF events with a given mark μ constitute a nonhomogeneous Poisson process for which the average arrival rate is given by

\frac{d N_{IIF}^{μ}}{d t} = J_{IIF}^{μ} (t) N_{u} (t)

(10)

Substituting the derivative of Eq. into the left-hand side of Eq. and re-arranging, one obtains

J_{IIF}^{μ} (t) = \frac{1}{N_{u} (t)} \sum_{k} δ (t - t_{k}^{μ})

(11)

where δ represents the Dirac delta function (which is the derivative of the Heaviside step function). Equation can now be substituted into Eq. to determine n_μ. If the function

N_{u}^{- 1} (t)

were continuous, then evaluation of the integral in Eq. would be straightforward, using the well-known sifting property of the Dirac delta function:

n_{μ} (t) = \sum_{k} \int_{0}^{t} \frac{1}{N_{u} (τ)} δ (τ - t_{k}^{μ}) d τ = \sum_{k} \frac{1}{N_{u} (t_{k}^{μ})}

(12)

However, because the function N_u(t) is discontinuous at the time points corresponding to each IIF event, its inverse

N_{u}^{- 1} (t)

must also be discontinuous at the event time points. Therefore, Eq. is ambiguous, and should not be used to evaluate Eq. . Instead, we have derived a modified version of the sifting property for discontinuous functions (see Eq. in the Appendix), allowing us to evaluate the integral in Eq. as follows:

n_{μ} (t) = \sum_{k} \frac{1}{2} (\frac{1}{N_{u} (t_{k}^{μ -})} + \frac{1}{N_{u} (t_{k}^{μ +})})

(13)

where $t_{k}^{μ -}$ and $t_{k}^{μ +}$ represent time points just before and just after the μ-type IIF event at $t_{k}^{μ}$ ⁠, respectively; the sum is taken over all values of k for which $t_{k}^{μ} \leq t$ ⁠. Thus, Eq. (13) represents the proposed new algorithm for quantifying IIF kinetics associated with a mechanism of interest: instead of computing a running sum of the number of IIF events of type μ (yielding the count $N_{IIF}^{μ}$ ⁠, which is of limited usefulness), one should maintain a running sum of the quantity n_μ. Each time a μ-event is observed, the number of unfrozen cells remaining (N_u) just before and just after the μ-event is to be recorded; the average of the inverse of these two values should then be computed and added to the quantity n_μ. The value thus evaluated can be substituted into Eq. (9) to yield a rigorous estimate of $P_{IIF}^{μ} (t)$ ⁠, or numerically differentiated to obtain a rigorous estimate of $J_{IIF}^{μ} (t)$ ⁠. The method embodied in Eq. (13) is a modified form of the Nelson-Aalen estimator [9,10], an alternative to the Kaplan–Meier estimator for survival analysis of life time data [11].

Equivalence of Conventional and New Methods for the Single-Mechanism Case

When there is only a single active mechanism of IIF (alternatively, if no attempt is made to categorize or otherwise assign distinct marks to the observed events in the IIF point process), the proposed method represented by Eq. must reduce to the conventional method. Noting that

P_{IIF} (t) = P_{IIF}^{μ} (t)

holds when there is only a single IIF mechanism, Eqs. and can be combined to evaluate the quantity n_μ using experimental data that have been analyzed using conventional methods:

n_{μ} (t) = \ln (\frac{N}{N_{u} (t)})

(14)

Because it is not immediately obvious that Eqs. (13) and (14) must yield equivalent results when μ is the sole mechanism of IIF, a proof of this assertion is presented below.

In the single-mechanism case, the number of unfrozen cells is

N_{u} (t) = N - N_{IIF}^{μ} (t)

⁠, and therefore it follows from Eq. that the right and left limits of N_u(t) at the IIF times

t_{k}^{μ}

are given by the following expression:

N_{u} (t_{k}^{μ +}) = N_{u} (t_{k + 1}^{μ -}) = N - k

(15)

Accordingly, Eq. can be simplified as follows:

n_{μ} (t) = \sum_{k = 1}^{N_{IIF}^{μ} (t)} \frac{1}{2} (\frac{1}{N - k + 1} + \frac{1}{N - k}) = \frac{1}{2 N} + \sum_{k = 2}^{N_{IIF}^{μ} (t)} (\frac{1}{N - k + 1}) + \frac{1}{2 (N - N_{IIF}^{μ} (t))}

(16)

Recognizing the summation above as a partial harmonic series, Eq. can be expressed in terms of harmonic numbers:

n_{μ} (t) = H (N - 1) - H (N_{u} (t)) + \frac{1}{2 N} + \frac{1}{2 N_{u} (t)}

(17)

where H(k) represents the kth harmonic number. For large k, the harmonic number can be approximated using an Euler–Maclaurin summation [12]:

H (k) = \ln k + γ + \frac{1}{2 k} - \frac{1}{12 k^{2}} + O {k^{- 4}}

(18)

where γ is the Euler–Mascheroni constant. Thus, as long as N_u(t) remains reasonably large, Eq. can be approximated by

n_{μ} (t) = \ln (N - 1) - \ln (N_{u} (t)) + \frac{1}{2 N} + \frac{1}{2 (N - 1)} - \frac{1}{12 {(N - 1)}^{2}} + \frac{1}{12 N_{u}^{2} (t)} + O {N_{u}^{- 4}}

(19)

Replacing all terms containing (N–1) by corresponding Taylor series approximations to order N⁻², Eq. becomes

n_{μ} (t) = \ln (\frac{N}{N_{u} (t)}) + \frac{1}{12 N_{u}^{2} (t)} - \frac{1}{12 N^{2}} + O {N_{u}^{- 4}} + O {N^{- 3}}

(20)

Comparison of Eqs. (14) and (20) indicates that for sufficiently large values of N and N_u, the new and conventional techniques for quantifying IIF kinetics yield equivalent results. For example, if N_u ≥ 3 and N ≥ 5, then the discrepancy between the two estimates of n_μ is expected to be of the order of 1% or less.

Comparison of Conventional and New Methods in a Simulated Experiment

To demonstrate the advantages of the new method for analyzing IIF data, a simulation of an experiment to measure IIF kinetics was generated for a hypothetical system with multiple independent IIF mechanisms, and the resulting data were then analyzed using both the conventional approach and the new technique proposed here. Although the new data analysis algorithm can accommodate any number of competing IIF mechanisms, the hypothetical experiment simulated in our example included only two distinct IIF mechanisms (labeled A and B, respectively), with characteristic rates given by $J_{IIF}^{A} = 1$ and $J_{IIF}^{B} = t$ (where a dimensionless time variable has been used, for simplicity). The theoretical cumulative probability distribution of IIF times for a system with J_IIF(t) = (1 + t) was determined from Eqs. (2) and (3), and is shown in Fig. 1. The kinetics of IIF were simulated in a population of N = 100 cells (a typical sample size for cryomicroscopy experiments), using a Monte Carlo algorithm previously described [4], and the arrival times of each A event and B event in the simulated experiment are shown in Fig. 1 (along with the cumulative IIF probability corresponding to the arrival time of each simulated IIF event, as estimated using Eq. (1)).

Fig. 1

View large Download slide

Cumulative incidence of IIF in a simulated experiment with 100 cells, as a function of dimensionless time. Individual IIF events were assumed to be distinguishable as associated with either hypothetical IIF mechanism A (open symbols) or hypothetical IIF mechanism B (closed symbols). The solid line shows the expected cumulative probability function for an infinite cell population, as theoretically predicted by Eqs. (2) and (3). See text for details.

Because the IIF events in the simulated experiment are marked (i.e., mechanisms A and B are experimentally distinguishable), it is possible to estimate the mechanism-specific kinetics of IIF. The theoretically expected cumulative probability functions for the two hypothetical IIF mechanisms were calculated using Eq. (4), and are shown in Fig. 2 along with estimated probabilities for mechanisms A and B, respectively. Data from the in silico experiment were first analyzed using the conventional approach, i.e., computing the cumulative incidence of A events and B events naïvely, using Eq. (7). As shown in Fig. 2(a), the resulting estimates of the probability of IIF settle at values less than 100%, because the conventional approach fails to account for the effect of competition between the two IIF mechanisms. In contrast, the new method of data analysis yielded more accurate estimates of the actual IIF kinetics. Figure 2(b) shows the cumulative probability distributions determined by Eq. (9), with n_μ(t) calculated from the simulated IIF data using Eq. (13). It is evident that calculation of IIF kinetics using the proposed new technique yields results that closely match the expected theoretical probability distributions, whereas the conventional approach yields results that are inaccurate except at the earliest stages of transformation (i.e., while the number of unfrozen cells remains sufficiently large that the confounding effects of competition can be neglected).

Fig. 2

View large Download slide

Estimated kinetics of IIF associated with mechanism A only (open symbols) or mechanism B only (closed symbols), compared to the theoretically expected results as predicted using Eq. (4) for mechanism A (solid line) and mechanism B (dashed line): (a) mechanism-specific IIF kinetics as determined with the conventional approach, using Eq. (7); (b) mechanism-specific IIF kinetics as determined with the proposed algorithm, using Eqs. (9) and (13).

It is also illustrative to compare calculations of the mean time of IIF using the two techniques. The conventional method (i.e., simply averaging the observed IIF time points associated with each mechanism) results in a dimensionless mean IIF time of 0.58 for mechanism A, and 1.01 for mechanism B. In comparison, when calculating means of the IIF probability distributions that have been determined using the new algorithm, the mean IIF time is estimated to be 0.77 and 1.44 for mechanisms A and B, respectively. Thus, in this example, the discrepancy between the two methods is of the order of ∼25% or larger. In general, the mean IIF times determined using the conventional approach will be lower than the corresponding mean IIF times determined using the new technique. This is a result of depletion of the pool of unfrozen cells due to competition between IIF mechanisms, which prematurely terminates the ability to observe IIF events in the experimental cell population (thus biasing the mean IIF time toward earlier events). The new technique for estimating the mean IIF time corrects for this effect by assigning greater weights to IIF events that occur later on in the freezing process.

The new method is not immune to inaccuracies that result from small sample sizes (although it is superior to the conventional approach in this regard). If IIF data were available for an infinite period of observation, then the expected mean dimensionless IIF time would be 1.00 for mechanism A and (π/2)^1/2 = 1.25 for mechanism B. For the example shown in Figs. 1 and 2, which represents a population of 100 cells, the magnitudes of the relative error in the estimated mean IIF time were 23% and 15% for mechanisms A and B, respectively, using the new algorithm (for the conventional algorithm, the corresponding errors were 42% and 19%, respectively). However, it is evident from Fig. 2(b) that the error in the estimated mean time of A-events is primarily due to the fact that IIF caused by mechanism A was only observable until t = 2, after which the pool of unfrozen cells was depleted by the competing mechanism B. The mean IIF time in the range t ≤ 2 can be shown to have a theoretical value of 0.70 for mechanism A; thus, the mean time of A-events occurring within this range was in fact estimated with a relative error of only 10% in the present example.

Application of New Method to Analysis of Experimental Intracellular Ice Formation Kinetics in Micropatterned Adherent Cells

In a 2009 study using novel high-speed video cryomicroscopy techniques, we reported evidence of previously unknown IIF pathways in adherent cells; for example, in adherent bovine pulmonary artery endothelial cells that had been allowed to spread into thin circular disks on a micropatterned substrate, the majority of IIF events observed during rapid cooling (130 °C/min) were found to originate at initiation sites along the periphery of the circular disk shape (Figs. 3(a)–3(c)), as opposed to sites inside the circle (Figs. 3(d)–3(f)) [7]. In Fig. 4, we present the IIF kinetics from our 2009 study in the form of a marked point process, in which each IIF event has been categorized by the location of its point of origin. As can be seen, the IIF events that had a peripheral-zone initiation site are interspersed with IIF events that originated from sites in the interior of the circular cross section of the spread cell (furthermore, there was also a small number of IIF events for which the location of the initiation site could not be determined [7]). Because these distinct IIF processes occurred concurrently, competing for the same pool of unfrozen cells, conventional methods of data analysis such as Eq. (7) cannot be used to estimate the kinetics associated with the individual IIF mechanisms.

Fig. 3

View large Download slide

Enhanced images from high-speed video recordings of IIF in two adherent endothelial cells that had been micropatterned in the form of thin circular disks [7]; dashed circles outline the perimeter of the spread cell (scale bar: 5 μm). The initiation sites of the IIF events (arrowheads) were frequently located at the periphery of the circle (a–c), and less frequently in the interior of the circle (d–f); the location of each IIF initiation site was determined by inspection, using frame-by-frame reverse playback to trace the solidification front back to its point of origin. To make possible visualization of the transparent intracellular ice solidification fronts in these still images, image processing has been applied to enhance changes in pixel values relative to the video frames just prior to IIF, using an algorithm previously described [7].

Fig. 4

View large Download slide

Cumulative incidence of IIF in adherent endothelial cells micropatterned in the form of thin circular disks. The data have been represented as a marked point process by classifying IIF events based on the observed location of the site of IIF initiation: the initiation site of each IIF event was categorized as being located either at the periphery of the circle (open symbols), or in the interior of the circle (closed symbols), or obscured by extracellular ice (crosses) [7]. The cumulative probability of IIF at each event temperature was calculated from experimental data using Eq. (1).

Therefore, we have used the formulas proposed in Eqs. (9) and (13) to analyze the marked IIF events shown in Fig. 4. Thus, we were able to estimate the cumulative probability of IIF caused by the mechanism responsible for intracellular crystal growth originating from the periphery of the cell disk, as well as the cumulative probability of IIF caused by the mechanism associated with intracellular crystallization originating from points in the disk interior (Fig. 5(a)). Based on the mechanism-specific cumulative probability distributions in Fig. 5(a), we can estimate that the median temperature of IIF for peripheral-zone initiation was approximately 20 K higher than the median IIF temperature for interior-zone initiation. Additional quantitative information about the kinetics of these two IIF mechanisms was obtained from the integrated IIF rates (mean cumulative number of IIF events per cell, n_μ) as determined directly from Eq. (13) and shown in Fig. 5(b); because Fig. 5(b) depicts n_μ(T), the slope of each curve has a magnitude that equals the product of the constant cooling rate and the temperature-dependent rate of IIF. For both mechanisms, the slope of the n_μ(T) curve was initially insignificant in magnitude, until the temperature decreased to approximately −15 °C; in this vicinity, there appeared to be a break in the kinetics as evidenced by the transition to a steeper slope (i.e., a faster rate of IIF). Because the slope of each curve remained relatively constant in the range −40 °C ≤ T ≤ 15 °C, a linear regression was performed in this temperature range to quantify the magnitude of the slope for each mechanism (yielding R² = 0.99 and R² = 0.93 for peripheral and interior IIF initiation, respectively). Based on the best-fit slopes, the rate of IIF for peripheral-zone initiation was 0.143±0.002 events/cell/s, while the rate of IIF for interior-zone initiation was only 0.035±0.003 events/cell/s; clearly, the peripheral-zone IIF initiation mechanism was dominant, with a rate quadruple that of the interior-zone IIF initiation mechanism. The best-fit lines for both mechanisms intercepted the temperature axis at approximately the same point (−16.5±0.3 °C for peripheral IIF initiation and −16.6±2.3 °C for interior IIF initiation), suggesting that both IIF mechanisms undergo a significant increase in activity in this temperature regime. The apparent deviations from linear (constant rate) kinetics at temperatures lower than −40 °C in Fig. 5(b) can in part be explained by the decreasing number of unfrozen cells (0 ≤ N_u ≤ 12) remaining in the population, which can adversely affect the estimated kinetics as shown in Eq. (20) above.

Fig. 5

View large Download slide

Estimated kinetics of IIF originating at initiation sites located either at the periphery (open symbols) or in the interior (closed symbols) of adherent endothelial cells micropatterned in the form of thin circular disks. (a) Cumulative probability distribution of mechanism-specific IIF as determined using Eqs.(9) and (13). (b) Mean cumulative number of mechanism-specific events per cell (n_μ), as calculated using Eq. (13); also shown are linear regressions to the kinetic data within the temperature range −40 °C ≤ T ≤ −15 °C for peripheral-zone IIF initiation (solid line) and interior-zone IIF initiation (dashed line).

Summary

When experimental conditions allow for discrimination between distinct mechanisms of IIF, conventional techniques for analysis of kinetic data cannot be used unless parametric mathematical models are available for every possible IIF mechanism, or individual IIF mechanisms can be observed in isolation (by experimentally suppressing all competing pathways for IIF). To address this deficiency, this work proposes a new algorithm for analysis of IIF kinetics, which is applicable to systems in which any number of competing mechanisms of IIF is concurrently active. The proposed method makes possible nonparametric estimation of IIF rates and probability distributions associated with each observed mechanism of IIF, without requiring experimental interventions to inhibit competing IIF mechanisms or the development of parametric theoretical models of the IIF mechanisms for curve fitting. In a hypothetical example problem with two independent mechanisms of IIF, the conventional technique for data analysis was shown to be inadequate; the probability distributions were unrealistic, and estimates of the mean times of IIF were inaccurate. In contrast, the proposed new method yielded reasonably accurate estimates of the probability distribution shapes as well as the means of the IIF times associated with each mechanism. Even though the two data analysis algorithms use dissimilar sets of mathematical operations to compute the probability of IIF, the two methods were shown to yield equivalent results for cases in which only a single mechanism of IIF is active, as would be expected. The main assumptions of the new nonparametric technique are that IIF is a Poisson process, that cell populations comprise cells with identical properties, and that the total number of cells in the population is reasonably large (ideally, on the order ∼10² or larger). These postulates are already accepted as valid in the cryobiology literature, because they are identical to the assumptions on which the conventional approach to analysis of IIF kinetics is based. By applying the proposed algorithm to real cryomicroscopy data, we were able to independently estimate the rates of IIF initiation for two competing mechanisms of IIF observed in adherent cell constructs, demonstrating the utility of the approach. In conclusion, the proposed new method represents an improved technique for estimating IIF kinetics when there are multiple distinguishable IIF mechanisms. Although the conventional method works well when all IIF events are due to a single mechanism, it should not be used when there is experimental evidence for competition among multiple mechanisms of IIF.

Acknowledgment

All endothelial cell experiments and cryomicroscopy video analyses were performed by Shannon L. Stott at the Georgia Institute of Technology.

Funding Data

Villanova University College of Engineering (Grant No. 421230; Funder ID: 10.13039/100007226).
Directorate for Engineering, National Science Foundation (Grant No. CBET-0541530; Funder ID: 10.13039/100000084).

Nomenclature

f =: discontinuous function
g =: continuous function
H =: harmonic number
J_IIF =: rate of intracellular ice formation
$J_{IIF}^{A}$ =: rate of intracellular ice formation by mechanism A
$J_{IIF}^{B}$ =: rate of intracellular ice formation by mechanism B
$J_{IIF}^{μ}$ =: rate of intracellular ice formation by mechanism μ
k =: index variable
N =: initial number of unfrozen cells
n_μ =: integrated rate of intracellular ice formation by mechanism μ
N_u =: number of unfrozen cells
$N_{IIF}^{μ}$ =: cumulative number of intracellular ice formation events by mechanism μ
P_IIF =: cumulative probability of intracellular ice formation
P_μ =: naïve estimate of cumulative probability of intracellular ice formation by mechanism μ
$P_{IIF}^{SCN}$ =: cumulative probability of intracellular ice formation by surface-catalyzed nucleation
$P_{IIF}^{VCN}$ =: cumulative probability of intracellular ice formation by volume-catalyzed nucleation
$P_{IIF}^{μ}$ =: cumulative probability of intracellular ice formation by mechanism μ
t =: time
$t_{k}^{μ}$ =: time of kth intracellular ice formation event by mechanism μ
u =: Heaviside step function

Greek Symbols

γ =: Euler–Mascheroni constant
δ =: Dirac delta function
Δf =: magnitude of discontinuity in function f(t)
τ =: integration variable

Subscripts or Superscripts

A =: hypothetical mechanism A
B =: hypothetical mechanism B
IIF =: intracellular ice formation
k =: index enumerating intracellular ice formation events
SCN =: surface-catalyzed nucleation
u =: unfrozen cells
VCN =: volume-catalyzed nucleation
μ =: index enumerating mechanisms of intracellular ice formation

Appendix

A function f(t) that is discontinuous at the origin can be decomposed into a linear combination of a Heaviside step function and a function g(t) that is continuous at the origin:

f (t) = g (t) + Δ f \cdot u (t)

(A1)

where Δf represents the magnitude of the discontinuity in f(t), which is the difference between the right and left limits of f(t):

Δ f = f (0^{+}) - f (0^{-})

(A2)

To derive a modified sifting property for the Dirac delta function, the integral of the product of the delta function with the discontinuous function f(t) must be evaluated:

\int_{- \infty}^{+ \infty} f (t) δ (t) d t = \int_{- \infty}^{+ \infty} g (t) δ (t) d t + Δ f \int_{- \infty}^{+ \infty} u (t) δ (t) d t

(A3)

where the right-hand side was obtained by substitution of Eq. . Because g(t) is continuous at the origin, the sifting property can be used to evaluate the first integral on the right-hand side of Eq. , yielding a value equal to g(0) = f(0⁻). To evaluate the second integral, one can take advantage of the fact that the Dirac delta function is the derivative of the Heaviside step function such that

\int_{- \infty}^{+ \infty} u (t) δ (t) d t = \int_{- \infty}^{+ \infty} u (t) \frac{d u}{d t} d t = \int_{0}^{1} udu = {\frac{u^{2}}{2} |}_{0}^{1} = \frac{1}{2}

(A4)

Substitution of these results into Eq. yields a modified sifting property that is appropriate when f(t) is discontinuous:

\int_{- \infty}^{+ \infty} f (t) δ (t) d t = f (0^{-}) + \frac{Δ f}{2} = \frac{f (0^{-}) + f (0^{+})}{2}

(A5)

Equation represents a generalized form of the sifting property that is applicable also when f(t) is continuous, because by definition, continuity requires equality of the left and right limits of f(t). The integral of the product of f(t) with an impulse train (i.e., a sequence of delta functions at multiple time-points t_k) can be determined from Eq. by taking advantage of the time invariance of the sifting property:

\int_{- \infty}^{+ \infty} f (t) \sum_{k} δ (t - t_{k}) d t = \sum_{k} \int_{- \infty}^{+ \infty} f (t) δ (t - t_{k}) d t = \sum_{k} \frac{f (t_{k}^{-}) + f (t_{k}^{+})}{2}

(A6)

Equation (A6) was used to evaluate the integral in Eq. (8), with the integrand specified by Eq. (11).

References

1.

Pitt

,

R. E.

,

Myers

,

S. P.

,

Lin

,

T.-T.

, and

Steponkus

,

P. L.

,

1991

, “

Subfreezing Volumetric Behavior and Stochastic Modeling of Intracellular Ice Formation in Drosophila melanogaster Embryos

,”

Cryobiology

,

28

(

1

), pp.

72

–

86

.10.1016/0011-2240(91)90009-D

Google Scholar

Crossref

PubMed

2.

Toner

,

M.

,

Cravalho

,

E. G.

, and

Karel

,

M.

,

1990

, “

Thermodynamics and Kinetics of Intracellular Ice Formation During Freezing of Biological Cells

,”

J. Appl. Phys.

,

67

(

3

), pp.

1582

–

1592

.10.1063/1.345670

Google Scholar

Crossref

3.

Karlsson

,

J. O. M.

,

Eroglu

,

A.

,

Toth

,

T. L.

,

Cravalho

,

E. G.

, and

Toner

,

M.

,

1996

, “

Fertilization and Development of Mouse Oocytes Cryopreserved Using a Theoretically Optimized Protocol

,”

Hum. Reprod.

,

11

(

6

), pp.

1296

–

1305

.10.1093/oxfordjournals.humrep.a019375

Google Scholar

Crossref

PubMed

4.

Irimia

,

D.

, and

Karlsson

,

J. O. M.

,

2005

, “

Kinetics of Intracellular Ice Formation in One-Dimensional Arrays of Interacting Biological Cells

,”

Biophys. J.

,

88

(

1

), pp.

647

–

660

.10.1529/biophysj.104.048355

Google Scholar

Crossref

PubMed

5.

Scheiwe

,

M. W.

, and

Körber

,

C.

,

1987

, “

Quantitative Cryomicroscopic Analysis of Intracellular Freezing of Granulocytes Without Cryoadditive

,”

Cryobiology

,

24

(

5

), pp.

473

–

483

.10.1016/0011-2240(87)90051-4

Google Scholar

Crossref

PubMed

6.

Berrada

,

M. S.

, and

Bischof

,

J. C.

,

2001

, “

Evaluation of Freezing Effects on Human Microvascular-Endothelial Cells (HMEC

),”

Cryo Lett.

,

22

(

6

), pp.

353

–

366

.https://www.jstor.org/stable/2958850

7.

Stott

,

S. L.

, and

Karlsson

,

J. O. M.

,

2009

, “

Visualization of Intracellular Ice Formation Using High-Speed Video Cryomicroscopy

,”

Cryobiology

,

58

(

1

), pp.

84

–

95

.10.1016/j.cryobiol.2008.11.003

Google Scholar

Crossref

PubMed

8.

Higgins

,

A. Z.

, and

Karlsson

,

J. O. M.

,

2013

, “

Effects of Intercellular Junction Protein Expression on Intracellular Ice Formation in Mouse Insulinoma Cells

,”

Biophys. J.

,

105

(

9

), pp.

2006

–

2015

.10.1016/j.bpj.2013.09.028

Google Scholar

Crossref

PubMed

9.

Nelson

,

W.

,

1972

, “

Theory and Application of Hazard Plotting for Censored Survival Data

,”

Technometrics

,

14

(

4

), pp.

945

–

966

.10.1080/00401706.1972.10488991

Google Scholar

Crossref

10.

Aalen

,

O.

,

1978

, “

Nonparametric Inference for a Family of Counting Processes

,”

Ann. Stat.

,

6

(

4

), pp.

701

–

726

.

Google Scholar

Crossref

11.

Colosimo

,

E. A.

,

Ferreira

,

F. F.

,

Oliveira

,

M. D.

, and

Sousa

,

C. B.

,

2002

, “

Empirical Comparisons Between Kaplan-Meier and Nelson-Aalen Survival Function Estimators

,”

J. Stat. Comput. Simul.

,

72

(

4

), pp.

299

–

308

.10.1080/00949650212847

Google Scholar

Crossref

12.

Havil

,

J.

,

2003

,

Gamma: Exploring Euler's Constant

,

Princeton University Press

,

Princeton, NJ

, Chap.

10

.

A Nonparametric Method for Estimating Rates of Intracellular Ice Formation Due to Independent Competing Mechanisms

Abstract

Introduction

Conventional Method for Analysis of Intracellular Ice Formation Kinetics

New Method for Analysis of Intracellular Ice Formation Kinetics

Equivalence of Conventional and New Methods for the Single-Mechanism Case

Comparison of Conventional and New Methods in a Simulated Experiment

Application of New Method to Analysis of Experimental Intracellular Ice Formation Kinetics in Micropatterned Adherent Cells

Summary

Acknowledgment

Funding Data

Nomenclature

Greek Symbols

Subscripts or Superscripts

Appendix

References

Get Email Alerts

Cited By

Journals

Conference Proceedings

eBooks

Resources

Opportunities

A Nonparametric Method for Estimating Rates of Intracellular Ice Formation Due to Independent Competing Mechanisms

Abstract

Introduction

Conventional Method for Analysis of Intracellular Ice Formation Kinetics

New Method for Analysis of Intracellular Ice Formation Kinetics

Equivalence of Conventional and New Methods for the Single-Mechanism Case

Comparison of Conventional and New Methods in a Simulated Experiment

Application of New Method to Analysis of Experimental Intracellular Ice Formation Kinetics in Micropatterned Adherent Cells

Summary

Acknowledgment

Funding Data

Nomenclature

Greek Symbols

Subscripts or Superscripts

Appendix

References

Get Email Alerts

Cited By

Related Articles

Related Proceedings Papers

Related Chapters

Journals

Conference Proceedings

eBooks

Resources

Opportunities

This Feature Is Available To Subscribers Only