Numerical Study of Heat and Mass Transfer during Cryopreservation Process with Application of Directed Interval Arithmetic

In the present paper, numerical modelling of heat and mass transfer proceeding in a two-dimensional axially symmetrical articular cartilage sample subjected to a cryopreservation process is presented. In the model under consideration, interval parameters were assumed. The heat transfer process is described using the Fourier interval equation, while the cryoprotectant transport (DMSO) across the cell membrane is analyzed using a two-parameter model taking into account the simulation of the water volume in the chondrocytes and the change in DMSO concentration over time. The liquidus tracking (LT) protocol introduced by Pegg et al. was used to model the cryopreservation process. This procedure divides the heating and cooling phases into eight and seven steps, respectively, allowing precise regulation of temperature and cryoprotectant (CPA) concentration of bathing solutions. This protocol protects chondrocytes from ice crystal, osmotic stress, and electrolyte damage. The obtained interval concentrations of cryoprotectant in chondrocytes were compared with previous simulations obtained using the deterministic model and they are mostly in agreement with the simulation data.


Introduction
Cryopreservation is the storage of cells, tissues, or other biological structure at low temperatures without injuring them. This approach has many practical applications, including usage in medicine. The process is applied to preserve stem cells (tissue engineering research) or to cryobank transported organs (in transplantology) [1][2][3].
Cryopreservation can be accomplished via two methods, slow freezing and vitrification. These techniques differ principally in cooling rate and in cryoprotectant (CPA) concentration. The slow freezing method is carried out at a rate of approx. −1 • C/min with a low CPA concentration of about 1-2 M. This procedure has a low risk of tissue damage from osmotic shock to cells or the chemical toxic effects of CPA. The disadvantage is the possibility of tissue injury caused by ice crystal formation. On the other hand, vitrification is the conversion from a liquid to a vitreous phase at a cooling rate of approx. −100 • C/min to avoid the ice crystal formation. This method uses a high concentration of CPA (about 4-8 M), which increases the risk of tissue damage through CPA toxicity and the possibility of osmotic shock to cells [2,3].
However, cells such as chondrocytes in articular cartilage are susceptible to injury due to the forming of ice crystals, which makes cryopreservation using conventional methods difficult [4]. To solve this problem, Pegg et al. [5] and Wang et al. [6] used another procedure called the liquidus tracking (LT) method. This technique was devised by Farrant (1965) [7] and further advanced by Elford et al. (1972) [8]. This approach precisely regulates the temperature and the CPA concentration of bathing solutions. The temperatures in the

