Analysis of crystallization during rewarming in suboptimal vitrification conditions: a semi-empirical approach.

Circumventing ice formation is critical to successful cryopreservation by vitrification of large organs. While ice formation during the cooling part of the cryogenic protocol is dictated by the evolving thermal conditions, ice formation during the rewarming part of the cryogenic protocol is also dependent on the history of cooling and storage conditions. Furthermore, while the exothermic effect of ice crystallization during cooling tends to adversely slow down the desired high cooling rates to ensure ice-free preservation, the same effect under some conditions tends to assist acceleration of rewarming during recovery of the specimen from cryogenic storage when limited crystallization does occur. The current study proposes a computational framework to study the thermal effects of crystallization during recovery from cryogenic storage, using a semi-empirical approach to account for the relationship between latent heat effects and the rewarming rate. This study adds another layer of computational capabilities to a recent study investigating similar effects during cooling. Results of this study demonstrate that the thermal effects of crystallization on the local cooling and rewarming rates cannot be neglected. It further explains how crystallization during rewarming helps in increasing the rewarming rate and, thereby, affects rewarming-phase crystallization. Counterintuitively, this study suggests that the fastest possible rewarming rate at the outer surface of the domain in an inwards rewarming problem is not always advantageous, while the proposed computational tool is essential to find an intermediate optimal rate.

Introduction

The global shortage of donated organs, combined with practical limitations on preservation time, have well been recognized as a major limiting factor to the widespread use of organ transplantation [42]. The short window for successful transplantation of vital organs, typically ranging between 4 and 12 h, creates a major public health challenge. By some accounts for example, if only one half of the currently discarded heart and lungs were utilized for transplantation, the waitlists for these organs could be eliminated in 2–3 years [25]. Cryopreservation— the preservation of biological material at very low temperatures, appears to be the only practical alternative for long-term organ preservation. While individual cells and small specimens can be preserved by classical cryopreservation [2,3,5,10,19,23,24,27,38,39], cryopreservation of large tissues and organs presents additional challenges associated with the underlying principles of heat transfer. The only viable approach to cryopreservation of large organs, such as the heart, appears to be by vitrification [1,12,13,32], where ice crystallization is suppressed and solutions loaded into the organ vitrify (i.e., turn into glass) [1,14,26].

Major advances in vitrification technology have recently been reported, and it is now possible to vitrify entire organs [13], but to do so with full recovery of viability and functionality after transplantation remains a challenge, partially due to rewarming phase crystallization (RPC) [18,22,29,36,37]. During cooling, ice formation is dependent on the nucleation rate and ice growth in the material. Upon rewarming, spontaneous nucleation takes place near the glass transition temperature, Tg, this is followed by ice growth and finally melting. The total ice formed during rewarming is the sum of two phenomena: ice nucleation and growth of those nuclei formed during rewarming, in addition to growth of preexisting nuclei and ice crystals formed during cooling. Consequently, the magnitude of the critical rate of temperature change required to suppress crystallization during rewarming is much higher than that during cooling. Notably, the temperature at which RPC reaches its maximum intensity shifts to a higher value with the increase in cooling and rewarming rates [28].

A variety of crystallization scenarios may hamper the ability to achieve and maintain an ice-free preservation protocol, based on the rate of temperature change history, for example: (i) the cryoprotective agent (CPA) solution in a completely vitrified specimen may crystallize during rewarming, (ii) the CPA solution which has partially crystallized during cooling may continue to crystallize during rewarming, or (iii) a completely crystallized CPA during cooling will be rewarmed as a solid. Furthermore, the rate of ice formation across the specimen may vary based on the local thermal history, making the process as a whole path dependent and spatially nonuniform.

Several theoretical models to calculate the amount of ice formation in CPAs during freezing exist in the literature [4,21,28]. For example, Boutron [4] proposed analytical expressions to describe single cell crystallization during constant rates of cooling and rewarming, where the freezing front is well defined. While these models provide a good insight into crystallization subject to the limitations of calorimetry, such as small sample size and constant rates of temperature change, scaling up the analysis to large tissues and organs requires additional work. In particular, the instantaneous temperature distribution across the organ is non-uniform during cooling and rewarming, and the subsequent thermal history of every given point in the domain is unique to that location. The objective for the current study is to provide a computational framework for that task and provide initial insight into various effects. In particular, a recent study presented a computational framework capable of analyzing partial crystallization during cooling [21], while the current study expands this line of research to include RPC effects.

