Real-Time Temperature Rise Estimation during Irreversible Electroporation Treatment through State-Space Modeling

To evaluate the feasibility of real-time temperature monitoring during an electroporation-based therapy procedure, a data-driven state-space model was developed. Agar phantoms mimicking low conductivity (LC) and high conductivity (HC) tissues were tested under the influences of high (HV) and low (LV) applied voltages. Real-time changes in impedance, measured by Fourier Analysis SpecTroscopy (FAST) along with the known tissue conductivity and applied voltages, were used to train the model. A theoretical finite element model was used for external validation of the model, producing model fits of 95.8, 88.4, 90.7, and 93.7% at 4 mm and 93.2, 58.9, 90.0, and 90.1% at 10 mm for the HV-HC, LV-LC, HV-LC, and LV-HC groups, respectively. The proposed model suggests that real-time temperature monitoring may be achieved with good accuracy through the use of real-time impedance monitoring.


Introduction
Irreversible electroporation (IRE) and second-generation high-frequency irreversible electroporation (H-FIRE) are techniques currently being developed to treat malignant tumors when other treatment methods, such as surgical resection or thermal ablation, are not appropriate [1]. It is also being developed for the treatment of arrhythmogenic cardiac tissues [2,3]. This technique uses short (1-100 µs) high-magnitude electric pulses (1-3 kV) to produce an electric field resulting in an increase in transmembrane potential (TMP), leading to cell death within the target zone [4][5][6]. When this TMP limit is reached (∼1 V), nanoscale defects, or "pores", form in the cellular membrane. When numerous pulses are administered, the formation of large, long-lived pores disrupts cellular homeostasis, resulting in cell death via different mechanisms [7,8]. Because cell death with IRE is primarily dependent on the generated TMP, IRE is classified as a non-thermal ablation method and may thus be used in a variety of settings. An electric field is produced between two or more needle electrodes placed directly into the tumor-containing area prior to pulse administration. In some cases, a single-insertion bipolar probe (two electrodes integrated into the same cylindrical shaft) is employed.
It has been demonstrated with IRE that ablation and thermal effects can be mutually exclusive phenomena [9][10][11]. Innumerable follow-up studies have found appropriate parameters for ablation in various contexts. However, unfavorable side effects caused by temperature increases continue to be the most important issue limiting the size of the IRE ablations [12]. Clinically, discrepancies in pulse delivery methods and clinical expertise might result in variable degrees of thermal tissue damage, which is most likely to blame for the wide range of complication rates and oncological outcomes documented in the literature [13][14][15][16]. Even at low voltages, the steep potential gradient at the electrode borders can generate quite large electric fields (and consequently temperatures), causing tissue The model is based on the development of a state-space model (SSM) that considers tissue as a system with inputs (applied voltage, tissue conductivity, and impedance changes) and outputs (tissue temperature at 4 mm and 10 mm). The state-space model would be easy to implement for estimating temperature from real-time impedance measurements during treatments while keeping all of the crucial aspects of more complicated models. Hence, it may be used for real-time assessment of tissue response or as a theoretical foundation for more sophisticated temperature controls. In this work, the data-driven model was built from experimental results in an agar tissue phantom and then validated using established theoretical models based on numerical modeling.

AGAR Model
The experimental design was based on a tissue-mimic agar phantom ( Figure 1). Four groups were evaluated to create a comprehensive model: high-voltage high-conductivity (HV-HC), low-voltage low-conductivity (LV-LC), high-voltage low-conductivity (HV-LC), and low-voltage high-conductivity (LV-HC). Data acquisition of n = 5 for each group was conducted (N = 20). A separate agar phantom was created for each sample. Deionized water combined with 1% agar (w/v) and 0.034% or 0.24% NaCl (w/v), to achieve 0.1220 or 0.526 S/m conductivities, respectively, were mixed and heated to 90 • C. Mixed agar solution was poured into a plastic cup mold with a top diameter of 71 mm, a bottom diameter of 46 mm, and a height of 66 mm and allowed to cool to room temperature. Conductivities of the agar were measured at 37 • C to ensure accurate conductive properties of tissues at physiological temperatures. The low conductivity model was set to mimic the conductivity of brain tissue while the high conductivity model was set to mimic initial conductivity of pancreatic tissue at a characteristic frequency of 1.8 kHz [32,33]. Schematic of the physical experimental setup for recording temperature and impedance measurements during treatment delivered from a high-voltage (HV) generator. A custom 3D printed electrode and fiber-optic temperature holder was fit on top of the agar tissue phantom. Two fiber-optic temperature probes were fed into the canals at a distance of 4 and 10 mm away from the electrode surface. Impedance measurements were calculated from voltage and current measurements collected on the oscilloscope. A thermal camera was placed at a distance 25 cm from the flat edge surface of the tissue phantom.

