A thermoelastic deformation model of tissue contraction during thermal ablation

Abstract

Purpose: Thermal ablation is an energy-based ablation technique widely used during minimally invasive cancer treatment. Simulations are used to predict the dead tissue post therapy. However, one difficulty with the simulations is accurately predicting the ablation zone in post-procedural images due to the contraction of tissue as a result of exposure to elevated temperatures.

Materials and methods: A mathematical model of the thermoelastic deformation for an elastic isotropic material was coupled with a three state thermal denaturation model to determine the contraction of tissue during thermal ablation. A finite difference method was considered to quantify the tissue contraction for a typical temperature distribution during thermal ablation.

Results: The simulations show that tissue displacement during thermal ablation was not bound to the tissue heated regions only. Both tissue expansion and contraction were observed at the different stages of the heating process. Tissue contraction of up to 42% was obtained with an applicator temperature of 90 °C. A recovery of around 2% was observed with heating removed as a result of unfolded state proteins returning back to its native state. Poisson’s ratio and the applicator temperature have both been shown to affect the tissue displacement significantly. The maximum tissue contraction was found to increase with both increasing Poisson’s ratio and temperature.

Conclusions: The model presented here will allow predictions of thermal ablation to be corrected for tissue contraction, which is an important effect, during comparison with post-procedural images, thus improving the accuracy of mathematical simulations for treatment planning.

Keywords:

• Thermal ablation
• thermal denaturation
• thermoelastic deformation model
• tissue contraction
• mathematical model

Introduction

Thermal ablation is an alternative to surgical resection increasingly used to treat cancer patients in a minimally invasive manner due to its advantages such as reduced morbidity and its ability to treat poor surgical candidates. Typical therapies such as radiofrequency ablation (RFA), microwave ablation (MWA) and laser ablation expose the cancer tissue to elevated temperatures causing death. For a safe clinical implementation of such techniques, it is important to accurately determine the volume of ablated tissue. Pre-procedural planning, such as experimental data and computer simulation tools, can provide such information and thus aid in determining the optimal patient specific treatment [1–3].

Studies have attempted to develop a computational framework to predict the performance of different ablation devices and procedures for clinical implementation [1–3], characterised from ex vivo and in vivo animal studies [4,5]. These models, however, have not considered the effects of tissue contraction during hyperthermic ablation therapies, which can be observed on immediate post-procedural images both in clinical and experimental settings [6,7]. As a result, they can underestimate the pre-treatment dimensions of tissue that has been killed. Assessing the degree of contraction associated with the different hyperthermic ablation techniques would therefore improve the accuracy of determining the ablated tissue.

Experimental studies on animals have reported tissue contraction as a result of MWA and RFA. Ex vivo studies have shown tissue volume contractions of around 40–70% during liver ablation [6–9], 50–60% during lung ablation [6] and 26–42% during kidney ablation [10]. Clinical studies have reported both contractions and expansions during ablation therapies. Treating hepatic haemangiomas and renal cell carcinomas using MWA reported tissue volume contractions of 62% (25% in diameter) [11] and 52% [12], respectively. Two independent studies on renal tumours reported tissue diameter contractions of –10 to 50% [13,14]. MWA and RFA on liver tumours have also reported both tissue contraction and expansion, the latter being observed more frequently during RFA [15].

Biological tissue subjected to a non-uniform temperature distribution experiences internal forces, leading to thermoelastic deformation. This deformation can be determined from the thermoelastic wave equation. Studies have attempted to determine the effects of tissue deformation as a result of ablation therapies on tissue samples by solving for the thermoelastic wave equation with simple pre-described ablation temperatures [16,17]. These models, however, only consider the mechanical deformation induced by thermal expansion.

In addition to this effect, several ex vivo experiments have associated contraction of tissue exposed to elevated temperatures with effects such as protein denaturation and dehydration [6,18,19]. Studies have developed mathematical models based on physiology to determine the behaviour of tissue contraction as a result of protein denaturation during hyperthermic ablation [20,21]. This was achieved using state models with first order reaction kinetics describing the transition of molecules between the different states. The change between the different states results in a change in length of the constituent molecules, which presents itself as a gross tissue contraction at the macroscale.