Mathematical formulation

This study represents an expansion of a previous study [21], which investigated partial crystallization during cooling. Some of the underlying concepts in those studies are similar and are presented here in brief only, for the completeness of presentation. Emphasis in this section is given to the new features of RPC modeling.

Geometrical model

A cylindrical container is analyzed in this study, consisting of two subdomains: CPA solution and the container walls (Fig. 1). It is assumed that the CPA-loaded biological specimen can be first order approximated as having physical properties similar to those of the pure CPA solution [7]. This is an assumption of convenience only for the purpose of the current proof-of-concept presentation, where no limitations exist on expanding the computation to account for spatially distributed physical properties, as would be the case with a real specimen.

Fig. 1.

Schematic illustration of (a) the CPA container analyzed in the current study and (b) its cross section, including applied boundary conditions for convective rewarming.

Two container sizes are analyzed in this study: (A) a diameter of 42 mm, a height of 70 mm, a CPA level height of 62.5 mm, and a wall thickness of 1 mm; and (B) a diameter of 21 mm, a height of 35 mm, a CPA level height of 31.25 mm, and the wall thickness of 1 mm. Container B is one eighth of the size of container A. The wall thickness is typical to ABS [41] containers, having physical properties which can reduce thermomechanical stress under typical conditions [31]. These sizes are chosen in order to analyze the different crystallization scenarios subject to marginal cooling and rewarming rates.

Heat transfer model

Heat transfer in the domain undergoing vitrification is assumed to be solely by conduction: where C is the volumetric specific heat, T is temperature, k is thermal conductivity, is volumetric heat generate rate, with nanowarming as an example. Continuity in temperature and heat flux is assumed on all interfaces.

For the case of nanowarming, the volumetric heating generation is calculated by: where SAR is the specific absorption rate of the nanoparticles-CPA cocktail, and C is the nanoparticle concentration.

The cryopreservation process is modeled to take place in the convective environment of a commercial controlled-rate cooling chamber [15]: where is a normal to the container outer surface, U is the overall heat transfer coefficient, and the indexes j and c refer to the container and cooling chamber, respectively. A previous experimental investigation under similar conditions suggested U = 350 W/m2-°C for a typical controlled-rate cooling chamber [15].

Enthalpy approach to solve the partial crystallization problem

The enthalpy approach can be conveniently presented as accounting for latent heat effects through the use of an effective specific heat property [33]: where, h12 is the enthalpy change within the temperature range of T1 and T2, C is the effective specific heat (also known as the apparent specific heat), C is intrinsic specific heat, and L12 is the latent heat released between the corresponding temperatures.

While the previous study applied Eq. (4) only during cooling [21], the current study expands its usage to the rewarming portion of the cryopreservation protocol. This expansion is non-trivial as the process is dependent on the thermal history, where the crystallization potential during rewarming is dependent on the crystallization events that took place during cooling. Examples of this path-dependency effect include: (i) completely vitrified material may partially crystallize during rewarming, when the rewarming rate is below the critical rewarming rate (CWR) and the temperature is below the heterogeneous nucleation point; (ii) partially crystallized material during cooling continues to crystallize during rewarming before melting; (iii) completely crystallized material during cooling is rewarmed as solid. Fig. 2 illustrates a generic effective specific heat model during cooling and rewarming, where the dashed line represents a simplified presentation of the effective specific heat, C. Recall that the area under the curve equals to the latent heat, L12, and this value is commonly measured.

Fig. 2.

Graphical illustration of an effective specific heat, combining the intrinsic specific heat and the latent heat, where a family of similar curves is used in the current study to correlate phase transition with (a) the cooling rate and (b) the rewarming rate. As the material cools down, crystallization may start at T and terminate at T when subjected to subcritical cooling rates, while rewarming-phase crystallization may take place between T and T when subjected to subcritical rewarming rates; either way, melting will follow between T and T.