Treatment Parameters
Biphasic pulsed electric fields were administered using a custom pulse generator (VoltMed Inc., Blacksburg, Virginia). Treatment voltages were applied at either 1250 V or 2500 V (achieving voltage to distance ratios of 834 and 1666 V/cm, respectively), and delivered at a rate of 1 burst per second for 300 s. A 5-5-5-5 µs waveform with 100 µs of energized time was delivered through a single insertion bipolar probe (diameter = 1.65 mm; electrode exposure = 7 mm; insulation = 8 mm) with voltage and current waveforms monitored with a WaveSurfer 3024 z oscilloscope (Teledyne LeCroy, Chestnut Ridge, NY, USA) equipped with a 1000× high voltage probe (BTX Enhancer 3000, Holliston, MA, USA) and a 10× current probe (3972, Pearson Electronics, Palo Alto, CA, USA).

Real-Time Impedance Acquisition
FAST impedance measurements were recorded using a 3024z Oscilloscope (Teledyne LeCroy, Chesnut Ridge, NY, USA) and a 1× current probe (2877, Pearson Electronics, Palo Alto, CA, USA). The recorded voltage and current waveforms from the custom FAST waveform ( Figure 2A) were converted into the frequency domain ( Figure 2B), and the complex tissue impedance was calculated as a ratio of the two following Ohm's law. For this study, impedance measurements were calculated from peak values at the 1.8 kHz frequency as the agar phantom is unaffected by tissue electroporation effects. The low-frequency impedance decreases as a function of temperature increase. After a considerable number of pulses, the low-frequency impedance mimics the stable high-frequency impedance. Further details on this technique can be found at [31]. The custom diagnostic FAST acquisition waveform delivered between therapeutic voltage bursts creates (B) a normalized frequency spectrum of the absolute FFT voltage magnitude, which was then used to collect real-time impedance measurements from our agar phantom. Delivered waveforms for the respective low-voltage (C) and high-voltage (D) groups are also shown.

Temperature Acquisition
To achieve reproducible conditions between electroporation treatments, a custommade electrode and temperature sensor holder was developed to precisely position the small fiber-optic temperature sensors (Luxtron m3300; LumaSense, Santa Clara, CA, USA) and single bipolar electrode. The holder was printed on a Form3 printer (FormLabs, USA) with a 25 µm print resolution using FormLabs Clear Resin. The holder consists of a 2.1 mm canal in the center for the electrode, two 1.1 mm canals for the fiber-optic temperature sensors with a center to electrode surface distance of 4 mm and 10 mm, and a depth of 30 mm, placing the temperature probes at a depth parallel with the active top electrode. Prior to treatment, a slice parallel to and at a 10 mm distance from the electrode surface was made to create a flat surface for generation of an immediate post-treatment thermal profile, recorded on a FLIR A325SC thermal camera (FLIR, Wilsonville, OR, USA). Shortly after treatment (within ∼5 s), the agar was sliced along the electrode surface to obtain a thermal profile of the highest temperature areas.

Data Processing
Voltage and current processing and impedance extractions were conducted in MAT-LAB v.R2022a (MathWorks Inc., Natick, MA, USA). Impedance measurements were transmitted to GraphPad Prism (GraphPad Software, San Diego, CA, USA) for fitting to a one-phase decay line to remove noise, and then sent back to MATLAB for difference calculations. Similarly, collected temperature and impedance measurements were loaded into GraphPad Prism for quadratic line fitting to remove noise and then imported back into MATLAB for calculation of temperature changes. All further data processing and statistical analysis were conducted in GraphPad Prism.