Mathematical Model
The change in water volume in the chondrocytes and the change in the intracellular CPA mole number over time is simulated using an interval two-parameter model [13,14,40], which, for a cylindrical coordinate system, is of the form: dV w (r, z, t) dt = −L p (r, z, t)AR g T(r, z, t) M e (r, z, t) − M i (r, z, t) , where V w (r, z, t) is the interval water volume in the chondrocyte, t is the time, r and z denote the cylindrical coordinates, N d is the CPA mole number in the chondrocyte, L p (r, z, t) is the interval hydraulic conductivity, P s (r, z, t) is the interval permeability of CPA through the cell membrane, A is the chondrocyte surface area, R g is the gas constant equal to 8.314 J·mol −1 ·K −1 , T(r, z, t) is the interval tissue temperature, and M(r, z, t) is the interval osmolarity, where the superscripts e and i refer to the extracellular and intracellular solutions, respectively. The interval hydraulic conductivity, which determines the water permeability of a cell membrane, and the interval permeability of CPA across the cell membrane are described by the following equations [13]: where A Lp and A Ps are pre-exponential factors, and E A,Lp and E A,Ps are activation energies for L p and P s , respectively. The osmolarity, which defines the number of moles of osmotically active substance in 1 L of solution, can be calculated on the basis of the osmolality of the solution. The osmolality is the number of moles of osmotically active substances dissolved in 1 kg of a solvent (e.g., water). For an undiluted solution of two solutes, the interval value of the osmolality can be expressed as [13,25]: where π f is the interval osmolality, k diss is the dissociation constant, m f is the interval molality, and B is the second osmotic virial coefficient, where the subscripts k and d represent the KCl and the CPA (exactly DMSO), respectively. Let us denote that the superscript f refers to the superscript e or i (extracellular or intracellular solution).
The molality of species j (e.g., k or d) is defined as the number of moles of the dissolved substance in 1 kg of solvent. When the interval mass fraction is known (the mass fraction is the ratio of the mass of a given component to the mass of the entire mixture), the interval molality can be determined: where ω f j is the interval mass fraction of species j and M at.j is the molar mass of the species j. The molar mass corresponds to the molecular mass: where M u,j is the molecular mass of species j, which is given by the sum of the relative atomic mass of the elements. The interval mass fraction can be calculated by the dependence [13]: where C f j is the interval concentration of species j and ρ j is the density of species j. It should be noted that concentration is also called molarity and is defined as the mole number per liter of the solution.
After defining the interval osmolality (see Equation (5)), the interval osmolarity can be determined [13]: where V f j is the interval volume per unit of mass of species j. It should be noted that subscript w represents H 2 O.
The interval volume per unit mass is calculated [13]: where v j is the partial molar volume of species j. It is also important to know the interval extracellular and intracellular concentrations, which are described by the following interval equations [13][14][15]18,22,27,40,41]: where D d is the interval diffusion coefficient, V 0 cell is the cell volume (the superscript 0 represents the initial moment of simulation), V b is the osmotically inactive volume of cells (for chondrocytes it is equal to 0.41 V 0 ), where V 0 is the isotonic volume of cells. The interval normalized cell volume is defined: The interval diffusion coefficient of CPA in extracellular matrix can be estimated by the Einstein-Stokes equation [41][42][43]: where k B is the Boltzmann constant, equal to 1.38 × 10 −23 J·K −1 , r s is the radius of the spherical particle molecule, and η is the dynamic viscosity. The interval temperature coupled to the interval mass transfer model (compare with Equations (1), (2) and (15)) is determined using the interval Fourier equation [13,41]: where λ is the interval thermal conductivity and c is the interval volumetric specific heat. Equations (11)- (13) and (16) need to be supplemented with boundary conditions (see Figure 1) and initial conditions of the following form: where T bulk is the temperature of the bathing solution, C bulk is the CPA concentration in the bathing solution, α is the natural convection heat transfer coefficient, T 0 is the initial temperature and C e,0 d the initial external DMSO concertation, C e,0 k is the initial external KCl concentration, and C i,0 is the initial internal DMSO and KCl concentration [13]. where kB is the Boltzmann constant, equal to 1.38 × 10 −23 J·K −1 , rs is the radius of the spherical particle molecule, and η is the dynamic viscosity. The interval temperature coupled to the interval mass transfer model (compare with Equations (1), (2) and (15)) is determined using the interval Fourier equation [13,41]: where  is the interval thermal conductivity and c is the interval volumetric specific heat.
Equations (11)- (13) and (16) need to be supplemented with boundary conditions (see Figure 1) and initial conditions of the following form: where bulk T is the temperature of the bathing solution, bulk C is the CPA concentration in the bathing solution, α is the natural convection heat transfer coefficient, 0 T is the initial temperature and  An axially symmetrical model of articular cartilage was considered, in particular the rectangle shown in Figure 1. Adiabatic boundary conditions were assumed on the symmetry axes (left and bottom sides), while on the other sides of the considered rectangle boundary conditions of 3rd type for temperature and 1st type for concentration were given. The presented mathematical model does not include the phenomenon of phase changes. This is due to the fact that the liquidus tracking (LT) method is used to model heat and mass transfer. The LT protocol regulates the temperature and concentration in such a way that the temperature of the sample is above or on the liquidus line, which eliminates the probability of ice crystallization in cells-see the calculated eutectic temperatures and the melting points of tissue in [5]. An axially symmetrical model of articular cartilage was considered, in particular the rectangle shown in Figure 1. Adiabatic boundary conditions were assumed on the symmetry axes (left and bottom sides), while on the other sides of the considered rectangle boundary conditions of 3rd type for temperature and 1st type for concentration were given. The presented mathematical model does not include the phenomenon of phase changes. This is due to the fact that the liquidus tracking (LT) method is used to model heat and mass transfer. The LT protocol regulates the temperature and concentration in such a way that the temperature of the sample is above or on the liquidus line, which eliminates the probability of ice crystallization in cells-see the calculated eutectic temperatures and the melting points of tissue in [5].