While the effective specific heat can be measured directly with a differential scanning calorimeter (DSC) and the outcome may be thermal-history dependent, it is inferred in this study from literature data for constant cooling and rewarming rates. Furthermore, the DSC measurements are performed on microsamples, which are assumed to have uniform temperature distribution at any given time, whereas nonuniform temperature distribution is evident in the analysis of large size vitrification. To account for that effect, a temperature-dependent and spatially distributed specific heat is used for the solution of Eq. (1), where the localized value is identical to the DSC measurements (or to the assumed values in the current study). Consistently, the extent of crystallization events along the cryopreservation protocol is calculated as:

Thermal properties and parameters

The objective of the current study is twofold: (i) to test the computational framework and compare crystallization during cooling to that during rewarming, and (ii) to investigate commonly used rewarming protocols. While the analysis approach presented here is material-independent, examples presented in this study relate to DP6 (3 M dimethyl sulfoxide and 3 M propylene glycol in a suitable vehicle solution) and its combination with 0.175 M sucrose as synthetic ice modulator (SIM), which is relevant to parallel experimental studies [40]. In particular, an idealized CPA (ICPA) physical behavior is selected as displayed in Table 1, where both the phase transition temperature range and the latent heat are temperature independent, for reasons that are further discussed in the Results and Discussion section.

Table 1

Material properties used for thermal analysis in the current study.

Property	ICPA	DP6 + 0.175 M Sucrose	ABS
Thermal conductivity, k (W/m – °C)	Eq. (6)	Eq. (6)	0.17 [41]
Thermal diffusivity, α (m2 /s)	–	–	1.11 × 10−7 [41]
Intrinsic specific heat, Cp (J /kg°C)	3000	Eq. (7)	–
Total latent heat, L (J /kg)	35,000	490 [40]	–
Density, ρ (kg /m3)	1000	1000	–
Critical cooling rate, CCR (°C /min)	5	<1 [40]	–
Critical rewarming rate, CWR (°C /min)	5	15 [40]	–
Onset of melting/Completion of RPC, T1 (°C)	−50	−45 [40]	–
Onset of RPC, T2 (°C)	−75	−57 [40]	–
Completion of Melting, T3 (°C)	−25	−33 [40]	–
Glass transition temperature,Tg (°C)	−100	−113 [40]	–

The thermal conductivity of DP6 within the phase-transition temperature range is calculated as a mixture of the amorphous and crystallized states [6]: where the temperature is given in Celsius degrees and the thermal conductivity in W/m-°C.

The specific heat for DP6+0.175 M sucrose was interpolated from pure DP6 [21] and DP6+0.6 M sucrose [30]: where the temperature is given in Celsius degrees and the specific heat in J/kg-°C.

In the absence of relevant data for nanowarming in DP6, representative values are taken from experiments on VS55 mixed with silica-coated nanoparticles (sIONP = EMG-308 Ferrotec) excited at a field strength of 20 kA/m and frequency of 360 kHz [17]. For this example, SAR is taken as 319 W/m3 for a C = 5 mg Fe/mL.

In order to test the proposed computational framework, three simplifying assumptions are made: (i) rewarming rate-independent onset and completion of crystallization, (ii) the probability of crystallization in an element is independent of the amount of crystallization in adjacent elements, and (iii) the mass transport effect on the heat transfer process is negligible. The latter assumption is associated with the high concentration of the CPA solution, where the viscosity is relatively high, and the magnitude of heat diffusivity is orders of magnitude higher than that of the mass diffusivity. In practice, these assumptions permit decoupling of the heat transfer problem from the mass transport problem, an approach that is consistent with the first-order analysis presented in this study. Nonetheless, a higher order solution is conceivable subject to additional computation efforts once a proof-of-concept is established.

Computational framework

The solution to the heat transfer problem is achieved with a commercial finite-elements analysis (FEA) code (ANSYS 19.1), based on the currently proposed framework. To account for the rewarming rate dependent phase transition process, an ANSYS parametric design language (APDL) script was composed using a multiframe restart technique, following the flowcharts illustrated in Figs. 3,4,5.

Fig. 3.

Computational framework for heat transfer simulations of CPA rewarming, which can account for partial vitrification.

Fig. 4.

Preprocessing detail in the computational framework displayed in Fig. 3.

Fig. 5.

Postprocessing detail in the computational framework displayed in Fig. 3.

