Thermomechanical stress analysis of rabbit kidney and human kidney during cryopreservation by vitrification with the application of radiofrequency heating.

This study presents a computational framework for thermomechanical stress analysis in a specimen undergoing cryopreservation, with emphasis on radiofrequency (RF) heating for recovering from cryogenic storage. In particular, this study addresses cryopreservation by vitrification, where the specimen is stored in the amorphous phase (vitreous means glassy). In broad terms, the relatively high cooling and rewarming rates necessary for vitrification result in differential thermal expansion in the specimen, which is the driving force for thermomechanical stress. Thermomechanical stress can lead to structural damage, such as fractures or plastic deformation, rendering the specimen useless. Not without technical difficulties, those hazardous effects during the rewarming phase of the protocol can be mitigated by applying volumetric heating, with RF heating as an attractive means. The proposed computational framework in this study addresses the coupled electromagnetic, thermal and solid mechanics fields, using commercially available solvers. This study advances from a spherical-case benchmark to realistic models of the rabbit kidney and the human kidney. Results of this study suggest that structural damage to the brittle material can be prevented when stress relaxation is facilitated around the glass transition temperature. Furthermore, this study suggests that volumetric heating is necessary to surpass the critical rewarming rate, while benefiting from lowering the overall thermomechanical stress during recovery from cryogenic storage. More broadly, the computational framework presented here can be used for the optimization of the RF heating parameters, chamber specifics, specimen container shape, and the thermal protocol in order to preserve structural integrity in the specimen.

Introduction

With direct implications to the future of biobanking and transplant medicine [24], cryopreservation via vitrification (vitreous in Latin means glassy) is considered to be the only promising technique for long-term preservation of large-sized tissues and organs [18,71]. With only 10% of the worldwide needs for organ transplantation being met [30], a supply of tissues and organs on demand—facilitated by cryopreservation via vitrification—can potentially create public health benefits on par with curing cancer [24]. Cryopreservation involves cooling tissues loaded with a cryoprotective agent (CPA) solution to cryogenic temperatures for indefinite storage [46]. Vitrification involves relatively rapid cooling of the CPA-permeated specimen, in order to trap the specimen in an amorphous solid, which is facilitated by the exponential increase in viscosity with the decreasing temperature [18]. The temperature below which the viscosity of the CPA is high enough for the specimen to be considered a solid for all practical purposes is known as the glass transition temperature (T). Indefinite cryogenic storage of a vitrified material is commonly targeted to take place below T.

While vitrification success for small tissues (typical size ranging from μm to mm) has been widely reported, such as for stem cells [5,26] and corneas [1,4], vitrification of large specimens (measured in cm or larger) remains a challenge, requiring a multi-parameter optimization process to reach its potential [31]. Key physical phenomena that affect successful cryopreservation include specimen preparation and handling, natural degradation of the biological specimen with temperature, ice crystallization [71,73], CPA toxicity [6], and thermomechanical stress [13,54,55,57,63,68].

Vitrification is the only viable alternative to eliminate the devastating effect of ice crystallization during cryopreservation. Vitrification involves rapid cooling of the CPA-permeated tissue specimen below T, such that the molecules of the CPA cocktail do not have enough time to rearrange and form crystals [18]. During rewarming, rapid heating rates are also necessary in order to prevent rewarming-phase crystallization (RPC) [40]. RPC is generally related to ice nucleation and growth during the rewarming phase of the cryogenic protocol (devitrification), or to crystal growth from nuclei and small crystals already trapped in the amorphous material since the cooling phase of the cryoprotocol [36]. In an absolute scale, the critical rewarming rates (CRR) required to avoid RPC are typically an order of magnitude higher than the critical cooling rates (CCR) required to avoid crystallization during cooling [40,66].

High cooling and rewarming rates, coupled with the typically low thermal conductivity and thermal diffusivity of the CPA and specimen [11], intensify the thermal gradients in the cryopreserved material, which can give rise to high thermomechanical stress (thermal stress in brief herein) [63]. In addition to thermal gradients, thermal stress may also intensify by thermal expansion associated with phase transition [62], thermal expansion mismatch between the specimen and its surrounding medium [67], the surrounding medium and the container [67], and by any combination of the above effects [13,54,55,57,68]. While significant stress typically develops below T, thermal stress may start building up already several degrees above the glass transition temperature depending on several factors [20,55]. The early stage of the rewarming phase has been identified as the most susceptible to structural damage resulting from thermal stress, such as fracturing [13,48,52, 56]. In a broader perspective, thermal stress is path-dependent, where the entire thermo-mechanical history must be considered when trying to interpret the causes for high-instantaneous levels of it [13].

Thermal stresses due to temperature gradients may be reduced by using volumetric heating techniques to rewarm the specimen from cryogenic storage [12,63]. Volumetric heating may also help achieving high rewarming rates to outrun RPC. Not without difficulties, electromagnetic heating in the radiofrequency (RF) range of the spectrum has been presented as an effective means for volumetric rewarming of cryopreserved specimen [60]. Inductive heating of a specimen permeated with magnetic nanoparticles—also known as nanowarming—has been recently demonstrated to moderate the temperature gradients during rewarming from cryogenic storage, thereby reducing thermal stresses [38,63]. While nanowarming utilizes the frequency range of 100–400 kHz [38], direct RF heating without intermediate foreign particle heaters requires electromagnetic frequencies which are three orders of magnitude higher [15]. The current study focusses on analysis of direct RF rewarming, without the addition of magnetic nanoparticles.

