Minimally Invasive Live Tissue High-fidelity Thermophysical Modeling using Real-time Thermography

We present a novel thermodynamic parameter estimation framework for energy-based surgery on live tissue, with direct applications to tissue characterization during electrosurgery. This framework addresses the problem of estimating tissue-specific thermodynamics in real-time, which would enable accurate prediction of thermal damage impact to the tissue and damage-conscious planning of electrosurgical procedures. Our approach provides basic thermodynamic information such as thermal diffusivity, and also allows for obtaining the thermal relaxation time and a model of the heat source, yielding in real-time a controlled hyperbolic thermodynamics model. The latter accounts for the finite thermal propagation time necessary for modeling of the electrosurgical action, in which the probe motion speed often surpasses the speed of thermal propagation in the tissue operated on. Our approach relies solely on thermographer feedback and a knowledge of the power level and position of the electrosurgical pencil, imposing only very minor adjustments to normal electrosurgery to obtain a high-fidelity model of the tissue-probe interaction. Our method is minimally invasive and can be performed in situ. We apply our method first to simulated data based on porcine muscle tissue to verify its accuracy and then to in vivo liver tissue, and compare the results with those from the literature. This comparison shows that parameterizing the Maxwell--Cattaneo model through the framework proposed yields a noticeably higher fidelity real-time adaptable representation of the thermodynamic tissue response to the electrosurgical impact than currently available. A discussion on the differences between the live and the dead tissue thermodynamics is also provided.


I. INTRODUCTION
I N this work, we evaluate the suitability of the hyperbolic Maxwell-Cattaneo model of heat propagation to support a high-fidelity thermodynamic representation of the live tissue response to electrosurgical impact, with applications to porcine liver tissue. As demonstrated in [1], [2], the movement speeds encountered in electrosurgical procedures (5-20 mm/s), when compared to the speed of heat propagation in tissue (around 2 mm/s), shows that wave-like phenomena due to the finite propagation of thermal information play a dominant role in the overall thermodynamics. An important effect of the slow thermal propagation speed of tissue is the fact that the electrosurgical heat source undergoes distortion from an ideally Gaussian shape to a tear-drop-like shape [2], causing significant discrepancies between simulations that use the Maxwell-Cattaneo formulation versus the classical Fourier heat equation in the case of a moving heat source [3]. Since accurately predicting the thermal spread in live tissues is of prime importance in enabling damage-conscious control of the electrosurgical process and preventing the high heat flux from violating safety margins [4], [5], we focus in this work on developing a method for estimating the relevant thermophysical parameters for this model based on thermographer readings.
Previous approaches to parameter estimation in service of modeling the electrosurgical process have been undertaken by [6], [7], focusing exclusively on estimating the thermal diffusivity, as well as [1], where the thermal relaxation time was the subject of study. These methods rely on direct solution of the underlying PDEs (partial differential equations), to find a thermal diffusivity that minimizes the error between the simulation result at specific locations that correspond to locations of embedded thermal sensors, and the temperatures recorded at the sensor locations. This method is invasive,  Fig. 1. Illustration of the proposed parameter estimation for controlled tissue thermodynamics (in this case, electrosurgery is used as an example). In step 1, a probing sequence is determined, consisting of (i) a fixed power stationary probing task, and (ii) a moving cutting task, possibly with variable velocity and power. These probing sequences are applied to the tissue in step 2, during which the tissue thermal response is recorded by a thermographer. These thermographer images are then processed in real-time in step 3, where attention-based noise-robust averaging (ANRA) is used to obtain: (A) the thermal diffusivity, (B) the heat source model, and (C) the thermal relaxation time. Subtasks (A) and (B) rely on the stationary probe data, while subtask (C) is based on the moving probe data. This yields an estimated Maxwell-Cattaneo hyperbolic heat transfer model. requiring embedding sensors into the tissue. In addition, obtaining parameter estimates requires a significant amount of time because the governing PDE must be explicitly solved multiple times to obtain an optimal parameter estimate, precluding the use of such an algorithm for real-time adaptation and safety region mapping during clinical procedures on live tissue. More importantly, multivariate optimization often fails to yield accurate parameter combinations owing to the existence of local minima encountered by optimization algorithms [8]. These local minima often necessitate the use of prior knowledge about the parameters to constrain the search space [6], [7], thereby inadvertently forcing preconceived biases on the obtained parameter estimates. Our proposed method does not require such assumptions, mainly by a proposed probing sequence that isolates various governing phenomena of the underlying heat transfer problem, allowing parameters to be obtained in series, rather than through the commonly observed error-prone method of simultaneous parameter fitting.