The selected ANSYS solver uses a nonlinear transient thermal solution based on the full Newton-Raphson method within each timestep. In order to increase the precision in rewarming rate calculations (Fig. 5), while eliminating the natural perturbations associated with the changing effective specific heat, a weighted moving average method was applied: where n is the number of the most recent rewarming rate data points used for a specific element, i is the data index, and f denotes the filtered rewarming-rate value. Note that Eq. (8) describes a linearly decreasing weight with time, where n = 10 was found adequate based on a convergence study for the specific thermal histories and material properties.

Specifically, at each time step, the ADPL script loops through every element in the mesh to sequentially calculate: (i) the rate of temperature change, (ii) the accumulated amount of crystallization, and (iii) all other material properties according to the state of the material and its temperature. Here, an element is considered to be undergoing crystallization when: (i) its temperature is found within the phase transition temperature range, (ii) the rewarming rate is below the critical rewarming rate for the specific CPA solution, and (iii) the accumulated amount of crystallization was less than 1 (i.e., less than 100% crystallization) in the previous time step. The amount of crystallization at each time step and in each element is calculated using Eq (5).

Rapid changes in thermal properties can result in periodic temperature instability and, in turn, in rewarming rate perturbations [21]. Hence, a convergence study was performed for the number of tetrahedral elements used (56,213 nodes and 15,456 elements), the simulation time step (Δt = 0.25 s), the time-interval for multiframe cycles (Δt = 0.5 s), and discretization of the effective specific heat (ten divisions of full scale) [21].

Results and Discussion

A thought experiment using a simplified protocol and an idealized CPA solution (ICPA)

It may be expected that the extent of crystal formation during cooling would be equal to the extent of crystal formation during rewarming of a completely vitrified material when the magnitude of the cooling rate in the first case equals the magnitude of the rewarming rate in the second case, providing that the critical cooling rate (CCR) and the critical warming rate (CWR) are also identical in magnitude and subject to all other assumptions listed in the mathematical formulation. This expectation can be tested by comparing simulations of two cryogenic protocols, one subject to a constant cooling rate from an initial uniform temperature down to a cryogenic storage temperature, and the other subject to a constant rewarming rate from the storage temperature back to the initial temperature. In order to make these two hypothetical cases similar for the purpose of comparison, the temperature difference between the initial temperature and onset of crystallization (between T0 and T1, respectively) in the cooling case is set to be equal to the temperature difference between storage temperature and the onset of crystallization during rewarming (between Ts and T2, respectively). Similarly, the phase transition temperature range, as well as the temperature difference between completion of crystallization and the final temperature is kept the same for both cases.

With the above thought experiment in mind, it is well appreciated that the CWR is typically an order of magnitude higher than the CCR for practical CPAs and that the cryogenic protocols are more complex than those outlined above. Nonetheless, testing an idealized CPA having the same CCR and CWR and subject to the above cooling and rewarming conditions seems like an insightful thought experiment in order to test the proposed computational framework in an increasing level of complexity, while developing a benchmark for the more complex and realistic scenarios discussed later in this study.

For the purpose of this thought experiment and inspired by the DP6 physical properties, Table 1 lists the properties of the selected ICPA. In addition, an initial temperature of −100 °C and final temperature of −5 °C were selected for the rewarming protocol, such that the difference between the initial temperature and onset of RPC (between Ts and T2, respectively) and the difference between end temperature of RPC and final temperature (between T1 and Tf, respectively) mirror the respective temperatures during the cooling protocol. These parameters are consistent with the cooling studies performed previously [21]. This part of the study was performed on the dimensions of Container A (see Geometric Model section).

Fig. 6(a)–6(d) display the accumulated relative crystallization throughout the cooling process at selected cooling rates. The cooling protocol had an initial temperature of −25 °C and a storage temperature of −120 °C. These protocol parameters were chosen to observe different partial crystallization scenarios. It can be observed from Fig. 6(a) that complete crystallization occurs in most of the domain for the case of subcritical cooling rate at the surface. However, complete crystallization is not observed at the center of the domain when the outer boundary is cooled at a subcritical rate, due to cooling acceleration associated with the heat transfer process in this inwards cooling process. While this effect may be counterintuitive, it has already been reported previously in different systems, and was termed the centerline effect [21,34]. It can be observed from Fig. 6(b)–6(d) that the size of the vitrified region close to the outer surface increases with the cooling rate in supra-CCR. Notably, crystallization is not completely avoided in any of those cases due to the decay of the cooling rate with the distance from the cooled surface, which is the result of an effect commonly referred to as thermal inertia [21]. Regardless, the centerline effect is still observed in those latter cases, whereas the maximum extent of crystallization is observed at some intermediate distance between centerline and the outer surface.