Successful experimental studies aimed at rewarming cryopreserved rabbit kidneys at electromagnetic frequencies of 2450 MHz [8] laid down the foundation for extensive research to better understand the role of RF rewarming in the field of cryopreservation. While electromagnetic frequencies of 2450 MHz may not be an issue for rewarming small rabbit kidneys (less than 2 cm thick), uniform heating is not achieved at these frequencies for larger organs such as canine kidneys [47]. Lower frequencies with larger penetration depths [15] are necessary for rewarming of larger organs. While lower frequencies offer larger penetration depths, the associated low power delivery calls for a compromise between the effective volume to be heated and the intensity of heating [15].

The choice of CPA and its concentration are not parameters that significantly affect the heating intensity penetration or its non-uniformity across the specimen. An idealized spherical phantom of a rabbit kidney model has been presented in literature as the basis for experimental [58] and theoretical [3] studies. For a cylindrical [16,60] or a rectangular [35] single-mode cavity used to study RF rewarming of cryopreserved materials, an electromagnetic frequency of 434 MHz was found suitable, when considering a spherical phantom of the kidney [16]. A single-mode RF cavity is excited by a single frequency, which provides a more uniform heating effect, for which the electromagnetic field is also simpler to compute [72]. By contrast, a multi-mode cavity incorporates complex components of the electromagnetic field, which may drive larger temperature non-uniformities in the specimen [69].

Represented by a complex number, the dielectric property of the material represents the extent of energy coupling (the real component) and energy dissipation in the form of heat into the material (the imaginary component) as a result of an external electromagnetic field [21]. Specifically, these components describe the magnitude of polarization and the ability to absorb electromagnetic energy, respectively. The dielectric properties of the CPA vary with its concentration, temperature and the frequency of electromagnetic excitation [17]. Given the sparsely available dielectric data for CPA-perfused tissue system [15], it is useful to note that the dielectric properties of the specimen are most-likely dominated by the dielectric properties of the perfusing CPA [39]. For that purpose, some dielectric data of CPAs are available between −80 °C and room temperature, for the specific electromagnetic frequencies of 27, 108, 434, and 2450 MHz [17,41,43]. Here, 27, 434 and 2450 MHz are a part of the approved industrial, scientific and medical (ISM) radiofrequency bands, while 108 MHz was included by Evans [15] as a practical mid-range frequency.

The uniformity of heating in the specimen correlates with the shape of the specimen [34,44,61], where excessive heating may occur on its boundaries, driven by the so-called effect of thermal runaway. Previous theoretical [15,34] and experimental [61] studies have demonstrated that the uniformity of heating can be related to the extent to which the specimen resembles a sphere, where a perfect sphere would minimize thermal runaway effects. To solve the effect of specimen-shape irregularities, it has been suggested to encase it with a CPA-based gel, such that the overall shape becomes a sphere [15].

Previous RF heating studies suggested the implementation of numerical methods such as the method of moments (MoM) [33] and the finite-difference time domain (FDTD) [25]. Those traditional numerical methods typically solve a less realistic problem, where the thermal and electrical fields are often decoupled but, nonetheless, are computationally expensive. More recently, commercial codes utilizing finite element methods (FEM) to solve the above coupled phenomena became available at affordable computation costs [44].

The current study presents a computational framework to analyze the thermal stress in a specimen as it undergoes a vitrification protocol, where RF heating assists in recovery from cryogenic storage, using the commercially available FEM codes COMSOL and ABAQUS. Given path-dependency effects, thermal stress modeling here is coupled with the thermal and electromagnetic fields. Results using the proposed framework are benchmarked against previous studies on a simplified spherical geometry [25]. Next, more realistic scenarios are taken into consideration by simulating a rabbit kidney and a human kidney model, when stored in cylindrical containers.

Mathematical formulation

Geometric model

First, an ideal spherical specimen is analyzed to benchmark the framework proposed in the current study with a previous study based on FDTD formulation [25], while focusing on RF heating during the rewarming phase of the cryopreservation protocol. Consistent with previous studies [25,34,60], the geometric model consists of a cylindrical RF resonance cavity with a monopole antenna from a coaxial transmission line as the power source, where the specimen is placed at the geometric center of the cavity. Fig. 1A displays a realistic kidney specimen for a similar case, which is substituted with a spherical model in the current benchmark study, having material properties similar to the kidney phantom model measured by Robinson et al. [60]. Since the benchmark study focused on rewarming of a frozen but not completely vitrified rabbit kidney phantom, phase change is included in this analysis using the enthalpy formulation [25], where latent heat effects are expressed in terms of an effective specific heat.

Fig. 1.

Kidney model used for the analysis of thermomechanical stress during cryopreservation by vitrification, involving RF rewarming: (A) a schematic illustration of the kidney model in a RF cavity (not drawn to scale); (B) a cross section of the geometric model displaying the FEA mesh; (C) the kidney model held in a cylindrical CPA container; and (D) a cross section of the kidney model immersed in CPA displaying the FEA mesh.

Next, stress is analyzed in a more realistic kidney model, when stored in a cylindrical CPA container. The geometry was developed by Ehrlich et al. [10], where the human kidney model is essentially a 21-fold scaled-up version of the rabbit kidney model. The cylindrical CPA container for the rabbit kidney is 33 mm in diameter and 46 mm in height. The human kidney CPA container is 90 mm in diameter and 127 mm in height. The edges of the cylindrical CPA container are curved (Fig. 1C) to avoid thermal runaway effects there [15,61], with a radius of 3.6 mm and 10 mm for the rabbit and human kidney containers, respectively. Due to prohibitive computation costs, the container material has not been modeled in this study, which, in practice, would yield similar results to the case of a highly compliant container [63]. For example, the latter scenario can be considered when selecting polymers for container material, such as fluorinated ethylene propylene (FEP) [63].