The purpose of this study was to estimate the degree of tissue contraction as a result of heating during typical ablation procedures, and thus to improve the accuracy of mathematical simulations for treatment planning. A previously developed three state model of protein denaturation [20] was coupled with the thermoelastic wave equation, which determines the deformation due to the internal forces. A typical temperature distribution during hyperthermic ablation therapies was considered to solve for the coupled equations using a finite difference method. The effects of the different parameters involved were finally quantified.

Methods

The thermoelastic deformation caused by the elevated temperature during thermal ablation procedures can be derived from the stress–strain equations for an elastic material with homogenous, isotropic property with thermal load. The effects of protein denaturation were introduced to the equations here by assuming that it had an additive effect to the strain with a linear relationship with the stress similar to the thermal expansion term. The general form in Cartesian co-ordinates can thus be expressed in standard tensor notation as: (1) where σ is the stress, E is Young’s modulus, ν is Poisson’s ratio, ɛ is the strain, δ is the Kronecker delta function, α is the coefficient of the volumetric thermal expansion, T is the temperature, β is the coefficient of the volumetric protein denaturation shrinkage and ξ is the relative shrinkage due to protein denaturation. The thermoelastic wave equation can then be expressed in general form as: (2) where ρ is the density, u is the displacement vector, t is time, G is the shear modulus, ∇ represents the nabla operator, K is the bulk modulus and f is a force density vector. The shear modulus and bulk modulus can be obtained from: (3)

A three state model was considered here for the thermal denaturation process with the thermal reaction rate being dependent on the Arrhenius equation. The equations and parameters to determine the strain can be found in Appendix 1, with the theory and the values of the parameters involved given in full in Park et al. [20].

To compare the ex vivo thermal ablation results, Equation (2)(2) was expressed in a spherical axi-symmetric co-ordinate system. Non-dimensionalising using the following characteristic values: (4) where R is the radial co-ordinate, and ω, τ, U and θ are the characteristic values for the length, time, displacement and temperature, respectively, and dropping the primes, Equation (2)(2) can be expressed as: (5)

Three non-dimensional quantities can be observed, the first of these being the mechanical stress term, the second being the thermal expansion term and the third being the protein denaturation term.

The tissue was assumed to have a uniform temperature and to be stress free before thermal ablation with the initial values of the temporal derivatives of temperature and tissue displacement taken to be zero. The spatial derivative of the tissue displacement at the origin is taken to be zero by symmetry. As r → ∞, the tissue was assumed to have zero temperature increase and stress. Taking the domain of the non-dimensional variables to be t ∈ [0,t_max] and r∈[0,a], where a is the non-dimensional radius of the sphere considered, the boundary and initial conditions can be expressed as u(0,t)=0, ∂u(0,t)/∂r=0, u(a,t)=0 and u(r,0)=0, ∂u(r,0)/∂t=0, respectively.

Tissue heating is considered to take place through an applicator, which is assumed to be inserted at r=0. A Gaussian function along the radial co-ordinate is considered for the temperature distribution with an exponential recovery in time so that: (6) where T₀ is the applicator temperature, T_b is the baseline temperature, H is the unit step function, t_m is the non-dimensional time shift of the unit step function and t_c is the non-dimensional time constant for the temperature change. Figure 1 shows the plot of the above temperature distribution. The temperature is at its maximum on the surface of the thermal applicator and decays with increasing distance. Tissue heating outside two non-dimensional characteristic lengths from the applicator (2 cm) is within one non-dimensional temperature unit (equivalent to 1 K) from baseline temperature.

Figure 1. Temperature distribution for a typical thermal ablation therapy.

The baseline value of β can be obtained from the isothermal experimental data of Chen et al. [22], where cylindrical tissue samples were heated under a zero stress condition and the tissue contraction measured in the axial co-ordinate system. Equation (1)(1) can thus be expressed after some re-arranging as: (7) where ɛ_zz represents the normal strain in the axial co-ordinate system. The axial tissue displacement at the macroscale can be expressed in dimensional form as: (8) where L is the length of the tissue sample. Under isothermal conditions, the relative shrinkage term is not dependent on the spatial co-ordinates and thus: (9)

