Mathematical models based on transfer functions to estimate tissue temperature during RF cardiac ablation in real time.

Radiofrequency cardiac ablation (RFCA) has been used to treat certain types of cardiac arrhythmias by producing a thermal lesion. Even though a tissue temperature higher than 50ºC is required to destroy the target, thermal mapping is not currently used during RFCA. Our aim was thus to develop mathematical models capable of estimating tissue temperature from tissue characteristics acquired or estimated at the beginning of the procedure (electrical conductivity, thermal conductivity, specific heat and density) and the applied voltage at any time. Biological tissue was considered as a system with an input (applied voltage) and output (tissue temperature), and so the mathematical models were based on transfer functions relating these variables. We used theoretical models based on finite element method to verify the mathematical models. Firstly, we solved finite element models to identify the transfer functions between the temperature at a depth of 4 mm and a constant applied voltage using a 7Fr and 4 mm electrode. The results showed that the relationships can be expressed as first-order transfer functions. Changes in electrical conductivity only affected the static gain of the system, while specific heat variations produced a change in the dynamic system response. In contrast, variations in thermal conductivity modified both the static gain and the dynamic system response. Finally, to assess the performance of the transfer functions obtained, we conducted a new set of computer simulations using a controlled temperature protocol and considering the temperature dependence of the thermal and electrical conductivities, i.e. conditions closer to those found in clinical use. The results showed that the difference between the values estimated from transfer functions and the temperatures obtained from finite element models was less than 4ºC, which suggests that the proposed method could be used to estimate tissue temperature in real time.

INTRODUCTION

Radiofrequency cardiac ablation (RFCA) is currently used to treat some types of cardiac arrhythmia. This technique uses radiofrequency (RF) current (≈500 kHz) to produce a thermal lesion and hence tissue necrosis in the target zone responsible for the arrhythmia. Electrical current is delivered to the tissue through a small active electrode placed at the tip of a percutaneous catheter and a large dispersive electrode located on the patient’s back. There are two modes of delivering radiofrequency energy: constant voltage and temperature-controlled. The constant voltage mode means that the applied voltage is constant; under this condition high temperatures (approximately 100ºC) can be reached in the tissue around the electrode tip in few seconds and lesion size is limited by the charring of tissue around the electrode. The temperature-controlled mode, which can make a larger lesion than the constant voltage approach [1], is the one most frequently used in clinical practice [1, 2] and consists of modulating the applied voltage to keep the temperature constant at the active electrode. The temperature sensor embedded in the electrode tip measures the approximate temperature at the tissue-electrode interface. Even though a tissue temperature higher than 50ºC is required to destroy the target, thermal mapping is not currently used during RFCA. This is an important issue, since tissue temperatures above 100ºC are associated with charring and can cause a thrombus [3]. Safer ablations could thus be achieved by knowing the temperature in the tissue in order to ensure the creation of the thermal lesion and to avoid overheating. Our aim was thus to develop mathematical models to estimate tissue temperature in real time from the tissue characteristics (electrical conductivity (σ), thermal conductivity (k), specific heat (c) and density (ρ)) measured or estimated at the beginning of the procedure and the applied voltage at any time.

Previous work has been conducted on estimating lesion progress on-line by studying impedance progress [4-6]. Numerical modeling (e.g. based on the Finite Element Method (FEM)) has also been proposed to estimate tissue temperature during RFCA [7]. This method calculates the temperature distribution by means of solving numerically partial differential equations problems [8] based on physical laws and hence the models using this method can be considered as structural models. We previously proposed and experimentally verified autoregressive moving average models with exogenous input (ARMAX) aimed at estimating tissue temperature during ablation [9]. These were clearly non-structural models since they were not based on physical laws. Our present proposal is based on mathematical models to estimate tissue temperature at a depth of 4 mm (T4) by measuring the applied voltage (U=V2) and initial value of the tissue characteristics, (see Fig. ). These models are neither directly based on physical laws (e.g. FEM based) nor are they totally non-structural (e.g. ARMAX models), but are based on transfer functions which consider the biological tissue as a system with inputs (applied voltage and initial value of tissue characteristics) and output (tissue temperature). They would be relatively simple to implement while retaining all the important characteristics of more complex models and hence could be used for real time estimation of tissue response or as a theoretical basis to obtain more sophisticated temperature controllers. In this study, we set up the models using theoretical models based on FEM, and then verified the feasibility of the proposal by other FEM models.

METHODS

Description of the Theoretical Model