Due to the underlying principles of RF heating, two different cylindrical heating chambers (i.e., heating cavities) were considered. For the small models—the spherical benchmark and the rabbit kidney—a cavity size of 0.5 m in dimeter and 0.58 m in height were adapted from previous studies [16,25,60]. This cavity resonates at about 434 MHz in the TE111 mode [25]. The geometry of the power source is adopted from the benchmark study [25]. Accordingly, the radii of the antenna and transmission line are 1.5 mm and 3.45 mm, respectively, where the penetration length of the antenna into the cavity is 60 mm. Although high electromagnetic frequency means better power delivery, a compromise must be made to avoid the associated lower penetration depth [15]. In order to achieve more uniform rewarming for the larger, human kidney model, an electromagnetic frequency of 108 MHz is selected. Here, frequency selection is informed by the sparsely available dielectric properties for Dimethyl Sulfoxide (DMSO), which is the key ingredient for commonly investigated CPAs [17]. At that frequency, a cavity having a diameter of 2 m and a height of 2.32 m was selected for heating excitation in the same TE111 resonance mode. The antenna penetration length is adjusted for the larger cavity to 0.24 m.

The geometric model is discretized using the meshing module of COMSOL Multiphysics with 10-noded quadratic tetrahedron elements (Fig. 1B and D). Identical mesh is used for solving the electromagnetic, thermal and stress fields. Using the COMSOL mesh refinement tool, finer meshing is performed close to the power source (Fig. 1B) and within the specimen. A maximum element size is selected while keeping at least six elements within a wavelength [32,75], resulting in 1.35 × 105 elements. Fig. 1D displays three virtual sensors placed at selected points of interest: Pa is selected at the center of the CPA container base, Pb is selected at mid-height of the CPA container wall, and Pc is selected at the volumetric center of the container. These points correspond to locations of maximum stresses at different points in time along the cryopreservation protocol [13,63].

Electromagnetic model

The spatial distribution of the electric field for a nonmagnetic dielectric material in a resonance cavity is determined by the Maxwell equations in frequency domain [2]: where is the electric field, ω = 2πf is the angular frequency of the electromagnetic (EM) field, f is the frequency, ε0 and μ0 are the electric permittivity of free space (8.854 × 10−12 F/m) and its magnetic permeability (1.257 × 10−7 H/m), respectively, and ε is a complex relative permittivity given by ε = ε′-jε”, where ε′ is the dielectric constant and ε” is the dielectric loss factor.

The walls of the resonance cavity are assumed to be perfect electric conductors [51]: where is a unit vector perpendicular to the resonance cavity wall.

The resonance frequency of the cylindrical cavity changes when the specimen is placed in it, depending on the dielectric properties of the specimen. These properties undergo a change of about one order of magnitude as the temperature increases during RF heating (Fig. 2). Consequently, the frequency of the power source must be adjusted continually during heating, in order to match the changing resonance frequency in the cylindrical cavity [35].

Fig. 2.

Dielectric properties of 7.05 M DMSO at 434 MHz and 108 MHz: (A) relative dielectric constant, ε’ (the real number component) and (B) relative dielectric loss factor, ε” (the imaginary number component).

In the current analysis, the change in resonance frequency is accounted for by adjusting the power source frequency f at every time step of the numerical simulation, such that the cavity is excited at TE111 resonance mode throughout the rewarming process. It has been experimentally demonstrated that such adjustment of power source frequency can improve temperature uniformity while increasing the rewarming rate of the cryopreserved specimen [35]. The numerical solvers applied in this study are used to calculate the resonance frequency of the system as a function of the specimen temperature. The results from the resonance frequency calculation are used as input for the RF heating simulation, where the operating frequency of the power source is defined as a function of the maximum temperature inside the specimen.

Heat transfer model

Following the assumption that heat generation due to viscous dissipation is negligible, the heat transfer problem is decoupled from the solid mechanics problem [12,13,20,63], where the thermal history is used as input for stress analysis. Heat transfer is assumed solely by conduction: where ρ is the density, C is the specific heat, T is the temperature, k is the thermal conductivity, P is the heat generation due to RF heating. No contact resistance to heat transfer is considered between the CPA and specimen. The heat generated within the specimen due to RF heating is given by [2]: where E is the root mean square value of the electric field.

The Einstein model for internal energy storage is used to describe the temperature-dependent specific heat [20]: where N is the number of oscillators, k is the Boltzmann constant (1.38 × 10−23 J/K), M is the molecular weight, T is the absolute temperature, θ is the Einstein temperature, ℏ is the reduced Planck’s constant (1.054 × 10−34 J s), and ω is the frequency of oscillation of the molecule (6.415 × 1013 Hz) [20]. Temperature-dependent thermal conductivity is assumed in this study (Table 1), which accounts for free convection effects in the CPA domain at temperatures above −60 °C [20].

Table 1

Representative material properties of the CPA (7.05 M DMSO) used in this study.

Property	Value
Viscosity, Pas	1.77 × 104	−106 °C ≤ T	[42]
	2.8190 × 10−27 e−0,6447 T	−143 °C ≤ T < −106 °C
	4.06 × 1014	T ≤ −143 °C
Glass transition temperature, °C	−132		[40]
Density, kg/m3	1100		[29]
Thermal conductivity, W/m °C	0.899 + 1.01 × 10−2 T	−60 °C ≤ T ≤ 25 °C	[20]
	0.312 + 2.54 × 10−4 T	−170 °C ≤ T ≤ −60 °C
Specific heat, J/kg °C	Eq. (5)		[20]
Thermal Expansion coefficient, 1/°C	1.1 × 10−4		[14,29,50,53]
Young’s modulus, GPa	0.8		[27,28]
Poisson’s ratio	0.25		[13]
Complex dielectric properties	Fig. 2		[17]