The experimental data suggest that the effects of thermal expansion on tissue displacement are much smaller than the protein denaturation effects, hence β can be approximated to equal three when applying the protein denaturation model [20].

A finite difference method was considered here to solve for the governing equations numerically. The resulting matrix form of Equation (5)(5) is shown in Appendix 2. The Courant–Friedrichs–Lewy condition for all simulations was set to be less than one by varying the time step. This is a requirement for convergence when solving hyperbolic partial differential equations numerically by the method of finite differences, such as the thermoelastic wave equation considered here.

Results

Simulations were run for various conditions. The applicator temperature, Young’s modulus and the Poisson’s ratio were varied one at a time to determine their individual effect on tissue contraction. Table 1 shows the values of the thermomechanical parameters of the tissue and the applicator temperature considered here. Both the baseline values and the range of values considered during the simulations are shown. The baseline applicator temperature, Young's modulus of the tissue and Poisson’s ratio were taken to be 90 °C, 1 MPa [23] and 0.45 [24], respectively. Table 2 shows the values of the additional parameters involved in the simulations. All values stated here onwards, including those in figures, are in non-dimensional form, unless followed by a unit.

Table 1. List of the tissue thermomechanical parameters and applicator temperature.

Parameter	Value	Unit	Range	Reference
ρ	1050	Kg m^−3		[32]
ν	0.45		0.25–0.45	[24]
E	1	MPa	O (10^–6) – O (1)	[23]
α	3 × 10^−4	K^−1		[17]
β	3
Tb	310
T0	363		313–368

ρ is the density, ν is Poisson's ratio, E is Young's modulus, α is the coefficient of the volumetric thermal expansion, β is the coefficient of the volumetric thermal denaturation, T_b is the baseline temperature and T₀ is the applicator temperature.

Table 2. List of simulation parameters.

Parameter	Value	Unit
ω	0.01	m
U	0.01	m
a	10
Δr	0.05
τ	1	s
tc	60
tm	600
tmax	2000
Δt	O (10^−3) – O (10^−5)
θ	1	K

ω, τ, U and θ are the characteristic values for the length, time, displacement and temperature, respectively, a is the non-dimensional radius of the sphere in the spatial domain, Δr and Δt are the grid spacing and time step in the discretisation, respectively, t_c and t_m  are non-dimensional time constants that determine the temperature profile and t_max is the maximum time in the temporal domain.

Figure 2 shows the displacement of tissue as a result of heating when the effects of protein denaturation were ignored. This was achieved by taking β to equal zero. There was a fast increase in the displacement close to the location of the applicator where the temperature was at its highest followed by a decay with increasing distance. However, the magnitude of the tissue expansion was of the order of 10⁻³ non-dimensional characteristic lengths (10⁻⁵ cm), thus not significant and justifying the approach taken in the methods section to obtain an estimate value for β.

Figure 2. Tissue displacement during a typical thermal ablation temperature distribution ignoring protein denaturation effects.

Figure 3 shows the variation of protein proportion in the (a) native and (b) unfolded states. A rapid drop in the proportion of proteins in the native state was observed in the region where the distance from the applicator was less than 1 cm, leading to a sharp increase in the proportion of proteins in the unfolded state in this region as observed. The unfolded protein changes to the denatured state with time, thus causing a drop in the value towards zero. There was no observable change in the proportion of protein in the native state where the distance from the applicator was greater than 1.5 cm as heating in this area did not lead to significant changes in the reaction rates of protein from the native to the denatured state. There is no apparent change in the proportion of protein states outside the heated region as expected.

Figure 3. Plot of protein proportion in its (A) native and (B) unfolded state during a typical thermal ablation temperature distribution.