Model Estimation
In contrast to the above methods, in [8] we presented a method that relies fully on non-invasive remote infrared thermography, in which a thermal camera is used to capture the surface temperature of the tissue. Our approach does not rely on explicitly solving the governing PDEs, but instead employs noise-robust gradient operators to directly obtain the PDE's governing parameters by balancing the underlying model differential equations based on the thermographic data. In this work, we build upon these results to achieve a method capable of automatic model estimation for a hyperbolic thermal equation, particularly the Maxwell-Cattaneo heat equation. Central to this endeavor is the adoption of a simple minimally invasive probing sequence, which eliminates ambiguity in the parameter estimation step by isolating the various driving phenomena, namely: diffusion, wave propagation, and external forcing. Our approach allows for in situ parameter estimation with a minimally invasive probing sequence (obviating biopsies and laboratory testing), with results being obtained in real-time. These thermophysical parameters are estimated such that our parsimonious thermodynamics model, the Maxwell-Cattaneo heat equation, best fits the observed thermodynamics. Any unmodeled dynamics are subsumed in this limited set of parameters, producing effective parameter estimates that can directly be used for prediction and control of the tissue thermodynamics. To the best of the authors' knowledge, this combination of probing sequence design and direct parameter estimation based on thermographer readings has hitherto not been attempted in the literature. Our method is outlined in Fig. 1.
This paper makes the following contributions to the existing literature: (i) development of a novel procedure for parameterizing the tissue model in real time under electrosurgical impact, based on thermographic data; (ii) using the latter, enabling a higher-fidelity than currently available modeling of the thermal tissue behavior though parametrizing the Maxwell-Cattaneo tissue model, directly applicable during regular operation without invasive sensor insertion; (iii) based on (ii), offering a high fidelity comparison of tissue thermodynamics under in vivo vs ex vivo electrosurgical procedures. This paper is organized as follows. In Sec. II, the Maxwell-Cattaneo heat transfer model is developed and the parameter estimation problem that will be considered is intro-duced. Sec. III presents the parameter estimation method, providing justification for our design choices and illustrating considerations in composing probing sequences to separate thermodynamic phenomena. In Sec. IV, simulation results demonstrating the efficacy of our method compared to a known ground truth are used for verification. In Sec. V, we apply our method to real-life thermographic data from electrosurgery on in vivo porcine liver tissue, discussing differences between the two, and showing improvements in the simulation results based on the updated parameters. Conclusions and future research directions are presented in Sec. VI.

II. PROBLEM FORMULATION
As mentioned in the introduction, in this work we seek to obtain a heat transfer model that captures the salient dynamics observed in energy-based surgery methods applied to live tissue. In our model, we wish to include the following phenomena: i) Diffusion (heat conduction); ii) Presence of a controlled external heat source; iii) Finite heat propagation time; The first of these phenomena, diffusion, is by far the most studied in the field of heat transfer in live tissue [6], [7], [9], [10] (the former two works studying porcine aorta, and the latter considering general animal tissues). The most commonly used formulation is the classical Fourier heat equation, which models thermal conductivity. When operating on a one-dimensional domain, the following expression is obtained based on the thermal diffusivity a: where x denotes the temperature, t the time, and η the spatial location.
After adding the heat source, we obtain the following formulation:ẋ In this expression, the heat source is position at location p(t), where T p(t) is an operator that serves to translate the center of the heat source to p(t). The applied heat source is modeled as a unit heat source field q scaled by the input power u(t).
A major drawback in the classical Fourier heat equation, is the fact that changes in the heat source propagate instantaneously throughout the entire domain [2], [3]. This assumption is unphysical, as heat waves propagate at a finite speed through any physical medium [1]- [3]. An important consequence of an infinite heat propagation speed is a typical assumption that the tissue immediately reacts to external thermal excitation, whereas in practice, when an external heat source is applied there is a finite time before the temperature starts rising [1]. In addition, in the case of electrosurgery, one does not observe a rise of temperature in front of the moving needle, even during regular surgical procedures [1]; this also points at the importance of a finite heat propagation time. Introduction of the latter is often achieved by introducing a second time derivative term, so as to induce wave-like behavior. The formulation thus obtained is commonly referred to as the Maxwell-Cattaneo equation, or the telegraph equation: Here, the second time derivative is scaled by τ , which is known as the thermal relaxation time. A so called second speed of sound is often obtained by considering both the thermal diffusivity and thermal relaxation time as follows [2]: We wish to obtain parameter estimates on a patient-bypatient basis, irrespective of underlying conditions. In light of this goal, we seek a method that is: 1) minimally invasive (i.e., no biopsies or unnecessary tissue resection is required); 2) applicable in situ, without the need for external specialized equipment or laboratory testing; 3) runnable in real-time, with results immediately available. This yields the following general problem statement, which can be classified as a parameter estimation problem: taken on subsets ofΩ with fixed time interval ∆t > 0, as well as control inputs (u(t), v(t)), estimate the following parameters of system (5): 1) Thermal diffusivity a; 2) Unit heat source q; 3) Thermal relaxation time τ . We now proceed by briefly introducing a well-posed mathematical formulation of our model on which our results in this work will be based.