Development of the Mathematical Model
A state-space model was chosen because of its capacity to decompose a complex higherorder differential equation into a sequence of first-order equations that can be easily solved, hence reducing computational burden. State-space control utilizes differential equations describing the time domain of the system using state variables in vector form. This makes it easy to evaluate the system using simple matrix algebra, which also enables the evaluation of complex multi-input multi-output (MIMO) systems. Other control system approaches require complicated Laplace transforms and Fourier transforms to transfer the system's time domain representation-supplied as a complex set of differential equations-into the frequency domain-given as a collection of algebraic equations.
The agar tissue mimic was seen as a dynamic system with the applied voltage (V), tissue identification by electrical conductivity (S) at the applied characteristic frequency, and intra-treatment impedance change measurements (∆Z) as the input signals, and the temperatures measured at depths of 4 mm (T 4 ) and 10 mm (T 10 ) as the output signals ( Figure 3). State-space modeling only requires unique state equations for values influencing the system that cannot be deduced from other model states already included; thus, a separate state for applied waveform is not necessary as its influence can be deduced from the tissue conductivity. We consider a system governed by the following continuous-time identified state- where x is the state vector, u is the input vector, y is the output vector, and K is the disturbance component. A, B, C and D are the state-space matrices which describe the system dynamics. By default, D is set to 0 for dynamic systems, indicating that the system has no feedthrough.
The state-space model was estimated by utilizing the numerical algorithms for subspace state-space identification (N4SID). Implementation of this method in MATLAB's System Identification toolbox uses input-output data to construct the model's weighting matrices. This method can be viewed as a type of multi-step, constrained prediction error method utilizing linear regression to solve the system matrix. The algorithm was configured to compute a fourth order SSM. Further details on the methodology for the N4SID method are outlined in Appendix A.

Validation of the Mathematical Model against FEM
After developing an SSM based on results from experimentally collected data, a finite element model was developed to create new test data to assess the validity of the state-space model. A three-dimensional Comsol Multiphysics 6.0 (Comsol, Stockholm, Sweden) model was designed to replicate the experimental conditions of the agar phantom.
The distribution of electrical fields and thermal effects were computed using standard methods [34,35]. The electric potential distribution at the end of a given pulse was calculated with a modified Laplace equation under the electroquasistatic approximation (Equation (4)), and the resulting electric field distribution was calculated with Equation (5).
where σ is the electrical conductivity, E is the local field magnitude, T is the temperature, and Φ is the local electric potential. Following the computation of the electric field distribution, tissue temperature was computed using the general heat conduction equation with the addition of a Joule heating term (Q J ), where ρ is the density of the medium; c p is the specific heat; k is the thermal conductivity; p is the burst on-time (100 × 10 −6 s) and τ is the period of the burst delivery (1 s). The initial conditions were set to match those of the experimental conditions with the outer boundaries of the geometry being assigned a convective heat flux due to air convective cooling with an external temperature of T ext = 18.8 • C. Due to the discrepancies between agar conductivity being measured at physiological temperatures (37 • C), set conductivities (σ 0 ) were scaled by the temperature of the agar at the time of pulsing (T agar ) to match the conductivity of the agar at the time of pulsing (σ f ) by using the following: producing initial conductivities of 0.331 S/m, 0.085 S/m, 0.093 S/m, and 0.331 S/m for the HV-HC, LV-LC, HV-LC and LV-HC groups, respectively. Material properties of the agar phantom were set to those of water. All other thermal and electrical conditions may be found in Table 1.

Experimental Data Results
The recorded temperature and impedance change data points were used as training data along with known parameters such as voltage and tissue conductivity ( Figure 4).

Model Validation with Training and Validation Data Sets
Selection of the final model was chosen through analysis of the accuracy of the training data being tested through the model, as well as a validation set in which specific data sets from each group were omitted from training the model.
The training and validation data set accuracies are shown in the heat map tiles ( Figure 5) for T 4 and T 10 . Specification of the grouping where a data set was omitted for use as the validation data set is specified by the column titles. The fit accuracy was calculated by: where y is the validation data point value and y pred is the predicted model from the output.

Model Implementation with Test Data
We validated our data-driven state-space model against a COMSOL finite element model, replicating our experimental setup. The calculated outputs of the COMSOL model were impedance and temperature at 4mm and 10mm for each voltage and conductivity condition. By importing the extracted impedance into MATLAB, we were then able to obtain predicted temperatures with our state-space model, trained on experimental data. These predicted temperatures were then used to calculate fit accuracy and maximum absolute error (Table 2).
We found that the mathematical model provided excellent estimations of temperature rise, with above 90% fit accuracy for all conditions, except for the LV-LC condition. The discrepancy for this point can be attributed to the low absolute change in temperature. At the 10 mm distance, the temperature increase is negligible, and the model was not able to fit such a small change. Though the LV-LC had the lowest fit accuracy, the absolute error in temperature prediction was under 0.5 • C. Further, all the test conditions had a maximum absolute error of under 0.5 • C, except for the HV-HC, which has a maximum absolute error of 1.45 and 0.52 • C for the 4 mm and 10 mm positions, respectively. The state-space model temperature predictions and finite element temperature simulations at both distances are given in Figure 6 for each voltage and conductivity condition.  The overall model quality was found to be 9.116e-10 as calculated by Akaike's Final Prediction Error (FPE): where d is the number of estimated parameters and the loss function, V, is represented by: where N represents the number of values in the estimation data set, (t) is the respective prediction error, and (θ) N is the estimated parameters.