Numerical Algorithm
Numerical model of thermal processes proceeding in domain of heating tissue is based on the finite difference method (FDM) in the version presented in [43,44].
A time grid with a constant step ∆t and a geometrical mesh are introduced at the beginning. The boundary nodes are located at the distance 0.5 h or 0.5 k with respect to the real boundary (h, k are the steps of regular mesh in directions r and z), respectively (see Figure 2). This approach gives a better approximation of the Neumann and Robin boundary conditions [38,41,43].

Numerical Algorithm
Numerical model of thermal processes proceeding in domain of heating tissue is based on the finite difference method (FDM) in the version presented in [43,44].
A time grid with a constant step t  and a geometrical mesh are introduced at the beginning. The boundary nodes are located at the distance 0.5 h or 0.5 k with respect to the real boundary (h, k are the steps of regular mesh in directions r and z), respectively (see Figure 2). This approach gives a better approximation of the Neumann and Robin boundary conditions [38,41,43]. The approximate form of the interval energy equation (Equation (16)) for the internal nodes (i, j) and transition t s−1 → t s is the following [41,43]: where , , , are the shape functions of differential mesh, while the thermal resistances are defined by the following formulas: The interval mass diffusion equation (Equation (11)) is transformed in a similar way The approximate form of the interval energy equation (Equation (16)) for the internal nodes (i, j) and transition t s−1 → t s is the following [41,43]: where are the shape functions of differential mesh, while the thermal resistances are defined by the following formulas: The interval mass diffusion equation (Equation (11)) is transformed in a similar way where the diffusion resistances between the central node and the adjoining ones are the following: Details of the above transformations are presented in [41]. The final form of the interval FDM equations (Equations (1) and (2)) of the twoparameter model for the internal nodes are The systems of Equations (18), (22), (25), and (26) have been solved using the assumption of the stability condition for the explicit differential scheme [43,44].
All mathematical operations leading to the determination of the temperatures, extracellular and intracellular concentrations and other sought quantities require the application of directed arithmetic rules-more information can be found in [41].

Results
The simulation of CPA transport across a cell membrane for the homogenous cylindrical articular cartilage sample was carried out. The dimensions of the sample were as follows: H = 1 mm and R = 3 mm. The thermophysical parameters of the tissue and the chemical properties of CPA, including DMSO, KCl and H 2 O, are given in Table 1. The chondrocytes' surface area and the chondrocytes' isotonic radius have been introduced, where A = 8.04 × 10 −10 m 2 and r cell = 8 × 10 −6 m, which allows to calculate the isotonic volume of the cells V 0 = 2.14 × 10 −15 m 3 . It is assumed that the initial cell volume is equal to the isotonic volume of the cells V 0 cell = V 0 . Additionally, the following input data were used: initial temperature T 0 = 22 • C, initial external DMSO concentration C e,0 d = 0% (w/w), initial external KCl concentration C e,0 k = 0.85% (w/w), initial internal DMSO and KCl concentrations C i,0 = 0% (w/w), pre-expositional factors, such as A Lp = 9.1 × 10 −6 m 2 ·s·kg −1 and A Ps = 1.2 × 10 12 m·s −1 and the activation energy, E A,Lp = 45.73 × 10 3 J mol −1 and E A,Ps = 107.40 × 10 3 J mol −1 [13,41]. Table 1. Thermophysical parameters of the tissue and chemical properties of the CPA [13,41,45].