A. Mathematical Heat Transfer Model
In this work, we consider the following thermodynamic model on a two-dimensional compact domain Ω ⊆ R 2 : where a > 0 denotes the thermal diffusivity, τ > 0 denotes the thermal relaxation time, and ∇ 2 is the Laplacian operator.
x(t) denotes the temperature field at time t, such that The heat source position is denoted by p(t) ∈ Ω, which is used as a parameter in the translation operator T , defined as where χ Ω is the indicator function on the set Ω. In our model, q ∈ H 1 (Ω) denotes a normalized heat source term, which is scaled by control input u(t) ∈ R + , which indicates the electrical power provided by the electrosurgical power supply.
We denote the full state by z ∈ H 1 (Ω) × L 2 (Ω) × C 1 (Ω), which includes the temperature, the temperature rate of change, and the probe position.
In this work, we assume having access to thermographer data comprised of temperatures that are sampled on a finite Cartesian grid with uniform spacing ∆η > 0, denoted byΩ ⊂ Ω. We consider dimΩ = M < ∞, such thatΩ : In the remainder of this work, we make the following assumptions: As can be seen in (5), we assume that there are no errors in the probe position, nor the heat source model, i.e., (5) is considered to be a nominal model of the system's thermodynamics, which can be used in a subsequent controller design or for predicting/assessing tissue damage. In the remainder of this work, we consider control input pair (u(t), v(t)) to be given for all t ∈ [0, ∞); we also assume that there is no error in the control inputs.
Existence and Uniqueness: To prove the existence and uniqueness of (5), we may apply the Picard iteration method on the system to obtain a unique solution (x(t, η)) for known input pair (u(t), v(t)) and with any initial condition (x 0 (η)) that satisfies assumption (A1) and (A2) [11, §2.3]; continuity of the ordinary differential equation (ODE) that defines p(t) follows directly from (A4). Thus, well-posedness of the heat equation under consideration is guaranteed.

B. Probing Sequence for Parameter Estimation
Rather than attempting to estimate all three parameters of system (5) at once, it bears mentioning that, in general, there is no unique solution to Problem 1 for an arbitrary electrosurgical cutting sequence. Indeed, the ratios a/τ and q/τ do not necessarily admit unique solutions without imposing further constraints, especially when considering standard methods that are prone to finding only locally optimal parameters, not to mention the case of unmodeled dynamics and noise when dealing with real-life data. To address this problem, a central idea is to precede the main cutting task by a series of minimally invasive probing tasks, to limit the number of phenomena simultaneously being excited to obtain a unique set of parameters.
Let us rewrite the dynamics of x(t) from (5): From observations in [2], it follows that a non-moving source with steady heat input results in τẍ(t) vanishing after a transient, whose duration is dictated by thermal relaxation time τ . Therefore, it is possible to study the following equation in isolation, i.e., the tissue's free response: To do so, one would need to sufficiently raise the temperature, e.g., by sustained application of a small amount of electrosurgical action without moving the electrosurgical probe's contact point. Starting from a low or steady-state temperature will lead toẋ(t) ≈ 0, which does not allow for observing a∇ 2 x(t). Having raised the temperature, after disengaging the heat source, one will need to wait for the thermal wave caused by the sudden lack of external forcing to have subsided as dictated by the second speed of sound.
It suffices to wait for the final heat wave to propagate outside of the field of view of the thermographer, after which the free response of the tissue can be observed. The latter has been the object of study in [8], in which a novel method known as attention-based noise-robust averaging (ANRA) has been developed to solve for the tissue's thermal diffusivity a given the free response.
We can then turn our attention to the stationary forced response, in which we study the time during which heat was actively supplied to the tissue. In this case, provided that a constant non-moving heat source is present, after the initial thermal wave has subsided, we observe a forced stationary response, governed by the following dynamics: where p 0 is the fixed position of the probe, and u 0 is the constant power supplied. In what follows, we shall discuss our model of q, and the method by which we obtain it. It must be noted that a is assumed to be constant throughout the entire cut. The latter assumption is justified, considering that during electrosurgery, undamaged tissue is constantly contacted by the moving probe. In contrast, damaged tissue is often immediately ablated, or is confined to a small region around the cut [10]. It also bears mentioning that obtaining a and q from the free and forced response can be achieved through a single stationary probing event, in which the tissue is heated (as dictated by (9)) and subsequently allowed to cool down (as governed by (8)). Finally, we may obtain the thermal relaxation time τ by studying any form of regular cutting action, as governed by (7). The exact method of estimating τ will be discussed in the following section.