Discussion
The goal of this study was to evaluate the feasibility of real-time temperature monitoring during an IRE or H-FIRE procedure without external temperature probes or devices. This feasibility study demonstrates a proof-of-concept construction of a black-box model to forecast tissue temperature rise during the ablation process based on observations of real-time impedance changes. In addition to not requiring additional sensors, one of the greatest potential advantages to implementing a state-space model is its ability to provide spatial data at several locations which would otherwise be challenging, if not impossible, in vivo. The model yielded a reasonably good estimate of temperature rise in most cases; however, in certain instances, considerable inaccuracies were detected (e.g., those in which the total treatment temperature rise was on average < 1.5 • C). Despite these discrepancies at low temperatures, temperature rises of this magnitude are not relevant in the context of clinical temperature monitoring. From Figure 7, it is evident that the conductivity of the agar plays the greatest role in temperature fluctuations and thus suggests that thermal monitoring may only be necessary in the treatment of tissues with high conductivities (e.g., pancreas, prostate, etc.).
Due to the fact that the electrical and thermal properties of the gel may differ from those of normal tissue-for instance, the gel was cooler and was not perfused-the absence of our model accounting for flow-rate perfusion, which in the clinical setting should permit a considerable degree of variability, means that these temperature curves cannot be directly applied to the clinical environment. Even though we have focused on the development of wholly data-driven prediction models, future computational modeling could be used to incorporate more biophysical information in the model's input parameters as a means for expanding the utility of the model without making it entirely data-driven. This strategy could improve the accuracy of predictions, as we believe the proposed technique could be enhanced by combining new data from a range of electrical and thermal tissue properties. Moreover, state-space models built on computationally driven data could make it simpler to create multi-output models with more locations for temperature predictions. Future studies may also consider developing a model which considers a multi-electrode setup. However, such geometries would introduce additional variables to be fit including electrode spacing and exposure. In addition, the process could be more complicated in the tissue domain due to the complex impedance spectrum in living tissue [38][39][40][41][42]. Because low-frequency impedance measurements are primarily confined to the extracellular region prior to tissue electropora-tion, membrane permeabilization during treatment has a substantial effect on impedance alterations. However, high-frequency impedance measurements, which correspond to currents that short the membrane reactance and penetrate the cell membrane, are less vulnerable to membrane pore formation. Thus, in the case of a tissue model, high-frequency impedance measurements may be used as a benchmark to identify impedance changes caused primarily by thermal effects.

Conclusions
State-space models can be created to predict tissue temperature increases during H-FIRE treatment by utilizing impedance changes calculated from rapid impedance spectroscopy, applied voltage, and the conductivity of the tissue as described at the characteristic frequency of the applied waveform. The best model produces a reasonably accurate prediction of tissue temperature with an overall model output to test data fit of 95.8% and 93.22% for T 4 and T 10 , respectively. The largest model discrepancies are observed in the scenario where the least amount of temperature rise is observed overall (temperature increases < 1.5 • C). Test data validation of the model suggests that the model is successful at reasonably predicting temperature increases in cases where temperature increases are of clinical relevance.

Abbreviations
The following abbreviations are used in this manuscript:

N4sid Oblique Projection by LQ Decomposition
The N4SID subspace method uses LQ decomposition, singular value decomposition, and applies an oblique projection approach for input and output Hankel blocks to produce subspace weighting matrices [43,44]. Past and future block Hankel matrices for inputs U p , U f and outputs Y p , Y f are constructed from experimental input and output data, The N4SID subspace method determines that our output matrix, Y k (t), can thus be found through the following relationship: . . .
with O k representing the observability matrix, X 0 representing the state initial state, Ψ k representing the Toeplitz matrix designating the shift-invariant properties of the system, and U k is the input matrix. Thus, the output matrix can be represented by In order to solve for the solution to our subspace weighting factors [A B C D], we must perform LQ decomposition of the data matrix (W p ) into the product of a lower triangular matrix (L) and a unitary matrix (Q): The resulting equations are as follows: Through substitution of terms, we end up with the following relationship: suggesting that the two terms U f and W p exist in two different sub-spaces with no overlapping bases. Thus from this realization, we can conclude that from Equation (A2), and use singular value decomposition of R 32 R # 22 W p to solve for O k and X f .