Heat transfer between the walls of the CPA and the cooling chamber environment is given by [20]: where h is the effective heat transfer coefficient which accounts for convection and radiation heat transfer [20], is the direction normal to the wall, and subscripts w and ∞ refer to the CPA wall and surrounding air temperature within the cooling chamber, respectively. Consistent with previous studies [20,64], the value of the overall heat transfer coefficient, h is kept constant at 346 W/m2K, which is typical to controlled-rate coolers.

Solid mechanics model

The Maxwell fluid model is used to determine the total strain rate in the rabbit kidney and human kidney models [28,42,55]: where creep, elastic, and thermal strain rates are calculated by: where σ is the stress tensor, is the deviatoric stress tensor, η is the viscosity, E is the elastic modulus, ν is the Poisson ratio, is the identity matrix, tr is the trace matrix, and β is the linear thermal expansion coefficient. Note that no scientific relationship exists between the Maxwell fluid model used in the thermal stress analysis, as referenced here, and the Maxwell model used to describe the electrical field in Eq. (1).

Materials and methods

Material properties

A successful vitrification protocol on a rabbit kidney specimen has been reported recently, using M22 as a CPA cocktail [19]. In the absence of a complete set of material properties for M22, and consistent with previous studies on thermomechanical stresses, a 7.05 M DMSO solution was selected as a representative CPA for all material properties in this study [10,28,50], Table 1. Note that previous thermal [10] and mechanical [28,29,50] studies have shown that the thermomechanical properties of the CPA-loaded tissues in cryogenic temperatures are essentially dominated by the thermomechanical properties of the vitrified CPA. The temperature dependency of the dielectric properties are presented in Fig. 2 [17]. Similarly to the mechanical and thermal properties, it has been shown that the dielectric properties of the CPA dominate those of the perfused tissue [39]. Due to unavailability of CPA dielectric material properties in the entire temperature range of interest for the current study, constant values are assumed outside of the available temperature range in Fig. 2, where the implication of this assumption is discussed below.

Thermal protocol

Two thermal protocols are investigated in this study. The first is related to the benchmark case, having an initial temperature of −60 °C, where only the rewarming phase of the cryopreservation protocol is considered [25]. The specimen is heated using a coaxial probe with either a 300V or a 600V electrical potential, combined with 50Ω impedance [25]. The simulation is terminated when a maximum temperature of 20 °C is reached anywhere in the specimen.

The second thermal protocol refers to a typical thermal history, as it pertains to vitrification of the kidney, Fig. 3, where compiled temperature curves are obtained from a previous study by Ehrlich et al. [10]. This protocol consists of four phases: cooling (“C”), storage (“S”), RF rewarming (“RF”), and convective rewarming (“CR”). Consistent with the application of M22 as a CPA for kidney cryopreservation [10,19], the kidney specimen is simulated to be cooled in three steps: (i) precooling to −22 °C while perfusing the organ with CPA, (ii) immersion into precooled CPA solution at a lower temperature of −50 °C, and (iii) cooling according to experimentally obtained thermal histories at the organ surface, down to a storage temperature of −135 °C which is 3 °C below glass transition, Fig. 3(A). Since the physical properties of 7.05 M DMSO is used for thermomechanical analysis in this study, the glass transition temperature of it was selected (−132 °C) rather than that of M22 (−123 °C). Nonetheless, the generality of the conclusion presented below are expected to be insensitive to the specific glass transition temperature selected, but only to temperature dependency of viscosity as it approaches that temperature. Once the outer surface of the specimen reaches −135 °C, it is held constant until mechanical stress simulation results indicate that the residual stress relaxes below 1% of the maximum stress at anytime and anywhere in the domain. Note that the viscosity value at −135 °C, although below glass transition, still permits stress relaxation.

Fig. 3.

Representative thermal history in the cooling chamber for three storage temperatures: (A) −135 °C, (B) −150 °C, and (C) −196°C. The thermal protocol consists of four phases: Cooling (“C”), Storage (“S”), RF rewarming (“RF”) and convective rewarming (“CR”). The duration of temperature hold for annealing during cooling (in Figs. B & C) and the duration of storage (in all cases) are varied such that thermal equilibrium is achieved at the end of each respective segment, and therefore the stress reached steady-state values.

The effect of storage temperature on thermomechanical stress was further investigated by incorporating an additional slow cooling step below −135 °C and after stress relaxation, at a rate of 0.27 °C/min down to storage temperatures of −150 °C and −196 °C [10] (Fig. 3B and C). The lower cooling rate below the glass transition temperature is required to eliminate substantiable thermal stress at temperatures where the CPA behaves primarily as an elastic solid [13,63]. Once the storage temperature is achieved in the cooling chamber, it is held constant thereafter until thermal and stress equilibrium are achieved.

A two-step RF heating protocol is considered for recovering the specimen from cryogenic storage. Initially, the specimen is low-power heated until its maximum temperature reaches −90 °C, and then it is high-power heated until its maximum temperature reaches 27 °C (±2 °C). An initial low heating power of 2 W and 40 W was selected for the rabbit kidney and the human kidney models, respectively, to achieve a rewarming rate below 2 °C/min. Such a low rewarming rate has shown to produce sufficiently low thermal stress in lower temperatures [13, 63]. The second heating step was designed to exceed 100 °C/min anywhere in the domain within the temperature range of −90 °C to −30 °C, where the rate of ice crystal growth is high. This requires a high heating power of 0.8 kW and 50 kW for the cases of rabbit kidney and human kidney, respectively. Note that the power rating refers to the overall power supplied to the system using the antenna power source (Fig. 1A). The power absorbed by the specimen, on the other hand, depends on the resulting spatial electric field in the RF cavity and the dielectric properties of the specimen according to Eq. (4). The power difference between the applied power and the absorbed power is reflected back due to impedance mismatch during rewarming. Such impedance mismatch has been previously described in experimental studies [58,60].

