CRYOPRESERVATION ANALYSIS CONSIDERING DEGREE OF CRYSTALLISATION USING FUZZY ARITHMETIC 1

This article presents numerical modelling of the heat transfer process in a sample during cryopreservation by vitriﬁcation in a microﬂuidic system. Single-phase ﬂow of the working ﬂuid in the microchannels during warming was considered, while two-phase ﬂow during cooling. The mathematical model is based on the Fourier equation with a source term that takes into account the degree of ice crystallisation. Fuzzy thermophysical parameters were assumed in the model. The problem was solved by the ﬁnite diﬀerence method and the fourth-order Runge-Kutta algorithm, using the concept of α -cuts. The results of numerical simulation were compared with the results from the literature.


Introduction
During modelling physical phenomena occurring in biological samples, uncertain parameters are often used.These variables are imprecise because they are determined experimentally and depend on the age, sex, and condition of the examined organism.As a consequence, physical processes in biological samples or in engineering systems are often simulated with deterministic models that introduce at the same time some assumptions and simplifications (Wang and Matthies, 2021).
On the other hand, uncertain quantities present in physical phenomena can be predicted by probabilistic or non-probabilistic approaches.Stochastic methods based on probabilistic algorithms are successfully suitable for describing uncertainties with a known database containing measurements of a given parameter.From these input data, a probability distribution function can be prepared to describe uncertainty characteristics.Unfortunately, for some engineering problems, probabilistic techniques are ineffective due to limited access to measurement data (Wang and Matthies, 2021).
Alternative approaches to the modelling of uncertain quantities are fuzzy set theory and interval set theory, which are qualified as non-probabilistic methods (Lü et al., 2017).The concept of fuzzy sets was suggested by Lofti and Zadeh in 1965 (Zadeh, 1965).According to this theory, a membership function is defined for each element of the set, which takes values in the range from 0 to 1.The membership function determines whether a given element belongs to the set completely, partially, or is external to it (Hanss, 2005;Skorupa, 2023).There are different types of membership functions, where the simplest include linear functions, such as triangular and trapezoidal ones, and the most complex functions are Gaussian or bell curves (Caniani et al., 2011).It is worth mentioning that fuzzy numbers are often performed using α-cuts concept (Giachetti and Young, 1997;Guerra and Stefanini, 2005).
An other effective method is the interval set theory introduced by Ramon and Moore in 1966 (Moore, 1966).In that case, uncertain variables are described by a set for which the upper and lower bounds are specified.The membership function is equal to 1 when the element belongs to the set, or 0 when it is not in the given interval (Lü et al., 2017;Skorupa, 2023).
The subject of this paper focusses on the use of fuzzy set theory to model phenomena that occur during cryopreservation.This is a process that involves reducing biological activity of a biological material and then storing it at a low temperature.The essence of cryopreservation is to prevent samples from damage after restoring them to physiological temperature (Jang et al., 2017;Zhao and Fu, 2017).
One of the phenomena present during cryopreservation is crystallisation, which involves the formation of ice crystals from the water contained in the biological sample.The crystallisation process is initiated by nucleation caused by the metastable homogeneous phase.Then around the nucleus, growth of ice crystals occurs (Tan et al., 2021).
The formation of ice crystals can damage biological samples.To prevent these problems, during cryopreservation, chemical compounds called cryoprotectants (CPAs) are used, and the cooling rate is regulated (Skorupa, 2023).On the basis of this, different cryopreservation methods can be distinguished.Vitrification, for example, involves overcooling the sample at a high cooling rate, thus vitrifying the solution without ice crystal formation (Jang et al., 2017).
To predict the potential damage to the cryopreserved sample, the cryopreservation process is modelled mathematically and numerically.For this, it is necessary to consider various transport phenomena, such as heat and mass transfer or osmotic transport (Skorupa, 2023).The governing equation for estimating the thermal distribution is the Fourier equation (Fourier, 1882).Furthermore, this relationship is coupled to the degree of crystallisation determined, for example, from the non-isothermal kinetic equation proposed by Boutron and Mehl (1990).For preparing a macroscopic model of the liquid solidification process, different approaches are applied, for example, the uncoupled method, the Stefan model (sharp interface method) and the zone model (Skorupa, 2023;Song et al., 2010;Zhou et al., 2013).
In the literature, it is possible to find examples of mathematical calculations of thermal processes during freezing coupled to the degree of crystallisation with deterministic models (Shi et al.,  The interesting position is the article (Piasecka-Belkhayat and Skorupa, 2023) published in 2023, which is devoted to modelling thermal processes, including the degree of crystallisation, using the interval set theory.Further research on applying non-probabilistic methods, including interval numbers and triangular fuzzy numbers, can be found in the thesis (Skorupa, 2023).This article presents the modelling of bioheat transfer and degree of crystallisation during cryopreservation in a microfluidic device with imprecise parameters.The uncertainties are considered by non-probabilistic methods, more specifically, algorithms for fuzzy numbers defined by a triangular and a trapezoidal membership function.The results obtained are compared with data from the literature (Piasecka-Belkhayat and Skorupa, 2023; Skorupa, 2023; Zhou et al., 2013).The problem is specified mathematically by the Fourier equation and the non-isothermal kinetic equation, while simulations are conducted with the finite difference method (FDM) and the Runge-Kutta algorithm.

Governing equations
The paper analyses the task of one-dimensional heat transfer, taking into account that the cell suspension is maintained as a thin layer on the chip.Figure 1 shows the microfluidic system model based on the concept presented in (Tuckerman and Pease, 1981;Zhou et al., 2013).The diagram also includes the marked nodes A, B and C where the corresponding boundary conditions are given.The energy equation describing the temperature distribution within a one-dimensional microfluidic chip can be formulated as the fuzzy Fourier equation (Fourier, 1882) where T is the fuzzy temperature, λ is the fuzzy thermal conductivity, c is the fuzzy specific heat, S h is the fuzzy heat source term, while ρ is the density.The fuzzy source component S h for the sample layer is expressed by the following relation (Shi et al., 2018) where χ is the fuzzy degree of ice crystallisation belonging to the interval ( 0, 1), ρ h is the density of water, and L h is the latent heat of water.In contrast, for the chip wall, this component is equal to zero.
In the above equation, values for water were implemented instead of the thermophysical parameters of the biological sample due to the fact that the thermal properties of biological cells are similar to those of their prevalent component, namely water (Shi et al., 2018).
The crystallisation process is described by the non-isothermal Mehl-Boutron kinetic equation (Boutron and Mehl, 1990) where ∂ χ/∂t means the growth rate of χ, while k a is the characteristic coefficient depending on the solution composition, T m is the freezing (melting) temperature, Q is the activation energy and R is the gas constant (R = 8.314 J mol −1 K −1 ).
The boundary conditions must be attached to the mathematical model defined in this way.Due to the use of thermal insulation (see Fig. 1), the heat flux between the devices and the environment is neglected.
Node A, shown in Fig. 1, lies at the boundary between the model and the working fluid, which is the main source of thermal changes.The fuzzy heat flux q at this node is described by the boundary condition of the 3rd type extended by the microchannel geometry (Zhou et al., 2013) where W f , W w and H f are the microchannel dimensions (see Fig. 1), T f is the temperature of the working fluid, α Γ is the external heat transfer coefficient, η is the fuzzy fin efficiency and subscripts w and f denote the chip wall and the working fluid, respectively.The fuzzy heat flux is defined as follows (Mochnacki and Suchy, 1993) where n is the normal vector.
The fuzzy fin efficiency is determined from the relationship (Zhou et al., 2013) where m is the fuzzy fin parameter (Zhou et al., 2013) where λ w is the fuzzy thermal conductivity of the chip wall.
Let us now turn to node B which is located at the interface between the chip wall and the sample layer.At this node, a boundary condition of the 4th type is assumed with ideal contact (Mochnacki and Suchy, 1993) where subscript s denotes the sample layer domain.
On the other hand, an adiabatic condition was assumed at node C due to symmetry of the system under consideration (Mochnacki and Suchy, 1993 To complete the mathematical description, it is necessary to take into account the initial conditions in which the temperature and the degree of crystallisation in the sample domain were determined at time t = 0 (Zhou et al., 2013) where T 0 and χ 0 are the initial values of temperature and degree of crystallisation, respectively.

Numerical model
The fuzzy finite difference method, which is a very good tool for solving nonlinear equations, was used to solve an unsteady heat transfer problem defined in this way, Eq. (2.1).The non-linearity present in Eq. (2.1) is related to both the varying values of thermophysical parameters c( T ), λ( T ) and the fuzzy source term S h ( T ).The numerical model of thermal processes occurring in the microfluidic system is based on the finite difference method as presented in (Mochnacki and Suchy, 1993) supplemented by the rules of fuzzy arithmetics, in particular, the concept of α-cuts dedicated for triangular and trapezoidal fuzzy numbers (Piasecka-Belkhayat and Korczak, 2020; Skorupa, 2023).
According to the fuzzy number theory, any fuzzy number can be written as the sum of all its α-cuts (Piasecka-Belkhayat and Korczak, 2020; Skorupa, 2023) where α-cuts are expressed as a set of closed intervals Considering the case of a triangular fuzzy number α = (a − , a 0 , a + ), we can define the α-cut as a set of closed intervals of the following form (Skorupa, 2023) where a 0 is the core of the number, while the values of a − and a + denote the left and right ends of the fuzzy number respectively; while in the case of the trapezoidal fuzzy number a = (x 0 , y 0 , σ, β), the α-cut is represented as follows (Piasecka-Belkhayat and Korczak, 2020) where x 0 and y 0 are the defuzzifiers from the left and right side, respectively; σ and β are the left and right fuzzinesses.In this way, complicated arithmetic operations on fuzzy numbers are avoided by performing calculations on closed intervals, which are interval numbers, using the rules of interval arithmetics (Skorupa, 2023).The calculations carried out involved different values of the parameter α.
Using the finite difference method, a time grid with a constant step ∆t and a geometric grid with a constant step h are introduced at the beginning.Figure 2 shows the three-point star idea, which is used to create the geometric mesh.The model assumes weak non-linearity of the specific heat.Therefore, the differential equation for internal nodes corresponding to the fuzzy Fourier equation written in the implicit scheme has the following form (Mochnacki and Suchy, 1993) or it can be written as where coefficients A i , B i and C i are calculated from the relationships while superscript f means a moment of time, the time step is where n el is the number of nodes, and where k denotes the node number (k = i − 1, i, i + 1).
The obtained system of fuzzy equations supplemented with boundary-initial conditions can be solved using the Thomas method (Skorupa, 2023).It is important to note that the advantage of the implicit scheme is its stability and the lack of restrictions on allowable values of the time step (Mochnacki and Suchy, 1993).
Let us now turn to linearisation of the fuzzy source component S h ( χ f i ), Eq. (2.2).The values of the fuzzy degree of ice crystallisation χ and ice growth rate χ ′ were calculated numerically using the fourth-order Runge-Kutta algorithm due to difficulty in calculating these relationships analytically (Piasecka-Belkhayat and Skorupa, 2023; Zhou et al., 2013) where The temperature dependence of thermal conductivity and specific heat of both the chip wall (made by silicon) and the sample layer (solution of ethylene glycol -EG and water) was assumed in the numerical model.These parameters were calculated as temperature-dependent polynomial functions using a linear regression method.In the case of silicon, the polynomial functions are of the form It should be noted that in the case of silicon, the measurement points taken from the literature (Desai, 1986;Glassbrenner and Slack, 1964) in the temperature range 20-300 K (−253 • C to 27 • C) and coefficients R 2 = 0.989 for thermal conductivity and R 2 = 0.999 for specific heat were adopted.On the other hand, for the EG solution, the relationships given the producer MeGlobal TM (MeGlobalTM, 2008; Zhou et al., 2013) were used.In addition, the fuzzy temperature appearing in Eqs.(3.12) and (3.13) should be expressed in • C.

Example of computations
The study analyses one-dimensional heat transfer in a microfluidic system (see Fig. 1) with the dimensions: The microchannel contains the working fluid.During the cooling process, liquid nitrogen is applied, and because of the presence of the liquid form and particles of evaporated nitrogen, its flow is considered two-phase.In contrast, during warming, water is introduced into the system, whose flow is assumed to be single-phase.The working fluid is characterised by the following variables: for cooling respectively (Piasecka-Belkhayat and Skorupa, 2023; Zhou et al., 2013).
In the numerical example, the parameters related to the time step and the given geometric grid are also specified: ∆t = 0.01 s, h 1 = 2.0202 • 10 −6 m, n el,1 = 100 (the number of elements is l 1 = 99).In addition, a numerical analysis was also performed for the modified grid, where h 2 = 1.005 • 10 −6 m, n el,2 = 200 (the number of elements is l 2 = 199).The initial condition determines the values: The interesting thing is the method of introducing trapezoidal fuzzy numbers into the model.At the time t = 0, the deterministic values of the uncertain thermal parameters are determined, and then the triangular and trapezoidal fuzzy numbers are defined according to the relations: λ w or s = (λ w or s −0.05λ w or s , λ w or s , λ w or s +0.05λ w or s ), c w or s = (c w or s −0.05c w or s , c w or s , c w or s + 0.05c w or s ) and λ w or s = (λ w or s − 0.025λ w or s , λ w or s + 0.025λ w or s , 0.025λ w or s , 0.025λ w or s ), c w or s = (c w or s − 0.025c w or s , c + 0.025c w or s , 0.025c w or s , 0.025c w or s ) -compare with the definition of trapezoidal fuzzy numbers provided in (Piasecka-Belkhayat and Korczak, 2020; Skorupa, 2023).As a consequence, the obtained results are triangular and trapezoidal fuzzy numbers as well.
Firstly, the results are reported for a mesh containing 100 nodes (n el,1 = 100).Figure 3 shows changes in the fuzzy temperature as a function of time for (a) cooling and (b) warming for α = 0.25.The results have been prepared for the point in the centre of the sample (point C in Fig. 1, when z = H w + H s ) as a solid line and for the point close to the contact between the sample layer and the chip wall (close to point B in Fig. 1, when z = 1.01 • 10 −4 m) as a dashed line.Because the obtained interval width is small, zoomed-in approximations of selected sections of the diagrams have been created.As can be observed, the minimum temperature was achieved after 14.1 s, confirming that a high cooling rate was used in the process (it is an average value of the fuzzy number).In addition, it can be concluded that point B responds more rapidly to temperature changes caused by the working fluid.The physiological temperature of the sample was restored after about 7.15 s during warming (it is the average value of the fuzzy number).Furthermore, from these plots and simulation data, it is possible to estimate the time required to pass through the "dangerous temperature region" (DTR), which typically occurs between −20 • C and −90 • C (Zhou et al., 2013).In our case, it is equal to 0.2 s and 0.26 s for cooling and warming, respectively.were produced, as well as an approximation of certain fragments of the curve.From Fig. 4, it can be deduced that the chart for cooling stabilises at a certain level.However, for warming, a sudden increase in the degree of crystallisation occurs at the time of passing through the DTR, and then goes to 0. This is caused by the phenomenon of recrystallisation.It should be noted that the recrystallisation phenomenon is also reflected in Fig. 3b, where at a certain moment the temperature drops despite the continuous warming process.
Based on the analysis of the change in the degree of crystallisation, its maximum values can be indicated (average of the given fuzzy number).For the cooling and warming process, χ = 1.075 • 10 −7 and χ = 0.999 were obtained, respectively.Figure 5 presents a comparison of the trapezoidal fuzzy temperature and the trapezoidal fuzzy degree of crystallisation with the results obtained for triangular fuzzy numbers (compare with (Skorupa, 2023)) for different values of the parameter α.The graphs depict the moment s of simulation when t = 0.3 for point C (compare with Fig. 1).It can be noted that the values for α = 0 coincide with each other, and that as the value of the parameter α increases, the intervals are narrower.One can see that for triangular fuzzy numbers for α = 1 the deterministic values are obtained.Next, an analysis was performed for a mesh that includes 200 nodes (n el,2 = 200).Very similar results were obtained, which in the graphical form would look close to those shown in Figs.3-5, so it was decided to compare the selected values in Table 1 for n el,1 and n el,2 .The results presented in Table 1 are expressed as trapezoidal fuzzy numbers for α = 0.25 at point C (z = H w + H s , cf.Fig. 1).

Discussion
Examining the results obtained, it is worth noting that vitrification was successfully modelled and the temperature given by the working fluid was reached relatively quickly in the sample domain.Furthermore, the study of the degree of crystallisation, which complements the thermal analysis, allows us to estimate the tendency of the solution to vitrify.It can also be treated as a marker of potential sample damage.As reported in the literature, this criterion assumes that χ < 10 −6 (Piasecka-Belkhayat and Skorupa, 2023; Skorupa, 2023; Zhou et al., 2013).In our simulation, this condition was fulfilled for freezing, but not for warming.This is due to the recrystallisation that occurs during the DTR transition.
The results can be compared with the data presented by Zhou et al. (2013).First, in the article published by Zhou et al. (2013), it was indicated that the DTR transition time was 0.042 s for cooling and 0.057 s for warming.It can be said that these values are more satisfactory than those calculated in our study.Similar observations appear in analysis of the maximum degree of crystallisation, which in Zhou et al. (2013) was χ = 2 • 10 −11 for cooling and χ = 2.4 • 10 −3 for warming.Those comparisons suggest that our mathematical and numerical models need some improvements.
On the other hand, it is worth exploring the results in terms of the uncertainties introduced.In this article, uncertain quantities are described by triangular and trapezoidal membership functions, which are non-probabilistic algorithms.Undoubtedly, fuzzy logic is a useful tool for modelling the inaccuracy of biological systems in the cryopreservation process.Comparing the results for triangular and trapezoidal fuzzy numbers, one can deduce that the calculated values are close to each other, where for α = 0 the same intervals were achieved.A regularity can also be noticed that the higher the values of parameter α, the narrower intervals are obtained.
When focusing on the fuzzy arithmetics, it is also necessary to refer to the works (Piasecka--Belkhayat and Skorupa, 2023; Skorupa, 2023), in which a similar problem was considered by applying the interval set theory.In this case, the data received coincide with the simulation results indicated in the literature (Piasecka-Belkhayat and Skorupa, 2023; Skorupa, 2023).Moreover, it can be observed that the results for α = 0 are the same as for interval numbers, for which a divergence was defined as 5% of the deterministic value of the thermophysical parameters.It is worth remembering the difference in defining the interval and fuzzy numbers, where for interval numbers the membership function is only equal to 1 or 0. For fuzzy numbers, the partial membership of a given quantity can be formulated (membership function in the range [0, 1]), and the modelled phenomenon can be described in more detail.Therefore, this approach allows for a more flexible representation of quantities and their analysis.
The convergence of the finite difference method is related to the mesh step.Reducing the grid step affects the results close to the real solution.In this paper, additional calculations were performed for a mesh step reduced by half.The results obtained were slightly different from the initial value of the grid step (see Table 1).However, it should be emphasized that a mesh step that is too small increases numerical simulation time.

Conclusion
The paper presents a model of the vitrification process performed in a microfluidic system.The mathematical and numerical model considers uncertain quantities by applying the fuzzy set theory and α-cuts concept.Comparison of the results with data from the literature confirms that fuzzy numbers are effective in the modelling of uncertainty.Therefore, imprecision of thermophysical parameters is included without drastically increasing the computation time.The differences in the results suggest that the prepared model should be improved.For example, it is possible to change the mesh parameters or reduce the time step.However, after analysing the results shown in Table 1, it can be deduced that by modifying the mesh by half, similar results are attained.Another interesting idea is to calculate the results using different numerical methods, for example, finite volume method (FVM) like Zhou et al. (2013) or finite element method (FEM).It is also worth considering extending the model to investigate energy activation and ice nucleation.By exploring the nucleation rate and ice growth (a change in the diameter of the ice crystals), one can further estimate the volume of ice crystals.
An interesting aspect of the paper is the discussion about differences between fuzzy numbers and interval numbers, which was used in the literature to model cryopreservation and crystallisation (Piasecka-Belkhayat and Skorupa, 2023;Skorupa, 2023).
In further research, it would also be worth analysing the nucleation process and ice crystal growth.
2018; Song et al., 2010; Zhang et al., 2017; Zhou et al., 2013).Song et al. (2010) presented an example of modelling vitrification process modelling performed by the droplet-based method.The same technique for cell aggregates was also examined by Shi et al. (2018).Zhang et al. (2017) studied vitrification, in which the tube with the sample was immersed directly into the working fluid.That numerical calculation was supplemented by a model of the probability of intracellular ice formation.Zhou et al. (2013) investigated the vitrification process, which was analysed in a microfluidic system.

Fig. 3 .
Fig. 3. Fuzzy temperature as a function of time during (a) cooling and (b) warming for α = 0.25 with zoomed fragments

Fig. 4 .Figure 4
Fig. 4. Fuzzy degree of crystallisation as a function of time during (a) cooling and (b) warming for α = 0.25 with zoomed fragments Figure 4 illustrates the change in the fuzzy degree of crystallisation as a function of time in (a) cooling and (b) warming phases for α = 0.25.Once again, functions for points B and Cwere produced, as well as an approximation of certain fragments of the curve.From Fig.4, it can be deduced that the chart for cooling stabilises at a certain level.However, for warming, a sudden increase in the degree of crystallisation occurs at the time of passing through the DTR, and then goes to 0. This is caused by the phenomenon of recrystallisation.It should be noted

Fig. 5 .
Fig. 5. Fuzzy temperature (a), (b), and fuzzy degree of crystallisation (c), (d) for t = 0.3 s during cooling for different values of the parameter α

Table 1 .
Results for α = 0.25 at point C for n el,1 and n el,2