We considered an active 7Fr and 4 mm long platinum-iridium electrode as used in clinical procedures. This electrode was positioned perpendicular to the tissue, which implied a rotational symmetry axis and allowed a two-dimensional model to be considered. Fig. () shows the geometry and dimensions of the model, which includes the active electrode with a fragment of the polyurethane probe, a thermistor embedded in the active electrode tip, a section of the thermistor coating material and a fragment of cardiac tissue. We considered the electrode to be inserted into the cardiac tissue to a depth of 1.25 mm [10] and the study focused on the temperature evolution at a depth of 4 mm (T) (see Fig. ). From a mathematical point of view, the model was based on a coupled electric-thermal problem which was solved using the Bioheat Equation [8]: where T is temperature, t is time, q is heat source produced by the RF power, Q is heat loss by blood perfusion and Q is metabolic heat generation. The terms Q and Q are insignificant in RFCA and were thus not considered. The heat source q (Joule losses) is given by q=JE where J is the current density and E is the electric field strength. The electrical problem was solved using the Laplace equation ∇(∇V)=0 in which V is voltage. The electric field was calculated by means of E and J using Ohm's law (J=σE), where σ is electrical conductivity. We used a quasi-static approach due to the frequencies used in RF ablation (≈500 kHz) and for the geometric area of interest the tissues can be considered as purely resistive [11]. The blood circulating inside the cardiac chamber was not considered in the model domain. However, the blood-tissue interface and blood-electrode interface were really model boundaries (see Fig. ).

The electrical boundary conditions were: voltage zero density in the transversal direction to the symmetry axis, blood-tissue interface and blood-electrode interface. A boundary condition of voltage zero was set at the dispersive electrode. The voltage at the active electrode tip was 15 V. This value was chosen in order to keep T below 100ºC. The thermal boundary conditions were: null thermal flux in the transversal direction to the symmetry axis and constant temperature at the dispersive electrode and at the outer end of the plastic probe. The circulating blood inside the cardiac chamber produces a cooling effect on the tissue and on the electrode surface, which was modeled by means of two thermal convection coefficients h and h which represent the cooling effect at the blood-tissue interface and at the blood-electrode interface respectively (see Fig. ). We considered a medium blood flow rate, with values for hof 3636 W/m2K and for h of 708 W/m2K [10] and an initial temperature of 37ºC. Table shows the values of the physical characteristics of the materials [10, 12]. Electrical and thermal conductivity variations with temperature were not considered in these FEM models in order to obtain a relatively simple mathematical model which could be used as a first step towards a more complex (non-linear model).

We used COMSOL Multiphysics software (COMSOL Inc., Burlington, MA, USA) to implement the numerical solution based on the FEM. In order to avoid boundary effects, the model dimensions R, Z and L (see Fig. ) were estimated by means of convergence tests in which the value of the maximum temperature reached in the tissue (T) after 120 s of RF heating was used as a control parameter. Parameter values were increased by 1 mm in each simulation. If there was a difference of less than 0.5% between T and the same parameter in the previous simulation, the former values were considered as appropriate. The mesh was initially heterogeneous and a finer mesh was used around the electrode-tissue interface, where the highest gradient was expected. The mesh was then refined to determine adequate spatial resolution and the time-step was reduced to determine adequate temporal resolution.

Identification of the Transfer Function

In the context of our study, the biological tissue was considered to be a dynamic system with applied voltage (V) as input signal and the temperature obtained at a depth of 4 mm (T) as an output signal (see Fig. ). It is assumed that the dynamic evolution of the system output (T) with respect to the system input (U) can be expressed mathematically in terms of a differential equation. If the differential equation is linear, the mathematical input/output relationship can be expressed as a transfer function.

The first step was to obtain an accurate identification of the system response of the modeled tissue, which involves obtaining the transfer function G(s). This was carried out firstly by applying a constant voltage of 15 V for 300 s. The The model is assumed to be U(t)=Vin which V was the applied voltage. We checked that the active electrode tip temperature stabilized around 55°C. We then obtained the evolution of T, and MATLAB was then used to estimate the model response.

A first-order transfer function was obtained from the simulated evolution. The estimated transfer functions G(s) in Laplace-transform quantities had the form: where ΔT is the temperature increase (above a base temperature of 37º), U(s) is the square of the voltage applied to the electrode, K is the steady state gain of the system expressed in ºC/V, and τ is the time constant. A first order transfer function was used since it is the simplest transfer function that leads to results reasonably similar to those obtained from finite element simulations.

From a practical point of view and using the above parameters for the first model obtained, the ablation final temperature could be calculated as follows: where ΔTbeing V the applied voltage and T= 37°C.