For benchmarking purposes, thermal stress development during RF rewarming was compared with the readily available method of convective rewarming. Similarly to RF rewarming, a two-step rewarming protocol is considered, starting with slow heating up to −90 °C, followed by fast heating thereafter. During the convective rewarming, the effective heat transfer coefficient (h, Eq. (6)) is kept constant at 346 W/m2K [20,64], while the average surface temperature from the RF rewarming counterpart is utilized as the chamber environment temperature (T∞, Eq. (6)).

Computation framework

The computational framework comprises a stepwise scheme, beginning with a bidirectional-coupled problem of the electromagnetic and thermal fields to solve for the transient temperature distribution within the specimen, while incorporating temperature-dependent dielectric properties [9,25,34,44,76]. This is followed by a unidirectional-coupled problem of the thermal and mechanical stress fields, where the specimen is modeled as a Maxwell fluid with temperature-dependent viscosity [12,13,63]. This unidirectional coupling implies that the thermal field may give rise to mechanical stress, but mechanical stress does not generate significant heat. Consistent with previous work, the bidirectional-coupled problem of the electromagnetic and thermal fields is solved using commercial FEA software package COMSOL Multiphysics [9,44] and the unidirectional-coupled problem of the thermal field and mechanical stress distribution is solved using commercial FEA software package ABAQUS [12,13,63,64]. The resulted transient temperature distribution from COMSOL serves as input to ABAQUS after conversion using a custom MATLAB function.

Results and discussions

Comparison with a benchmark solution

Fig. 4 displays comparison of the maximum temperature for a 36 mm in diameter spherical model and the benchmark cases [25], where data in the benchmark study is reported at one location only. Very good agreement can be observed below −20 °C, in a temperature range which is relevant to preventing crystallization during rewarming. By contrast, poor agreement is observed above −20 °C, where phase transition takes place, an effect which is irrelevant to the study of cryopreservation by vitrification. At least three reasons may contribute to this disagreement: (i) the previous study was based on a finite difference simulation method and a proprietary code, while the current study is based on FEM and a robust commercial code; (ii) the electric field was updated periodically in the benchmark study, when the dielectric properties of the specimen changed by more than 5% between two consecutive time steps, while the electric field is being updated continually in the current study; and (iii) differences in numerical discretization. It is difficult to determine which solution approach and code is more accurate in the absence of experimental data, but the authors of the current study express a higher confidence in the robust commercial code and the optimized runtime parameters, and especially in lieu of convergence studies. Nonetheless, similarity in trends between the two studies below −20 °C is encouraging for the purpose of vitrification.

Fig. 4.

Benchmarking study of the thermal history in a simplified spherical geometry [25].

Storage temperature effect

The current study investigates the effect of the storage temperatures on the developing stress during rewarming. For this purpose, the following storage temperatures are considered: −135 °C (3 °C below T), −150 °C (an intermediate temperature for cryogenic coolers), and −196 °C (a simple storage by immersion in liquid nitrogen). The maximum principal stress results at key points are listed in Tables 2 and 3, while the complete thermal and stress histories for a storage temperature of −196 °C are displayed in Figs. 5 and 6. In general, the maximum stress during rewarming increases with the decreasing storage temperature, where storage at −135 °C can facilitate complete stress relaxation. The latter observation is consistent with previous studies [10,63].

Table 2

Summary of simulation results for a rabbit kidney model, for the selected points illustrated in Fig. 1, where the average rewarming rate is calculated between −90 °C and −30 °C.

Property			Value
Storage Temperature			−135°C		−150°C		−196°C
Type of Rewarming			RF	Convective	RF	Convective	RF	Convective
Max. Stress (MPa)	Cooling	Pa	0.17
		Pb	0.16
		Pc	0.32
	Storage	Pa	0 (Stress Dissipation)		1.04 × 10−3		1.05 × 10−3
		Pb			1.73 × 10−3		1.74 × 10−3
		Pc			2.04 × 10−2		1.98 × 10−2
	Rewarming	Pa	4.36×10−3	7.8×10−3	2.63×10−3	1.88×10−2	4.70×10−3	2.02×10−2
		Pb	6.39×10−3	1.11×10−2	2.03×10−3	2.59×10−2	2.17×10−3	2.81×10−2
		Pc	7.07×10−3	0.18	6.34×10−2	0.46	7.58×10−2	0.48
Average Rewarming Rate (°C/min)		Pa	608	134	414	136	593	132
		Pb	157	132	141	134	160	130
		Pc	214	11.1	187	11.1	218	10.9

Table 3

Summary of simulation results for a human kidney model, for the selected points illustrated in Fig. 1, where the average rewarming rate is calculated between −90 °C and −30 °C.

Property			Value
Storage Temperature			−135°C		−150°C		−196°C
Type of Rewarming			RF	Convective	RF	Convective	RF	Convective
Max. Stress (MPa)	Cooling	Pa	0.11		0.42		0.42
		Pb	9.24 × 10−2		0.43		0.45
		Pc	6.04 × 10−2
	Storage	Pa	0 (Stress Dissipation)		−2.71 × 10−4		−2.18 × 10−4
		Pb			5.92 × 10−4		6.87 × 10−4
		Pc			0.12		0.12
	Rewarming	Pa	5.66×10−3	1.74×10−2	1.09×10−5	5.83×10−2	1.46×10−3	0.10
		Pb	1.71×10−2	6.58×10−2	2.79×10−2	0.11	6.48×10−2	0.18
		Pc	2.87×10−3	0.15	0.18	0.95	0.31	2.04