III. PARAMETER ESTIMATION METHOD
In this section, we present our parameter estimation method to solve Problem 1. We identify three distinct steps: A) Determining the thermal diffusivity a; B) Determining the heat source q; C) Determining the thermal relaxation time τ . It will become clear that these steps must be performed in order using both the stationary and moving thermal response, as illustrated in Fig. 1.
Our method relies solely on surface thermographer readings, and thus serves to obtain a model of the apparent thermodynamics of the tissue, as observed on a 2D surface. Extensions to a 3D domain based on 2D observations have been discussed in related work [8], where an observer for 3D tissue temperature has been developed; such an extension is, however, beyond the scope of the current work.
To enable direct estimation of parameters, we leverage a set of noise-robust gradient operators introduced by P. Holoborodko [12], [13]; these combine the capabilities of finite difference operators with those of a low-pass filter, allowing for noise-robust gradient computation [8]. In essence, our method reduces the parameter estimation problem from an optimization problem over the solution of (5) to that of balancing differential equations, as will be shown next.
Since it was previously hard to obtain gradient field from thermography data on account of high-frequency noise, such gradient-based balancing techniques have not existed before the concept was introduced by the authors in [8].
In the following, we shall denote noise-robust gradient operators by∇ and∂ t . We now proceed with the parameter estimation method for the thermal diffusivity.

A. Thermal Diffusivity Estimation
We require that (8) holds, with a sufficiently large magnitude ofẋ to obtain the thermal diffusivity. As mentioned in Sec. II-B, this is achieved by raising the temperature of the tissue by application of a small amount of sustained electrosurgical action, after which the free response is studied. Focusing on (8), it is apparent that solving for at each η ∈Ω, will yield a fieldâ as opposed to a single scalar. Many of these parameter estimates will be erroneous, safe for those estimates obtained in regions with a high rate of temperature change [8]. This observation has previously led to the development of attention-based noise-robust averaging (ANRA) [8], where an attention layer is constructed to allow for an average parameter to be obtained. We refer the reader to [8] for a full explanation of this method.

B. Heat Source Estimation
We assume that the unit heat source is modeled as a Gaussian heat source of the form [8]: for η ∈ Ω, where q 0 is such that Ω q(η) dη = 1, and ρ controls the spread of heat. Also, the heat source scales linearly with the power applied, denoted by u in (5). We also introduce some additional notation: we shall callẍ the 'temperature acceleration,' and ...
x the 'temperature jerk.' When applied to real-life data, these repeated derivatives are well defined.
Given a diffusivity estimateâ as obtained from the free response in Sec. III-A, we can now consider the forced response of the tissue, i.e., the time when the tissue was subject to stationary and constant electrosurgical action. In other words, we stipulate that the forced part of the tissue response be obtained by subjecting the tissue to a constant power, with a constant electrosurgical probe position.
To ensure that the forced stationary response model (9) is valid, it is necessary to wait for the initial thermal waves to subside; the modeling errors accrued when not doing so will be demonstrated in due detail in Sec. IV. Since this time period is dependent on the thermal relaxation time τ , which has not been obtained yet, we elect to study the second time derivative of the observed temperature field instead. In particular, we considerx max [k] := max i,j |∂ 2 txi,j [k]|; the rationale for studying this value lies in the fact that a small magnitude ofẍ results in (9) being an accurate approximation of (7). In this work, we wait until the maximum temperature acceleration x max [k] is below 10% of the all-time maximum temperature accelerationx max, total = sup k=1,...,Nxmax [k], before assuming that (9) is sufficiently valid; this is illustrated in Sec. IV.
Based on (9), we defineŜ Given that the heat source is located at position η 0 ∈Ω and has power u 0 > 0, we choose to apply nonlinear least-squares fitting to obtain (q 0 , ρ) by considering the residual: where (i, j) are chosen in a small radius around η 0 ; in this work, a region around η 0 where S assumes positive values is considered. This yields estimatesq 0 =Q/u 0 andρ, after taking an unweighted time-average over the period following subsidence of the initial thermal wave, as described previously.

