Computational modelling has received significant attention in the research community over the past few decades due to its pronounced contribution in better understanding of the nature, as well as in the development of advanced technologies. The modelling of more and more complex transport physical systems helps the community to address important issues, like identification of environmental problems, improvement of technological processes, development of biomedical applications, etc. However, the physical modelling is only one part of the entire problem. In most realistic scientific computing applications, a physical model cannot be solved in a closed form, thus adequate numerical approach is required. The continuous physical domain is replaced by its discrete form. In the majority of the simulations, the classical numerical methods, such as the Finite Volume Method (FVM) (Ferziger and Perić, 2002), Finite Difference Method (FDM) (Őzisik, 2000, Shashkov, 1996), Boundary Element Method (BEM) (Białecki, et al., 2006, Wrobel, 2002) or the Finite Element Method (FEM) (Rappaz, et al., 2003), are used to solve these problems.
An important part of the numerical approach is the effective implementation of the solution procedure on modern computer architectures. Realistic computations might take vast amount of time to complete, especially in 3-D domains with high number of computational nodes. The computing performance can be improved by numerous approaches, such as network computing, grid computing, cloud computing and other variants of "distributed computing" (Trobec, et al., 2009).
In recent decades, computer simulations have proved a great help in understanding and solving a variety of problems in science. Advances in computer technology enable simulation of natural phenomena that cannot be subject to experiment in reality since experiments would be ecologically problematic, dangerous to humans, or would cost vast amounts of money (Martino, et al., 1994). Especially in medicine, experiments are often difficult to perform because human subjects are involved. Measurements during clinical procedures are time consuming and often not as accurate as desired, because many parameters are difficult to control (McMaster, et al., 1978). Many examples can be found in medicine, where in vivo experiments and measurements are often difficult, dangerous or even impossible (Trunk, et al., 2003), especially if deep tissues or vital organs are in question (Tikuisis and Ducharme, 1991). However, simulation can provide the only safe and inexpensive insight into physiological processes. With the use of computer simulation, it is possible to calculate, analyse, and visualize both stationary temperature fields and the changes of temperature in time in living biological tissues (Grana, 1994). The temperature of the human tissue is an important factor in many fields of physiology, sports (Olin and Huljebrant, 1993), cryotherapy (Trobec, et al., 1998), etc. It has been recognized that the temperature field is influenced by the environmental conditions, the temperatures of neighbouring tissues, the muscle metabolism and the blood circulation. Different tissues have different physical and thermodynamic properties and respond diversely to temperature changes (Bernardi, et al., 2003, Pennes, 1948, Trunk, et al., 2003). The temperature field varies in space and time in different parts of the investigated domain.
In this chapter we present numerical solution procedure for the treatment of heat transfer in biological tissues. Different body parts are modelled and given as numerical examples. All cases are focused on realistic three dimensional (3-D) models and realistic process parameters. In all numerical examples we use 3-D closed cavity as a computational domain. Within the domain, high resolution models of body parts are used. The models are obtained through expert assisted segmentation of human body parts from the Visible Human Dataset (VHD) (Ackerman, 1998). The body parts are composed of different tissues (muscles, bones, cartilage, fat, vessels and nerves). To respect the complexity of the modelled body parts, a one-domain approach is employed through spatial dependent thermo-physical properties. The final geometrical domain is composed of small cubic voxels, resulting in millions of discretization points.
The Pennes’ Bio-heat equation (Pennes, 1948) represents the basis for the physical model. It describes the energy continuity within biological tissues by incorporating heat convection, heat transfer between blood and tissues, and heat production by metabolism. Some important extensions have been introduced, in particular, an inhomogeneous spatial model composed of tissues with different characteristics and modelling of the heat contributions from arteries, blood perfusion, and metabolism as functions of the surrounding tissue temperature. This improved Bio-heat equation was evaluated in terms of sensitivity and accuracy and solved numerically on single and parallel computers (Akl, 1997, Martino, et al., 1994). Outside biological tissues, the Bio-heat equation is reduced to diffusion equation as both additional source terms (heat transfer between blood and tissues, and metabolism) vanish. The flow outside tissues is modelled as the Navier Stokes fluid (Ferziger and Perić, 2002), while in the tissue no flow is allowed. We consider only incompressible fluids.
A substantial amount of work on closed form and numerical solutions of the Bio-heat equation has been published (Bernardi, et al., 2003, Bowman, et al., 1975, Charny, 1992, Pennes, 1948). In our examples, the governing system of partial differential equations (PDEs) is numerically solved through spatial discretization with the FDM and the Euler time stepping (Őzisik, 2000). We present three different numerical studies regarding different topics in medicine, where results are computed using the proposed physical model and numerical solution procedure.
The first numerical example is dedicated to a study of the temperature changes in a resting proximal human forearm. By computer simulation, we explained the variations in the Pennes’ well-known in vivo measurements of the steady-state temperatures along the transverse axis of the proximal forearm (Pennes, 1948). Such in vivo measurements are no longer practiced and are very valuable, as are their additional analyses. The domain of the simulation is a region of the forearm 5 cm – 15 cm below the elbow. We use Robin boundary conditions on the contact between the skin and the ambient air, while on the domain boundaries the heat flow is set to zero. The heat production by tissue metabolism is modelled using the Q10 rule (Tikuisis and Ducharme, 1991), while the heat exchange between the blood and the tissue is modelled as a function of the local temperature and the regional blood flow. We confirm the stability and accuracy of the method by varying the simulation parameters, the initial and boundary values, and the model dimensions, with subsequent analysis of the results. Our simulations indicate that the fluctuations of the measured steady-state temperatures are a natural consequence of a complex interplay between the position of the measuring probes, the anatomical position of the main arteries, the dimensions of the forearm, the blood flow, the inhomogeneity of tissues, and the environmental temperature.
The second illustrative example is simulation of a human knee exposed to topical cooling after anterior cruciate ligament (ACL) surgery (Martin, et al., 2001, McMaster, et al., 1978, Tikuisis and Ducharme, 1991, Trobec, et al., 2008). The computational domain is composed of the knee model and surrounding layers: protective bandage, cooling layer, isolating blanket and ambient air. The ambient air is held at constant initial temperature as it is assumed to be well mixed. The heat flux from the first and last slices is kept constant in order to imitate the influence of the leg not exposed to the cooling. Simulations are performed on the described model using realistic process parameters gathered from measurements. Two different methods of topical cooling are compared; first, the use of a gel-pack filled with refrigerated gel, which is exposed to ambient temperature and therefore becomes less and less effective as the gel is warming; and second, the use of a cryo-cuff cooled by a liquid at constant temperature maintained by an external cooling device. The results are presented in terms of temporal temperature development in critical points of the knee. The temperature dynamics is also validated by in vivo temperature measurements during cryotherapy for different regimes of cooling. We show that the second cooling approach gives much better results regarding the desired temperature distribution within the knee.
The cooling of a human heart during surgery is taken as a final, most complex example. In the last example full model is incorporated. Besides the Bio-heat equation Navier Stokes fluid flow is solved as well. The topical cooling (Olin and Huljebrant, 1993) is simulated through consideration of partially submerged heart in the cooling liquid. The effect of natural convection in the cooling liquid is tackled. We show that considering solely diffusive heat transport, as done in present (Šterk and Trobec, 2005, Trobec, et al., 1998, Trunk, et al., 2003), simulation does not supply adequate results.
The presentation of the numerical examples is organized as follows. First, a short historical overview of the tackled subject is introduced. Next, the spatial 3-D geometrical model of the designated body part is described. Then, different simulation modalities and parameters are stated. Afterwards, analysis of the simulated results and their comparison with the measured values, if available, is presented. Finally, conclusions are presented.
Before presenting the numerical examples, the considered physical model and solution procedure are introduced.
2. Physical model
The present chapter is focused on the formulation of the physical model suitable for describing biological tissues. The main part of the modelling is heat continuity, where biological response is treated as volumetric heat source. The power of such sources depends on the temperature of the tissue. Besides bio-heat transfer within tissues, the surrounding is modelled as a fluid in order to simulate the heat exchange with the ambient. The surrounding is treated as a Newtonian fluid described by Navier-Stokes equation. The fluid is assumed to be incompressible.
The model in the present chapter is developed from the classical continuum conservation equations for energy, momentum and mass. The differential forms of the conservation equations are:
for energy (1), momentum (2) and mass (3), where stand for density, time, velocity, pressure, deviatory stress tensor, body force, volumetric heat sources, specific heat, temperature and heat flux, respectively.
2.1. Model assumptions
To construct an accurate and stable model, on the one hand, and still solvable in a reasonable way, on the other hand, some basic assumptions are made:
all material properties are assumed to be temperature independent,
no flow is allowed in the tissues,
Newtonian fluid is considered,
the biological response of the tissues on temperature change is modelled with the Pennes’ Bio-heat equation,
the metabolic heat production per unit mass is assumed to obey the Q10 rule,
all materials are isotropic.
2.2. Bio-heat transport
The local mass heat flux is modelled by the Fourier law
with standing for the thermal conductivity tensor. Considering the model assumption for material isotropy, the thermal conductivity diffusion tensor is rewritten in the following form:
where stand for identity matrix and thermal conductivity.
where and are the density and the specific heat of the blood, and stands for the heat generated from the biological response. The dimensionless parameters and are used to control the spatial/material dependency of the biological response and were not included in the original paper. The functions characterize the behaviour of heat generation with respect to the temperature of the tissue and are defined latter in this chapter. The temperature stands for the reference tissue temperature.
Equation (7) models the heat transfer between the arterial blood and the tissues, and the heat production arising from tissue metabolism, in addition to the classical heat transport.
To close the energy transport model, both biological response functions have to be defined.
Published measurements of the blood flow rate show that it increases with the tissue temperature (Ducharme and Frim, 1988). Earlier measurements of the forearm blood flow were limited to plethysmography that measures the mean total blood flow through a given region of the forearm (Barcroft and Edholm, 1942). Later, more data became available also for regional blood flow of specific tissues (Cooper, et al., 1955, McElfresh and Kelly, 1974, Sugimoto and Monafo, 1987). Several studies have established that the muscle blood flow does not increase with temperature; most of the increase takes place in the skin (Detry, et al., 1972). The value of the mean blood flow rate V as a function of the skin temperature Ts has been reported in several publications (Barcroft and Edholm, 1946, Cooper, et al., 1955, Ducharme and Tikuisis, 1994, Wenger, et al., 1986, Yue, et al., 2007). We use approximate exponential function:
with newly introduced case dependant dimensionless parameters. The reference value is set to V0 = 1.667 10-4 s-1. Although a more detailed blood flow model could be used in our simulation, for example, incorporating the autonomic thermoregulation of the body, its development is beyond the scope of this chapter.
The metabolic heat production in the human body can be divided into unregulated heat production from voluntary muscle contraction and normal metabolic pathways, and into regulated heat production for maintaining temperature homeostasis at lower ambient temperatures (Bullock and Rosendahl, 1992, Proulx, et al., 2003). The rate of metabolic heat production per unit mass is assumed to obey the Q10 rule (Bullock and Rosendahl, 1992), which is expressed as a function of the tissue temperature:
where stands for the reference metabolic heat production of a tissue at reference temperature Tr = 35 °C. An example of the biological response against temperature change is presented in Figure 1, where the parameters are taken as in the following first numerical example.
Outside the tissue, the Bio-heat equation is reduced to the classical heat transfer equation (Amimul, 2011), as both biological sources vanish:
The model assumes no flow in the tissues and, as a consequence, the advection term vanishes. In future discussions, a tissue is referred as a solid material.
2.3. Fluid flow
In order to specify the momentum equation (2), the stress tensor has to be specified. It is constituted as a Newtonian fluid:
where stands for the dynamic viscosity. Further, the body force in the momentum conservation equation is modelled by the Boussinesq approximation for temperature dependent density:
To close the system, the mass continuity (3) is considered. The model assumes incompressible fluids, thus the velocity is divergent free:
2.4. Thermo-physical propertiesThe thermo-physical properties of living tissues are still not uniformly defined. For present simulation purposes, the average of published data is used (Bowman, et al., 1975, Moses and Geddes, 1990, Ponder, 1962, Poppendiek, et al., 1967)]. For muscles at rest, hr was taken to be about half the human basal metabolic rate, which equals 0.6 W/kg (Harris and Benedict, 1918, Ho, et al., 1994, Silverthorn, 2001). The impact of the metabolic rate production is almost negligible for the temperature field in the forearm; even so, it can be adjusted in the model for each tissue independently.
3. Solution procedure
3.1. Geometric model preparation
The computer model domain of the forearm, the knee, and the heart is derived from the axial anatomic images of a male human body that are available from the Visible Human Dataset (VHD) (Ackerman, 1998). A similar model of the whole body at a resolution of 1 mm3 has been developed within the VHD project, but was not available for public use at full extent when this chapter was written. The axial anatomic images are 24-bit colour cross-sections photographs of a frozen male human cadaver, with resolution of 0.33 mm 0.33 mm, and taken at spatial interval of 1 mm.
To construct the three-dimensional model, several consecutive slices (cross-sections) (Figure 2) of VHD are used. First, exact overlapping of all slices is required to ensure the integrity of the model. Therefore, a reference point on every slice is needed to align the slices. The white circle in the top right-hand corner of the VHD cross-section is selected for this purpose. A simple automatic procedure is used to find the centre of this circle and thus the absolute position of each slice. The advanced elastic matching (Dann, et al., 1989) is not applied at this stage. Second, several impurities and other errors in the original photographs have to be discarded manually. Last, different tissues have to be identified on the slices. This is done manually, under the supervision of an expert for anatomy, with image editing software, by painting over the photographs with predefined colours, always the same colour for the same tissue. Only the most prominent tissues are marked, the rest are approximated with the surrounding tissue. An automated procedure is then applied to scale down the slices to 1 mm1 mm resolution, convert them from colour images to a format recognized by the computer simulator, and stack them together into a 3-D model. Finally, the 3-D model is checked and edited in custom modelling software to smooth out the jaggedness in z axis appearing as a result of the previous processing done only on x-y cross-sections.
Each voxel of the geometric model is characterized by thermo-physical properties.
The described physical model does not have a closed form solution and has to be treated by a numerical approach. The present work deals with the transient problems described by parabolic partial differential equations. Such problems are space and time dependent and therefore discretization in space and time is required.
In the present solution procedure, one of the simplest and most intuitive numerical techniques for time dependent simulation is used. Temporal discretization is performed through two-level explicit discretization. The most commonly and widely used unconditionally stable schemes for temporal stepping are implicit and Crank-Nicolson approaches. Besides two-level methods, also multilevel methods can be used to achieve higher-order approximations. The disadvantage of the multilevel methods is the need to store data from several previous time steps. The multilevel data storage issue can be avoided by using Runge‐Kutta methods.
For the spatial discretization, the Finite Difference Method (FDM) (Őzisik, 2000, Shashkov, 1996) is used, as it has low calculation complexity, is sufficiently accurate, and is simple to implement. In addition, FDM can directly use the presented discretized domain, where uniform discretization elements are used. There are also several numerical techniques with higher order of accuracy and better stability properties. For example, spatial discretization can be performed through weak form of Finite Element Method (FEM) (Rappaz, et al., 2003) or Boundary Element Method (BEM) (Białecki, et al., 2006, Wrobel, 2002). Furthermore, there are several meshless techniques (Atluri and Shen, 2002, Atluri, 2004, Chen, 2002, Kansa, 1990, Liu, 2003, Šarler, 2007) with convenient properties.
However, the complication in the temporal and/or spatial discretization of the domain gains benefits in accuracy and stability, but loses much in terms of computational time. In the light of this work, the simplest two-level method is used with the lowest CPU complexity.
In general, a partial differential equation
is solved. In equation (15) stand for an arbitrary field/quantity, a second order partial differential operator and a position vector, respectively. In the present work, a transport equation is treated, which is governed by the Laplace, the gradient and the divergenceoperators. First and second spatial derivatives are needed to construct those operators. Both derivatives are evaluated through Taylor’s expansion, known as finite differences:
The left side of equation (15) is discretised through two-level time discretization with the finite difference representation of the first temporal derivative. It is rewritten in an explicit form as:
For the finite difference spatial discretization, the Neumann (or Fourier) stability analysis shows that both implicit and Crank-Nicolson are unconditionally stable, while the explicit scheme suffers from stability problems (Smith, 2003). For diffusive-convective problems, as the ones considered in the present work, the differential operator is generally written as:
and for the FDM, the stability criteria derived by Fourier analysis is stated as:
wherestands for the component of the velocity.
Despite the stability issues, in the present work the explicit method is used because of its straightforward local applicability and simplicity. The main drawback of the implicit and semi-implicit forms is the need for solving the global problem, where systems (usually large) of linear equations have to be considered in order to proceed to the next simulation time-step. On the other hand, the explicit scheme uses series of pre-solved linear systems solutions. In other words, the explicit scheme uses the local information from the current time-step () to compute spatial derivatives, while the implicit scheme uses the current local time-step information only for the evaluation of time derivative. Finally, the Crank-Nicholson scheme uses linear combination of the current and the next time-step to improve the accuracy of the method.
One of the most convenient properties of the Explicit FDM (EFDM) is its straightforward parallel implementation, as there are no requirements of any kind for global communication.
3.3. Pressure velocity coupling
An important part of the solution procedure represents the solution of the momentum equation. An additional problem lies in the momentum-mass conservation coupling. The pressure is not explicitly included in the continuity equation and therefore, a special treatment is needed to solve the pressure-velocity coupling problem. The most widely known solution for the problem is the Semi-Implicit Method for Pressure Linked Equations SIMPLE (Ferziger and Perić, 2002), where the problem is translated to the Poisson equation problem. However, the SIMPLE algorithm demands the solution of the global system, which is complex from parallel implementation point of view. There are local alternatives, like the Local Pressure Velocity Coupling (LVPC) used in context with local meshless methods (Kosec and Šarler, 2008a), the artificial compressibility method (ACM), which has been recently under intense research (Malan and Lewis, 2011, Rahman and Siikonen, 2008, Traivivatana, et al., 2007), and the framework for the Finite Difference Method in the SOLA approach (Hong, 2004). In this chapter, the more standard and well accepted pressure velocity coupling, based on the SIMPLE algorithm, is used.
The problem domain is discretized to voxels that can be either liquid or solid. For solid voxels only the velocity is assumed to be zero and no pressure computations take place. Otherwise, standard staggered grid approach (Ferziger and Perić, 2002, Shashkov, 1996) is used, with the temperature and the pressure defined in the middle of each voxel, and the velocity components defined on the respective perpendicular voxel edges.
The Navier-Stokes equation is discretized in time by using the explicit time discretization:
where stands for the intermediate velocity and the superscript 0 denotes that the current time-step values are used. As the explicit time discretization is employed, the new time-step velocity is computed from the velocity values in the previous time-step. However, instead of a new time-step value, the intermediate velocity is introduced, which does not take into account the mass continuity (3). This problem is solved by introducing the velocity and pressure corrections. The added velocity correction drives the intermediate velocity towards divergent free field:
where and stand for the velocity and the pressure correction terms, and the superscript 1 denotes the next time-step. With the corrected velocity and pressure, the momentum equation can be rewritten as:
From equation (24), directly follows:
In order to solve the Poisson equation, additional boundary conditions are required. The boundary conditions for pressure field can be formulated by multiplying the momentum equation with the normal boundary vector. The boundary conditions are obtained with a similar approach used to construct the equation for pressure correction, stated as:
where stands for the velocity boundary condition and denotes the normal boundary vector. It is common to use the normal pressure boundary condition when no-slip and no-permeable velocity boundary conditions are used:
With the pressure corrected according to (23), the velocity can be corrected as:
where stands for the next time-step velocity. With computed velocity and pressure, the simulation can proceed to the next time-step. The computation cycle is closed.
3.4. Implementation of parallel program
The simulator is written in C++ and compiled with GCC. The parallelization is done by domain decomposition. The domain decomposition is one-dimensional: the domain is split in the z axis into subdomains of approximately equal size. The number of subdomains equals the number of processors participating in the solver execution. The sizes of the subdomains are determined at the program initialization and are constant for the whole execution; therefore, the load-balancing of the processors is static. One-dimensional decomposition was chosen because of its simplicity. To calculate a point on the border of a subdomain, data from one or more points from neighbouring subdomains might be required. In such cases, point-to-point communication between processors of those neighbouring subdomains is performed. MPI library is used for communication. It enables the program to execute in parallel on multicore computers, computer clusters or other distributed computer systems.
The multigrid method, used in our solver, is parallelized using similar principles. Calculations on fine grids are done in parallel using the same domain decomposition, while a single processor does the computation on the coarsest grid, which is negligible compared to the computation on the finest grid.
The developed physical model and solution procedure have been evaluated through comparing some of the simulated results against the measurements. The well-known Pennes (Pennes, 1948) in vivo measurements have been reproduced with the proposed model. The model results have been also compared with real case measurements of temperature distribution within the knee after the ACL surgery. The steady state profiles and the temporal transients have been analysed. For all cases, it is shown that an appropriate modelling of the bio-heat energy transport results in a good agreement with experimental data. The fluid flow solution procedure, however, is designed as a proven approach and tested on several standard benchmark tests, e.g., driven cavity, de Wahl Davis test, etc. The presentation of benchmark tests are beyond the scope of this chapter and are as such omitted. All results related to the bio-heat transfer are presented and elaborated in details in the following sections.
The numerical scheme used to solve the proposed physical model is already well accepted and extensively tested in the past. A detailed convergence of several numerical techniques applied the diffusion problem was published in (Trobec, et al., 2012). It was shown that the FDM has an adequate convergence rate. Based on the obtained results and taking into account that we use the same approach we do not insert the convergence analysis in this chapter.
4. Heat transfer in resting human forearm
The invaluable measurements of the steady-state temperature fields in forearms of subjects that are not anesthetized, published sixty years ago by Pennes (Pennes, 1948), have still not been explored in all details (Wissler, 1998). Pennes’ used a thermocouple probe on a wire that was pulled through the forearm in a straight line (referred to as the measurement path in the rest of the section) in steps of less than 1 cm. At each step, the measurement was read after the probe temperature has stabilized. In Figure 7 the original measurements and the author’s interpolated temperature curves are reproduced.
Other experimental data are available for the temperature of the forearm during an immersion in water at various temperatures, either as it evolves in time (Barcroft and Edholm, 1946) or near its steady-state. Since then, these measurements have been elaborated by others (Ducharme, 1988). In some recent studies, measured lower resolution temperature profiles in the forearm cross-sections are reported, together with estimated values for the in-vivo thermal conductivities of muscle and subcutaneous tissues (Ducharme, 1988). We recreated the conditions similar to those of the Pennes’ measurements and simulated the steady-state temperatures of the proximal forearm. In this section, we consider the effects of blood perfusion, metabolism, position of the arteries, and inhomogeneity of tissues on the steady-state temperature profile of the forearm.
4.1. Geometric model of the forearm
Cross-sections of the VHD male thorax (numbered 1601 – 1800) are used for the forearm model. The right forearm is cropped from those cross-sections as rectangles of 506 (width) 632 (height) pixels. The VHD right forearm is inclined by 32.5° in the x-z plane and by 43.1° in the y-z plane. It is more suitable for the simulation if the z axis of the forearm is parallel to the z axis of the model; therefore, the elements of the model are sampled out of the forearm rectangles at an angle, making the z axes parallel. Finally, a regular grid of size 140116100 elements with spatial step of 1 mm is created. For the 2-D simulations, “model slices” are used – individual x-y planes of the grid, numbered by their position on the z axis.
The 3-D domain represents a section of the human forearm, 50 mm distal to the tip of the elbow and 100 mm long (50 mm from the probe insertion point in both directions). For better visibility of the bones and arteries in the model, muscles, skin, and subcutaneous tissues are shown semi-transparent. In Figure 3 model example slice is shown.
The Pennes’ measurements (Pennes, 1948) were performed exactly “8 cm distal to tip of ulna olecranon” – or 8 cm from the tip of the elbow, with no respect to the length of the forearm. While the lengths of the forearms are not reported, the diameters of the measured forearms differed greatly (7 – 10 cm). Therefore, the measurements were likely performed at different anatomical positions. In case that the 8 cm were strictly respected, the variation in the position of the measuring probe, regarding the forearm anatomy, could be significant − up to 25% of the distance to the elbow.
The diameter of the described geometric model of the forearm at slice 50, which is in the middle of the model, is 12.4 cm. To obtain a closer match to the measured forearms, which were significantly smaller, we scale down the domain size by linearly scaling down the domain elements. In the rest of the section, we use forearm diameter of 8.3 cm, except when explicitly specifying otherwise.
4.2. Numerical setup
We simulate the steady-state temperatures of a resting forearm near the thermo-neutral conditions. It is known that the thermo-neutral conditions of an unclothed arm with skin temperature of approximately 33 °C are obtained at ambient air temperature of 25.2 ± 1.1 °C. The initial conditions of the simulation comprise the temperatures of the modelled tissues except arteries, and are set to 36 °C. The arteries acts as heat generators, and ambient air acts as heat sink. They are set, in accordance to previously measured values, to constant values of 36.8 °C and 26 °C, respectively. During the simulation, the tissue temperatures converge with time towards their stationary values. After the simulating three hours of a resting forearm, the maximal change of temperatures during the last hour is less than 0.01 °C, which confirms that the near-steady state is reached.
Moving air is not simulated because of its significant contribution to the calculation complexity. We simulate the heat flow between the skin surface of the forearm and the ambient air using Robin boundary conditions, based on the continuity of the heat flow perpendicular to the surface of the simulated body part:
where stands for the ambient temperature. The convection coefficient H = 16 W m2/°C is determined by numerical experiments, based on simulation of the thermo-neutral conditions ().
The influence of the parts of the forearm that were not simulated is considered by setting the heat flow to zero at both “ends” of the modelled arm (slice 1 and 100), indicating a thermal equilibrium between the modelled and excluded parts of the forearm.
The blood flow is modelled with the parameters set to:
The simulation setup is presented in Figure 4.
4.3. Results of numerical integration
The simulated steady-state temperature field of the model slice 50 with diameter of 8 cm, after three hours of simulation, is shown in Figure 5a. The impact of the two main arteries, with constant blood temperature, is visible as peaks in the temperature profile at the locations of the arteries. The peaks are marked by two vertical lines. We can also notice that the skin temperature depends on the position because its distance from the arteries varies. The steady-state temperature profile taken from the same slice along the path of the thermocouple measurements are shown in Figure 5a by a black curve and Figure 5b by a thick black curve. The dotted curve shows the analytical solution of the Bio-heat equation (Yamazaki and Sone, 2000) for a cylindrical homogeneous model with blood flow V = 4.7 10-4 1 s-1, hm = 0.62 W/kg, skin temperature TS = 32.5 °C, and with the other constants set the
same as in the simulation. We see that the closed form solution of the homogenous model can reproduce the general trends in the temperature gradients; however, it cannot reproduce the measured results because it is symmetric and uniform throughout the model.
By disabling, one at a time, the metabolism, main arteries, and blood perfusion in the simulation, we obtain the last three temperature profiles shown in Figure 5b. The blue curve shows the simulated temperature without metabolism. We can see that the metabolism minimally contributes to the final temperature. The red curve was obtained by not keeping the arteries at a constant temperature in the simulation, thus cancelling their heat generating role. The impact of the arteries not heating up the surrounding tissue makes this temperature profile stand out among the simulated profiles. It shows lower deep muscle temperature than normal, which is more symmetric. The green curve, obtained by the simulation with the arteries at a constant temperature and disabled blood perfusion, confirms that the contribution of the blood perfusion is essential. The average temperature of this profile is significantly lower with pronounced peaks corresponding to the two main arteries.
4.4. Sensitivity analysis
The simulated steady-state solution is evaluated regarding its sensitivity. The input parameters are varied within selected ranges and the variations in the resulting temperature profiles are analysed. In Figure 6 the temperatures of the forearm after three hours of simulated resting are shown.
The temperatures are taken from the model slice 50 along the measurement path. All results are calculated based on the same geometric model with varying forearm diameter, thermal diffusivity, blood flow, and metabolism, all for +20% (dotted line) and -20% (dashed line) of their nominal values (solid line). The nominal values are taken from Table 1. The nominal value for the forearm diameter was determined to be 10 cm.
The most significant impact on the temperature profiles comes from varying the forearm dimension (Figure 6a). Larger arms result in a plateau slightly above 36.3 °C, while smaller arms are cooled more intensively, reducing the maximal temperature to 35.6 °C and with larger temperature gradients in the superficial regions of the forearm. Changes in the thermal diffusivity (Figure 6b) have a much smaller impact on the temperature profiles than the forearm dimension. The effect is also non uniform across the profile, but is rather equalizing the superficial and the deep muscle temperatures. The impact of the blood flow (Figure 6c) is again uniform across the profile, with higher forearm temperatures when the blood flow is larger. The impact of the metabolism (Figure 6d) is analogous to the one of the blood flow − greater metabolism results in higher forearm internal temperatures − but significantly smaller in scale.
4.5. Reproducibility of measurements
We found out that the most important parameter that influences the simulation results is the dimension of the simulated forearm. The impact of the diameter is evident from Figure 6a. Other thermal parameters, like blood perfusion and structural details, have much smaller
impact on the variations in the simulated temperatures and have been therefore considered in the same way in all individual cases. Only two main vessels have been considered, which are anatomically quite uniform at most subjects. The remaining smaller vessels are encompassed through the blood perfusion as an average temperature regulator. The temperature gradients are very pronounced and therefore the location of the temperature probe provides the second most important contribution in the analysis of the temperature profiles. Because the anatomical probe positions are not known exactly, we try to find them to see if they can explain the unexpected variations in the measured temperature profiles.
We can deduce from Pennes work further sources of measuring errors. The positions of the measuring points and the reading of the temperature values have been with limited accuracy, which could result in errors of 0.1 °C. The dimensions of the thermocouples were up to 0.2 mm and hence, in the areas with steepest temperature gradient, near the forearm surface, an additional error could be introduced.
To reproduce the measurements, we first simulate the resting forearms with the same diameter and under the same conditions as in the real measurements. Using a matching 3-D simulation and a temperature profile measurement, we perform a search for the probe trajectory on which the simulated temperature profile fits the measured temperature profile the most. We do so by defining a fitness function as the root mean square (RMS) error between the simulated temperature profile and the measured referential profile, weighted with Gaussian function to diminish the effect of greater inaccuracies in the measurements near the forearm surface. Then, we use the fitness function in Differential Evolution (Price and Storn, 2008) – a stochastic optimization procedure – to find a probe trajectory that fits the real measurement well. The similarity between the measured and the simulated profiles for all analysed measurements is given by the Pearson correlation coefficient in Table 2. The simulated paths of the measuring probes, which provide the highest correlation between the simulated and the measured profiles, are also determined in Table 2 by two offsets and two angles to the measurement path. They can be obtained by sequentially applying translations in y and x directions and rotations in x-y and x-z planes to the trajectory of the thermocouple probe.
|y offset [mm]||-2.6||2.7||4.5||1.8||5.7||-0.4||2.4||1.2|
|z offset [mm]||-21.4||-32.7||-30.3||-9.3||-28.8||-9.0||-36.3||-4.4|
|Angle in x-y plane [°]||0||-0.4||-1.6||-0.4||0||0.4||0||3.5|
|Angle in x-z plane [°]||5.5||5.9||0||0.8||7.8||4.3||3.9||4.3|
Two typical simulated temperature profiles for measurements 4 and 9, together with points from the original measurements, are shown in Figure 7. The central peak in the profile of measurement 4 is expected and can also be predicted by the theory. The temperature profile of measurement 9 with dual maxima can be explained by the probe moving close to and about the same distance from both arteries. The left and right shifted maxima visible in some other measurements can be explained by the path of the temperature probe being closer to one of the arteries. The difference in the forearm dimensions results in significantly different temperatures in the central part of the forearm. As visible from Figure 7, larger forearms induce higher temperatures. We are able to reproduce the original measurements using the same geometric model for all subjects (just scaled to the measured diameter), and by introducing minor changes in the position of the measuring paths.
Instead of first simulating and then evaluating the simulated results, our study is designed in reverse direction. First, the well-known in vivo measurements of the steady-state temperatures of the proximal forearm are taken from the famous Pennes' publication; then, the experimental environment is reproduced; next, the bio heat transfer is modelled, and lastly, the measurements are simulated. The primary goal was to explain, by computer simulation, the variations in the Pennes’ in vivo measurements. We have known that the measurements vary significantly between subjects and between forearm cross-sections. Several attempts have been made in the past to average the measured steady-state temperatures in order to smooth fluctuations and obtain a standardized response.
The anatomical positioning of the measuring probes is assumed to vary because of the different forearm sizes of the subjects; therefore, their possible paths have been reconstructed by searching for the simulated paths that result in the best agreement between the simulated and the measured temperature fields. Our simulated results suggest that the smoothing approach is inappropriate, and that the fluctuations are natural result and the consequence of a complex interplay between the position of the measuring probe, the anatomical position of the main arteries, the dimensions of the forearm, the blood flow, the inhomogeneity of tissues, and the environmental temperature.
Several innovative approaches have been introduced in the simulation method. A procedure for spatial model generation based on digitized slice data has been developed. A mathematical model and a 3-D computer simulation program have been implemented for the simulation of steady-state temperatures. The heat transfer in the model of non-homogenous tissue was modelled with a well-known Bio-heat equation. The heat production by tissue metabolism was modelled using the Q10 rule. The heat exchange between the blood and the tissue was modelled as a function of the local temperature and the regional blood flow that was modelled as an exponential function fitted to experimental data. The regional blood flow was assumed to depend on the regional temperature in the same way as the mean blood flow rate depends on the skin temperature. Such an approach produced reasonable results. Nevertheless, it needs to be supported by further physiological research. The developed simulation, however, is sufficiently robust to incorporate arbitrary model functions for blood flow and metabolism. The stability of the method was confirmed by varying the simulation parameters, the initial and boundary values, and the model shape, with subsequent analysis of the results.
We have confirmed that blood flow has a significant and complex impact on the steady-state temperatures of the forearm. Lower blood flow results in lower temperature of the central part of the forearm. When the effect of blood flow on the temperature profile is increased by disabling the direct effect of the main arteries, the resulting temperature profiles match very well to the analytically produced profiles from the literature, but not well with the experimental measurements. We have shown that only when the effect of both sources, the heat distributed via blood flow and the heat emitted directly from the arteries, are considered, the measured temperature profiles can be fully explained. Furthermore, we have revealed the dependence of the temperature profile shape on the dimensions of the forearm. Smaller forearms are more likely to feature off-centre peaks in the vicinity of the main arteries than larger forearms. In larger forearms, a central plateau may be formed, masking all other potential peaks. The thermodynamic constants of the tissue and the tissue metabolism were confirmed to have a relatively minor impact on the temperature field.
Our simulation methodology has several limitations. The anatomy of the forearm differs between individuals and changes in time. Consequently, the simulation on one generalized model can differ from the experimental measurements. Minor errors in tissue segmentation and inaccurate thermodynamic constants could produce errors in the simulated results. The possible influence of different mechanisms of thermoregulation on blood perfusion has not been simulated because we assumed that their effects are stable when the forearm temperature is in its steady-state. We confirmed that the above essential findings are in close agreement with in vivo measurements, particularly if the inaccuracy of the position of the measuring path is accepted.
5. Simulation of a human knee exposed to topical cooling after ACL surgery
Topical cooling after surgery or after knee injuries is often performed using gel-packs filled with gel that can be cooled well below 0 °C, yet still remaining soft and flexible; or alternatively, with cryo-cuffs, which provide a constant temperature maintained by circulating cooled liquid (Hubbard and Denegar, 2004). Desired cooling conditions are obtained by differently shaped gel-packs or cryo-cuffs, often influenced by environmental conditions and the temperatures of internal tissues (Barber, 2000). However, different tissues around the knee joint have different physical and thermodynamic properties, and respond diversely to the cooling. The knee temperature is influenced by the muscle metabolism and the circulating blood, both of which depend on the tissue temperature. Therefore, the real cooling efficacy varies in space and time in different parts of the knee. The effects of cooling on the tissue metabolism and the response of inflamed or traumatized tissue have been researched and have shown beneficial effects (Ho, et al., 1994, McGuire and Hendricks, 2006). There is no generally accepted doctrine about topical cooling. Therefore, further investigation of cooling effects and different cooling regimes are needed, but these are difficult to perform by experiments alone.
In this section, we are interested in simulating possible modes of knee cooling after injury or surgery (Grana, 1994). We focused on the comparison of two different methods of topical cooling. First is the cooling with gel-pack filled with refrigerated gel (McMaster, et al., 1978). The gel-pack is exposed to ambient temperature and therefore becomes less and less effective with time. Second is the cooling with a cryo-cuff filled with a liquid. The temperature of the liquid is maintained constant by an external cooling device (Shelbourne and Nitz, 1990).
Lowering the tissue temperature reduces the need for medicaments and shortens the rehabilitation period (Knight, 1985). It is not clear, however, whether different tissue regions in the knee should be cooled uniformly. Some recent findings put the benefit of topical cooling under question (Hubbard and Denegar, 2004, Knight, 1985, McGuire and Hendricks, 2006). The aim of our work is to simulate topical cooling of the knee after injury or surgery, and to calculate and show the development of temperature distribution in all tissues of the knee region.
5.1. Geometric Model of the Knee
The knee area is cropped from the VHD male cross-section photographs number 2200 to 2390 as rectangles of 550 (width) 610 (height) pixels. Cropped VHD slice 2301 is shown on the left in Figure 8.
The 3-D model is shown in Figure 9. Skin, joint liquid and subcutaneous tissues are not shown. The artero-lateral quadrant is removed for better visibility into the central knee region.
Some surrounding space for the protective bandage, the gel-pack, the blanket, and the ambient air is added around the 3-D knee model. The simulation environment is imitated by an isolated cube composed of resulting in roughly voxels independently characterized by thermodynamic properties and initial temperatures. Each voxel has a volume of 1 mm3. The knee is covered by a 2 mm thick protective bandage and embraced by a 12 mm thick cooling layer (gel-pack or cryo-cuff). Additionally, the cooling layer is covered with a 5 mm thick isolating blanket to reduce the convection and slow down the warming from outside. The protective bandage, the cooling layer, and the isolating blanket are inserted into the model automatically by a computer program.
5.2. Numerical setup
The boundary layers of the simulated box are held at constant initial temperatures to mimic the ambient air. The influence of the rest of the leg, which is not simulated, is managed by setting the temperature flux at the boundaries. The simulation can be carried out in 2-D with the heat flux in the axial direction set to zero. In this way, an infinitely long "knee", with homogeneous structure in the axial dimension, can be simulated. Simulation in 3-D is performed on the described model with all boundaries set to values that are similar to those of the measuring conditions. The simulation setup is presented in Figure 4.
The blood perfusion in bones, cartilage, and ligaments from the knee joint are assumed to be very small. Thus, for these tissues, ten times smaller numerical values for are taken than for the other tissues. The blood flow is modelled using the following parameters:
As in the case of the blood flow, the metabolism of bones, cartilage and ligaments is assumed to be very small.
|subcutaneous tissue||35.6||venous blood||36|
5.3. Results of numerical integration
First, we simulate the steady-state temperatures of a resting knee in thermo-neutral conditions (Ducharme, et al., 1991) at ambient air temperature of 27, which is equivalent to the ambient water, e.g., a water bath at 33. A steady-state is reached after three hours of simulation, with the maximal change in the temperature near the end of the third hour being less than 0.01.
At that state, the simulated temperatures can be recorded at arbitrary positions. For example, the temperature field over the model slice is shown with a 2-D surface in Figure 11. For better visibility of the temperature field, the point is on the right and the frontal part of the knee with the patella on the left. The temperature along the transverse axis for is used in a later analysis. The obtained steady-state is used as an initial condition in all our subsequent simulations. The tissues nearer the patella, toward the surface, are seen to be colder than the internal part. The impact of the main knee artery with its constant blood temperature is visible as a peak in the temperature field. The location of the bones is also evident as shallow depressions in the temperature field, mainly as a result of the lower blood flow in the bones.
In Figure 12 the solid curve shows the simulated steady-state temperature obtained with the 3-D simulation. The dotted curve is obtained under the same conditions with a 2-D simulation and zero flux in the third dimension (infinite knee). The difference between the simulations is in the range of 0.5. It can be explained by underestimating the amount of muscle tissue in the whole knee region. Slice 102 is taken, namely, from the central part of the knee with less muscle, which could contribute additional heat through the blood flow. However, the shapes of the profiles are similar. Therefore, for the initial analysis we use 2-D simulation to reduce the simulation time.
We can observe and analyse arbitrary simulated points or regions of the domain. For example, in Figure 13, the steady-state temperature profiles of the knee are shown for temperatures on the transverse knee axis and model slice for, and. It is clear that there can be significant differences in the temperature, even for analysed points as close as 10 mm. The effect of altering the observed position in other directions is similar to that seen in Figure 13.
The stability of the simulated steady-state solution is evaluated by varying the simulation parameters. Input parameters are varied within selected ranges and the variations in the solution were analysed. In Figure 14 the simulated steady-state temperatures from model slice 102 on the transversal axis for are shown for various knee dimensions, diffusion constants, blood flow and metabolism, all differing for +20 % and -20 % of their nominal values.
The most important impact on the temperature profiles seen in Figure 14 arises from varying the knee dimension (Figure 14a). Larger knees result in a temperature plateau slightly above 36.7, while smaller knees are cooled more intensively to the central temperature of about 36, with larger temperature gradients in the superficial regions of the knee.
Changes in thermal diffusivity (Figure 14b) have a smaller impact on the temperature profiles than changes in the knee dimensions, and in the opposite direction, i.e., larger diffusivity constants result in lower central temperatures.
The impact of blood flow (Figure 14c) is similar to that of knee dimensions, with higher internal temperatures arising from greater blood flow. However, an important difference is that the shapes of the temperature profiles remain unchanged with variations of the blood flow.
The impact of metabolism (Figure 14d) is analogous to that of the blood flow – increased metabolism results in higher internal temperatures, but its impact is so small that it can be neglected in our experiments.
Simulated steady-state temperatures of a resting knee under thermo-neutral conditions (TN) are used as an initial condition for all further simulations. First, we simulated a naked knee in steady-state at ambient air temperature of 25. Then, we simulate a two-hour period of an arthroscopic operation, during which the knee joint is washed out by sterilized water at 22, and therefore cooled. During the following two-hour period we simulate the temperature evolution in the operated knee while it is resting and is covered by a blanket, and therefore warming. Finally, we simulate a subsequent two-hour postoperative topical cooling.
5.4. Washing out during arthroscopy
During arthroscopic reconstruction of ligaments, the central part of the knee is washed by sterilized water at 22. The water circulates in the space around the femoral intercondylar notch, normally filled by the joint liquid. The initial temperatures are taken from the steady-state of a naked knee at an ambient temperature of 25. The temperature of the joint liquid is fixed at 22. The temperature profile after two hours of washing out are shown in Figure 15 for the same plane as before; that is, along the axis for, but for the 7 mm higher slice at, because this is nearer the actual position of our measuring probes for validation of the simulated results. The internal knee temperature decreases significantly, maximally by more than 14in places that have direct contact with the washing water, as can be seen from the temperature profiles in Figure 15.
5.5. Resting after arthroscopy
Next, we simulate the two hour period immediately after surgery. The knee is resting and is covered with a blanket in a room temperature of 25. The evolution of the temperatures in voxel, which is at the level of the femoral intercondylar notch in the central part of the knee, and in voxel, which is nearer to the knee surface, 1 cm below the skin in the subcutaneous tissue, are shown in Figure 16. Note that arbitrary voxels can be selected for the analysis. The knee is initially colder in the central region because of the previous washing out with cold water. During resting, its temperature increases and approaches the steady-state, with colder regions nearer the skin.
5.6. Postoperative topical cooling
Finally, we simulate postoperative topical cooling with two different cooling methods, gel-pack and cryo-cuff. In both cases, the knee is bound with a protective blanket surrounded by ambient air at 25.
In Figure 17 the temperatures of voxels and are shown for the two-hour simulated period, for cooling with a gel-pack (initial temperature 0) and a cryo-cuff (water with constant temperature of 15). The effects of the two methods on the knee temperatures are quite different.
When cooling with gel-pack, the temperature of the voxel in the central knee region initially slightly increases because of the arterial blood perfusion and metabolism, and the weak influence of the initial cooling. After 5 minutes, the temperature of the voxel starts to decrease. However, after 40 minutes, the gel-pack has received enough heat from
the knee surface and the ambient air, allowing the inner knee temperature to increase during the second part of the cooling period. For the voxel in the subcutaneous tissue, 10 mm below the skin, the temperature first sharply decreases; the effectiveness of the gel-pack then becomes weaker, the voxel starts warming up, and after 120 minutes reaches almost 36.0.
The cooling with a cryo-cuff is found to be more effective. It induces lower tissue temperatures, even if the temperature of the cooling liquid is as high as 15. In the initial phase, both voxels experience the same cooling rate as the one with the gel-pack. However, there is no subsequent increase in temperature because the cryo-cuff is a constant heat sink which gradually cools the knee. After two hours of cooling, the near surface voxel reaches a temperature of 27. In the same way, but with smaller intensity, the inner voxel is cooled to 33.
In the case of topical cooling with gel-pack, the heat of fusion should also be simulated, which is necessary for the transition between aggregate states. Obviously, for a crushed gel-pack, a significant part of the heat is needed for such a transition, which would prolong the effective cooling time. This phenomenon has not been incorporated in our mathematical model. Instead, we increase the heat capacity of the gel-pack to recognize and account for this behaviour. The simulated results show that the topical cooling with a cryo-cuff provides more constant lowering of the temperatures in the whole region of the knee. Cooling with gel-packs is less stable; consequently, they should be changed every half an hour in order to be effective.
In Figure 18 temperature profiles after one hour of simulated cooling with gel-pack and cryo-cuff are shown for a cross-section from the patella to the lateral side of the knee on our standard axis, i.e., along the (anteroposterior) axis for on the model slice.
After one hour of simulated cooling, the gradients in the temperature profiles were much more pronounced than in the initial state. When cooling with gel-pack, the temperature of the outer knee layers at the skin level remained cooled to 32, and in the centre of the knee to 36. The peaks in the tissue temperatures around result from the simulated heat conduction from the middle popliteal artery. Significantly lower temperatures are observed for the cooling with cryo-cuff, even though its constant temperature was as high as 15.
Given the above results, it would be interesting to test the effectiveness of a simple cooling with ambient air. The knee would remain uncovered and exposed to the ambient air temperature of approximately 20. We expect, from the results under thermo-neutral conditions, that the skin temperature would be about 25, which could lower the temperatures inside the knee. This could be tested preliminarily using the proposed simulation method.
5.7. Validation of results
To evaluate the simulated results, we made two control measurements of the knee temperature following surgery in a room with constant ambient temperature of 25. We measured the temperature of the knee covered with a blanket during the two-hour resting period immediately after surgery; then, during the next two-hour period in which topical cooling with a gel-pack was applied.
Two small thermistors were placed into the knee in thin sterile tubes (Foley-catheter with temperature sensor, 3 mm, Ch8-thermistor; Curity, Degania Silicone Ltd., Degania Bet, Israel). Similar tubes, without thermistors, are ordinarily inserted for wound drainage following surgery. The thermistors were connected to a registration device for continuous measurement with a sample rate of 0.1 Hz and resolution of 0.01.
The first thermistor was placed in the centre of the knee near the voxel, and the second approximately 1 cm below the skin in the subcutaneous tissue, near the voxel. The measurements were approved by the Slovenian State Medical Ethics Committee and the patient gave written informed consent prior to participation.
In Figure 19 the simulated and measured temperatures are shown for the two-hour resting period after washing out during arthroscopic surgery. The knee was wrapped in a protective blanket at ambient temperature of 25under the same conditions as in the simulation.
The simulated temperature evolution for the resting period shows very good agreement with the measured values for the test point in the central part of the knee. However, the simulated rate of cooling in point in the subcutaneous tissue was much smaller in the initial phase than the one obtained by the measurements. One of the possible reasons is the fact that we did not simulate the cold washing out inlet that also cooled the surrounding tissue from the skin to the central part of the knee. The subcutaneous thermistor was placed in such cooled environment which could then exhibit faster warming than in our simulation. In Figure 20 the simulated and measured temperatures are shown for two hours of cooling with a gel-pack, which follows immediately after the resting period. The initial temperature of the gel-pack was 0and the ambient temperature 25. The knee was bound with elastic bandages approximately 2 mm thick, surrounded with fixed gel-packs and wrapped in a protective blanket.
The simulated evolution of temperature for the cooling period shows good agreement with the measured values for both test points. The simulated rates of warming in the second hour are slightly greater than by measurement. One of the possible reasons is the incomplete mathematical model that does not include the heat of fusion for the gel-pack. In fact, all the thermodynamic characteristics of the gel were not available, and we just took some approximate values provided by the supplier. Another possible reason is inaccurate measurement. We did not collect detailed data during the measurements, e.g., the wetness of the protective bandage, which could have a significant impact on the cooling intensity
We demonstrated a practical application of the simulation program for two different methods of topical postoperative cooling of a knee. We simulated topical cooling of a knee with a gel-pack. The inner knee tissues reached their lowest temperature in 40 minutes. For continued effective cooling, the gel-pack has to be replaced. The topical cooling with a cryo-cuff was more effective. We simulated situations with relatively small cooling rates to analyse the small influences of the blood flow and metabolism. The thickness of a protective bandage and isolating blanket, or their thermal conductivity, together with the cooling temperature, can be used to regulate the cooling intensity. The results have been validated by experimental measurements of knee temperatures. The simulation results confirm that the model and methodology used are appropriate for the thermal simulation of bio-tissues.
We have shown that the blood flow has a significant and complex impact on the stationary state temperatures, and the gradient of temperature change in the subcutaneous tissue and in the tissues nearer the central part of the knee. Lower blood flow results in linear temperature profiles with larger gradients, for example in bones. The dimensions of the knee are also important factors influencing the temperature distributions and gradients. Temperatures in smaller knees will differ more from the arterial blood temperature than those in larger knees. At the same time, the temperature gradient will be much greater in smaller knees. Thermal constants and metabolism have a relatively minor impact on the temperature field.
The simulated results were visualized and compared with measured values. Good agreement was obtained, leading to the conclusion that the model and method used in our simulation are appropriate for such medical simulations. Although there are not many studies of knee temperatures measured in vivo, and the measuring conditions are often not described in sufficient detail, we have also run our simulation software with the initial conditions described in some published measurements. In (Martin, et al., 2001) the temperature in the lateral gutter of the knee decreased by 4, one hour after knee arthroscopy. Similar values have been obtained by our simulation program.
The method described has several limitations. The 3-D knee model used in our simulation is not complete, as only a small part of the leg above and below the knee was included. The remaining part of the leg that was not included was compensated by a constant flux in boundary conditions. Spatial models differ with different persons and with time, and consequently, the simulated results can differ. Minor errors in tissue segmentation and inaccurate thermodynamic constants could also produce small errors in the simulated results. Moving air was compensated by an artificial thin layer of air with non-constant temperature. Incorporation of a fluid-flow model would be needed for even more accurate results. The possible influence of blood perfusion by different regulatory mechanisms has not been simulated. Personal regulatory mechanisms have not been included in the simulation model, but could easily be incorporated. Listed limitations could have some impact on the simulated temperatures, but the essential findings agree remarkably well with existing experimental in vivo measurements.
6. Simulation of human heart cooling during surgery
During longer surgery the human body and the heart have to be managed appropriately in order to slow down their vital functions. One of the options for lowering metabolic requirements is the cooling of the body and the heart, e.g., by whole body hypothermia or infusion of a cold solution in coronary arteries (cardioplegia). For even better cardiac cooling, a method of topical cooling of the heart is sometimes applied (Olin and Huljebrant, 1993) by using packs filled with cooled gel or ice slush.
The efficiency of cooling is reflected in the temperature distribution within the heart, which is hard to determine. In vivo temperature measurements are invasive and, if applicable, limited to a few test points. These issues can be circumvented to some extent by physical modelling and computer simulation. Appropriate simulation tools provide insight in the temperature distribution through the whole heart volume and its neighbourhood, and therefore, various preliminary evaluations and analyses of different cooling regimes can be assessed.
Initial results on the simulation of heart cooling based only on diffusion have been reported in (Trunk, et al., 2003). In this section, we present simulation results obtained through much more complex physical model, incorporating, beside diffusion, also advection, bio-heat sources and fluid flow. Such a model stands for more realistic variant of already developed simulation tools (Šterk and Trobec, 2005, Trobec, et al., 1998).
Our final goal is to implement an accurate and reliable simulation framework that enables simulations of various diagnostics or therapeutic medical procedures. Simulation of different protocols of myocardial protection is taken as a test-case because they incorporate most of the physical phenomena during surgery that are related to the bio-heat transport. In this final section we focus on topical cooling as an adjunct to cardioplegia in classical non-beating heart operation (Allen, et al., 1992). During cardiac surgery, local hypothermia is reported to be successful in slowing down the tissue metabolism (Birdi, et al., 1999). This is achieved by infusing a cold cardioplegic solution into the coronary vessel system. Alternatively, topical cooling is applied with a cold saline solution, sometimes mixed with ice slush, which is distributed around the heart to further lower the temperature of the heart tissue. In non-perfused tissues, ischemia could provide deleterious effects on the cells, finally producing necrosis. Lowering of tissue temperature lengthens the time in which the ischaemic damage to the cells is reversible, which can prolong the safe time available for cardiac surgery. To provide adequate protection to the whole heart tissues, the heart has to be cooled down as uniformly as possible (Gonzales, et al., 1985). The anterior part of the heart, especially the right ventricular wall is likely to remain outside the topical cooling solution during the operation, because it represents the uppermost part of the heart in supine patient. The importance of uniform temperature distribution in heart tissues comes also from the fact that areas with different temperatures could change myocardial contraction velocity (Riishede and Nielsen-Kudsk, 1990) and trigger cardiac arrhythmias. Some previous computer simulations have already shown irregular temperature distribution in the human heart during topical cooling (Geršak, et al., 1997). Opinions about the cryo-therapeutic effects on the heart during cardiac surgery and during post-surgery rehabilitation period are still controversial. Furthermore, recent findings put the real benefit of topical cooling under question (Cordell, 1995). The optimal temperature at which the heart tissues should be cooled down and the time intervals of optimal cooling are still not clearly defined, probably because of the lack of in vivo measurements and difficult detailed post-surgery patients screening.
The aim of our work is to accurately simulate topical cooling during induced cardiac arrest in a typical heart operation. We want to use 3-D geometric model with all details, including coronary vessels. The desired result is the temperature distribution in the heart tissues as evolving in the time of surgery. The topical cooling is mimicked by submerging the heart in a cooling liquid. We are particularly interested in the influence of the liquid dynamics on the cooling performance.
Detailed quantitative validation of the results is not the main goal of our investigation. We are more concerned about the adequacy analysis of the presented physical model for simulation of heat transfer phenomenon during open heart surgery. For the sake of simplicity and transparency of results, we assume some simplifications, which can be declared as limitations of our study.
First, we simulate only the topical cooling and do not take into account the heart cardioplegia and the whole body hypothermia. It allows us to be concentrated only on a single cooling modality. The impacts of the remaining cooling methods remain to be studied in the future. The whole body hypothermia could be simulated by setting the domain boundary conditions to lower temperatures. The cardioplegia could be approximated by setting the temperature of blood in the coronary vessels to a constant value with Dirichlet boundary conditions on the vessel walls, assuming well mixed blood with a significant flow. Next, we simplify the domain boundary conditions to the Dirichlet, although the cooling would affect the body as well; however, the body temperature is assumed to remain constant through physiological regulatory loops. Next, because we assume well mixed air we neglect the natural convection of air in the operating room. The ambient temperature is set to constant value. Finally, blood in heart chambers is assumed to be steady.
6.1. Geometric model of a heart
VHD male cross-section photographs, numbered 1350 to 1505, are selected for the construction of the geometric model. Z axis of the VHD is parallel to the heart axis, pointing approximately from the apex towards the atria and the base of the heart. The selected cross-sections are cropped to the size of 512 (width)512 (height) pixels. Different tissues are classified (segmented) manually with an assistance of expert anatomist. The space occupied by voxels not belonging to the heart tissues or major vessels (aorta, vena cava, pulmonary vessels) are replaced either by voxels of air, water or the box wall. The final 3-D model used in the simulation is shown in Figure 22.
During the tissue segmentation, a few problems that could deteriorate the desired spatial accuracy of the heart model have to be resolved. The details of the heart structures have to be preserved as much as possible. The accuracy of the simulations relies on the resolution of the geometric model, which serves as a substrate for the simulation. For example, it is impossible to distinguish between the epicardium and pericardium on the VHD photographs; therefore, these two tissues are treated as a single tissue. The wall of the large vessels originating in the heart chambers is also not distinguishable from the neighbouring pericardium. The myocardial tissue contains several layers of muscular fibres that run in different directions. This could influence the heat propagation through the tissue, if different directions are examined; however, these layers could not be distinguished on the original cross-sections. Thus, the myocardium is classified as a homogenous mass of muscular tissue. The left ventricle, usually much larger in a heart in function, is quite small in the VHD dataset, evidently contracted and being in a systole like phase. On the contrary, the right ventricle is huge and elliptically shaped with an extremely thin wall. The endocardium, which normally lines the entire interior of the heart chambers, is microscopically thin, so it could not be differentiated in our model.
The border between two tissues was sometimes difficult to be determined from the VHD photographs because of the similar pixel colours. In such cases, the borders were drawn subjectively according to the anatomical knowledge. The majority of problems are detected in the area of the pulmonary and other large vessels at the base of the heart. These vessels can be identified only with concurrent examination of neighbouring slices and subsequent merging of multiple-slice information. Also, the pericardium is discontinued in some places, particularly on the posterior surface of the heart, because it is not distinguishable from the surrounding tissues on the posterior side of the heart. In such cases, the borders are drawn subjectively with some variations, taking into account a few consecutive neighbouring cross-sections. The small coronary vessels are also hard to be identified. The left and right coronary arteries, the circumflex artery and the largest veins are visible and one can follow such vessels from their origin to the periphery. On the periphery, coronary arteries and veins sometimes lie in close proximity and look like a single large vessel. Some vessels do not have a uniform course in the 3-D view. There are discontinuations that occurred because the vessels are not recognized correctly on every cross-section. Some further artefacts are also detected, like spots painted on some cross-sections as vessels, but actually not corresponding to vessel tissues. This results in a single pixel or a group of pixels apart from the coronary vessels. Therefore, the coronary vessel system was re-segmented again in a separate phase.
The segmentation errors of the generated spatial heart model have to be corrected, either manually or automatically by a computer program. A spatial editor program is implemented to allow the user to easily spot and correct the smaller failures in the generated 3-D models. The spatial editor is also used for several automatic modifications of the 3D-model, e.g., a uniform reduction of the model size by a user supplied factor n. Another function is designed for finding small openings or discontinuities in a selected tissue (exposing other tissues) and repairing them by adding small portions of missing tissue onto adequate places of the cross-sections. This function was mostly used to fill in small holes in the pericardium.
6.2. Numerical setup
The thermal conditions during cardiac surgery are imitated with an isolated cube large enough to contain the simulated heart model. The computational domain is shown Figure 22. The lower part of the cube represents a cooling container filled with a cooling liquid that is initially set to a temperature of 0.2 °C, assumed to be cooled via ice slush. The container walls are kept at a constant temperature of 36.5 °C, simulating the body temperature. The initial temperature of the heart is also set to 36.5 °C. The upper part of the cube represents the open thorax surrounded by air at room temperature. The surrounding air is assumed to be well mixed and, as such, is kept on a constant temperature of 20 °C. The cube is composed of voxels of 1 mm3 representing a 3-D box with size x = 163 mm, y = 196 mm, z = 163 mm independently characterized by specific thermodynamic constants, defined earlier in Table 1 (Section 2) and Table 3 (Section 5).
The blood flow in a Bio-heat equation (7) is modelled using the following parameters:
The cooling liquid is assumed to be water with thermal expansion coefficient and viscosity. The heat transfer in water is simulated without bio heat sources (10).
The present analysis is focused on the detailed assessment of the impact of fluid flow on the cooling process of the heart. Note that in cases described in previous sections the cooling is modelled only as a diffusive process.
The heat flux from the boundaries and the heart warms the surrounding water. The temperature dependent density governed by the Boussinesq approximation induces non-uniformly distributed volumetric body force in the Navier Stokes momentum equation, thus, a velocity field develops. The phenomenon is also known as natural convection (Amimul, 2011). The moving fluid extensively affects the cooling process through the advection (second term in the right side of equation (7)) term. In this section, we show the importance of heat advection. Simulation time of 10 minutes is considered, as the characteristics relevant to our discussions already develop in that time. The results are presented in terms of velocity profiles and temperature contour plots on x-y cross-section at z = 80 mm in four characteristic time instants t = [60, 240, 420, 600] s, and in terms of cross-lines in x-y plane for four different cross-sections z = [30, 60, 90, 120] mm at the end of simulation (600 s).
We compare two different scenarios; first, with fluid flow and second, without fluid flow. The main difference is in the treatment of the cooling liquid. The first case represents complex simulation with diffusion, advection and bio-heat sources, while the second case is similar to previously presented simulations of arm and knee. Such a cooling method with the "still" water can be approximated by cold gel packs or ice slash packs, where the flow is not possible.
In Figure 23 the temporal development of the setup with fluid flow is presented. The effects of advection, due to the natural convection, are clearly visible. Near the boundaries, the water is warmer as it is constantly heated through the domain walls and the heart. The heat flux generates upward fluid flow that eventually turns inwards as it reaches the boundary of the free liquid media. The consequent fluid flow significantly increases the heat transport from the domain boundaries, as well as from the heart into the cooling water. Note that in the bottom of the domain, local flow instabilities can be identified. The local bursts of upward flow result from instable buoyant force distribution. The warmer fluid pressures upward and only a small perturbation is needed to start a local upward burst of fluid, resulting in a local vortex. The induced circulation significantly increases the heat transfer between the walls and the heart.
In Figure 24 the simulation without fluid flow is presented. There is no advection and the heat extracted from the heart is transported only by means of diffusion, which is much slower in comparison with the advection. There is significantly less heat transported from the walls and, therefore, more heat is extracted from the heart before water reaches steady state temperature. The steady state, although not simulated in the present work, is also reached considerably later if the advection is not in place, as the heat flow is much slower (see water temperature on the last figure panel - 600 s). In the case with diffusion only, the water plays an additional role of an insulator between the heart and the boundary walls.
Finally, the whole domain is examined by considering temperature cross-section plots after the simulated 600 s. Three cross lines in x-y plane at z = [30, 60, 90, 120] mm for both simulated cases (with and without natural convection) are shown in Figure 25. These additional representations confirm the earlier conclusions. The water retains its initial temperature for much longer period of time in the case without flow. The immediate effect, visible after ten minutes of simulation, is lower heart temperature in the case without advection (curves marked by circles).
We confirmed that the fluid flow plays a crucial role in the process of topical cooling of the heart. The natural convection in the cooling liquid increases the dynamics of the heat transport. The advection increases the flux from the domain boundaries, as well as from the
heart walls. In the presented cases, the walls are kept at a constant temperature. Hence, the increased flux primarily affects the heating of the water. Consequently, the cooling process is less effective. If the domain boundaries would be kept at lower temperatures or kept isolated, e.g., by heart transplantation, the advection induced cooling would be much more pronounced. Alternatively, the water could be externally cooled or completely replaced after certain period of time to keep the cooling process more effective; however. The above cooling options are beyond the scope of our work and remain to be further investigated. In the case of omitted advection, the impact of cold "still water" lasts longer, resulting in lower simulated temperatures of the heart. We expect that new standards for the cooling protocols have to be developed in the future, supported also with the results and methodology presented in our work.
7. Discussion and conclusions
We developed and implemented a physical model for the solution of the Bio-heat equation described in Section 2. The model incorporates heat diffusion, heat advection, heat transfer through blood perfusion and heat generation by tissue metabolism or other bio sources. A corresponding mathematical model and a computer simulation program have been implemented. Some important extensions have been introduced, in particular, an inhomogeneous spatial model composed of tissues with different characteristics, and modelling of the blood perfusion and heat sources as functions of the surrounding tissue temperature.
The heat produced by tissue metabolism was modelled using the Q10 rule, while the heat exchange between the blood and the tissue was modelled as a function of the regional blood flow and the local temperature. The regional blood flow was modelled as an exponential function fitted to experimental data. It was assumed to depend on the regional temperature in the same way as the mean blood flow rate depends on the skin temperature. Such an approach produced reasonable results. Nevertheless, it needs to be supported by further physiological research. The developed simulation, however, is sufficiently robust to incorporate arbitrary model functions for blood flow and metabolism. The stability of the method was confirmed by varying the simulation parameters, the initial and boundary values, and the model shape, with subsequent analysis of the results.
The construction of geometrical model, described in Section 3, was based on VHD slices. The models were composed of equal voxels, which results in a regular spatial discretization. Consequently, an explicit finite difference scheme has been developed and optimized for this purpose. A general methodology for the development of 3-D models, based on digitalized slice data, is described and applied for generating our spatial geometric models.
We have shown that the computational complexity of the simulation is proportional to the geometric resolution. The long run times of the presented 3-D simulation, which was in the range of several hours per run, constitute a technical limitation. The initial simulations were done in 2-D on a single model slice, which gives first merits about the temperature of deep tissue. Although the 2-D simulation times were shorter by several orders of magnitude, there were only minor differences in the steady-state temperature profiles. The final results presented in this chapter have been obtained by the 3-D simulation. Simpler models with lower geometric resolution could not reproduce vessels, which significantly deteriorate the accuracy of the simulation results. Regarding the available measuring equipment and environmental clinical conditions (e.g., variations in room temperature), the simulation accuracy should be approximately 0.1 °C.
We also implemented a parallel version of the proposed method, which runs efficiently on 16 or more connected computers in a computing cluster. In this way, the computation time can be shortened significantly. This approach enables the solution of several millions of equations for each time-step. Besides the simulated results, the computer simulations also support the development of new ideas and theories because a lot of "unexpected" situations can be explained. The proposed method can be applied in investigations of a variety of living tissues.
The solution of the Bio-heat equation was used in three different cases in the field of medicine. In the first case, described in Section 4, our primary goal was to explain, by computer simulation, the variations in the Pennes’ in vivo measurements of the steady-state temperatures along the transverse axis of the proximal forearm. Assuming that the anatomical positioning of the measuring probes varied because of the different forearm sizes of the subjects, we have reconstructed their possible paths by searching for the simulated paths that result in the best agreement between the simulated and the measured temperature fields. Several attempts have been made in the past to average the measured steady-state temperatures in order to smooth fluctuations and obtain a standardized response. Our simulated results suggest that this approach is inappropriate, and that the fluctuations are natural result and the consequence of a complex interplay between the position of the measuring probe, the anatomical position of the main arteries, the dimensions of the forearm, the blood flow, the inhomogeneity of tissues, and the environmental temperature.
We have confirmed that blood flow has a significant and complex impact on the steady-state temperatures of the forearm. Lower blood flow results in lower temperature of the central part of the forearm. When the effect of blood flow on the temperature profile is increased, by disabling the direct effect of the main arteries, the resulting temperature profiles match very well to the analytically produced profiles from the literature, but not well with the experimental measurements. We have shown that only when the effect of both sources, the heat distributed via blood flow and the heat emitted directly from the arteries, are considered, the measured temperature profiles can be fully explained. Furthermore, we have revealed the dependence of the temperature profile shape on the dimensions of the forearm. Smaller forearms are more likely to feature off-centre peaks in the vicinity of the main arteries than larger forearms. In larger forearms, a central plateau may be formed, masking all other potential peaks. The thermodynamic constants of the tissue and the tissue metabolism were confirmed to have a relatively minor impact on the temperature field.
In the second test case, described in Section 5, we demonstrated a practical application of our simulation program in topical postoperative cooling of a knee with two different cryotherapeutic methods. First, we simulated the topical cooling of a knee with a gel-pack. The inner knee tissues reached their lowest temperature in 40 minutes. For continued effective cooling, the gel-pack has to be replaced. The topical cooling with a cryo-cuff was more effective. We simulated situations with relatively small cooling rates in order to be able to analyse the small influences of the blood flow and metabolism. We observed that various thicknesses of a protective bandage and isolating blanket, or variations in their thermal conductivity and cooling temperature, can be used to regulate the cooling intensity.
The results have been validated by experimental measurements of knee temperatures. The simulation results confirm that the model and methodology used are appropriate for the thermal simulation of bio-tissues. We have shown that the blood flow has a significant and complex impact on the stationary state temperatures, and on the gradient of temperature change in the subcutaneous tissue and in the tissues nearer the central part of the knee. Lower blood flow results in linear temperature profiles with larger gradients, for example, in bones. The dimensions of the knee are a very important factor influencing its temperature distribution and gradient. Temperatures in smaller knees will differ more from the arterial blood temperature than those in larger knees. At the same time, the temperature gradient will be much greater in smaller knees. Thermal constants and metabolism have a relatively minor impact on the temperature field.
We confirmed that fluctuations in temperature profiles are natural results and a consequence of the complex interplay between the positions of the measuring probes, the dimensions of the investigated body part, the anatomical positions of the main arteries and bones, and the environmental temperature. The fluctuations and augmented peaks in the temperature profiles can be explained by our simulation results, since the temperature changes in the centre of the knee are influenced by the vicinity of the main leg artery with constant blood temperature. The steepest, almost linear, gradients in the area of the subcutaneous tissues are observed in the outer parts of the simulated temperature profiles, particularly in cases with lower cooling temperatures. The temperature plateaus, measured earlier by other investigators, in the inner knee region with increased dimensions of the knee, were also demonstrated by our simulation.
The simulated results were visualized and compared with measured values. Good agreement was obtained, leading to the conclusion that the model and method used in our simulation are appropriate for such medical simulations. Although there are not many studies of knee temperatures measured in vivo, and the measuring conditions are often not described in sufficient detail, we have also run our simulation software with the initial conditions described in some published measurements. In (Martin, et al., 2001) the temperature in the lateral gutter of the knee decreased by 4, one hour after knee arthroscopy. Similar values have been obtained by our simulation program. More detailed control studies should be done to compare experimentally measured and simulated results in order to fine-tune the simulation method.
The last test case, described in Section 6, stands for the most complex computation in the present chapter. It incorporates all physical model elements, including incompressible fluid flow. We confirmed that the fluid flow plays a crucial role in the process of topical cooling of the heart. We showed that, under simulation assumptions, the cooling process is less effective if natural convection is considered. In the case with no advection, the impact of rigid cold substance lasts longer, resulting in lower simulated temperatures of the heart. However, if the cooling liquid would be externally regulated, the results might be different. The advection increases the heat transfer and if appropriate cooling water management would be chosen, the cooling efficiency might be improved.
We are aware that our simulation methodology still has several limitations. A separate study should be devoted to the validation of simulated results; however, in many cases, in vivo measurements are not applicable or significantly limited. The anatomy of the human body differs between individuals and changes in time; consequently, the simulation on a single generalized model can differ from the experimental measurements. Some kind of personalized approach is foreseen in the future, also because of possible influences of different mechanisms of thermoregulated blood perfusion. When known, these mechanisms could be easily incorporated into the simulation model. Finally, the moving air and the blood in heart chambers were approximated without fluid flow in order to simplify the simulation. Although all these limitations could have some impact on the simulated temperatures, the essential findings are in close agreement with in vivo measurements.
Future work includes the simulation of the temperature field evolution in time and analysis of the possible impact of the measuring probes, surgical instruments and other medical equipment on the temperature fields. Using an approach analogous to the presented one, other parts of the human body could also be simulated and their temperature fields could be analysed. The proposed methodology could be useful for researchers in developing methods for accurate prediction of the temperature distribution in inhomogeneous tissues. Computer simulations could help in the study of various medical applications. It has been proved that the techniques described can be used to predict the temperature distribution in deep tissues, at any body point and time of interest, and for different environmental conditions.
For future development of the numerical approach we will consider the promising local meshless methods, like Radial Basis Function Collocation Method (LRBFCM) (Šarler, 2007) or Diffuse Approximate Method (DAM) (Prax, et al., 1996). The meshless methods are not limited to the equidistant uniform nodal distribution as the FDM is. The meshless approach would improve the accuracy near irregular geometric model boundaries. The DAM shows good stability and accuracy behaviour on non-regular nodal distributions (Trobec, et al., 2012), and consequently, is convenient for treatment of the presented bio-medical computations. The adaptive nodal distribution (Kosec and Šarler, 2011) would also improve the treatment of areas with high gradients, e.g., the interfaces between different materials with discontinues in thermo-physical properties. The meshless methods have been recently under intense research in fluid flow context, as well. The local pressure velocity coupling procedure in combination with LRBFCM has been proposed (Kosec and Šarler, 2008a, Kosec and Šarler, 2008b). The local fluid flow solution procedure is simpler to paralelize as it does not require anykind of global treatment.
The spatial models and the simulation programs, presented in this work, are available from the authors for research purposes that could contribute to any kind of knowledge in the area of biomedical simulation.
|S||heat source||specific heat|
|j||heat flux||thermal conductivity tensor|
|identity matrix||density of the blood|
|specific heat of the blood||spatial/material indicators|
|blood flow rate||metabolism response|
|reference temperature||reference blood flow rate|
|blood flow rate parameters||reference metabolism heat production|
|tissue reference temperature||dynamic viscosity|
|Boussinesq reference density||Boussinesq reference temperature|
|differential operator||position vector|
|spatial step indirection||component of the position vector|
|component of the velocity||arbitrary field|
|current time step||next time step|
|intermediate velocity||velocity correction|
|pressure correction||velocity boundary condition|
|normal boundary vector||temperature of skin|
|ambient temperature||convection coefficient|
We thank the staff of the University Medical Centre Ljubljana, Department of Cardiovascular Surgery and Department of Traumatology, Slovenia. We also acknowledge the financial support from the state budget by the Slovenian Research Agency under the grant P2-0095.