The effects of protein denaturation on tissue contraction were examined. Figure 4 shows the displacement of tissue as a result of protein denaturation. During heating, an initial expansion in tissue is observed up to just under 50 s followed by a tissue contraction with increasing time. Contractions reach a maximum at around 500 s followed by a slight recovery in contraction and approaching a steady value post heating. This can be observed more clearly in Figure 5 where the plots of the tissue displacement along the distance from the applicator at several time points are shown. A maximum radial contraction of around 42% can be observed in the simulations at a distance of 0.85 cm from the centre of the applicator followed by a decay with increasing distance.

Figure 4. Tissue displacement during a typical thermal ablation temperature distribution.

Figure 5. Tissue displacement at various time points during a typical thermal ablation temperature distribution.

The effects of tissue properties and applicator temperature were next examined. Changes in Young’s modulus within the physiological range did not show visible changes in the magnitude of the contractions. Changes in Poisson’s ratio, however, had a significant impact on the tissue displacement. Figure 6 shows the tissue displacement for a Poisson’s ratio of (a) 0.35 and (b) 0.25. The plots of the tissue displacement for the three different Poisson’s ratio values considered differed only in the magnitude of the peak tissue contraction, decreasing with decreasing Poisson’s ratio. The peak tissue contraction was observed at the same distance from the applicator for the different Poisson’s ratio values.

Figure 6. Tissue displacement during a typical thermal ablation temperature distribution with a Poisson's ratio of (A) 0.35 and (B) 0.25.

A change in applicator temperature also affected the tissue displacement. Figure 7 shows the tissue displacement for an applicator temperature of (a) 80 °C and (b) 70 °C. A drop in temperature led to a drop in the peak tissue contraction. Furthermore, the distance of the peak values from the applicator decreased with decreasing applicator temperatures. Figure 8 shows the value of the maximum tissue contraction for varying applicator temperature. The maximum tissue contraction varies in a nonlinear manner with temperature. It has a sigmoidal curve shape with a low absolute slope value initially, increasing to an inflection point, followed by a levelling off. The effects of tissue contraction become large at applicator temperatures above 333 K (60 °C).

Figure 7. Tissue displacement during a typical thermal ablation temperature distribution with an applicator temperatures of (A) 80 °C and (B) 70 °C.

Figure 8. Maximum tissue contraction for varying applicator temperature during a typical thermal ablation temperature distribution.

Discussion

The degree of tissue contraction was found to be dependent on the applicator temperature and Poisson’s ratio, each having a positive correlation with tissue contraction; an increase in any of these values increased the degree of tissue contraction. The tissue contraction results obtained from the simulations were found to be within the range observed in ex vivo experiments, where ablation zone contractions were up to 15–50% in diameter and 25–75% in volume. The maximum contraction value obtained from the simulation here is around three times greater than values reported in ex vivo bovine liver experiments [8,25], in accordance with experimental observations [22,26]. Results suggested that changes in Young’s modulus within the physiological range did not have a significant impact on the tissue contraction. Recovery of tissue contraction was also observed post-heating as noted in experiments [22].

Regions of expansion and contraction were observed along the tissue; an initial expansion was observed as a result of the thermal expansion effect followed by a contraction due to the effects of protein denaturation. This was due to the response of the thermal expansion being immediate to the initial change in temperature whilst changes in protein states required temperatures above a certain threshold for them to become significant. Since the mechanical stress and the thermal denaturation terms are the dominant effects for the given parameter baseline values the effects of thermal expansion are only appreciable during the initial stages of the heating process.

Tissue heating leads to an increase in the forward and backward reaction rates between the native and unfolded states. The reaction rates increase at different amounts with increasing temperature. An overall forward reaction rate occurs above a threshold temperature leading to protein changing from the native to the unfolded state. Since protein length decreases along the denaturation process, this leads to an initial tissue contraction. This is followed by proteins changing from the unfolded to the denatured state causing further contraction. Returning the temperature back to baseline causes an overall backward reaction rate between the native and unfolded state, which is significantly greater than the forward reaction rate from the unfolded to the denatured state. As a result, any unfolded protein left will return back to its native state, causing a recovery in tissue contraction. The simulation results, however, found this to be small (around 2% depending on applicator temperature) as most of the protein had already progressed to the denatured state. Tissue contraction was also observed in areas outside the heated region, thus further highlighting the importance of tissue displacements during hyperthermal ablation therapies. Note that the largest contractions were not observed where the tissue temperature was at its highest.