C. Thermal Relaxation Time Estimation
Having obtainedâ,q 0 , andρ, we may now obtain an estimate for the thermal relaxation time τ . The latter is particularly difficult to observe, since its effects manifest only in the case of transient behavior, such as initial application of electrosurgical power, or when disabling the electrosurgical probe. In steady-state cutting, as is often observed in surgical practice, wave-like effects are rarely observed, given that the process is in a quasi-steady-state whereẍ is of very small magnitude [2]. Instead, a diffuse shock front can be observed when cutting speeds are sufficiently high, as is commonly the case [2]; in practice, the resulting Mach cone is so diffuse thaẗ x is not of a significant magnitude to be useful in estimation. This latter shock front will be the subject of study in future work.
The equation we wish to solve for is similar to (10), but is instead governed by (7): where we have used the heat source model obtained in Sec. III-B.
Our approach is mostly identical to that of Sec. III-A (see (10)), with the exception that the variables are replaced by those appearing in (13). To obtain a time average, instead of considering the intensity I[k], we elect to impose an eventdriven method. Any time we have a sufficiently large change in u(t) (i.e., ∂ t |u(t)| >u threshold > 0), we consider the time starting at that change, up until the time when the third timederivative of the temperature,.
.. Thermal diffusivity running estimate as jerk, a commonly encountered concept in mechanics. The latter condition based on the maximum absolute jerk. ..
x max [k] provides a clear means of delineating the time during which the effects of a sudden change are apparent; this is illustrated in Sec. IV.

IV. SIMULATION
To verify the efficacy of the approach outlined in Sec. III, we simulated (5) based on the properties of porcine liver tissue. These parameters have been obtained from the literature [1], [2], [9]. They are as follows: a thermal diffusivity a of 0.144 mm 2 /s, a thermal relaxation time τ of 0.1 s, and a heat source spread ρ of 0.5 mm. In all of the simulations, the heat source has a power of 30 watts when engaged.
The problem domain is a 4 by 6 cm grid, with a grid resolution of 0.2 mm; a discretization time step of 2 ms is employed. We consider two data sets: one in which the heat source is stationary (to obtain a and q), and the other in which the heat source is moving (to obtain τ ). The movement speed is set at 20 mm/s. In all cases, the simulation time is 2 seconds, and the probe is active for 1 second in the stationary case, and 1.5 seconds in the moving case. The simulation was carried out using a custom simulation code based on a finite difference scheme, as described in [2].

A. Thermal Diffusivity
Employing the method introduced in Sec. III-A, we focus on the last 1 second portion of the stationary simulation run. A cumulative estimate of the diffusivity is shown in Fig. 2, expressed as the ratioâ/a, where the estimate is based on the samples obtained until that instant. The estimated thermal diffusivityâ is 0.102 mm 2 /s, which deviates 11% from the ground truth. As one may observe, the estimate quickly converges to the ground truth within 100 samples (0.2 s), showing little response to perturbations over time.

B. Heat Source
As mentioned in Sec. III-B, the heat source must only be estimated after the initial heat wave subsides, as illustrated bŷ x max . We find that the sampling index at which magnitudes below 10% of the maximum acceleration are maintained is 250 (0.5 s).  We apply our heat source fitting method at index 250, after the thermal wave has subsided, yielding the results shown in Fig. 3. Here, we have shown the heat source fit with respect to the available data based on (9) (i.e., a Fourier heat equation assumption); we have compared this data to the actual case of the Maxwell-Cattaneo model, in which we have considered (7) using the ground truth value of τ = 0.1 seconds. As can be seen in Fig. 3(b), we have obtained an adequate fit with respect to the unknown Maxwell-Cattaneo-based heat source. We may assess the accuracy of our model by considering the spread ρ in (11), which is obtained to be 0.56 mm, as opposed to the ground truth 0.5 mm.
To demonstrate what the results would be if we considered the earlier stages of the forced response when the thermal acceleration is still high (i.e., the thermal wave has not subsided), we have developed the same estimates as in Fig. 3 at index 50 (0.1 s), as shown in Fig. 4.