Thermophysical Tissue Parameters
Parameter Value The process was simulated using the interval finite difference method. The time step and the mesh steps were assumed as ∆t = 0.001 s, h = 0.0001 m and k = 0.00005 m.
During modelling of the cryopreservation process, the LT protocol introduced by Pegg et al. in [5] was used. The procedure includes 8 steps in a cooling phase and 7 steps in a heating phase. During each step, the temperature and concentration of the bathing solution were controlled by computer. This ensured that no ice crystals were formed and that no cells were damaged by toxicity. It is worth mentioning that the impact of CPA causes a change in the melting point, and then the temperature of the cryopreserved sample was close to the liquidus track. Table 2 and Figure 3 present the assumption of the LT approach. The given values of the temperature and concentration of the bathing solution corresponded to the boundary conditions of the mathematical model. Table 2. The assumption of the LT protocol [5,13].

Phase
Step in a heating phase. During each step, the temperature and concentration of the bathing solution were controlled by computer. This ensured that no ice crystals were formed and that no cells were damaged by toxicity. It is worth mentioning that the impact of CPA causes a change in the melting point, and then the temperature of the cryopreserved sample was close to the liquidus track.   Firstly, simulation of the CPA concentration in the intracellular membrane using directed interval arithmetic was executed. Figure 4 illustrates the distribution of the concentration in chondrocytes located at the point with the coordinates r = 0.1 mm and z = 0.45 mm. This refers to the first 20 s in step 3 of the cooling phase ( Figure 4a) and step 7 of heating phase (Figure 4b). It should be pointed out that, for each node of the domain considered, there are two curves representing the beginning and end of concentration intervals. The upper limit of the concentration interval is marked with a red line and the lower limit with a blue line. Figure 5 represents the CPA concentration in the entire analyzed area after 10 s in step 3 in the cooling phase. The assumption of the LT protocol [5,13].

Phase
Step Time Firstly, simulation of the CPA concentration in the intracellular membrane using directed interval arithmetic was executed. Figure 4 illustrates the distribution of the concentration in chondrocytes located at the point with the coordinates r = 0.1 mm and z = 0.45 mm. This refers to the first 20 s in step 3 of the cooling phase ( Figure 4a) and step 7 of heating phase (Figure 4b). It should be pointed out that, for each node of the domain considered, there are two curves representing the beginning and end of concentration intervals. The upper limit of the concentration interval is marked with a red line and the lower limit with a blue line. Figure 5 represents the CPA concentration in the entire analyzed area after 10 s in step 3 in the cooling phase.   Afterwards, to better understand the phenomena occurring in chondrocytes, the simulation of CPA and water volume changing in the cell was performed without using interval arithmetic. In this model, the thermophysical parameters were defined as deterministic values (see deterministic values of λ and c in Table 1). Figure 6 shows a diagram of the number of moles of CPA over time during the whole process. Figure 7 presents changes in the intracellular water volume during cryopreservation. This diagram expresses the phenomenon of dehydration, which consists of water transport from the inside cell to the extracellular matrix. Afterwards, to better understand the phenomena occurring in chondrocytes, the simulation of CPA and water volume changing in the cell was performed without using interval arithmetic. In this model, the thermophysical parameters were defined as deterministic values (see deterministic values of λ and c in Table 1). Figure 6 shows a diagram of the number of moles of CPA over time during the whole process. Figure 7 presents changes in the intracellular water volume during cryopreservation. This diagram expresses the phenomenon of dehydration, which consists of water transport from the inside cell to the extracellular matrix.  Afterwards, to better understand the phenomena occurring in chondro ulation of CPA and water volume changing in the cell was performed wi terval arithmetic. In this model, the thermophysical parameters were define istic values (see deterministic values of λ and c in Table 1). Figure 6 show the number of moles of CPA over time during the whole process. Fig  changes in the intracellular water volume during cryopreservation. Thi presses the phenomenon of dehydration, which consists of water transport cell to the extracellular matrix.  Then, the intracellular volume of water and CPA was also analyzed using interval arithmetic rules. Figure 8 demonstrates the CPA moles' number in the cell for the first 20 s in step 3 of the cooling phase and step 7 of the heating phase. These diagrams correspond to the intracellular CPA concertation (compare with Figure 4). Figure 9 shows a diagram of water volume over time for transition from 22 • C to −5 • C. It should be noted that the obtained intervals are narrow. Then, the intracellular volume of water and CPA was also analyzed arithmetic rules. Figure 8 demonstrates the CPA moles' number in the cell s in step 3 of the cooling phase and step 7 of the heating phase. These diagra to the intracellular CPA concertation (compare with Figure 4).  Then, the intracellular volume of water and CPA was also analyzed using interval arithmetic rules. Figure 8 demonstrates the CPA moles' number in the cell for the first 20 s in step 3 of the cooling phase and step 7 of the heating phase. These diagrams correspond to the intracellular CPA concertation (compare with Figure 4). Figure 9 shows a diagram of water volume over time for transition from 22 C to −5 C. It should be noted that the obtained intervals are narrow.   Table 3 includes the values of intracellular quantities at the end of th The exact results of the extracellular temperature and CPA concentration c previous work [41]. The obtained interval concentrations were also comp simulation data received by Yu et al. [13]. In individual steps, these valu each other. Significant discrepancies occur for steps where the temperature Table 3. Values of intracellular quantities at the end of each step.