Average Rewarming Rate (°C/min)		Pa	964	152	1.15×103	150	1.78×103	159
		Pb	478	109	629	109	551	112
		Pc	495	1.77	490	1.77	500	1.79

Fig. 5.

Thermal history and the maximum principal stress history for a rabbit kidney model undergoing vitrification, at selected points illustrated in Fig. 1D, and for a storage temperature of −196 °C: (A) rewarmed with a two-step RF heating protocol; and (B) rewarmed with a two-step convective heating protocol.Thermal protocol phases: Cooling (“C”), Storage (“S”), RF rewarming (“RF”) and convective rewarming (“CR”).

Fig. 6.

Thermal history and the maximum principal stress history for a human kidney model undergoing vitrification, at selected points illustrated in Fig. 1D, and for a storage temperature of −196 °C: (A) rewarmed with a two-step RF heating protocol; and (B) rewarmed with a two-step convective heating protocol. Thermal protocol phases: Cooling (“C”), Storage (“S”), RF rewarming (“RF”) and convective rewarming (“CR”).

The significance of presenting the maximum principal stress in this study is that the vitrifying (glassy) specimen is brittle and is most vulnerable to uniaxial tensile stress, where the strength of brittle materials in compression is typically several folds higher than in tension [23]. The strength of the vitrified CPAs in tension (that is stress to fracture) is typically of the order of 3.2 MPa [55].

It can be observed from Figs. 5 and 6 that no significant stress develops in the kidney during the initial cooling portion of the protocol. This is attributed to the low viscosity of the material at higher temperatures, which allows it to flow as a liquid and prevent from stress build up. Significant stresses start to develop only several degrees above the glass transition temperature, after the viscosity has increased by 9–10 orders of magnitude from its initial value, but the specific temperature for that effect to become significant is also dependent on the cooling rate [55]. Either way, significant stress starts to develop on the outer region of the specimen, as it gets colder earlier. The stress reaches a maximum value as the temperature at the point of interest reaches the intermediate hold temperature of −135 °C, followed by stress relaxation down to negligible values.

During cooling but above T, higher stress develops in the rabbit kidney compared to the human kidney, which is counter intuitive. Furthermore, the maximum stress in the human kidney is found at the outer surface while the maximum stress in the rabbit kidney is found at the center of the specimen. The counterintuitive observation that a higher maximum thermal stress may develop in a smaller specimen is further discussed in a recent study by Solanki and Rabin [65]. In broad terms, this observation is associated with the path-dependency nature of the developed thermomechanical stress and variation in the temperature gradient history across the specimen. Either way, the stress magnitude above T under the cooling conditions studied here do not exceed 0.32 MPa and do not reach hazardous levels.

As the cooling process progresses below T, compressive stress develops at the center of the specimen, while the outer surface experiences tensile stress for both kidney sizes. The maximum tensile stress at Pb on the outer surface of the specimen is 6.03 × 10−2 MPa and 0.45 MPa for the rabbit and human kidney models respectively. As the specimen temperature equilibrates at the storage temperature, tensile stress develops at the center of the specimen while the stress at the outer surface decays to negligible values. The maximum residual tensile stress at Pc is 0.02 MPa and 0.12 MPa for the rabbit and human kidney specimens respectively. Note that without the stress relaxation stage around T, the stress level at cryogenic storage of −196 °C would easily exceed the strength of the material [13,63], where that intermediate hold temperature is crucial for cryopreservation success.

Rewarming strategy

The change in material behavior around T calls for different rewarming strategies below and above this threshold when recovering the specimen from cryogenic storage. Fig. 5A displays results from rewarming a rabbit kidney subject to a two-step RF heating, starting with a power of 2W for 3250 s (54.2 min) up to −90 °C, followed by 0.8 kW for 40 s, until the maximum temperature in the domain reaches 24.7 °C. These parameters are selected such that the initial rewarming rate would not exceed 2 °C/min and the maximum temperature difference across the specimen will not exceed 3.8 °C below −90 °C, in order to minimize stress development when there is no risk for RPC [13,63]. These parameters are also selected such that rewarming rate will exceed 100 °C/min everywhere in the specimen when it passes the temperature range of −90 °C and −30 °C, where the CPA is most susceptible to RPC but when significant stress levels are prohibited due to low viscosity values.

The specifically selected parameters yielded a maximum temperature difference of 34.9 °C within the rabbit kidney specimen at the end of RF heating, and average rewarming rates between −90 °C and −30 °C of 593 °C/min, 160 °C/min and 218 °C/min at Pa, Pb and Pc, respectively. The results displayed in Fig. 5A indicate lower maximum stress during RF heating than during cooling by 23.6%. This is very significant as previous studies suggested a much higher stress at the onset of rewarming by other means when compared with stress during cooling. The change in rewarming rate during recovery from cryogenic storage of −196 °C using two-step RF heating is displayed in Fig. 7. It can be observed from Fig. 7 that RF heating can lead to rewarming rates above the CRR in the temperature range of interest for rapid crystal formation.

Fig. 7.

Rewarming rate history for vitrified rabbit kidney and human kidney models, as they are rewarmed using a two-step RF heating from a cryogenic storage temperature of −196 °C.

Recall that the high-power rewarming in the kidney is set to start when the maximum temperature reaches −90 °C. The temperature variations in the rabbit kidney and the human kidney at this stage are 3.8 °C and 21.7 °C, respectively, which is reflected in Fig. 7 through the onset of the maximum rewarming curves below −90 °C. Additionally, the maximum and minimum cooling rates in Fig. 7 are not at a particular location in the specimen but are the maxima and minima of cooling rates at any given temperature, at a location which may change as the process progresses. As such, correlating the individual trends of the maximum or the minimum rewarming rates per se with the RF power may prove difficult. The significance of Fig. 7 is in the global range of rewarming rates, and how they relate to the threshold values to maintain vitrification.