C. Thermal Relaxation Time
To obtain the thermal relaxation time τ , we apply the method introduced in Sec. III-C on the moving probe simulation. At each time index, we have given the estimated thermal relaxation constant, as shown in Fig. 5. Clearly, many of the values are erroneous, except for the initial few values, and more accurately, the values that appear shortly after the heat source is disengaged (as indicated by the dotted vertical line). Following the averaging criteria given in Sec. III-C, we obtain an estimate ofτ = 0.09 s, deviating 10% from the ground truth τ = 0.1 s.

V. APPLICATION
In this section, we apply our proposed method to the thermal response of in vivo porcine liver tissue to electrosurgical excitation. The trial was performed at the University of Illinois at Chicago's (UIC) Surgical Innovation Training Laboratory (SITL), in accordance with NIH Vertebrate Animal Section standards. The specimen was provided by the Biological Resources Laboratory of the UIC. After exposing the liver tissue, a stationary point probe was performed using a custom robotic platform based on a Franka Emika Panda robotic arm, equipped with visual and thermal sensors (four Optris Xi 400 thermographers (Optris GmbH, Berlin, Germany)). The electrosurgical generator used was a Valleylab Force FX (Valleylab Inc., CO, USA), applying the monopolar pure cut mode at a power of 5 watts. For the stationary cut, a Covidien Valleylab 2.8" needle electrode (Medtronic Inc., MN, USA) Fig. 6. Maximum temperature recorded during the stationary cut trial on live porcine liver tissue. At 2.8 seconds, the power source is engaged, with disengagement taking place at 5.6 seconds, as indicated by the dotted black lines.
was inserted at 5 mm depth.
To illustrate the various segments of the stationary probe experiment, the maximum temperature recording during the trial over time is shown in Fig. 6. During the interval 2.8-5.6 seconds, the electrosurgical unit was engaged. Similar to the procedures in Sec. IV, the period during which the diffusivity was estimated was chosen to be 6.2-12.4 seconds based on the second derivative of temperature. Applying the same method as in Sec. III-A, we obtain at thermal diffusivity of 0.222 mm 2 /s, which is about 53% higher than the nonperfused ex vivo value reported in [9, Tab. 2.3, p. 16] [14]. This difference likely stems from the presence of blood perfusion, which causes an increase in effective diffusivity of 1.4-1.5 [9, Sec. 2.1.5.2] , consistent with our result of 1.53 here.
Having obtained the thermal diffusivity, we may now estimate the heat source model. We delineate the period during which the heat source model is estimated based on the maximum absolute temperature acceleration criterion (Sec. III-B), obtaining a time period of interest between 3.8-4.8 seconds. In this case, the heat source estimate field is much more spurious when compared to Sec. IV as shown in Fig. 7, but most of the noise has been suppressed by the low-pass properties of the noise-robust gradient operators. In the line corresponding to the 'up' direction, however, one may notice an increase in heat flux around 2.75 mm; this anomaly can be attributed to thermal sensor noise, as well as the presence of particulates (smoke), which inevitably introduce discrepancies. Hence, we have based our heat source model on the average heat flux across multiple directions.
To assess the quality of the heat source, we note that our estimated spread isρ = 2.26 mm, which is reasonable considering the diameter of the probe (1 mm). The correctness of this result will further be discussed in Sec. V-A.
We now consider obtaining the thermal relaxation time. Due to issues of image alignment as a result of tissue motion caused by breathing and uncertainties in the thermographer location, we have chosen to estimate the thermal relaxation time based on the stationary probe response, rather than a moving probe response. In an ideal situation, the latter method would yield more conclusive results, since the change in heating may be much more substantial when compared to what can be achieved by stationary probing, while minimizing unnecessary tissue damage. While tissue motion rejection will  be the subject of future study, in this work we show that an adequate estimate may also be obtained based on the stationary response. Applying the method of Sec. III-C, we find relaxation times as shown in Fig. 8. In the latter figure, we find that the initial engagement of the probe results in two spikes in the thermal relaxation constant estimate. This is a result of the probe temporarily losing contact with the tissue after initial ablation, after which it reconnects with the remaining tissue. This effects can also be observed in the maximum temperature upon close inspection of Fig. 6, where there is a short period of stagnation at the 3.5 second mark. After averaging the relaxation time as described in Sec. III-C, we find a relaxation time ofτ = 0.16 s. This yields a second speed of sound of 1.2 mm/s by (4), which is close to the results reported in [1], [2] of around 2 mm/s. The finite heat propagation time can also be observed in the delay between the times where usable estimates of τ are available in Fig. 8  (around 2τ seconds).
We may observe the effect of neglecting the thermal relaxation constant based on the error between the Fourier-based and Maxwell-Cattaneo-based heat source approximations with the thermal relaxation time being taken as 0.16 s, as shown in Fig. 9. For ease of viewing, we have smoothened the raw heat flux curves using a third-order Savitzky-Golay filter, with a window size of 31 samples. Clearly, the Fourier heat approximation exhibits a large offset close to the probe location, where an error of 10 W/m 2 is observed in Fig. 9(b) at the point of application (0 mm distance); such an error  will yield poorly agreeing simulation results across the entire domain. In the next subsection, we will show that our heat model produces an adequate response when compared to the real-life thermal measurements.