Phase
Step   Table 3 includes the values of intracellular quantities at the end of the given steps. The exact results of the extracellular temperature and CPA concentration can be found in previous work [41]. The obtained interval concentrations were also compared with the simulation data received by Yu et al. [13]. In individual steps, these values are close to each other. Significant discrepancies occur for steps where the temperature is low.

Discussion
In the present study, numerical analysis of heat and mass transfer occurring in axially symmetrical articular cartilage sample subjected to the cryopreservation process is presented. The mathematical model is based on heat and mass diffusion interval equations, in conjunction with a model of transmembrane transport of CPA. A two-parameter model was used to analyze the transport of CPA (solution of DMSO, KCl and water) across the cell membrane, considering the simulated volume of water in chondrocytes and the change in the cryoprotectant concentration over time. Interval values of volume specific heat and thermal conductivity of articular cartilage were assumed for the modelled process. The interval FDM was used to solve the problem, applying the rules of directed interval arithmetic.
Since thermophysical parameters are often determined experimentally, as well as the phenomena that biological structures are unpredictable, it is reasonable to include in the mathematical model imprecise values instead of deterministic values. The application of the stochastic and randomized model is time-consuming and complicated, therefore applying interval arithmetic rules is proposed in the research. Using interval values, solutions in the form of intervals were obtained, which better reflect the analyzed process.
In the paper, the history of the intracellular concentration, the history of the intracellular water volume in the chondrocytes and the history of the intracellular CPA moles number in the cell have been presented. The results depict the phenomena that occur in cells during the process, such as dehydration. It can be also seen that, for low temperatures, wider intervals were obtained, while for positive temperatures the intervals were narrow. The obtained interval concentrations were compared with simulation data calculated using a deterministic model [13]. The interval results mostly coincide with Yu et al.'s simulation data. Slight differences are observed in those process steps in which the concentration is high (e.g., step 8 in the cooling phase or step 1 in the heating phase). The comparison of the deterministic model with the interval model allows us to conclude that the model modified by interval arithmetic is justified because the results are compatible, and the received ranges include the correct results, whereas the computation time is only slightly extended.
In a further stage of the study, the macroscale mass transfer equation will be supplemented with a velocity vector determined from the Navier-Stokes equation.