For comparison with the RF heating cases discussed above, Fig. 5B displays the thermal and maximum principal stress histories for the rabbit kidney model when rewarmed by convective heating at the container boundaries, starting from the same cryogenic storage of −196 °C. As with the RF heating strategy, a two-step convective rewarming protocol is applied above and below T, in order to maintain low rewarming rates when the material behaves like solid. This effort resulted in a maximum temperature difference of 108 °C in the specimen during convective rewarming, but the average rewarming rate at higher temperature did not exceed 10.9 °C/min at the center of the specimen, which does not satisfy the vitrification requirements. Under these conditions, a maximum stress of 0.48 MPa was calculated at the center of the specimen, which is consistent with previous reports [13,63].

Fig. 6 displays results from a two-step RF rewarming protocol and a two-step convective rewarming protocol for the human kidney. Due to the human kidney size and the required parameters listed above, an initial RF power of 40W was applied for the first 2800 s (46.7 min), followed by a power of 50 kW for 50 s. Due to its larger size, a maximum temperature difference of 21.7 °C and 28.1 °C was calculated at the end of the two stages of RF heating, respectively. The average initial rewarming rate was again limited to 2 °C/min, by selecting a lower initial heating power. The average rewarming rate between −90 °C and 30 °C was calculated to be 1780 °C/min, 551 °C/min and 500 °C/min, at points Pa, Pb and Pc respectively.

The effect of convective heating on the thermal stress is even more prominent during convective heating of the human kidney model, Fig. 6B, with stresses exceeding 2 MPa. The average rewarming rate between −90 °C and −30 °C at the center of the specimen is 1.79 °C/min, which is significantly below the CRR required for successful vitrification.

In general, due to the underlying principles of the electromagnetic field, a high power density may develop at the tip of the antenna (Fig. 1A), which can cause the dielectric breakdown of the air [74]. In such a case, the electric field exceeds the dielectric strength of the material (i.e., the air), and the dielectric material changes from an insulator to a conductor [59]. Dielectric breakdown can have a variety of undesirable effects such as noise, excessive heating of conductors, arcing between metal surfaces, and disintegration of the antenna material. The exact value of the dielectric strength of air is dependent on various factors including electromagnetic frequency, geometry of the antenna, pressure, humidity, and temperature [45]. An approximate value of 3 × 106 V/m is considered herein for the dielectric strength of air at atmospheric pressure conditions [74].

For the human kidney case, a maximum electric field strength of the order of 104 V/m is calculated during high power rewarming within the specimen, which is below the dielectric strength of the material. While the electric field in the air surrounding the specimens is also below the dielectric breakdown threshold, it can reach that level in close vicinity to the tip of the antenna. This phenomenon has not received attention in the previous benchmark study [25], although the studied case there also leads to conditions surpassing the dielectric strength of air around the antenna.

In general, the critically high electric field around the antenna may be avoided by a careful design of its dimensions and a coating layer. For rabbit kidney case as an example, the maximum electrical field at the antenna tip can be reduced by 42%, dropping below the dielectric strength of air, when the antenna radius is increased from 1.5 mm to 3 mm, the transmission line radius is increased from 3.45 mm to 6.9 mm, and the antenna penetration depth is reduced from 60 mm to 51 mm. This alternate design would not affect the temperature field and its history within the specimen, while avoiding the dielectric breakdown of the air.

The human kidney case would require a more involved design to overcome the dielectric breakdown challenge, with the addition of a Teflon (Polytetrafluoroethylene) cap as another example. Here the antenna radius is increased from 1.5 mm to 8 mm, its penetration length is decreased from 240 mm to 166 mm, the transmission line radius is increased from 3.45 mm to 18.4 mm, and a Teflon cap of 8 mm thickness and 50 mm length covers the antenna tip. The Teflon cap covering the tip of the antenna benefits from low relative permittivity, low dielectric loss, high dielectric strength and high melting temperature making it ideal for this purpose [70]. The use of a Teflon cap covering the antenna to avoid corrosion of the antenna due to high electric fields is demonstrated by Gao and Bilén [22]. This alternative antenna design would reduce the maximum field in the free space around the tip of the antenna by 92%, while unaffecting the thermal and thermomechanical stress histories in the human kidney.

The above alternative antenna designs serve as examples only, but do not represent the optimum designs with respect to maximizing heating efficiency, which is a challenge left for future studies.

Rewarming and RF heat generation

A notable observation is that while the human kidney requires a much higher RF power input, the maximum temperature difference at the end of RF heating, when the warmest point in the specimen reaches 25 °C (±1 °C), is larger in the rabbit kidney case due to the underlying principles of RF heating. In this context, the stability factor K [15,49] has been described as a convenient parameter for evaluating RF heating as a rewarming technique.

The stability factor K is a material property derived from the dielectric properties of the material, which can be used to estimate the extent of thermal runaway [15,49]. In broad terms, a thermal runaway occurs when an increase in temperature changes the local conditions in a way that promotes further heating, which may result in a destructive outcome to the specimen [7,37]. While thermal runaway is an inherent risk with RF heating, it can be mitigated by designing a setup and a protocol of operation such that more power would be absorbed in the colder region of the specimen. Such favorable conditions can be achieved when K decreases with increasing temperature (negative dK/dT) as discussed previously [15,49]. For example, RF heating of 7.05 M DMSO at 434 MHz displays favorable condition to prevent thermal runaway in the range of −40 °C to 30 °C while adverse condition is observed in the temperature range of −196 °C to −40 °C [15]. For the same material at 108 MHz, the cutoff between favorable and adverse conditions shits from −40 °C to −55 °C.