Fig. 6.

Representative results for cumulative ice formation in ICPA contained in the vial illustrated in Fig. 1 during cooling (a)–(d), when the convective surrounding cools at the specified rates from an initial temperature of −25 °C down to −120 °C and during rewarming (e)–(h), when the convective surrounding warms at the specified rates from a storage temperature of −100 °C up to −5 °C.

Fig. 6(e)–6(h) display the maximum portion of crystallization that existed at any point in time across the domain during rewarming. Recall that here the analysis is provided for a previously vitrified material, and that the crystallization presented occurred during rewarming only. Similar trends are observed between the cooling cases and the respective warming cases, although the extent of crystallization is much lower during rewarming. For example, the maximum amount of crystallization at any point in the domain during rewarming at 30 °C/min (Fig. 6(h)) is about 30% while the maximum amount of crystallization during cooling at the same rate resulted in 100% crystallization (Fig. 6(d)). Moreover, the region size of complete crystallization in sub-CWR is always much smaller than the respective region size during sub-CCR.

The reason for the lower rate of crystallization during rewarming can be explained in association with the exothermic effect of crystallization and the endothermic effect of melting. While crystallization during cooling always slows down the cooling process and, thereby, further promotes crystallization, crystallization during rewarming initially accelerates rewarming but then slows it down during melting; these effects are also reflected in Fig. 2(a) and (b), respectively. The combined effects of crystallization-melting during rewarming yields a net acceleration of the process in comparison with cooling. To put that in context, recall that the ICPA and the protocol analyzed here are imaginary only, which were specifically designed in order to isolate a specific effect—the contribution of crystallization during rewarming to the local rewarming rate. Evidently, this effect is potentially significant.

Recall that the CWR is an order of magnitude higher than the CCR and that a higher ratio of crystallization is expected in a real CPA under the same cooling rates but, nonetheless, the thermal effects of crystallization on the local cooling and rewarming rates cannot be neglected. In other words, the local cooling and rewarming rate is not only dependent on the cooling and rewarming rates at the boundaries, but also on a combination of physical properties, including the latent heat, the specific heat, and the thermal conductivity.

Analysis of partial crystallization during rewarming

This part of the study focuses on DP6 as a representative CPA, which is investigated experimentally in parallel and recent studies [1,9,40]. Additionally, this part of the study focuses on Container B (see Geometric Model section) and the addition of 0.175 M Sucrose to the CPA as a synthetic ice modulator (SIM) in order to explore marginal conditions of vitrification, Table 1 [40]. The range of thermal protocols investigated here is displayed in Fig. 7. Since this study focuses on rewarming, the cooling history was selected such that complete vitrification will be achieved in all protocols: (a) initial rapid cooling of 40 °C/min from −20 °C to −120 °C; (b) temperature hold at −120 °C until the specimen reaches thermal equilibrium below the T [35]; (c) further cooling to a cryogenic storage temperature of −150 °C; and (d) indefinite temperature hold at cryogenic storage. While the selection of segments b and c has no bearing on the thermal analysis, their variation may affect the likelihood to fracture, which is beyond the scope of the current study [8, 35]. During rewarming from cryogenic storage, the solution is either continually rewarmed from −150 °C to 0 °C (Protocols I-III), or slowly rewarmed to an intermediate temperature T (−120 °C) and held there until thermal equilibrium, before subsequent rapid rewarming to 0 °C (Protocol IV). Once the outer surface of the container reaches 0 °C, it is held constant thereafter.

Fig. 7.

Thermal history in the cooling chambers (a convective environment) used to investigated RPC in DP6+0.175 M Sucrose.