A. Comparison with Simulation
We rerun the simulation using the parameter estimates obtained in the foregoing; we have chosen the initial temperature to be 28 • C, which is the average temperature prior to engaging the electrosurgical unit. To illustrate the agreement between the observed real-life thermal response, and the simulated response, the maximum temperature time series plotted in Fig. 10 provides a clear overview. We note that given the parameters obtained by application of our method, the simulated maximum temperature discrepancy is under 5 • C (6%) with respect to real-life data across the entire simulation run. To provide further evidence of the parameter fidelity, we compare the temperature fields obtained at 2.5 seconds, when the probe is disengaged. In Fig. 11 the observed and simulated temperature fields are shown. Most of the error can be attributed to misalignment of the camera and the liver tissue's curvature, but despite these effects, the maximum temperature error is less than 12 • C (14%). We have previously addressed mismatches between simulation and sensor measurements in detail, proving that an observer structure with output feedback can eliminate these discrepancies [8].  (b) Sim. temperature field. Fig. 11. Comparison between the temperature field at the point of probe disengagement, for both (a) the measured real-life data, as well as (b) the simulated data.

VI. CONCLUSION
We have presented a novel approach to real-time parameter estimation for a hyperbolic heat transfer model using thermographic data, as applied to in vivo porcine liver tissue subjected to electrosurgical action. Unlike past methods, which require the solution of complex PDE systems, our approach, based on our previous work on attention-based noise robust averaging (ANRA) [8], operates directly on thermographic data and requires only the PDE structure to be given. We have applied our approach to a theoretical test case involving a Maxwell-Cattaneo model with parameters based on porcine liver tissue, as well as real-life data from an in vivo porcine liver subject to electrosurgical action. In all cases, our method produces parameter estimates close to the ground truth or parameters reported in the literature. Using the results obtained by our method, we have also presented a comparison between live and dead tissue thermal behavior. In addition, a simulation was run using these new parameter estimates, achieving clear correspondence with the real-life porcine liver thermal response. With this validation, we have shown that our method is capable of in situ real-time minimally invasive high-fidelity thermophysical parameter estimation.
Some of the limitations of this work, planned to be addressed in our future research, arise due to the challenging live tissue geometry and breathing motion, and electrosurgical smoke and particles causing temperature reading distortion. The future work may extend this approach to other thermodynamic models, such as the reaction-diffusion equation. In addition, we plan to develop methods to discern tissue types in real-time based on their thermal response. Related to this is planned work on comparing tissue properties in ex vivo and in vivo tissue, with the goal of deriving robust relationships between the two. We also plan to further investigate thermal shock waves and their significance to energy-based surgery and the tissue thermal response.

ANIMAL USE APPROVAL DETAILS
The animal use protocols (ACC NO: 19-151 and 22-097) were reviewed in accordance with the Animal Care Policies and Procedures of the University of Illinois at Chicago. Original protocol approval: 9/24/2019 (3 year approval with annual continuation required). Current approval period: 8/9/2022 to 7/19/2025.