Description of the Analyzed Cases

In order to study the influence of the tissue characteristics on the mathematical behavior of the ablation process we carried out several computer simulations with varying values of σ, k and c. We considered a variation from −75% to 100%, in steps of 25%, beginning with the value shown in Table . We show an example: where k=0.531 W/mK and k-25%= 0.398 W/mK is the thermal conductivity variation. A first-order transfer function (2) was identified for each case, with corresponding values of K and τ. This group of transfer functions was obtained in order to study the relationship between the transfer function gain (K) and time constant (τ) with the variations in tissue specific heat and electrical and thermal conductivity. From the series of the K and τ values, a general mathematical relationship for the variation of K and τ for different values of σ, k and c was obtained.

Building the Mathematical Model

We divided the study into two parts: firstly, we assessed the influence of the values of σ, k and c independently, i.e. for each variable separately in order to study the relationship between gain and time constant in each tissue characteristic. From these relationships we were able to obtain a mathematical model.

Secondly, in order to obtain better-fitting parameters in the first-order transfer function of the mathematical model, we considered simultaneous variations around the ends of σ, k and c and then used the linear minimum square fit to obtain more accurate parameters.

Assessment of the Performance of the Models Under Different Conditions

Having built the mathematical model based on a first-order transfer function with fixed applied voltage and temperature non dependent on the electrical and thermal conductivity, in order to assess the accuracy of the models we carried out new computer simulations using a temperature-controlled protocol (variable applied voltage to maintain target temperature in the temperature sensor) and considering the temperature to be dependent on electrical and thermal conductivity, mimicking conditions similar to those found in clinical use.

In the models, we used random values around the center of variations (see Table ). In addition, for a better approximation to clinical reality we considered the simulations to have a controlled temperature protocol and also that temperature was dependent on tissue thermal and electrical conductivity, with the aim of keeping the temperature in the active electrode approximately constant. Therefore, we used a PI controller with a target temperature of 55ºC. The values of the Kp and Ki of the PI controller were 4.78 and 3.39, respectively [13]. The change in cardiac tissue electrical conductivity was taken as +1.5%/ºC [10] and in thermal conductivity as +0.001195K-1 [14]. The temperature evolution at a depth of 4 mm was considered as the real temperature, in order to be comparable to that estimated by the mathematical models. From the transfer functions temperature progress was estimated by means of: where T4(t) is the inverse Laplace transform of G(s)V2(s), and is obtained by numerical integration. The error was calculated as the difference between the estimated temperature (obtained by the mathematical model) and the real temperature (obtained by COMSOL).

RESULTS

Study of Temperature Evolution at a Depth of 4 mm

Fig. () shows the temperature evolution at a depth of 4 mm (T) when electrical conductivity was varied. No changes were observed in the time needed to reach Tduring these variations. The final temperature was seen to be directly proportional to the variation of electrical conductivity. Fig. () shows temperature evolution with variations in thermal conductivity. The value of T was found to be inversely proportional to the variations of thermal conductivity, i.e. high values for T were obtained with low thermal conductivity values and vice versa. Changes were observed in the time needed to reach T, high thermal conductivity values produced faster heating than low values. Finally, Fig. () shows temperature evolution for changes in specific heat. In this case, T was also inversely proportional to the variation of specific heat. Changes were observed in the time needed to reach T, low specific heat values produced faster heating than high values.

Static and Dynamic Characterization

Fig. () shows the values of the gain (K) and time constant (τ) against variations in tissue characteristics. Fig. () shows a linear variation of the gain with respect to electrical conductivity. Gain was practically constant when specific heat was varied, and was a hyperbolic function of thermal conductivity. Fig. () shows the evolution of τ against variations in tissue characteristics. Note that the value of the time constant was not affected by changes in electric conductivity. Conversely, the time constant was linearly dependent on specific heat variations and inversely proportional to thermal conductivity.

From this behavior we tried to obtain transfer functions for each tissue characteristic:

Fit for Variations in Electrical Conductivity

In this case, the relation was constant with time constant and linear with the gain was linear. We therefore proposed the following relationship: where the gain can be represented by Kv=A and the time constant τ=τwhere As=0.07632 ºCm/s obtained by using a linear minimum square fit on the obtained data and t is at 52.7 s.

Fit for Variations in Thermal Conductivity

Gain and the time constant depended on thermal conductivity as a hyperbolic + constant function. Thermal conductivity variations were therefore considered as: where A=0.01900 ºCW/mK, B=0.005059ºC, D=25.2171 and E=5.423 were obtained by using a linear minimum square fit on the obtained data.