It can be observed from Fig. 8, that despite the supra-CWR by convection at the boundary, partial crystallization did occur in all constant-rate rewarming cases (Protocols I through III). Regardless of the variation in rewarming rate at the boundary between cases, more than one quarter of the domain experienced some level of crystallization. Fig. 8(d) displays the volume portion of the domain that experienced crystallization above selected thresholds. The significance of these thresholds is twofold: (i) some mammalian cells can tolerate marginal volumes of ice under specific conditions [1], and (ii) there exists some numerical uncertainty in the assessment of the crystallization ratio, which is inherent to FEA simulations as well as to the application of the enthalpy approach. With the numerical considerations in mind, while the 0% threshold is unrealistically strict, it is reasonable to expect that the confidence in simulated crystallization ratio is in the 0.5%–1% range. The uncertainty evaluation is also related to the numerical convergence study. Either way, the crystallization ratios are weakly dependent on the rewarming rate on the outer surface.

Fig. 8.

Representative results for the rewarming rate study on DP6+0.175 M sucrose for the thermal protocols I, II and III displayed in Fig. 7, where (a)–(c) display the amount of crystallization formed throughout the rewarming process under various boundary conditions, starting from a completely vitrified solution at cryogenic storage; (d) displays the volume ratio with RPC for various crystallization threshold values.

While the results in Fig. 8 refer to the accumulated crystallization, the analysis of the instantaneous rewarming rates provides additional insight into the relationship between rewarming conditions at the surface and the thermal response within the domain. For this purpose, Fig. 9 displays the localized rewarming rate as a function of temperature along the CPA-air interface at the top of the container for Protocols I through III. This interface was selected for the analysis, as it experiences the slowest rewarming rates in the domain. In particular, two temperature ranges are of interest for the analysis, near glass transition, where nucleation rate is maximal, and within the phase transition temperature range, where ice growth is maximal.

Fig. 9.

Rewarming rate as a function of boundary temperature and distance from the outer wall surface along the CPA-air interface (isolated upper boundary in Fig. 1), for the thermal protocols I, II and III displayed in Fig. 7.

It can be seen from Fig. 9 that the rewarming rate within the domain around glass transition of DP6+0.175 M Sucrose (T = 113 °C) increases with the increasing rewarming rate at the boundary. This means that the material is exposed for shorter periods to these temperatures at higher rewarming rates, which reduces the potential for nucleation, although the analysis of nucleation is not included in this study. Note that the thermal effect associated with ice nucleation is negligible, and that the current computational framework ignores it. Further note that the rewarming rate does not change monotonically with temperature, where its slow down around glass transition can be inferred from Fig. 9. This observation is associated with the threefold increase in the specific heat as the rewarmed material passes the glass transition temperature, Eq. (7).

It can be observed from Fig. 9 that the rewarming rate within the phase transition temperature range is only moderately dependent on the rewarming rate at the outer surface, an effect that diminishes with the increasing distance from that surface. This observation suggests that minimizing crystal growth by increasing the rewarming rate at the surface may present an ineffective approach, to an extent that is dependent on the specific container size and the CPA physical properties.

At 100 °C/min in Protocol III, the chamber reached 0 °C four times as fast as in Protocol I (25 °C/min). In turn, the outer boundary was held at 0 °C thereafter for a longer time in Protocol III until the domain reached thermal equilibrium. Consequently, the rewarming rate at the center of the domain is lower in Protocol III as it approaches the phase transition temperature range. This effect can explain the counterintuitive observation that a higher ratio of crystallization occurs in protocols having higher rewarming rates at the outer surface in Fig. 8(d).

A potentially adverse effect of rapid rewarming, and especially in lower temperatures is thermomechanical stress, which may results in structural damage [8,35]. To overcome this risk, a three-step rewarming protocol has been suggested previously, with Protocol IV in Fig. 7 as an example. The total volume of solution that underwent RPC in Protocol IV was 31% as compared to 32% in Protocol III, for the case of 0% crystallization threshold (Fig. 8). This shows that starting with rapid rewarming at higher temperatures near glass transition temperature did not adversely affect the likelihood for vitrification in the domain.

Nanowarming as an alternative for surface rewarming