The focus of this study was to couple the body forces and protein denaturation to determine biological tissue contraction as a result of heating and to determine the degree of tissue contraction during typical thermal ablation therapies. A simplified model of each of the involved properties was considered as a result. The temperature distribution was based on a Gaussian function and the physical properties of tissue were taken to be elastic and isotropic. The effects of protein denaturation were introduced to the stress–strain equations here by assuming that it had an additive effect to the strain with a linear relationship with the stress similar to the thermal expansion term. Furthermore, the protein thermal denaturation process was simplified to a three state model with the reaction rates being dependent on the Arrhenius equation.

A Gaussian function was determined to be a good representation of a typical temperature distribution during hyperthermic ablation therapies, having been considered in previous studies [16,27], as it matches the results obtained from simulations where thermal energy was directly applied as heating source points [17]. A decrease in temperature with increasing distance from the heating source was observed. Furthermore, there was a rapid increase in temperature during the initial stages of the heating procedure. The rate of change in temperature reduced with increasing duration of applied heating. It also met the requirement of the boundary conditions set for this study.

Accurate quantitative representation of tissue contractions during thermal ablation therapies requires solving for the appropriate governing equations, which is dependent on the ablation technique considered [5]. Common methods such as RFA and MWA use different methods to heat the target tissue, thus requiring knowledge of different tissue properties. Further, the tissue exposure temperatures differ between the two techniques; RFA usually being below 100 °C whilst MWA can reach to temperatures of up to 180 °C. As a result, variations in tissue thermophysical properties are also expected at the higher temperatures experienced by the tissue during MWA, driving greater motion of water [18]. Water content has been found to play an important role in cell swelling effects post-treatment [28] and thus careful consideration will be required for clinical translation.

An assumption made here was that tissue was homogeneous and had an elastic behaviour, thus having equal contractions along the different orientations. Tissue consists of a complex matrix of elastin, collagen and ground substances, each having different physical properties. Given these observations, it is likely that the tissue will have a non-isotropic layout and will undergo a complex non-elastic deformation due to the induced stress. In fact, an experimental study showed different contraction degrees across the different longitudinal and transverse directions for each tissue sample during ex vivo MWA [7]. Furthermore, different organs have different tissue compositions, and hence different mechanical properties. As a result, the degree of tissue contraction from heating has been found to vary across different organs.

Another important parameter based on the physical property of tissue is the protein thermal denaturation procedure and the reaction rate parameters between the different states. Simplifications have been used as the exact protein thermal denaturation process is not known. The reaction rate parameters were obtained from bovine pericardium tissue samples in this study, which have shown tissue contraction rates around three times greater than bovine liver tissue samples. Thus, the tissue contraction reported here is around three times greater than results obtained from liver tissue. These parameters will vary across organs as the degree of tissue contraction is known to vary significantly. Thus for clinical translation, these parameters will have to be obtained for the specific organ being treated. Furthermore, β was obtained here by assuming an isotropic material property such that the tissue displacement as a result of protein denaturation was equal in all co-ordinate directions. This is an assumption that needs further validation as the value of β was found to have a significant effect on the tissue displacement (results not shown); the degree of tissue contraction increased with increasing β.

The parameter values for the physical property of the tissue considered here were obtained from experimental studies on healthy tissue rather than tumour tissue. A recent study on tumour contraction immediately after RFA, however, found the results to approximate the response of normal tissue [13]. This observation combined with the simulation results suggest that Poisson’s ratio and the protein denaturation reaction rates do not differ significantly between healthy and tumour tissue.