The implication of the stability factor K on rewarming in the first stage of the proposed cooling protocol in this study (from storage to 90 °C) is that, while the specimen is prone to a thermal runaway for both frequencies and, hence, kidney models (108 MHz and 434 MHz for the human and rabbit kidney models, respectively), the larger kidney model requires a higher power, which results in a larger maximum temperature difference across the specimen: 21.7 °C for the human kidney versus 3.8 °C for the rabbit kidney when they each get to −90 °C. However, in the second stage of the RF rewarming protocol, when the stability factor K changes trend earlier for the lower frequency, the maximum temperature difference for the smaller specimen of rabbit kidney is found higher than that of the human kidney: 34.9 °C as opposed to 28.1 °C, respectively, thus supporting the above observation.

While the effects of CPA dielectric properties on the thermal and stress analyses are clear from the above discussion, they are generally not available in the entire temperature range of interest for the current study. As presented in the material properties section above, the dielectric properties in the current study are taken as temperature-dependent only within the −80 °C to 20 °C range displayed in Fig. 2. Outside of that range and for the purpose of the current study, those properties are assumed constant at their respective range boundary values. It is within the above temperature range that the material is most susceptible to significantly rapid crystal growth rate, where rapid and uniform rewarming are desired. Crystallization and thermal stress are not feasible above that temperature range and the assumption of constant dielectric properties at high temperatures has no implications on the objectives of the current study. Below that temperature range on the other hand, the assumption of constant dielectric properties may lead to lower temperature gradients and thermal stresses. This means that extrapolating those decaying properties to lower temperatures would result in higher temperature gradients and thermal stresses.

More specifically to lower temperatures, it appears that the real component of relative permittivity approaches a constant value with decreasing temperature near −80 °C (Fig. 2A), while its imaginary component significantly decays with the decreasing temperatures below that temperature (Fig. 2B), based on the available data. Such decreasing of ε” with the decreasing temperatures, could result in a positive gradient of stability factor K [15,49]. In turn, this means that the risk to a thermal runaway is lower with the currently assumed constant properties, compared to the anticipated thermal runaway with extrapolation of the available data to lower temperatures. However, accurate estimation of thermal runaway will only be possible with expansion of the material properties database. Measurement of data for improving the accuracy of the thermomechanical stress analysis is therefore a subject for future studies.

An additional effect at higher temperatures is the contribution of ionic components of the CPA solution to power absorption. While the current study considers 7.05 M DMSO as a reference solution [28,42, 50], the actual selection of CPA cocktail is governed by various biophysical aspects, possibly affecting the effect of the ionic component [19]. Such components may increase the dielectric loss factor (ε”) in the specimens at higher temperatures (say, above −25 °C), while having no effect on the relative permittivity (ε’) [15]. The contribution of the ionic component to the dielectric loss factor is more significant in frequencies lower than 100 MHz, an effect which decreases with the increasing frequency [15]. A higher value of the dielectric loss at higher temperatures would result in a larger power dissipation in the specimen, while the positive gradient of stability factor K results in a larger thermal runaway. It follows that the current study presents conservative results in the sense that the ionic components in the solution would only decrease the needed electromagnetic power required to rewarm the specimen at the reported rate. Alternatively, using the calculated electromagnetic power would result in higher temperature gradients and rewarming rates.

Recall that the goal in the current study is to present a computation framework for calculating thermomechanical stresses in a specimen undergoing an RF-assisted vitrification protocol, as well as to provide insight on the thermal process. While meeting that goal, it might seem that an unrealistically large cavity size and input RF power is required for the human kidney case. However, practical approaches to cavity design do exist, which can affect the RF heating parameters and bring its size to practical dimensions, including altering the cavity shape [16,35], addition of an antenna [49], and modifications to the specimen container geometry. Note that those changes would also affect the resulting thermomechanical stress, and the coupled process must be optimized as a whole. Nonetheless, the computational framework presented in this study could be used as tool in the process of RF heating optimization, which is the subject matter of ongoing studies.

Summary and conclusions

RF heating has been presented as an effective means for uniformly rewarming a cryopreserved specimen, with the ultimate goal of reducing the thermomechanical stress during rewarming. This study proposes a computation framework to solve the coupled fields of heat transfer, solid mechanics, and electromagnetic radiation, with the goal of providing means to avoid thermomechanical stress damage.

Comparison of the proposed FEA framework with a benchmark study, using a proprietary numerical solution to solve a simplified spherical problem, suggest good agreement in lower cryogenic temperatures relevant to cryopreservation by vitrification.

Next, this study moves to a more realistic scenario of a thermo-mechanical stress problem in rabbit kidney and human kidney models, in cylindrical CPA containers as they undergo cryopreservation by vitrification. While the general problem has been studied previously, the unique contribution of the current investigation is in integrating thermomechanical stress analysis in the kidney as it undergoes vitrification and subsequent RF heating. This study includes a two-phase cooling protocol, a two-stage RF rewarming, intermediate stages for stress relaxation, and a parametric study for the storage temperature.

Results suggest that it is feasible to rewarm the specimen from cryogenic storage with favorably high rewarming rates, while maintaining unharmful levels of stress, provided that special attention is given to RF chamber design, optimization of the rewarming protocol, and specimen container design. The developed knowledge from this study can serve as the basis for optimization of cryopreservation protocols, with the goal of reducing thermomechanical stress and preserving structural integrity.

Inductive heating using magnetic nanoparticles with CPA has been recently presented as a potential approach for uniform rewarming of the specimen from cryogenic storage [38]. Notably, the computation framework presented in this study can be straightforwardly applied to thermomechanical analysis of vitrification involving nanowarming, provided that the dielectric properties of CPA-nanoparticle solutions are measured.