While container B was carefully selected to study partial crystallization in marginal conditions, and thereby the applicability of the proposed computational framework, it is clear that a size limit exists above which RPC cannot be avoided during an inwards rewarming process. To demonstrate that, the discussion now returns to the larger container A and two additional cases, the first with a very rapid rewarming rate at the boundary of 100 °C/min and a practically infinite heat transfer coefficient by convection, and the second with volumetric rewarming, as would be expected during the application of nanowarming [11,17,26]. Fig. 10 displays the results from those two cases assuming DP6+0.175 M Sucrose properties and a uniform concentration of 5 mg Fe/mL of silica coated nanoparticles [17]. It is clear from Fig. 10 that uniformly distributed nanoparticles can lead to the conditions that prevent RPC. Nonetheless, nanowarming may not prevent RPC in some cases when nanoparticles are distributed nonuniformly. Parallel studies are now focused on the effect of nonuniform nanoparticles distribution, considering the specific organ detail and nanoparticles loading limitations, while taking advantage of the computational framework presented above [16,20].

Fig. 10.

Cumulative crystallization subject to various rewarming methods using DP6+0.175 M Sucrose for the container illustrated in Fig. 1, where the initial temperature is −150 °C the final temperature is 0 °C, and subject to the following conditions: (a) convective rewarming at 100 °C/min, (b) maximum possible convective heating (infinite heat transfer coefficient) at the boundary at a rate of 100 °C/min, and (c) volumetric heating (nanowarming).

Summary and conclusions

A computational framework is proposed in this study to simulate vitrification processes with emphasis on crystallization formed during rewarming. The computational framework uses the enthalpy approach to model phase change effects, while assuming the availability of empirical correlations of the thermal history with the likelihood of crystallization. Such correlations must be obtained from DSC measurements, which is beyond the scope of this study. This semi-empirical approach simplifies the analysis, which accommodates spatial distribution of the key variables including the thermophysical properties, the thermal history, the instantaneous extent of crystallization, and the dependency of crystallization on the local rewarming rate. To demonstrate its applicability, the computational framework is applied using the FEA commercial software ANSYS, while examples are given for CPAs contained in a cylindrical polymeric container.

In order to reduce evaluation complexity of the multi-variable problem, an idealized CPA is investigated first, with a simplified temperature-independent latent heat expression upon rewarming. The selection of the ICPA properties is inspired by parallel studies on DP6 as a CPA. Results of this part of the study demonstrate that the exothermic effect of crystallization increases the local rewarming rate, with a trend that slows down further localized crystallization. On the other hand, while the endothermic effect of the subsequent melting slows down the localized rewarming rate, the net effect of melting is smaller than the net effect of crystallization on the localized rewarming rate. Either way, the latent heat effects during rewarming are significant and cannot be neglected.

The second part of the analysis in this study focuses on DP6 combined with the SIM 0.175 M Sucrose. The rewarming protocol consists of two key segments: constant rate rewarming at the boundary from below glass transition to a set temperature above the heterogeneous nucleation point, and a constant outer surface temperature thereafter. While faster rewarming rates at the surface were found beneficial to suppress crystallization in areas closer to the wall, their effect diminishes towards the center of the domain. Notably, the rewarming rate accelerates towards the center of the domain in comparison with some intermediate radial location due to the underlying principles of heat transfer, and the effect has been reported previously as the centerline effect.

From a holistic perspective on the rewarming thermal protocol, the outer surface of the domain needs to be held constant longer in the second segment of the protocol until the risk of crystallization is eliminated, when the rewarming boundary condition is initially faster. This means that there may be more opportunity for crystallization in inner regions in that case. To mitigate this adverse effect, an intermediate rewarming rate is required, which is not known a priori but can be determined by using the computational framework proposed in this study. With this conclusion in mind, it must be emphasized that the execution of the proposed computational framework relies on significant simplification and the results should be carefully considered. Furthermore, refining the underlying empirical correlations by means of experimental studies can improve the certainty of modeling predictions. Regardless, the counterintuitive observation that the amount of overall crystallization in the domain changes trend around some intermediate rewarming rate at the outer surface is independent of the quality of the empirical correlations and is rooted in the underlying principles of heat transfer.

Finally, the computational framework developed here is not limited to inward heating cases and can handle volumetric heating with nanowarming as an example. While the current study demonstrates the superiority of nanowarming over surface rewarming, nonuniform nanoparticle distribution due to the loading process into biological materials may affect the uniformity of rewarming rate with yet to be presented thermal effects on RPC.