Jorge Aguilar Barrera
Ingeniero Superior Industrial (Especialidad Químico)
e-mail: jorgeab77@hotmail.com
| |
Jorge Aguilar BarreraIngeniero Superior Industrial (Especialidad Químico) e-mail: jorgeab77@hotmail.com |
| Degree Thesis |
|
Degree Thesis:
A NEW TOOL FOR COMPUTING THE REAL TRASIENT TEMPERATURE AND MASS FIELDS OF PERFUSED TISSUES AND ORGANS FROM A DIGITAL IMAGE AND FINITE-DIFFERENCE METHOD
J. Aguilar1, S. Álvarez2, J.M. Salmerón2 and R. Risco1
1CryoBioTech; Escuela Superior de Ingenieros, Universidad de Sevilla, 41092 Sevilla (Spain)
2Depto. de Ingeniería Energética y Mecánica de Fluidos; Escuela Superior de Ingenieros, Universidad de Sevilla, 41092 Sevilla (Spain)
ABSTRACT:
A calculation of the temperature field during the transient of a slice of perfused heart is presented. The novelties of the method are several: the way of introducing the boundary condition (directly from a digital picture), the sub-structuring approach and the computing technique, a variant of finite-difference. It can be of great utility in tissue engineering and cryobiology.
1.- INTRODUCTION:
The problem of heat and mass transport in perfused tissue and organ is crucial for a good understanding of several physiological processes as well as the application of some techniques in the fields of diagnosis, treatment, tissue engineering, cryosurgery and cryopreservation among others.
Several attempts of computing the temperature field of perfused biological samples has been published in the past. A general analytical theory was settled by Klinger (1). Later on, the use of Pennes (Pennes, H.H. “Analysis of Tissue and Arterial Blood in the Resting Human Forearm” Journal of Applied Physiology, 1, 93-22 (1948)) and Weinbaum-Jiji ( Weinbaum, S. and and Jiji, L.M. “ A New Simplified Bioheat Equation for the Effects of Blood Flow on Local Average Tissue Temparature”, ASME Journal of Biomechanical Engineering 107, 131-139 (1985)) bio-heat transfer equations, and their interesting generalizations (Zhu, M., Weinbaum, S., Jiji, L.M. and Lemons D.E. “On the generalization of the Weinbaum-Jiji Bio-heat Equation to Microvessels of Unequal Siza; The Relation Between the Near Field and Local Average Tissue Temperatures”, ASME Journal of Biomechanical Engineering 110, 74-80 (1988), Song, W.J., Weinbaum, S., Jiji, L.M. and Lemons D.E. “A Combined Macro and Mocrovascular Model for Whole Limb Heat Transfer”, ASME Journal of Biomechanical Enginnering 110, 259-268 (1988), Wisser. E.H., “An Analytical Solution Countercurrent Heat Transfer Between Parallel Vessels Wih a Linear Axial Temperature Gradient” ASME Journal of Biomechanical Enginnering 110, 254-258 (1988) and Zhu, M., Weinbaum, S., and Lemons, D.E., “A Three Dimensional Variable Geometry Countercurrent Model for the Limb Heat Transfer” ASME Journal of Biomechanical Engineering 114, 366-376 (1992)) allowed an average calculation of the perfused tissue. A new approach was given by the use of statistical models (Baish J.W. “Formulation of Statistical Model of Heat Transfer in Perfused Tissue”, Journal of Biomechanical Engineering 116, 521-527 (1994)), showing information on the most probable temperature. With the arrival of computers, these analytical and statistical methods have been complemented with numerical solutions. However, three important limitations of these methods are present: i) in some cases they do not deal with the transient (Blanchard C.H., Gutierrez, G., White, J.A. Roemer, R. B. “ Hybrid Finite Element-Finite Difference Method for Thermal Analysis of Blood Vessels”, Int. J. Hyperthermia 16, 341-353 (2000); Huang, H.W., chen, Z. P., and Roemer, R.B., “ A Counter Current Vascular Network Model fo Heat Transfer in Tissues, ASME Journal of Biomechanical Engineering, 118, 120-129 (1996)); ii) some of them only deal with a 1D approach (Zhang, Y.T., Liu, J., “Numerical study on three-region thawing problem during cryosurgical re-warming”, Medical Engineering and Physics 24, 265-277 (2002); Cui, Z. F., Dykhuizen, R. C., Nerem, R. M., and Sembanis, A., “Modeling of Cryopreservation of Engineered Tissues with One-Dimensional Geometry”, Biotechnol. 18, 354-361 (2002)) iii) or the perfusion through the vascular system is not considered (Cui, Z. F., Dykhuizen, R. C., Nerem, R. M., and Sembanis, A., “Modeling of Cryopreservation of Engineered Tissues with One-Dimensional Geometry”, Biotechnol. 18, 354-361 (2002)) iv) and finally, all the geometries used are not real, they are only approximations generated in different ways (all of them).
In this paper we have developed a new tool that tries to overcome some of theses limitations. It is based on a digital treatment of imagines, finite differences and sub-structuring.
Our main interest is tissues and organs cryopreservation. In this case a good knowledge of the temperature and the cryoprotector fields are necessary in order to guarantee that at any place and at any time, the concentration of cryoprotector is high enough to avoid ice formation. Several steps have been done in this direction (Zhang, Y.T., Liu, J., “Numerical study on three-region thawing problem during cryosurgical re-warming”, Medical Engineering and Physics 24, 265-277 (2002); Cui, Z. F., Dykhuizen, R. C., Nerem, R. M., and Sembanis, A., “Modeling of Cryopreservation of Engineered Tissues with One-Dimensional Geometry”, Biotechnol. 18, 354-361 (2002)) by others. In the case of cryopreservation, the dehydratation of the cell is frequently very important for practical purposes, due to the fact that a reduction in the amount of water implies a reduction in the amount of the potential ice, and besides, a depression of the melting point as a colligative effect. In this paper we do not deal with dehydratation. However we set the basis for a computation of it, because in order to do that it is necessary a previous knowledge of the transient temperature and mass fields, as our model does.
To understand the way in which our tool works, let us fix our attention into a perfused heart through its coronary system during a cryopreservation protocol. Then, let us suppose that we are interested on the temperature conditions of a certain two dimensional region of this heart. A transversal slice of this part is taken and a digital picture is made from it. With the picture in BMP format, a first algorithm, CRYOJAB-Colour, is run in order to asses three different codes to the three relevant parts of the picture: capillary, cardiomyocyte and intercellular tissue. From the file generated, a second algorithm, CRYOJAB, is started. This algorithm has two inputs: the matrix generated from CRYOJAB-Colour and the cooling profile of the perfusion of the cryoprotector agent. The output of this program is a set of files with the temperature field of the 2D sample a any time during the transient. These files can be graphically displayed for example with Matlab.
Although in this paper we only tackle the 2D problem and we compute the temperature field as an example, a natural extension of this work will be to study the three dimensional problem as 3D=2D+1D, and the mass field. These works are underway.
The structure of this paper is the following. In section 2 we shall describe the geometry that we are going to study, we shall justify the 2D approach, and show some consideration about the sizes of the pixels taken as optimal in the image treatment. In section 3 we explain how to deal with the digital picture of the tissue in order to generate the corresponding tem BMP file, and the . In section 4 we explain and clarify all the hypothesis that we have made in order to simplify the problem. Some of this simplifications can be easily relaxed in future extensions of this work. Section 5 is devoted to develop the mathematical equations of the model, based on finite differences. In order to deal with tractable matrices we have implemented a technique called sub-structuring, allowing a faster inversion of the present matrices; this is explained in section 6. Section 7 is an explanation in detail of the algorithms that we have develop. Finally, the solution of a real example is given in section 8, displaying the files generated and the final result.
2.- GEOMETRY
In order to introduce our tool, we shall fix our attention on a concrete case: a perfused heart through its coronaries with a cardioplegic solution from 36ºC to –10ºC. In a real case, the solution should contain some cryoprotector agent in a proper concentration (DMSO, PEG, etc) in order to avoid ice formation at this low temperature. Our tool allows to decide the cooling profile that we are going to apply. In this case we shall suppose a constant cooling rate of 0.1ºC/s; that is: T(t)=36-0.1 t (ºC). Of course, we can modify this profile at will. In particular, the importance of the cooling and thawing rates have been reported elsewhere (Zhu, Q., Layne JR, J.R., Claydon, M., Hicks JR, G.L., and Wang T., Freezing preservation of the mammalian cardiac explant VI. Effect of the thawing on functional recovery". Cryobiology 29, 478-484 (1992); Pegg, D.E., Jacobsen, I.A., Diaper, M. and Foreman J. The effect of cooling and warmaing rate on cortical cell function of glycerolized kidneys. Cryobiology 21, 529-535 (1984).; Jacobsen, I.A., Pegg, D.E., Starklink, H., Chemnitz, J., Hunt, C., Barfort, P. and Diaper M.P. Effect of cooling and warming rate on glycerolized rabbit kidneys. Cryobiology 21, 637-653 (1984)). The tronco-cilindrical form of the cardiomiocyte and the disposal of the vessels along them, make our 2D approach specially indicated for this case, as we shall see in section 4.
Suppose we are interested on the temperature field of a particular part of the heart. The first step is to take a digital picture of a transversal slice of it (figure 1).
figure
As we can see, these cardiomyocytes and capillaries exist in a ratio of 1:1 in any typical part of the heart absent of big vessels. The transversal size of the capillaries are around #microns and that of the cell of #microns. Although this is in principle irrelevant for the working of our algorithms, it was necessary to have some estimations in order to optimise it.
3.- ALGORITHMS
In this section we shall describe the four algorithm in which our tool is based on. They are called GIMP, CRYOJAB-Colour and CRYOJAB.
GIMP: Is a commercial utility that allows to change from a bmp file to a c-source file, an operation that is needed due to the fact that CRYOJAB and CRYOJAB-Colour are in C. Furthermore, GIMP allows to make sharper the contrast in the digital picture of the tissue. This operation is needed in some cases in order to better define the different regions of the picture.
CRYOJAB-Colour: It has two inputs: the file with the digital picture of the tissue and the number of nodes and substructures that we want to make from it. A large number of nodes increase the precision of the result, but the computer spends more time to get the solution. So, this number should be a compromise between these two aspects. The importance of the substructures will be discussed in the section 6. CRYOJAB-Colour assign a code to the different material of the sample cell=0, interconnetive=1 or capillary=2, depending on the predominant material of the node. This information is store in the output set of files that CRYOJAB takes as inputs.
CRYOJAB: Is the core of the tool. It has three inputs: the file generated by CRYOJAB-Colour, the physical properties of the three materials (specific heat, thermal conductivity and density), and the cooling profile (starting temperature, end temperature and cooling rate). From these date, CRYOJAB generates a set of files that gives the experimenter the temperature field of the tissue at any time during the cryopreservation protocol.
MATLAB: This commercial software was used in order to display a graphical representation of our results. It takes the set of files produced by CRYOJAB and shows two or three dimensional picture of the temperature field.
4.- HYPOTHEISIS OF THE MODEL
In order to simplify the treatment, and taking into account that we are mainly interested on the heat transfer during a heart cryopreservation protocol, the algorithm is based on a series of hypothesis that we should clarify.
The first of this is geometry selected: 2D. We have already said that the cardiomyocytes have a tronco-cilindrical shape, and that capillaries are along them. This translates into parallel lines (z exe) of heat transfer. Thus, a transversal slide (xy plane) of this lines gives a good approximation of the real situation. Of course, in every z=const. plane there is a different temperature field. A natural extension of this work would be to treat the 3D case as 2D (xy) +1D (z). Nevertheless, if the speed of the cooling fluid is high enough, the 2D approach represents a good approximation.
The second approximation consist on considering only three elements: cell, interconective and capillary. This situation represents quite well a piece of myocardium. However, there are other tissues that could require a more detailed a description. This would represent a straightforward generalization of our tool. It is also true that, due to the similar values of the physical properties of the interconetive and cardiomyocytes, in principle, only two material could be considered. However, as we said before, the dehydratation is something very important in during a cryopreservation protocol. This means that, from this point of view, the behaviour of the cell and interconetive, are very different. In order to take dehydration into account in future works, we have preferred to consider them as different.
The fourth approximation consists on supposing that the physical properties of the material do not change with time. This assumption is based two considerations: First, in the present case there is no change of phase. This can be achieve by using a proper cryoprotector agent. And second, we are studying a range of temperatures smaller than that in which a significant change of these properties takes place. Nevertheless, a slight modification of this tool can also deal with variable physical properties as a function of the temperature.
The fifth is to suppose that there is no mass transfer during the process. So, the heat is only transferred by thermal conductivity. Probably this hypothesis is the stronger one. Besides, diffusion is an essential ingredient of almost any cryopreservation protocol. However, an important fact is the similar treatment of diffusion and conduction from a mathematical point of view. They are represented by the same equations, by replacing the thermal conductivity by the diffusion coefficient. This means that our tool can also deal with diffusion. By coupling diffusion with conduction we could give a more realistic description of the temperature field. This work is underway.
6.- SUBSTRUCTURING
The resolution method described in previous section is valid for any kind of problem. Nevertheless it has an important limitation, the computational time. As the number of nodes used in order to model correctly the problem is very high the matrix dimensions are huge. The limit will be imposed as a consequence of the elevated calculation time required. In fact this could make inefficient the described procedure and, in addition, the user could lose his interest in the software tool.
Therefore, it is strongly recommended to reduce the computational time by means of reduce the size of the systems of equations and the number of eigenvalues and eigenvectors to assess, without decrease the density of the mesh necessary to obtain accurate results.
Numeric models in which a structure is divided into multiple substructures are extensively used in Structures theory.
This technique consists on modelling the substructures individually and then to establish constraints between the degrees of freedom of every one. As a result the system could be divided in N sub-zones where each zone is linked in some way to the other, so all together conform the whole system. This technique is known as sub-structuring. Next figure shows a schematic representation of the idea of sub-structuring
The main advantage of sub-structuring procedures is the saving in computational time that they produce as a consequence of lower dimensions of the matrix- it is easier to invert various matrixes of reduced dimensions that only one inversion of a matrix which dimension is the addition of those.
Nevertheless, these procedures have drawbacks. For example, it could be mentioned that as a result of our decision of avoid assess the whole system in one step the eigenvalues can not be obtained directly, in this situation we must to use a traditional method out of the infinite array in the case that we need this additional information. This is not going to be our case because we are not solving the problem using transfer functions but finite differences.
The strategy is the folloing:
There are a high number of ways to carry out the link-up. They can be classified as a function of:
The more interesting strategy for the present work is those that use as link-up variable the temperatures in the adjacent nodes.
3. Iterative sub-structuring method with coupling through equivalent temperature.
The assessment procedure will be iterative in order to avoid working with systems of equations of very high dimensions. Starting from an initial solution of temperature field we could solve, one by one, all the sub-structures; when a sub-structure is resolved its temperature field in the proximities of contact with another sub-structure will be used in order to assess the next sub-structure. The temperature field in the proximities of contact with the next sub-structure is called equivalent temperature.
The sub-structure could have homogeneous or heterogeneous properties. In this study we have used sub-structures with heterogeneous properties; however, the mesh has been refined reducing the cell size to be able to assign to each node a single property.
The order in which the sub-structures are solved can be changed in order to improve the convergency by means of accelerate the propagation of boundary conditions across the whole structure. In this case, as the solutions of multidimensional conduction problems are very complex and it is necessary to process all the cases in general, it has been chosen an isotropic order of resolution without preferential sub-structures.
4. Interpolation.
The previous method carries out the coupling between two sub-structures node by node due to this the method has the next limitation, if two adjacent sub-structures have different cell size of the mesh, the nodes of both in their boundaries will not be each one opposite to the other. To avoid this problem, the coupling could be carried out through interpolation.
The interpolation will consist of establish a linear approach between the temperatures of two nodes in the same substructure and then calculate the corresponding temperature of the opposite node in this line.
Next graphics shows the described case:
 |