Fit for Variations in Specific Heat

The time constant showed a linear dependence on thermal conductivity and gain was almost constant. Specific heat variations were therefore considered as: where the gain can be represented by Kv=0.041ºC and D=0.01478 and E=5.009 are obtained by using a linear minimum square fit.

General Mathematical Model

We know that the time constants in thermal systems can be expressed as τ=c/k where τ is the time constant, c is the specific heat and k is thermal conductivity. By taking this into account, we combined equations (6), (7) and (8) to propose the following mathematical model: where A, B, D and E could be obtained by using a linear minimum square fit on the obtained data. However, to obtain a better fit we considered simultaneous variations of σ, k and c. The resulting values were A=0.0342ºCW/sK, B=0.0114ºCm/s, D=0.0083 and E=3.8908.

Assessment of the performance of the models under different conditions

Fig. () shows the temperature evolutions obtained from FEM and from the mathematical model, both of which correspond to tissue temperature at a depth of 4 mm. Ten different tissue characteristic values were considered (see Table ). The mathematical models underestimated the temperature in all cases. Fig. () shows the error evolution in each case, which was obtained from the difference between the estimated tissue temperature and the real temperature (obtained from FEM-based models). It can be seen that the differences in temperature are less than −4ºC.

DISCUSSION

With the aim of developing mathematical models to estimate in real time tissue temperature from the tissue characteristics measured or estimated at the beginning of the procedure (electrical conductivity, thermal conductivity, specific heat and density) and applied voltage at a given time, we built a mathematical model based on a first-order transfer function. The temperatures estimated from the mathematical models proved to be reasonably similar to those obtained from numerical models, since they always kept lower than 5ºC, which is around the range of lethal isotherm found in experiments using cardiac tissue [15]. The maximum error of 4ºC was obtained from the maximum difference between the tissue temperature estimated by the mathematical model and tissue temperature obtained from COMSOL Multiphysics. In all the cases the temperature was underestimated. The temperature dependence of the thermal and electric conductivity was not initially considered in the model and may be the likely reason for the consistent underestimation of tissue temperature.

In order to obtain the input parameters for the mathematical models the initial values of the tissue characteristics (σ, k and c) must be known and can be obtained from experimental measurements, as described in previous studies. For example, techniques based on automatic swept-frequency network and impedance analysers have been proposed to measure electrical conductivity [16], and the empirical relation between the power delivered to a heated thermistor by a heating pulse and the temperature rise recorded by the sensing thermistor placed at a distance from the sample could be used to measure thermal conductivity [14]. Further studies will be undertaken in order to obtain a mathematical model that takes into account temperature-dependent parameters and new methods of estimating initial tissue characteristics.

Clearly, tissue temperature estimation is currently a challenge. Although previous studies estimated temperature by means of impedance measurements [4-6] and numerical modeling [7], here we present a new method based on simple mathematical models. Our findings suggest that if the applied voltage (U=V2) and the initial value of the tissue characteristics (σ, k and c) are known, mathematical models based on transfer functions could be used to estimate tissue temperature at specific points with reasonable accuracy. This study could lead to the further development of more detailed mathematical models by considering additional factors, such as insertion depth and circulation blood flow.

CONCLUSIONS

The results suggest that it is possible to estimate temperature at a specific point (at a depth of 4 mm in this case) from the evolution of the applied voltage when the initial value of the tissue characteristics is known. The proposed mathematical model agreed reasonably accurately (around 4ºC) with the results given by FEM.

Table 1.

Characteristics of the Materials Used in the Model [10, 12]

Material	Region	σ (S/m)	ρ (kg/m3)	c (J/kgK)	k (W/mK)	Reference
Cardiac Tissue	Myocardial	0.541	1060	3111	0.531	[10]
Pt-Ir	Electrode	4×106	21.5×103	132	71	[10]
Insulation	Coating	10-5	32	835	0.038	[10]
Glass fiber	Thermistor	10-5	32	835	0.038	[12]
Polyurethane	Catheter body	10-5	70	1045	0.026	[10]

Note: σ: electrical conductivity, ρ: density, c: specific heat, k: thermal conductivity. Tissue characteristics evaluated at 37°C.

Table 2.

Cases Considered in Validating the Mathematical Model

Cases	k	σ	c
1	-25	0	25
2	25	0	-25
3	-50	-25	0
4	0	25	50
5	50	0	-50
6	-50	-25	-50
7	50	-25	25
8	50	-50	-25
9	25	25	25
10	-25	50	25