- Email: [email protected]
Prediction of local cooling rates and cell survival during the freezing of a cylindrical specimen
CRYOBIOLOGY
25, 67-82 (1988)
Prediction
of Local Cooling Rates and Cell Survival Freezing of a Cylindrical Specimen’
LINDA J. HAYES*V
KENNETH R. DILLER,t HAK S. LEES
during
the
HSIEN-JAN CHANG,* AND
*Texas Institute for Computational Mechanics, Department of Aerospace Engineering and Engineering Mechanics; and ?Biomedical Engineering Program, Department of Mechanical Engineering, The University of Texas at Austin, Austin, Texas 78712; and SGeneral Electric Company, Schenectady, New York 12309
A finite element numerical model was implemented to simulate the freezing process of an aqueous salt solution in a cylindrical container. Local cooling rates within the container were computed for several defined cooling protocols applied at the boundary. Characteristic cell survival signatures were used to predict the associated local survival rates throughout the system. These calculations show that there are two definite time domains during a typical freezing process: (1) while the surface temperature is changing and (2) after the surface temperature reaches a constant storage value. The calculations also show significant spatial variations in the local cooling rates within the container and considerable local deviation from the volumetric average survival for various simulated freezing protocols. 0 1988 Academic Press, Inc.
measured post-thaw survival is plotted as a function of cooling rate during freezing. The characteristic shape of the signature includes a maximum value at an intermediate cooling rate, above and below which the survival decreases monotonically. To describe an actual cryopreservation procedure, the freezing process is typically characterized in terms of a single value of the cooling rate, and this single value is used to correlate and interpret cell survival as indicated in Fig. 1. While this approach lends itself to convenience and simplicity, it will be physically accurate only under very special conditions of freezing in which a nearly constant localized cooling rate can be effected throughout the specimen in both time and position. Two factors are primary deterrents to achieving an idealized constant cooling process throughout a Received May 8, 1986; accepted September 14, 1987. specimen. First, it is often not possible to i This research was sponsored in part by the National Science Foundation Grant No. CBT-8713600 manipulate the environmental temperature so as to create the transient boundary conand a University Research Institute Grant at the Uniditions required to maintain constant surversity of Texas at Austin. * To whom correspondence should be addressed at face cooling, and second, because all real Aerospace Engineering/Engineering Mechanics, systems of interest have finite dimensions, WRW 305, The University of Texas, Austin, TX thermal transport properties, and latent 78712. The long-term frozen preservation of living tissue for transplantation has been developed as a clinically efficacious technique during the past several decades. The survival of living cells after cryopreservation has been clearly established to be a strong function of the temperature history of the specimen (8). Numerous thermal parameters are required to fully specify a freeze-thaw protocol. These include the cooling and warming rates, nucleation and storage temperatures, and intermediate isothermal holding periods (2). Extensive research has identified the cooling rate as a crucial factor in determining the ultimate fate of a frozen cell. This phenomenon is illustrated by the set of simulated survival signatures presented in Fig. 1, where the
67 001l-2240/88 $3.00 Copyright 0 1988 by Academic Press, Inc. All rights of reproduction in any form reserved.
68
HAYES ET AL. A Modified
Experimental
Data
Normal
Dist.,
S.D.=l
‘Urnin
o Normal
Dist.,
S.D.=4
OC/min
q
Cooling
Rate
(MED)
(“Clmin)
FIG. 1. Three survival signatures, all normalized to 100% survival at S”C/min. The two signatures with normal distributions are hypothetical and are introduced to illustrate varying degrees of sensitivity of cell survival to thermal history. The third signature is adapted from data reported in the literature [12] and scaled to the same optimal cooling rate.
tion of cell survival throughout the container for a specific freezing protocol. The local survival rates can then be volumetrically averaged to give an overall survival rate for the bulk sample. This is a concept that has also been explored by prior investigators (1). We present here a numerical model by which the local, transient temperature field which develops within a freezing system may be simulated. Cooling rates are calculated locally and are a function of boundary conditions, of system transport properties, and of geometry, and the corresponding local values of cell survival are determined from an appropriate survival signature. A volumetrically averaged cell survival rate, a, is obtained for the bulk sample by integration of these local data. The effect which manipulation of the exterior cooling rate and holding temperature has upon net cell survival is illustrated for several example survival signatures. MATHEMATICAL
heats, constant surface conditions are not propagated uniformly to interior locations. These limitations can result in a nonuniform distribution of temperature histories and cooling rates throughout the spatial domain. Such local variations in thermal history have been clearly demonstrated experimentally (7, 12) and appreciated in the context of their influence on the response to freezing in finite dimensioned systems (10) for many years. Nonetheless, the question remains as to how to characterize and quantify the effect of these inhomogeneities. The present study addresses this problem. Given the improbability of achieving a cooling process which is uniform in either space or time, it is necessary to define a local value for the cooling rate which is averaged in time during the solidification process (6). This local value of cooling rate can be used with a survival signature in order to determine the localized distribu-
MODEL
Heat transfer problems involving phase change are very difftcult to model both analytically and numerically. The position of the freezing front is dependent upon the unknown temperature field; the locations of the phase boundaries are not known a priori and must be carefully determined during the numerical procedure. This situation, coupled with the possibility of rapidly changing thermal properties such as specific heat and thermal conductivities near the phase front, makes the numerical modeling more difficult than in the case of simple heat diffusion. The volumetric latent heat of fusion, A (kJ/m3), which is liberated during freezing depends on the mass fraction of water undergoing solidification and may be expressed as A = LWf,
ill
where L is the mass specific latent heat of fusion for water (333 kJ/kg) and wf is the
PREDICTION
OF LOCAL
mass fraction of water available for freezing of the medium (kg/m3). In systems of impure chemical composition, such as tissue, the phase change process is distributed over a range of temperatures rather than occurring at a specific temperature. The numerical instability problems inherent in modeling phase change processes increase as the energy release associated with latent heat occurs over a smaller temperature range. This effect may be interpreted in terms of the ratio of latent heat and sensible heat effects within the system as defined by a modified Stefan number, Ste’ (4) Ste’ = PC, 7
,
PI
where AT, is the temperature range across which latent heat is released, C, is specific heat, and p is density. The modified Stefan number describes the ratio of sensible to latent energy changes during solidification. Increased values of A or decreased values of AT, result in lower values of Ste’, indicating that the energy release due to latent heat dominates that of sensible heat and that numerical instabilities will be more pronounced. The freezing process can be described by the time-dependent partial differential equation
CELL
69
SURVIVAL
tainer and T, is the transient surface temperature. In this study the range of temperatures for liberation of latent heat was taken to extend from -2.6”C to -21.3”C, (AT, = 18.7”C). The nonlinear pattern of latent heat release as a function of temperature shown in Fig. 2 was derived by simple application of the lever rule to the phase diagram in Fig. 3 for a binary aqueous solution with an initial composition of five times isotonic electrolyte concentration. This derivation is valid if one assumes that the system passes continuously through a sequence of quasistationary states during the freezing process. The starting composition was chosen since it illustrates effectively the balance between latent and sensible effects over a range of temperatures, although it may not be optimal for cryobiological applications. The influence of different initial compositions, including isotonic, on thermal history has been addressed in a companion paper (5). NUMERICAL
METHOD
The numerical scheme used in the present investigation is the modified apparent heat capacity method (5, 11) and was implemented in a transient two-dimensional finite element code (3). Although the model is general, in this study the code was applied to a simple one-dimensional cylinit (pC,T + A) = V - (kVZ’) in fI, 131 drically symmetric system, which was nonetheless adequate to illustrate the spawhere t is time, T is temperature (“C), k is tial variations in thermal histories and thermal conductivity (kW/m-“C), and R is cooling rates and the resulting effects on the region of interest. The initial thermal local cell survival. By assuming the extent condition in the system is specified as of latent heat release to be proportional to the remaining liquid fraction during the T(x,O) = T,(x) for x in a, 141 freezing process, the energy state of the where x is the position vector and To is an system can be expressed as a function of initial temperature field, which is taken to temperature. The apparent heat capacity is be uniform. The temperature on the outer defined as surface of the container is specified as c* = per + g T(x,t) = T,(t) for x in aR, [A where dKI is the outer surface of the con-
and the finite element method is used to
70
HAYES ET AL. loo-
zo-
0. 0
-10
-15
Temperature
-20
PC)
FIG. 2. Nonlinear latent heat release pattern during freezing for an initial liquid solution composition of 5 x isotonic. Zero supercooling is assumed.
solve the equation
tions. For elements whose nodal point temperatures are in the phase transition zone, c* g - V - (kVT) = 0 in fl [71 i.e., - 2.6 to -21.3”C, the apparent heat capacity c* was evaluated directly at the with prescribed boundary and initial condi- Gauss integration points based on the local
MOLE
FRACTION
OF
SALT
FIG. 3. Phase diagram for a binary, aqueous solution of sodium chloride. Various prefreezing liquid phase compositions are indicated.
71
PREDICTION OF LOCAL CELL SURVIVAL
freezing processes were effected by lowering the surface temperature at a constant rate from 1OS”C to holding temperatures of -80, - 120, or -240°C. Exterior surface cooling rates of - 2, -5, -7, - 10 and - 20”C/min were considered. The time averaged cooling rate at each nodal position in the finite element grid was calculated between the temperatures -2XC, which is 0.2”C below the equilibrium phase change temperature for the initial system composition, and -3O.O”C. This average cooling rate was then applied to obtain a single value of the survival rate, SR(r), at each radial node position in the grid for each of the three cell survival signatures in Fig. 1. It was assumed that time variations in cooling rate had no effect on the cell survival signature and that temperatures above -2XC and below - 3O.O”C did not influence survival. For each of the cases studied a volumetric average survival rate, SR, was obtained for the entire specimen from the local survival rates.
TABLE 1 Material Properties Used in the Simulation J+OpefiY
System
A TO PC,
Values
Units J/cm’ “C J/cm3-“C
328.0
10.5 Liquid
4.08
Solution Solid Container
k
1.88 2.08
J/cm3-‘T
Liquid
0.0059
W/cm-T
Solid
0.0222 0.00328
W/cm-T
Solution Container
temperatures. The model produces values of temperature as functions of position and time, which are then used to calculate local temperature histories and cooling rates. A long cylindrical container was modeled with the objective of calculating both local and volumetrically averaged cell survival rates for specific simulation protocols. The container was assumed to be made of high density polyethylene with a 0.1~cm wall thickness and an internal radius of 4.23 cm and was filled with a 5 x isotonic salt-water solution. End effects were neglected. Although specifying particular dimensions for the system limits the general applicability of the analysis, it does afford an explicit example of the phenomenon we wish to illustrate. Values of the thermal properties of the model system are given in Table 1. All property values were assumed constant for each phase throughout the simulation and were varied linearly during phase change according to the fraction of the solution solidified. The finite element grid used to represent this system is shown in Fig. 4. A variety of SOLUTION
SR=-
1
R
SR(r)27~rdr, dW2 i 0 where R = 4.23 cm.
Three different simulated survival vs cooling rate curves, shown in Fig. 1, were used in the model: (1) a normal distribution centered at a cooling rate of - S”C/min with a standard deviation (SD) of l.O”C/min, (2) a normal distribution centered at - S”C/min with a standard deviation of 4.0”C/min, and (3) a curve scaled from experimental data (9) to indicate maximum survival at a cooling rate of - SUmin. The original survival signature was obtained for freezing
_ CONTAINER
iilllill I#, : s f 0 ( 0.u 0.14
” / 1.03
*I” 1.c.s
33 I 2.51
ml
39 A0
I 43 I 8.23
“I’ ST 4.094 !,
NODE # RADIUS
(cm)
FIG. 4. Finite element grid for the model cylindrical container. Representative dimensions and the corresponding node numbers are indicated. The container wall is denoted by the crosshatching.
72
HAYES
mouse marrow stem cells suspended in Tyrode’s solution with a 1.25 M concentration of glycerol and was cooled to - 196°C. However, the specific conditions under which this signature was derived are not crucial to the present study, since these data are being used as a generic example to illustrate how local variations in cell survival may compare with and affect bulk measurements. All three of the curves were normalized to a maximum survival of 100% at - 5°C to simplify interpretation of the results. The container surface temperature is lowered at a constant rate to one of three final values, - 80, - 120, or -240°C and then held constant for the duration of the simulation. If when the final boundary temperature is reached while there are interior nodes within the system still at a temperature within the range for which the cooling rate is calculated for determination of cell survival, in general, there will be a significant effect on the computed thermal histories and survivals at these nodes. For example, at the slowest cooling rate of - 2”C/ min, almost all of the specimen is below - 30°C before the surface is held at - 80°C. However, at the fastest cooling rate of -2O”C/min, approximately 10% of the volume is below -30°C when the surface reaches - 80. RESULTS
Figures Sa-5c show the transient thermal histories at selected interior points for cooling rates of -5, - 10, and -20°C with a holding temperature of - 80°C. The node numbers given in the legends correspond to the radial positions identified on the grid in Fig. 4. On each graph the band of temperatures between which average cooling rates were calculated (- 2.8 to -3O.O”C) is indicated by two horizontal lines. A vertical line indicates when the surface reaches the holding temperature. The shapes of the thermal history curves in
ET AL.
Fig. 5 for the three different cooling rates at the surface are qualitatively similar to each other. At locations near the surface the temperature-time curve is strongly influenced by the cooling rate at the boundary; at interior points the history is dominated by diffusion through and capacitance of the outlying solution so that variations in boundary conditions are propagated less strongly. For example, consider any nodal point near the center of the container (nodes l-25). The time required for interior points to reach the temperature at which phase change is initiated is roughly inversely proportional to the boundary cooling rate in each of the live cases. However, once local phase change begins (at T = - 2.8”C), the temperature histories for a given point are very similar for all three boundary cooling rates. To some degree this effect can be attributed to the decrease in the unfrozen volume near the center of the container and to the increase in thermal conductivity and volume of the frozen phase as the protocol proceeds. However, it will be shown that the exterior holding temperature is the dominant factor in determining interior cooling rates. Figure 6a shows an overlay of the -5 and -2O”Umin cooling rates with the surface temperature held at -80°C. These thermal histories have been displaced so that the cooling rates at the center of the container lie on top of each other. A vertical line is used to indicate the time after which both surfaces have reached the holding temperature of -80°C. These results show that when the boundary arrives at the holding temperature, the thermal histories for the region wherein survival is determined are approximately the same, independent of the surface cooling rate. This effect can be attributed to the outer frozen phase having a relatively high thermal diffusivity which quickly propagates boundary conditions into the interior of the container. Figure 6b shows a similar
73
PREDICTION OF LOCAL CELL SURVIVAL
” 0.0 -1s.o e5 --3o.o z ; -45.0 ii +5-a-60.0
- .._...57 5, -----. .-.--.-.. 43 39 - * *33 ._..... 5 --+--. ,, *
NODE
,
-75.0
FIG. 5. (a) Temperature histories at selected points for a surface cooling rate of - S”C/min with a holding temperature of -80°C. (b) Temperature histories at selected points for a surface cooling rate of - lO”C/min with a holding temperature of - 80°C. (c) Temperature histories at selected points for a surface cooling rate of - 20Wmin with a holding temperature of - 80°C.
overlay for the -5 and -20°C cooling rates with a holding temperature of - 120°C. Again, after both protocols reach the holding temperature - 12O”C, the interior temperatures are approximately the same. Figure 6c shows a corresponding overlay with - 5 and - 20°C cooling rates when the surface temperature is lowered to -240°C. In this case -240°C is reached at 3000 set for the -5Wmin protocol, at which time the entire specimen is below the -30°C limit where cooling rate calculations are
terminated. It is clear in this overlay that the thermal histories in the interior of the container are substantially different for the two cooling rates. This condition occurred because the outside holding temperature was not reached until after the entire interior had solidified. Figure 6d shows an overlay of the -5 and -20°C cooling rates with the surface held at -80°C. Here overlays have been positioned so that the -80°C holding temperature is reached at the same time, t,. In this graph curves corresponding to the
74
HAYES
ET AL.
-90.0 1 Time ---5’C --------20T
Cooling rate * -24OT Cooling rate@ -24O’C
c 15.0 0.0 -15.0 0 e s -30.0 : B j -45.0 -60.0 -75.0 -90.0
FIG. 6. (a) - S”C/min and - 20Wmin cooling rates with a holding temperature of - 80°C. The time scales are identical to those in Fig. 5, but have been offset so that interior histories are aligned. (b) The -5Wmin and -2WWmin cooling rates with a holding temperature of - 120°C. The time scales are identical to those in Fig. 5, but have been offset so that interior histories are aligned. (c) The -WI min and - 20”C/min cooling rates with a holding temperature of - 240°C. The time scales are identical to those in Fig. 5, but have been offset so that interior histories are aligned. (d) The -5Wmin and - 20YYmin cooling rates with a holding temperature of - 80°C. The time scales are identical to those in Fig. 5, but have been offset so that holding temperature is reached at time to.
same node are connected by arrows with the node number. Figures 6a and 6d are identical except for the different time offsets. Figure 6d clearly illustrates that thermal histories do not correspond immediately upon reaching the boundary holding temperature. This set of data can be considered as the solution to the solidification process at a common final temperature of -80°C with two different sets of initial
temperature distributions, one corresponding to the - 5°C cooling rate and one corresponding to the -20°C cooling rate. The value of the minimum surface temperature is a dominant factor in determining the thermal histories within the region which has not solidified when the final holding temperature is reached. Figure 6d illustrates that the cooling rate has a direct effeet on the internal temperature distribu-
PREDICTION
OF LOCAL
tion that exists within the system at the completion of the surface cooling process. The constitutive law for this material indicates that after solidification the thermal capacitance drops and the thermal conductivity increases, both by a factor of approximately four, thereby giving thermal diffusivity in the solid phase that is about 16 times higher than that in the unfrozen portion. This phenomenon causes boundary temperature alterations to propagate very quickly through the frozen shell of the system into the interior region where latent heat is being removed. Figures 7a-7c show the local averaged cooling rate as a function of position in the container for each of the five boundary cooling protocols at - 2, - 5, - 7, - 10 and - 20°C and for each of the holding temperatures of - 80, - 120, and - 240°C. Cooling rates for the - 80 and - 120°C holding temperatures display a very distinctive behavior. Near the exterior of the container the curves exhibit a local maximum cooling rate (minimum in curve) that exceeds the value at nodes further interior until the centerline is approached. All of the curves unify at the center and have very high cooling rates independent of the surface cooling rate. Therefore the curve can be described in terms of three distinct portions, including a decrease from the boundary cooling rate to a point of minimum value, followed by an increase in cooling rate at locations more interior leading to a merging with other curves, and finally, as the midline is approached, a unified curve for all protocols having the common minimum surface temperature. For each cooling rate, the outermost point on the curve corresponds to the location in the container which is just beginning phase change when the boundary reaches the holding temperature. This state of the system also identifies the point where the curves merge, which corresponds to the location in the container which has just reached -30°C. For instance, the two
CELL
SURVIVAL
75
points in Fig. 7a (which are indicated by the symbols A and a) correspond to the two points in Fig. 5a marked with the same symbol. The symbols A and $r identify the range of temperatures across which the cooling rate is measured, representing the minimum and maximum values, respectively. The cooling rates for a holding temperature of -240°C do not exhibit this behavior because all of the sample has cooled below -30°C when the minimum surface temperature is reached. Figure 7d is a composite graph illustrating the influence of the minimum final surface temperature of a protocol on the internal distribution of cooling rates within the container. The figure shows that at the point in time when the surface temperature begins to be held constant there is a distinct alteration in the interior thermal pattern that is nearly independent of the prior cooling rate. The more rapid the surface cooling rate, the greater is the portion of the system affected by this phenomenon. After this transition state is achieved, the cooling rates at interior locations are a strong function of the constant surface temperature . Tables 2a-2c show cooling rate and survival rate data for each of the holding temperatures: - 80, - 120, and - 240°C. These tables give the values of the local cooling rates at selected nodes in the container for each of the surface cooling rates simulated: -2, -5, -7, - 10, and -20°C. These data were also shown graphically in Figs. 7a-7c, and they illustrate that there was little variation in cooling rates at positions toward the interior of the container. The interior nodes (9 to 25) appear to be essentially insensitive to boundary cooling rates. However, the data in Tables 2b and 2c and in Figs. 7b and 7c show that the cooling rates in the interior are a strong function of the minimum surface temperature. When the surface is lowered to -8O”C, this temperature is reached before the interior has a chance to cool substantially; thus, for all
a
0.
-5.
O-
-10.
o-
-15.
O-
.E m -20.
O-
F ._ E ; e aa : Ly
7j 8
O-
-25.
i
-
BCR
-_-
__-_-_
-----_
BCRz-7
-.-.__.___._.__.-. -
-30.1 0
-35.1
D-
c-2
BCRS_-J
BCflC_,
-
I
-
0
BCR=--20
I
I
I
Container
f 1
I
5 9
I
I
17
I
25
33
II
39
43
51
Node
g
Radius
Container
Ill
II
I
I
159
17
25
33
“ES 66;;
z
z
3 82 &tin’
II
39 43
51
Node#
g 4
Radius(cm1
# (cm)
i .O
7
Position
(cm)
r
8’
d
0.0 0.0 ,
1.0
2.0
I
1
3.0
I
4.0
I
/
i
Container i
59
17
25
33
39 43
51
Node#
2k
3$
39 4$
5’1’ Node
FIG. 7. Local cooling rate as a function of interior radial position for -2, -5, - 7, - 10 and -2OWmin cooling rates with a holding temperature of (a) -8O”C, (b) - 12O”C,and (c) -240°C. (d) Composite map of local cooling rates for cooling to surface holding temperatures of - 80, - 120, and - 240°C. 76
#
77
PREDICTION OF LOCAL CELL SURVIVAL TABLE 2a Cooling and Survival Rates for Holding Temperature - 80°C Local cooling rates at selected points (Wmin)
Boundary CR Coordinate
Node
0.48 09 1.58 25 2.905 39 3.23 43 4.03 51 4.33 57 Volumetric average CR Range of 80% cooling rates
-2.0
-5.0
-7.0
- 10.0
- 20.0
-8.43 -3.39 -2.40 -2.18 - 1.62 -2.00
-8.48 -3.53 -4.31 -4.84 -3.66 -5.00
-8.63 -3.51 -4.41 -5.36 - 4.95 -7.00
-8.36 -3.54 -4.45 -5.45 - 6.75 - 10.00
- 8.33 -3.54 -4.41 -5.47 - 11.81 - 20.00
-2.69
-4.23
-4.80
-5.33
-6.30
- 1.4 to -3.1
-3.1 to -4.7
-3.4 to -5.1
-3.4 to -7.0
-3.4 to -9.8
MED SD = 1 SD = 4
Volumetric average ??i?(%) based on local CR 57.4 87.5 90.1 11.8 61.2 65.7 80.2 95.4 96.0
MED SD = 1 SD = 4
65.3 7.08 84.2
Survival Rate z (%) based only on CR 90.4 97.5 73.6 96.7 98.4 99.8
85.5 46.0 92.4
76.5 38.3 79.4
97.9 93.6 99.7
90.0
45.5 94.8
TABLE 2b Cooling and Survival Rates for Holding Temperature - 120°C Local cooling rate at selected points (Wmin)
Boundary CR Coordinate
Node
0.48 09 1.58 25 2.905 39 3.23 43 4.03 51 4.33 57 Volumetric average CR Range of 80% cooling rates MED SD = 1 SD = 4 MED SD = 1 SD = 4
-2.0
-5.0
-7.0
- 10.0
- 20.0
- 10.21 -3.70 -2.40 -2.18 - 1.62 - 2.00
- 14.63 -7.37 -5.62 - 5.20 -3.66 -5.00
- 14.67 -7.59 -7.66 - 7.07 - 4.95 -7.00
- 14.92 -7.64 -9.77 -9.66 -6.74 - 10.00
- 15.12 -7.63 - 10.89 - 13.70 - 12.59 -20.00
-2.85
-5.94
-7.28
-8.65
-11.54
1.4 to -3.3
-3.1 to -7.1
-4.1 to -8.0
-5.4 to -9.6
- 7.6 to - 14.6
Volumetric average m (%) based on local ??i? 57.7 86.0 78.7 13.6 52.5 28.0 80.1 90.4 83.1 Survival Rate z (%) based only on CR 70.3 94.1 77.8 9.1 65.8 7.36 86.1 97.4 84.6
67.4 6.8 67.6
50.0 0.66 31.2
64.7 0.16 65.8
40.4 0 26.3
HAYES ET AL. TABLE 2c Cooling and Survival Rates for Holding Temperature - 240°C Boundary CR Coordinate
Node
0.48 09 1.58 25 2.905 39 3.23 43 4.03 51 4.33 57 Volumetric average CR Range of 80% cooling rates MED SD = 1 SD = 4 MED SD = 1 SD = 4
Local cooling rate at selected points (“Urnin) -2.0
-5.0
-7.0
- 10.0
- 20.0
- 10.21 -3.70 -2.40 -2.18 - 1.62 -2.00
- 19.45 - 7.89 -5.62 -5.20 -3.66 -5.00
- 23.50 - 10.33 - 7.66 - 7.07 -4.95 -7.00
- 30.73 - 13.74 - 10.24 -9.66 -6.74 - 10.00
- 36.97 -21.96 - 18.37 - 17.41 - 12.48 - 20.00
-2.85
-6.27
- 8.34
-11.11
- 18.73
- 1.4 to -3.3
-3.1 to -7.2
-5.4 to - 12.6
-9.3 to -9.8
-4.1 to -9.5
Volumetric average a (%) based on local CR 57.7 84.4 73.5 13.6 52.0 27.2 80.1 87.7 73.2 Survival Rate E (%) based only on CR 70.3 90.4 66.6 9.7 47.2 0.57 86.1 95.0 70.4
surface cooling rates the temperatures at internal locations are about the same. When the surface is lowered to - 12O”C, approximately one-half of the container has cooled significantly prior to reaching the holding temperature, and this condition produces a different pattern in the interior thermal histories. This effect is illustrated most dramatically for a minimum boundary temperature of -240°C. Here all of the internal nodes have cooled below -30°C so that the effective thermal history of the entire container is defined before the holding temperature is reached; as a result, the cooling rates in the system are directly dependent on the boundary thermal history. Tables 2a-2c also list the volumetrically averaged cooling rate, CR, for each one of these trials. For lower holding temperatures, CR is more strongly a function of the exterior cooling rate. This table also includes values for the median 80% of cooling rates achieved within the specimen. For example Table 2a shows that for the container surface cooled at - lO”C/min
55.0 6.1 45.7
33.3 0 4.8
43.1 0 31.4
32.2 0 0.32
to a holding temperature of - 80°C 80% of the sample experiences cooling rates between - 3.4 and -7”C, 10% experiences rates lower than -3.4”C, and 10% experiences rates higher than -7°C. This parameter is thus a measure of the breadth of distribution of thermal histories within the specimen. The volumetric average survival, SR, for the entire container, as calculated by Eq. [6], is presented in Table 2. SR is the overall survival rate that would be identified after thawing and remixing of the specimen. Alternatively, CR was calculated and applied as a look-up table in conjunction with Fig. 1 to determine a single value approximation for the survival rate, @. This procedure is similar to monitoring the temperature history at a single location and using that one cooling rate value to predict survival for the entire container from one of the curves in Fig. 1. Table 2 shows that the single value E differs by as much as 100% from the integrated SR obtained from local thermal history data distributed
PREDICTION OF LOCAL CELL SURVIVAL
over which cooling occurs, and the shape of the specific survival signature must be considered in predicting total survival for cells frozen in bulk containers. Figure 8 shows the percentage of the total volume which was cooled at or below any specified cooling rate. This data can be used to investigate the volumetric distribution of thermal histories within the container. The curves provide an indication of how broad a spectrum of cooling rates occurs for a given boundary cooling rate and holding temperature. This figure illustrates that with more rapid cooling at the boundary, the internal thermal histories become spread over a much broader range, and it is less appropriate to represent the cooling rate within the cylinder by a single value. Local cell survival was calculated as a function of position for each of the three survival signatures for the - 120°C holding temperature. These data are presented in Figs. 9a-9c. As expected the survival is highest at points in the container where the
throughout the system. A survival rate based on a single value for the cooling rate does not reflect the range of local variations in thermal history which may occur within a container. For example, with a holding temperature of -80°C and a - lO”C/min cooling rate, the volumetric average cooling rate was -5.33”C/min, which is very near the optimum value for survival defined in Fig. 1. However, considerably lower local survivals are expected when the actual distribution of cooling rates between - 3.4 and - 7Wmin is considered. The predicted & valu3 which are obtained from the single CR value, agree moderately well for broad survival signatures when a relatively narrow range of local cooling r&s is experienced. However, a and SR may vary substantially if the survival signature is not broad and if the container experiences a signifcaniange of local cooling rates. In general, SR will not have a good correlation with the actual survival rate, %. Thus, the average cooling rate, the range
-
Holding
Temperature
-24OOC
---
Holding
Temperature
-180°C
.--.---..
Holding
Temperature
-gO°C
Cooling
79
Rate
PCImin)
FIG. 8. Volumetric distribution of cooling rates within the container as a function of the surface cooling rate and minimum temperature of exposure. The labels on individual curves represent the cooling rate in (“Clmin).
80
HAYES ET AL.
Node Radius
c (cm)
Porkion
1 5 e
17
25
(cm)
33
39 43
FIG. 9. (a) Survival rates determined for the modified experimental data from Fig. 1 (holding temperature - 120°C). (b) Survival rates using a normalized curve with a standard deviation of l”C/min (holding temperature - 120°C). (c) Survival rates using a normalized curve with a standard deviation of 4Wmin (holding temperature - 120°C). (d) Distribution of local survival rates within the container as a function of surface cooling rate and minimum temperature. Survival was obtained from local cooling rates in terms of the survival signature in Fig. 1 derived from modified experimental data.
cooling rate is close to -5Wmin. However, the survival patterns for a given cooling protocol differ greatly for each of the three signatures. One can use the average cooling rate and the range of cooling rates given in Table 2 to predict the variation in survival rates for each of the three signatures in Fig. 1. The predicted
survival distributions from Fig. 1 for the modified experimental data and for the normal distribution with a standard deviation of 4YYmin are shown in Figs. 9a and 9c. For these two cases, the local survival rates are relatively insensitive to position for most cooling rates effected at the boundary. In contrast, the predicted local
PREDICTION OF LOCAL CELL SURVIVAL
survival for a narrower signature corresponding to a normal distribution with a standard deviation of l”C/min (Fig. 9b) varies greatly. For this signature, the survival decreases rapidly as cooling rates are displaced from the optimal value of -SC/ min, and the resulting survival rates are extremely sensitive to the local cooling rate. These graphs illustrate the effect which nonuniform cooling of bulk specimens can have upon total survival. Figure 9d is a composite illustration of survival as a function of holding temperature for - 5, - 10, and - 20Wmin surface cooling protocols, corroborating that survival is a strong combined function of the surface cooling rate and the minimum temperature of exposure.
81
single probe to monitor cooling rates to try to ensure optimal survival of a specific cell type frozen in a finite dimensioned container. The fastest cooling occurs at the center of the container, so that a probe placed in this position would represent a poor choice for purposes of monitoring, even though the geometric symmetry may appear appealing. Quantitative results described in this paper are dependent on the assumed container size and geometry, solution composition, and cell survival signatures. Application to other systems would be dependent upon implementing values of these parameters consistent with the specific system being considered. However, the general trends illustrated by the present study serve to identify the interaction beCONCLUSIONS tween physical and physiological pheThe finite element freezing model pre- nomena in determining the response of dicts that thermal boundary conditions are cells to freezing in finite dimension not propagated uniformly from the exterior systems. into the interior of a cylindrical system, ACKNOWLEDGMENT thus resulting in a nonuniform distribution The authors wish to thank one of the reviewers for of temperature histories and cooling rates throughout the spatial domain. The cooling helpful comments regarding physical interpretation of rates at interior points in the system are de- the analysis. pendent upon both the surface cooling rate REFERENCES and the exterior holding temperature. As 1. Chmiel, H., Stroh, N., and Walitza, E. Vefahrenthe minimum surface holding temperature stechnische aspecte der gefrierkonsevierung decreases, a larger portion of the interior is von blubestandteilen. Biomed. Tech. 25, 32-57 (1980). sensitive to the boundary cooling condi2. Farrant, J. Genera1 observations on cell preservations. This effect occurs because local tion. In “Low Temperature Preservation in cooling rates are calculated from the Medicine and Biology” (M. J. Ashwood-Smith thermal history between the range of -2.8 and J. Farrant, Eds.), pp. l-18. Univ. Park and - 3O”C, and for lower holding temperaPress, Baltimore, 1980. 3. Hayes, L. J. “A users guide to PARAB-A twotures more of the container is cooled to this dimensional linear time-dependent finite elerange of temperatures before the isoment code. TICOM Report 80-10.” Univ. thermal boundary conditions are invoked. Texas, Austin, 1980. As the exterior cooling rate is increased, 4. Hayes, L. J., and Diller, K. R. Implementation of the interior volumetric distribution of phase change in numerical models of heat cooling rates in the interior broadens. transfer. Trans. ASME, J. Energy Resour. Technol. 105, 431-435 (1983). A single average value of cooling rate 5. Hayes, L. J., Diller, K. R., and Chang, H. J. A within a container is in general of limited robust numerical method for latent heat release use in describing survival rates during during phase change. In “Numerical Methods cryopreservation. Neither is it possible to in Heat Transfer” (J. L. S. Chen and K. Vafai, Eds.), pp. 63-70. ASME, New York, 1986. predict a best location for positioning a
82
HAYES ET AL.
6. Hayes, L. J., Diller, K. R., and Lee, H. S. On the definition of an average cooling rate during cell freezing. Cryo-Letr. 5, 97- 110 (1984). 7. Luyet, B. J. An attempt at a systematic analysis of the notion of freezing rates and at an evaluation of the main contributory factors. Cryobiology 2, 198-205 (1966). 8. Mazur, l? Cryobiology, the freezing of biological systems. Science 168, 939-949 (1970). 9. Mazur, P., Leibo, S. P., Farrant, J., Chu, E. H. Y., Hanna, M. G., and Smith, L. H. Interactions of cooling rate, warming rate and protective additive on the survival of frozen
mammalian cells. In “The Frozen Cell” (G. E. W. Wolstenholme and M. O’Connor, Eds.), pp. 69-84. Churchill, London, 1970. 10. Meryman, H. T. The interpretation of freezing rates in biological materials. Cryobiology 2, 165-170 (1966). 11. Morgan, K., Lewis, R. W., and Zienkiewicz, 0. C. An improved algorithm for heat conduction problems with phase change. Znt. .I. Numer. Merhods Eng. 12, 1191-1195 (1978). 12. Rinfret, A. P. Thermal history. Cryobiology 2, 171-180 (1966).
25, 67-82 (1988)
Prediction
of Local Cooling Rates and Cell Survival Freezing of a Cylindrical Specimen’
LINDA J. HAYES*V
KENNETH R. DILLER,t HAK S. LEES
during
the
HSIEN-JAN CHANG,* AND
*Texas Institute for Computational Mechanics, Department of Aerospace Engineering and Engineering Mechanics; and ?Biomedical Engineering Program, Department of Mechanical Engineering, The University of Texas at Austin, Austin, Texas 78712; and SGeneral Electric Company, Schenectady, New York 12309
A finite element numerical model was implemented to simulate the freezing process of an aqueous salt solution in a cylindrical container. Local cooling rates within the container were computed for several defined cooling protocols applied at the boundary. Characteristic cell survival signatures were used to predict the associated local survival rates throughout the system. These calculations show that there are two definite time domains during a typical freezing process: (1) while the surface temperature is changing and (2) after the surface temperature reaches a constant storage value. The calculations also show significant spatial variations in the local cooling rates within the container and considerable local deviation from the volumetric average survival for various simulated freezing protocols. 0 1988 Academic Press, Inc.
measured post-thaw survival is plotted as a function of cooling rate during freezing. The characteristic shape of the signature includes a maximum value at an intermediate cooling rate, above and below which the survival decreases monotonically. To describe an actual cryopreservation procedure, the freezing process is typically characterized in terms of a single value of the cooling rate, and this single value is used to correlate and interpret cell survival as indicated in Fig. 1. While this approach lends itself to convenience and simplicity, it will be physically accurate only under very special conditions of freezing in which a nearly constant localized cooling rate can be effected throughout the specimen in both time and position. Two factors are primary deterrents to achieving an idealized constant cooling process throughout a Received May 8, 1986; accepted September 14, 1987. specimen. First, it is often not possible to i This research was sponsored in part by the National Science Foundation Grant No. CBT-8713600 manipulate the environmental temperature so as to create the transient boundary conand a University Research Institute Grant at the Uniditions required to maintain constant surversity of Texas at Austin. * To whom correspondence should be addressed at face cooling, and second, because all real Aerospace Engineering/Engineering Mechanics, systems of interest have finite dimensions, WRW 305, The University of Texas, Austin, TX thermal transport properties, and latent 78712. The long-term frozen preservation of living tissue for transplantation has been developed as a clinically efficacious technique during the past several decades. The survival of living cells after cryopreservation has been clearly established to be a strong function of the temperature history of the specimen (8). Numerous thermal parameters are required to fully specify a freeze-thaw protocol. These include the cooling and warming rates, nucleation and storage temperatures, and intermediate isothermal holding periods (2). Extensive research has identified the cooling rate as a crucial factor in determining the ultimate fate of a frozen cell. This phenomenon is illustrated by the set of simulated survival signatures presented in Fig. 1, where the
67 001l-2240/88 $3.00 Copyright 0 1988 by Academic Press, Inc. All rights of reproduction in any form reserved.
68
HAYES ET AL. A Modified
Experimental
Data
Normal
Dist.,
S.D.=l
‘Urnin
o Normal
Dist.,
S.D.=4
OC/min
q
Cooling
Rate
(MED)
(“Clmin)
FIG. 1. Three survival signatures, all normalized to 100% survival at S”C/min. The two signatures with normal distributions are hypothetical and are introduced to illustrate varying degrees of sensitivity of cell survival to thermal history. The third signature is adapted from data reported in the literature [12] and scaled to the same optimal cooling rate.
tion of cell survival throughout the container for a specific freezing protocol. The local survival rates can then be volumetrically averaged to give an overall survival rate for the bulk sample. This is a concept that has also been explored by prior investigators (1). We present here a numerical model by which the local, transient temperature field which develops within a freezing system may be simulated. Cooling rates are calculated locally and are a function of boundary conditions, of system transport properties, and of geometry, and the corresponding local values of cell survival are determined from an appropriate survival signature. A volumetrically averaged cell survival rate, a, is obtained for the bulk sample by integration of these local data. The effect which manipulation of the exterior cooling rate and holding temperature has upon net cell survival is illustrated for several example survival signatures. MATHEMATICAL
heats, constant surface conditions are not propagated uniformly to interior locations. These limitations can result in a nonuniform distribution of temperature histories and cooling rates throughout the spatial domain. Such local variations in thermal history have been clearly demonstrated experimentally (7, 12) and appreciated in the context of their influence on the response to freezing in finite dimensioned systems (10) for many years. Nonetheless, the question remains as to how to characterize and quantify the effect of these inhomogeneities. The present study addresses this problem. Given the improbability of achieving a cooling process which is uniform in either space or time, it is necessary to define a local value for the cooling rate which is averaged in time during the solidification process (6). This local value of cooling rate can be used with a survival signature in order to determine the localized distribu-
MODEL
Heat transfer problems involving phase change are very difftcult to model both analytically and numerically. The position of the freezing front is dependent upon the unknown temperature field; the locations of the phase boundaries are not known a priori and must be carefully determined during the numerical procedure. This situation, coupled with the possibility of rapidly changing thermal properties such as specific heat and thermal conductivities near the phase front, makes the numerical modeling more difficult than in the case of simple heat diffusion. The volumetric latent heat of fusion, A (kJ/m3), which is liberated during freezing depends on the mass fraction of water undergoing solidification and may be expressed as A = LWf,
ill
where L is the mass specific latent heat of fusion for water (333 kJ/kg) and wf is the
PREDICTION
OF LOCAL
mass fraction of water available for freezing of the medium (kg/m3). In systems of impure chemical composition, such as tissue, the phase change process is distributed over a range of temperatures rather than occurring at a specific temperature. The numerical instability problems inherent in modeling phase change processes increase as the energy release associated with latent heat occurs over a smaller temperature range. This effect may be interpreted in terms of the ratio of latent heat and sensible heat effects within the system as defined by a modified Stefan number, Ste’ (4) Ste’ = PC, 7
,
PI
where AT, is the temperature range across which latent heat is released, C, is specific heat, and p is density. The modified Stefan number describes the ratio of sensible to latent energy changes during solidification. Increased values of A or decreased values of AT, result in lower values of Ste’, indicating that the energy release due to latent heat dominates that of sensible heat and that numerical instabilities will be more pronounced. The freezing process can be described by the time-dependent partial differential equation
CELL
69
SURVIVAL
tainer and T, is the transient surface temperature. In this study the range of temperatures for liberation of latent heat was taken to extend from -2.6”C to -21.3”C, (AT, = 18.7”C). The nonlinear pattern of latent heat release as a function of temperature shown in Fig. 2 was derived by simple application of the lever rule to the phase diagram in Fig. 3 for a binary aqueous solution with an initial composition of five times isotonic electrolyte concentration. This derivation is valid if one assumes that the system passes continuously through a sequence of quasistationary states during the freezing process. The starting composition was chosen since it illustrates effectively the balance between latent and sensible effects over a range of temperatures, although it may not be optimal for cryobiological applications. The influence of different initial compositions, including isotonic, on thermal history has been addressed in a companion paper (5). NUMERICAL
METHOD
The numerical scheme used in the present investigation is the modified apparent heat capacity method (5, 11) and was implemented in a transient two-dimensional finite element code (3). Although the model is general, in this study the code was applied to a simple one-dimensional cylinit (pC,T + A) = V - (kVZ’) in fI, 131 drically symmetric system, which was nonetheless adequate to illustrate the spawhere t is time, T is temperature (“C), k is tial variations in thermal histories and thermal conductivity (kW/m-“C), and R is cooling rates and the resulting effects on the region of interest. The initial thermal local cell survival. By assuming the extent condition in the system is specified as of latent heat release to be proportional to the remaining liquid fraction during the T(x,O) = T,(x) for x in a, 141 freezing process, the energy state of the where x is the position vector and To is an system can be expressed as a function of initial temperature field, which is taken to temperature. The apparent heat capacity is be uniform. The temperature on the outer defined as surface of the container is specified as c* = per + g T(x,t) = T,(t) for x in aR, [A where dKI is the outer surface of the con-
and the finite element method is used to
70
HAYES ET AL. loo-
zo-
0. 0
-10
-15
Temperature
-20
PC)
FIG. 2. Nonlinear latent heat release pattern during freezing for an initial liquid solution composition of 5 x isotonic. Zero supercooling is assumed.
solve the equation
tions. For elements whose nodal point temperatures are in the phase transition zone, c* g - V - (kVT) = 0 in fl [71 i.e., - 2.6 to -21.3”C, the apparent heat capacity c* was evaluated directly at the with prescribed boundary and initial condi- Gauss integration points based on the local
MOLE
FRACTION
OF
SALT
FIG. 3. Phase diagram for a binary, aqueous solution of sodium chloride. Various prefreezing liquid phase compositions are indicated.
71
PREDICTION OF LOCAL CELL SURVIVAL
freezing processes were effected by lowering the surface temperature at a constant rate from 1OS”C to holding temperatures of -80, - 120, or -240°C. Exterior surface cooling rates of - 2, -5, -7, - 10 and - 20”C/min were considered. The time averaged cooling rate at each nodal position in the finite element grid was calculated between the temperatures -2XC, which is 0.2”C below the equilibrium phase change temperature for the initial system composition, and -3O.O”C. This average cooling rate was then applied to obtain a single value of the survival rate, SR(r), at each radial node position in the grid for each of the three cell survival signatures in Fig. 1. It was assumed that time variations in cooling rate had no effect on the cell survival signature and that temperatures above -2XC and below - 3O.O”C did not influence survival. For each of the cases studied a volumetric average survival rate, SR, was obtained for the entire specimen from the local survival rates.
TABLE 1 Material Properties Used in the Simulation J+OpefiY
System
A TO PC,
Values
Units J/cm’ “C J/cm3-“C
328.0
10.5 Liquid
4.08
Solution Solid Container
k
1.88 2.08
J/cm3-‘T
Liquid
0.0059
W/cm-T
Solid
0.0222 0.00328
W/cm-T
Solution Container
temperatures. The model produces values of temperature as functions of position and time, which are then used to calculate local temperature histories and cooling rates. A long cylindrical container was modeled with the objective of calculating both local and volumetrically averaged cell survival rates for specific simulation protocols. The container was assumed to be made of high density polyethylene with a 0.1~cm wall thickness and an internal radius of 4.23 cm and was filled with a 5 x isotonic salt-water solution. End effects were neglected. Although specifying particular dimensions for the system limits the general applicability of the analysis, it does afford an explicit example of the phenomenon we wish to illustrate. Values of the thermal properties of the model system are given in Table 1. All property values were assumed constant for each phase throughout the simulation and were varied linearly during phase change according to the fraction of the solution solidified. The finite element grid used to represent this system is shown in Fig. 4. A variety of SOLUTION
SR=-
1
R
SR(r)27~rdr, dW2 i 0 where R = 4.23 cm.
Three different simulated survival vs cooling rate curves, shown in Fig. 1, were used in the model: (1) a normal distribution centered at a cooling rate of - S”C/min with a standard deviation (SD) of l.O”C/min, (2) a normal distribution centered at - S”C/min with a standard deviation of 4.0”C/min, and (3) a curve scaled from experimental data (9) to indicate maximum survival at a cooling rate of - SUmin. The original survival signature was obtained for freezing
_ CONTAINER
iilllill I#, : s f 0 ( 0.u 0.14
” / 1.03
*I” 1.c.s
33 I 2.51
ml
39 A0
I 43 I 8.23
“I’ ST 4.094 !,
NODE # RADIUS
(cm)
FIG. 4. Finite element grid for the model cylindrical container. Representative dimensions and the corresponding node numbers are indicated. The container wall is denoted by the crosshatching.
72
HAYES
mouse marrow stem cells suspended in Tyrode’s solution with a 1.25 M concentration of glycerol and was cooled to - 196°C. However, the specific conditions under which this signature was derived are not crucial to the present study, since these data are being used as a generic example to illustrate how local variations in cell survival may compare with and affect bulk measurements. All three of the curves were normalized to a maximum survival of 100% at - 5°C to simplify interpretation of the results. The container surface temperature is lowered at a constant rate to one of three final values, - 80, - 120, or -240°C and then held constant for the duration of the simulation. If when the final boundary temperature is reached while there are interior nodes within the system still at a temperature within the range for which the cooling rate is calculated for determination of cell survival, in general, there will be a significant effect on the computed thermal histories and survivals at these nodes. For example, at the slowest cooling rate of - 2”C/ min, almost all of the specimen is below - 30°C before the surface is held at - 80°C. However, at the fastest cooling rate of -2O”C/min, approximately 10% of the volume is below -30°C when the surface reaches - 80. RESULTS
Figures Sa-5c show the transient thermal histories at selected interior points for cooling rates of -5, - 10, and -20°C with a holding temperature of - 80°C. The node numbers given in the legends correspond to the radial positions identified on the grid in Fig. 4. On each graph the band of temperatures between which average cooling rates were calculated (- 2.8 to -3O.O”C) is indicated by two horizontal lines. A vertical line indicates when the surface reaches the holding temperature. The shapes of the thermal history curves in
ET AL.
Fig. 5 for the three different cooling rates at the surface are qualitatively similar to each other. At locations near the surface the temperature-time curve is strongly influenced by the cooling rate at the boundary; at interior points the history is dominated by diffusion through and capacitance of the outlying solution so that variations in boundary conditions are propagated less strongly. For example, consider any nodal point near the center of the container (nodes l-25). The time required for interior points to reach the temperature at which phase change is initiated is roughly inversely proportional to the boundary cooling rate in each of the live cases. However, once local phase change begins (at T = - 2.8”C), the temperature histories for a given point are very similar for all three boundary cooling rates. To some degree this effect can be attributed to the decrease in the unfrozen volume near the center of the container and to the increase in thermal conductivity and volume of the frozen phase as the protocol proceeds. However, it will be shown that the exterior holding temperature is the dominant factor in determining interior cooling rates. Figure 6a shows an overlay of the -5 and -2O”Umin cooling rates with the surface temperature held at -80°C. These thermal histories have been displaced so that the cooling rates at the center of the container lie on top of each other. A vertical line is used to indicate the time after which both surfaces have reached the holding temperature of -80°C. These results show that when the boundary arrives at the holding temperature, the thermal histories for the region wherein survival is determined are approximately the same, independent of the surface cooling rate. This effect can be attributed to the outer frozen phase having a relatively high thermal diffusivity which quickly propagates boundary conditions into the interior of the container. Figure 6b shows a similar
73
PREDICTION OF LOCAL CELL SURVIVAL
” 0.0 -1s.o e5 --3o.o z ; -45.0 ii +5-a-60.0
- .._...57 5, -----. .-.--.-.. 43 39 - * *33 ._..... 5 --+--. ,, *
NODE
,
-75.0
FIG. 5. (a) Temperature histories at selected points for a surface cooling rate of - S”C/min with a holding temperature of -80°C. (b) Temperature histories at selected points for a surface cooling rate of - lO”C/min with a holding temperature of - 80°C. (c) Temperature histories at selected points for a surface cooling rate of - 20Wmin with a holding temperature of - 80°C.
overlay for the -5 and -20°C cooling rates with a holding temperature of - 120°C. Again, after both protocols reach the holding temperature - 12O”C, the interior temperatures are approximately the same. Figure 6c shows a corresponding overlay with - 5 and - 20°C cooling rates when the surface temperature is lowered to -240°C. In this case -240°C is reached at 3000 set for the -5Wmin protocol, at which time the entire specimen is below the -30°C limit where cooling rate calculations are
terminated. It is clear in this overlay that the thermal histories in the interior of the container are substantially different for the two cooling rates. This condition occurred because the outside holding temperature was not reached until after the entire interior had solidified. Figure 6d shows an overlay of the -5 and -20°C cooling rates with the surface held at -80°C. Here overlays have been positioned so that the -80°C holding temperature is reached at the same time, t,. In this graph curves corresponding to the
74
HAYES
ET AL.
-90.0 1 Time ---5’C --------20T
Cooling rate * -24OT Cooling rate@ -24O’C
c 15.0 0.0 -15.0 0 e s -30.0 : B j -45.0 -60.0 -75.0 -90.0
FIG. 6. (a) - S”C/min and - 20Wmin cooling rates with a holding temperature of - 80°C. The time scales are identical to those in Fig. 5, but have been offset so that interior histories are aligned. (b) The -5Wmin and -2WWmin cooling rates with a holding temperature of - 120°C. The time scales are identical to those in Fig. 5, but have been offset so that interior histories are aligned. (c) The -WI min and - 20”C/min cooling rates with a holding temperature of - 240°C. The time scales are identical to those in Fig. 5, but have been offset so that interior histories are aligned. (d) The -5Wmin and - 20YYmin cooling rates with a holding temperature of - 80°C. The time scales are identical to those in Fig. 5, but have been offset so that holding temperature is reached at time to.
same node are connected by arrows with the node number. Figures 6a and 6d are identical except for the different time offsets. Figure 6d clearly illustrates that thermal histories do not correspond immediately upon reaching the boundary holding temperature. This set of data can be considered as the solution to the solidification process at a common final temperature of -80°C with two different sets of initial
temperature distributions, one corresponding to the - 5°C cooling rate and one corresponding to the -20°C cooling rate. The value of the minimum surface temperature is a dominant factor in determining the thermal histories within the region which has not solidified when the final holding temperature is reached. Figure 6d illustrates that the cooling rate has a direct effeet on the internal temperature distribu-
PREDICTION
OF LOCAL
tion that exists within the system at the completion of the surface cooling process. The constitutive law for this material indicates that after solidification the thermal capacitance drops and the thermal conductivity increases, both by a factor of approximately four, thereby giving thermal diffusivity in the solid phase that is about 16 times higher than that in the unfrozen portion. This phenomenon causes boundary temperature alterations to propagate very quickly through the frozen shell of the system into the interior region where latent heat is being removed. Figures 7a-7c show the local averaged cooling rate as a function of position in the container for each of the five boundary cooling protocols at - 2, - 5, - 7, - 10 and - 20°C and for each of the holding temperatures of - 80, - 120, and - 240°C. Cooling rates for the - 80 and - 120°C holding temperatures display a very distinctive behavior. Near the exterior of the container the curves exhibit a local maximum cooling rate (minimum in curve) that exceeds the value at nodes further interior until the centerline is approached. All of the curves unify at the center and have very high cooling rates independent of the surface cooling rate. Therefore the curve can be described in terms of three distinct portions, including a decrease from the boundary cooling rate to a point of minimum value, followed by an increase in cooling rate at locations more interior leading to a merging with other curves, and finally, as the midline is approached, a unified curve for all protocols having the common minimum surface temperature. For each cooling rate, the outermost point on the curve corresponds to the location in the container which is just beginning phase change when the boundary reaches the holding temperature. This state of the system also identifies the point where the curves merge, which corresponds to the location in the container which has just reached -30°C. For instance, the two
CELL
SURVIVAL
75
points in Fig. 7a (which are indicated by the symbols A and a) correspond to the two points in Fig. 5a marked with the same symbol. The symbols A and $r identify the range of temperatures across which the cooling rate is measured, representing the minimum and maximum values, respectively. The cooling rates for a holding temperature of -240°C do not exhibit this behavior because all of the sample has cooled below -30°C when the minimum surface temperature is reached. Figure 7d is a composite graph illustrating the influence of the minimum final surface temperature of a protocol on the internal distribution of cooling rates within the container. The figure shows that at the point in time when the surface temperature begins to be held constant there is a distinct alteration in the interior thermal pattern that is nearly independent of the prior cooling rate. The more rapid the surface cooling rate, the greater is the portion of the system affected by this phenomenon. After this transition state is achieved, the cooling rates at interior locations are a strong function of the constant surface temperature . Tables 2a-2c show cooling rate and survival rate data for each of the holding temperatures: - 80, - 120, and - 240°C. These tables give the values of the local cooling rates at selected nodes in the container for each of the surface cooling rates simulated: -2, -5, -7, - 10, and -20°C. These data were also shown graphically in Figs. 7a-7c, and they illustrate that there was little variation in cooling rates at positions toward the interior of the container. The interior nodes (9 to 25) appear to be essentially insensitive to boundary cooling rates. However, the data in Tables 2b and 2c and in Figs. 7b and 7c show that the cooling rates in the interior are a strong function of the minimum surface temperature. When the surface is lowered to -8O”C, this temperature is reached before the interior has a chance to cool substantially; thus, for all
a
0.
-5.
O-
-10.
o-
-15.
O-
.E m -20.
O-
F ._ E ; e aa : Ly
7j 8
O-
-25.
i
-
BCR
-_-
__-_-_
-----_
BCRz-7
-.-.__.___._.__.-. -
-30.1 0
-35.1
D-
c-2
BCRS_-J
BCflC_,
-
I
-
0
BCR=--20
I
I
I
Container
f 1
I
5 9
I
I
17
I
25
33
II
39
43
51
Node
g
Radius
Container
Ill
II
I
I
159
17
25
33
“ES 66;;
z
z
3 82 &tin’
II
39 43
51
Node#
g 4
Radius(cm1
# (cm)
i .O
7
Position
(cm)
r
8’
d
0.0 0.0 ,
1.0
2.0
I
1
3.0
I
4.0
I
/
i
Container i
59
17
25
33
39 43
51
Node#
2k
3$
39 4$
5’1’ Node
FIG. 7. Local cooling rate as a function of interior radial position for -2, -5, - 7, - 10 and -2OWmin cooling rates with a holding temperature of (a) -8O”C, (b) - 12O”C,and (c) -240°C. (d) Composite map of local cooling rates for cooling to surface holding temperatures of - 80, - 120, and - 240°C. 76
#
77
PREDICTION OF LOCAL CELL SURVIVAL TABLE 2a Cooling and Survival Rates for Holding Temperature - 80°C Local cooling rates at selected points (Wmin)
Boundary CR Coordinate
Node
0.48 09 1.58 25 2.905 39 3.23 43 4.03 51 4.33 57 Volumetric average CR Range of 80% cooling rates
-2.0
-5.0
-7.0
- 10.0
- 20.0
-8.43 -3.39 -2.40 -2.18 - 1.62 -2.00
-8.48 -3.53 -4.31 -4.84 -3.66 -5.00
-8.63 -3.51 -4.41 -5.36 - 4.95 -7.00
-8.36 -3.54 -4.45 -5.45 - 6.75 - 10.00
- 8.33 -3.54 -4.41 -5.47 - 11.81 - 20.00
-2.69
-4.23
-4.80
-5.33
-6.30
- 1.4 to -3.1
-3.1 to -4.7
-3.4 to -5.1
-3.4 to -7.0
-3.4 to -9.8
MED SD = 1 SD = 4
Volumetric average ??i?(%) based on local CR 57.4 87.5 90.1 11.8 61.2 65.7 80.2 95.4 96.0
MED SD = 1 SD = 4
65.3 7.08 84.2
Survival Rate z (%) based only on CR 90.4 97.5 73.6 96.7 98.4 99.8
85.5 46.0 92.4
76.5 38.3 79.4
97.9 93.6 99.7
90.0
45.5 94.8
TABLE 2b Cooling and Survival Rates for Holding Temperature - 120°C Local cooling rate at selected points (Wmin)
Boundary CR Coordinate
Node
0.48 09 1.58 25 2.905 39 3.23 43 4.03 51 4.33 57 Volumetric average CR Range of 80% cooling rates MED SD = 1 SD = 4 MED SD = 1 SD = 4
-2.0
-5.0
-7.0
- 10.0
- 20.0
- 10.21 -3.70 -2.40 -2.18 - 1.62 - 2.00
- 14.63 -7.37 -5.62 - 5.20 -3.66 -5.00
- 14.67 -7.59 -7.66 - 7.07 - 4.95 -7.00
- 14.92 -7.64 -9.77 -9.66 -6.74 - 10.00
- 15.12 -7.63 - 10.89 - 13.70 - 12.59 -20.00
-2.85
-5.94
-7.28
-8.65
-11.54
1.4 to -3.3
-3.1 to -7.1
-4.1 to -8.0
-5.4 to -9.6
- 7.6 to - 14.6
Volumetric average m (%) based on local ??i? 57.7 86.0 78.7 13.6 52.5 28.0 80.1 90.4 83.1 Survival Rate z (%) based only on CR 70.3 94.1 77.8 9.1 65.8 7.36 86.1 97.4 84.6
67.4 6.8 67.6
50.0 0.66 31.2
64.7 0.16 65.8
40.4 0 26.3
HAYES ET AL. TABLE 2c Cooling and Survival Rates for Holding Temperature - 240°C Boundary CR Coordinate
Node
0.48 09 1.58 25 2.905 39 3.23 43 4.03 51 4.33 57 Volumetric average CR Range of 80% cooling rates MED SD = 1 SD = 4 MED SD = 1 SD = 4
Local cooling rate at selected points (“Urnin) -2.0
-5.0
-7.0
- 10.0
- 20.0
- 10.21 -3.70 -2.40 -2.18 - 1.62 -2.00
- 19.45 - 7.89 -5.62 -5.20 -3.66 -5.00
- 23.50 - 10.33 - 7.66 - 7.07 -4.95 -7.00
- 30.73 - 13.74 - 10.24 -9.66 -6.74 - 10.00
- 36.97 -21.96 - 18.37 - 17.41 - 12.48 - 20.00
-2.85
-6.27
- 8.34
-11.11
- 18.73
- 1.4 to -3.3
-3.1 to -7.2
-5.4 to - 12.6
-9.3 to -9.8
-4.1 to -9.5
Volumetric average a (%) based on local CR 57.7 84.4 73.5 13.6 52.0 27.2 80.1 87.7 73.2 Survival Rate E (%) based only on CR 70.3 90.4 66.6 9.7 47.2 0.57 86.1 95.0 70.4
surface cooling rates the temperatures at internal locations are about the same. When the surface is lowered to - 12O”C, approximately one-half of the container has cooled significantly prior to reaching the holding temperature, and this condition produces a different pattern in the interior thermal histories. This effect is illustrated most dramatically for a minimum boundary temperature of -240°C. Here all of the internal nodes have cooled below -30°C so that the effective thermal history of the entire container is defined before the holding temperature is reached; as a result, the cooling rates in the system are directly dependent on the boundary thermal history. Tables 2a-2c also list the volumetrically averaged cooling rate, CR, for each one of these trials. For lower holding temperatures, CR is more strongly a function of the exterior cooling rate. This table also includes values for the median 80% of cooling rates achieved within the specimen. For example Table 2a shows that for the container surface cooled at - lO”C/min
55.0 6.1 45.7
33.3 0 4.8
43.1 0 31.4
32.2 0 0.32
to a holding temperature of - 80°C 80% of the sample experiences cooling rates between - 3.4 and -7”C, 10% experiences rates lower than -3.4”C, and 10% experiences rates higher than -7°C. This parameter is thus a measure of the breadth of distribution of thermal histories within the specimen. The volumetric average survival, SR, for the entire container, as calculated by Eq. [6], is presented in Table 2. SR is the overall survival rate that would be identified after thawing and remixing of the specimen. Alternatively, CR was calculated and applied as a look-up table in conjunction with Fig. 1 to determine a single value approximation for the survival rate, @. This procedure is similar to monitoring the temperature history at a single location and using that one cooling rate value to predict survival for the entire container from one of the curves in Fig. 1. Table 2 shows that the single value E differs by as much as 100% from the integrated SR obtained from local thermal history data distributed
PREDICTION OF LOCAL CELL SURVIVAL
over which cooling occurs, and the shape of the specific survival signature must be considered in predicting total survival for cells frozen in bulk containers. Figure 8 shows the percentage of the total volume which was cooled at or below any specified cooling rate. This data can be used to investigate the volumetric distribution of thermal histories within the container. The curves provide an indication of how broad a spectrum of cooling rates occurs for a given boundary cooling rate and holding temperature. This figure illustrates that with more rapid cooling at the boundary, the internal thermal histories become spread over a much broader range, and it is less appropriate to represent the cooling rate within the cylinder by a single value. Local cell survival was calculated as a function of position for each of the three survival signatures for the - 120°C holding temperature. These data are presented in Figs. 9a-9c. As expected the survival is highest at points in the container where the
throughout the system. A survival rate based on a single value for the cooling rate does not reflect the range of local variations in thermal history which may occur within a container. For example, with a holding temperature of -80°C and a - lO”C/min cooling rate, the volumetric average cooling rate was -5.33”C/min, which is very near the optimum value for survival defined in Fig. 1. However, considerably lower local survivals are expected when the actual distribution of cooling rates between - 3.4 and - 7Wmin is considered. The predicted & valu3 which are obtained from the single CR value, agree moderately well for broad survival signatures when a relatively narrow range of local cooling r&s is experienced. However, a and SR may vary substantially if the survival signature is not broad and if the container experiences a signifcaniange of local cooling rates. In general, SR will not have a good correlation with the actual survival rate, %. Thus, the average cooling rate, the range
-
Holding
Temperature
-24OOC
---
Holding
Temperature
-180°C
.--.---..
Holding
Temperature
-gO°C
Cooling
79
Rate
PCImin)
FIG. 8. Volumetric distribution of cooling rates within the container as a function of the surface cooling rate and minimum temperature of exposure. The labels on individual curves represent the cooling rate in (“Clmin).
80
HAYES ET AL.
Node Radius
c (cm)
Porkion
1 5 e
17
25
(cm)
33
39 43
FIG. 9. (a) Survival rates determined for the modified experimental data from Fig. 1 (holding temperature - 120°C). (b) Survival rates using a normalized curve with a standard deviation of l”C/min (holding temperature - 120°C). (c) Survival rates using a normalized curve with a standard deviation of 4Wmin (holding temperature - 120°C). (d) Distribution of local survival rates within the container as a function of surface cooling rate and minimum temperature. Survival was obtained from local cooling rates in terms of the survival signature in Fig. 1 derived from modified experimental data.
cooling rate is close to -5Wmin. However, the survival patterns for a given cooling protocol differ greatly for each of the three signatures. One can use the average cooling rate and the range of cooling rates given in Table 2 to predict the variation in survival rates for each of the three signatures in Fig. 1. The predicted
survival distributions from Fig. 1 for the modified experimental data and for the normal distribution with a standard deviation of 4YYmin are shown in Figs. 9a and 9c. For these two cases, the local survival rates are relatively insensitive to position for most cooling rates effected at the boundary. In contrast, the predicted local
PREDICTION OF LOCAL CELL SURVIVAL
survival for a narrower signature corresponding to a normal distribution with a standard deviation of l”C/min (Fig. 9b) varies greatly. For this signature, the survival decreases rapidly as cooling rates are displaced from the optimal value of -SC/ min, and the resulting survival rates are extremely sensitive to the local cooling rate. These graphs illustrate the effect which nonuniform cooling of bulk specimens can have upon total survival. Figure 9d is a composite illustration of survival as a function of holding temperature for - 5, - 10, and - 20Wmin surface cooling protocols, corroborating that survival is a strong combined function of the surface cooling rate and the minimum temperature of exposure.
81
single probe to monitor cooling rates to try to ensure optimal survival of a specific cell type frozen in a finite dimensioned container. The fastest cooling occurs at the center of the container, so that a probe placed in this position would represent a poor choice for purposes of monitoring, even though the geometric symmetry may appear appealing. Quantitative results described in this paper are dependent on the assumed container size and geometry, solution composition, and cell survival signatures. Application to other systems would be dependent upon implementing values of these parameters consistent with the specific system being considered. However, the general trends illustrated by the present study serve to identify the interaction beCONCLUSIONS tween physical and physiological pheThe finite element freezing model pre- nomena in determining the response of dicts that thermal boundary conditions are cells to freezing in finite dimension not propagated uniformly from the exterior systems. into the interior of a cylindrical system, ACKNOWLEDGMENT thus resulting in a nonuniform distribution The authors wish to thank one of the reviewers for of temperature histories and cooling rates throughout the spatial domain. The cooling helpful comments regarding physical interpretation of rates at interior points in the system are de- the analysis. pendent upon both the surface cooling rate REFERENCES and the exterior holding temperature. As 1. Chmiel, H., Stroh, N., and Walitza, E. Vefahrenthe minimum surface holding temperature stechnische aspecte der gefrierkonsevierung decreases, a larger portion of the interior is von blubestandteilen. Biomed. Tech. 25, 32-57 (1980). sensitive to the boundary cooling condi2. Farrant, J. Genera1 observations on cell preservations. This effect occurs because local tion. In “Low Temperature Preservation in cooling rates are calculated from the Medicine and Biology” (M. J. Ashwood-Smith thermal history between the range of -2.8 and J. Farrant, Eds.), pp. l-18. Univ. Park and - 3O”C, and for lower holding temperaPress, Baltimore, 1980. 3. Hayes, L. J. “A users guide to PARAB-A twotures more of the container is cooled to this dimensional linear time-dependent finite elerange of temperatures before the isoment code. TICOM Report 80-10.” Univ. thermal boundary conditions are invoked. Texas, Austin, 1980. As the exterior cooling rate is increased, 4. Hayes, L. J., and Diller, K. R. Implementation of the interior volumetric distribution of phase change in numerical models of heat cooling rates in the interior broadens. transfer. Trans. ASME, J. Energy Resour. Technol. 105, 431-435 (1983). A single average value of cooling rate 5. Hayes, L. J., Diller, K. R., and Chang, H. J. A within a container is in general of limited robust numerical method for latent heat release use in describing survival rates during during phase change. In “Numerical Methods cryopreservation. Neither is it possible to in Heat Transfer” (J. L. S. Chen and K. Vafai, Eds.), pp. 63-70. ASME, New York, 1986. predict a best location for positioning a
82
HAYES ET AL.
6. Hayes, L. J., Diller, K. R., and Lee, H. S. On the definition of an average cooling rate during cell freezing. Cryo-Letr. 5, 97- 110 (1984). 7. Luyet, B. J. An attempt at a systematic analysis of the notion of freezing rates and at an evaluation of the main contributory factors. Cryobiology 2, 198-205 (1966). 8. Mazur, l? Cryobiology, the freezing of biological systems. Science 168, 939-949 (1970). 9. Mazur, P., Leibo, S. P., Farrant, J., Chu, E. H. Y., Hanna, M. G., and Smith, L. H. Interactions of cooling rate, warming rate and protective additive on the survival of frozen
mammalian cells. In “The Frozen Cell” (G. E. W. Wolstenholme and M. O’Connor, Eds.), pp. 69-84. Churchill, London, 1970. 10. Meryman, H. T. The interpretation of freezing rates in biological materials. Cryobiology 2, 165-170 (1966). 11. Morgan, K., Lewis, R. W., and Zienkiewicz, 0. C. An improved algorithm for heat conduction problems with phase change. Znt. .I. Numer. Merhods Eng. 12, 1191-1195 (1978). 12. Rinfret, A. P. Thermal history. Cryobiology 2, 171-180 (1966).