The temperatures considered in this study were below 100 °C, thus being more suitable for RFA procedures. To expand this model to procedures such as MWA, where tissue is exposed to temperatures of up to 180 °C, the incorporation of effects such as dehydration associated with the contraction of tissue at elevated temperatures will be necessary. Clinical translation will also require implementing a suitable blood-perfusion model as heat dissipation is known to occur close to blood vessels [29,30]. Studies have shown that the degree of cooling of the tissue varies strongly with vessel size and thus with generations of the vasculature [31]. Hence, knowledge of vasculature structure will also be required.

As well as applicator temperature, studies showed that heating time was also an important factor that affected the degree of tissue contraction [8,9], with a recent study using thermal dose, obtained by taking the integral of the temperature distribution with time, to find the relationship between heating time and degree of tissue contraction [25]. The mathematical model presented here is able to investigate such features as temporal and spatial variations in temperature can be captured by the thermal expansion and protein denaturation terms in the governing equations. This is an aspect that will need to be considered in a future study.

The assumption that the tissue is homogeneous and isotropic is a limitation to this model. This approach was taken due to the limited availability of experimental data as well as to simplify the coupling of the protein denaturation into the stress–strain equation. Loosening of this assumption would require representing the governing equations in a 3-D co-ordinate system form with the respective mechanical and denaturation parameters for each orientation. This is also a possible future study that will need to be considered once experimental data becomes available.

Future work will require coupling the thermoelastic equation, including the effects of protein denaturation, with the bioheat equation and consider accurate parameter values for physical processes during ablation. The latter has been shown to be important in predicting the ablation zone size during RFA [32]. Note that the temperature profile obtained from the heat equation will cause tissue contraction, which will in turn affect the temperature profile, thus requiring the use of moving meshes. This will allow the model to obtain relationships between different ablation therapies and the degree of tissue contraction, thus providing useful information during intervention planning. It will be especially useful to investigate tissue expansion observed in post clinical intervention of tumour organs [15].

Conclusions

The mathematical model presented here, which coupled the effects of protein denaturation with the thermoelastic wave equation, was able to quantify the effects of tissue contraction for typical temperature distributions observed during hyperthermal ablation therapies. Both tissue expansion and contraction were observed at the different stages of the heating process. Tissue contraction observed outside the heated region has emphasised once more the need to incorporate such models during predictions of thermal ablation and comparison with post-procedural images. Varying applicator temperature, β and Poisson’s ratio showed significant changes in the degree of tissue contraction.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This research has been funded by the European Commission under Grant Agreement No. 600641, FP7 Project Go-Smart.

Appendices

Appendix 1. Protein denaturation

The protein denaturation is modelled using three states: native, N, unfolded, U, and denatured, D, thus grouping all the different intermediate states into one overall state. The reaction equation is thus: (A.1) where k_i is the reaction rate. Assuming these to be governed by the Arrhenius equation: (A.2) where E_A is the activation energy, R is the universal gas constant, T is the absolute temperature measured in Kelvin (K) and a is the frequency factor. The equations describing this system are thus: (A.3)

Taking N, U and D to represent the proportion of the different states such that the sum of all the three states is equal to one and assuming that all the protein is in the native state initially, the relative shrinkage can be expressed as: (A.4) where L is the length of the proteins with the subscripts representing the different states.

Appendix 2. Finite difference method

Using a central difference finite difference method, the temporal and spatial derivatives of the tissue displacement can be expressed as: (B.1) and as: (B.2) respectively. Substituting Equations (B.1(B.1) ) and (B.2(B.2) ) into Equation (5)(5) , the tissue displacement can be expressed as: (B.3) for i=0, 1, …, M and j=0, 1, …, N, where:

Taking the boundary conditions u(0,t)=0, ∂u(0,t)/∂r=0 and u(a,t)=0, and the initial conditions u(r,0)=0 and ∂u(r,0)/∂t=0, Equation (B.3(B.3) ) can be expressed in matrix form as: (B.4) where u_j=(u_m)^j, T_j=(T_m)^j, ξ_j=(ξ_m)^j and: for m,n=1, 2, …, M-1.

For the protein state, a forward difference approach is considered so that: (B.5) where X is the state of the protein.