Calibration of a three-state cell death model for cardiomyocytes and its application in radiofrequency ablation

Objective. Thermal cellular injury follows complex dynamics and subcellular processes can heal the inflicted damage if insufficient heat is administered during the procedure. This work aims to the identification of irreversible cardiac tissue damage for predicting the success of thermal treatments. Approach. Several approaches exist in the literature, but they are unable to capture the healing process and the variable energy absorption rate that several cells display. Moreover, none of the existing models is calibrated for cardiomyocytes. We consider a three-state cell death model capable of capturing the reversible damage of a cell, we modify it to include a variable energy absorption rate and we calibrate it for cardiac myocytes. Main results. We show how the thermal damage predicted by the model response is in accordance with available data in the literature on myocytes for different temperature distributions. When coupled with a computational model of radiofrequency catheter ablation, the model predicts lesions in agreement with experimental measurements. We also present additional experiments (repeated ablations and catheter movement) to further illustrate the potential of the model. Significance. We calibrated a three-state cell death model to provide physiological results for cardiac myocytes. The model can be coupled with ablation models and reliably predict lesion sizes comparable to experimental measurements. Such approach is robust for repeated ablations and dynamic catheter-cardiac wall interaction, and allows for tissue remodelling in the predicted damaged area, leading to more accurate in-silico predictions of ablation outcomes.


Introduction
Thermal treatments are common procedures performed on patients with cardiac arrhythmias. For severe cases, thermal energy is typically delivered to the target arrhythmogenic tissue via radiofrequency (Huang and Wood 2014), inducing thermal damage. At a cellular level, prolonged exposure to high temperatures results in protein denaturation and irreversible damage (Pearce 2013, Liu et al 2018. Moreover, a delayed cellular death has also been observed for cells that suffered injuries beyond recovery, through intrinsic mechanisms such as, among others, apoptosis or necroptosis (O'Neill et al 2011, Pearce 2013. While Penne's Bioheat equation is the standard to describe the power delivery to the tissue in computational models of radiofrequency ablation (RFA, see, e.g. Singh and Melnik 2020), several models have been introduced to assess the irreversible cellular damage and estimate the resulting lesion size. As irreversible damage is expected for temperatures above 50°C (Haines 2011), the simplest approach is the use of the 50°C isotherm contour to identify the lesion boundary (Petras et al 2019, Ahn and Kim 2021, Gómez-Barea et al 2022, Gu et al 2022. This approach is typically overestimating the myocardial lesion size (Haines 2011, Liu et al 2018 and does not allow a physiologically accurate assessment of the lesion, particularly in the case of dynamic catheter-tissue interactions. Another approach uses a two-state Arrhenius model (Pearce 2013, Qin et al 2014 to model the transition of the Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. cells from native to denaturated state. Despite being widely used, in particular in RFA (González-Suárez et al 2022), it is well known that this model greatly overpredicts the shoulder region at low temperatures, where the dynamics of the thermal damage are slower (O'Neill et al 2011, Pearce 2013, Liu et al 2018. This shortcoming is particularly relevant for thermal damage of myocytes, since it is known that at 48°C the damage is reversible for all treatment times (Zaltieri et al 2021), while irreversible damage is expected between 50°C and 56°C for up to 60 s of ablation (Haines 2011). The recovery of the tissue damage, correlated to the heat shock proteins HSP70 expression, among other intrinsic protective mechanisms (Rylander et al 2005), cannot be captured by the twostate Arrhenius model, thus leading inevitably to inaccurate injury assessment. In this direction, the three-state thermal cell damage model initially presented in Park et al (2016) is capable of capturing the reversible damage of a cell, which makes it appealing for lesion estimation in radiofrequency catheter ablation models (Qadri et al 2017, Singh and Melnik 2019, Molinari et al 2022a, 2022b, although the model itself is not specific to myocytes. As a matter of fact, not only was the three-state model calibrated to bovine chordae tendineae tissue, but also all the existing two-state Arrhenius models were calibrated for a variety of different cells (Hep G2 human liver hepatocellular carcinoma and MRC-5 human lung fibroblast (O'Neill et al 2011), AT-1 prostate cancer cells, HDF human dermal fibroblasts, LNCaP prostate tumor cell and HuH-7 liver cancer cell (Qin et al 2014) and CHO Chinese hamster ovary cells, among others (Pearce 2013)), none of which were myocytes: hence the direct application of these models in the assessment of lesion formation is, at best, not advisable. Moreover, many cell types exhibit higher energy absorption when exposed to low temperatures and lower energy absorption when exposed to higher temperatures (Dewey et al 1977, Pearce 2013. None of the existing Arrhenius-type models of the three-state from Park et al (2016) is taking this aspect into account, as they all consider a unique absorption rate instead.
In this work, in order to surpass all the above limitations, we consider a modification of the three-state model from Park et al (2016) that accounts for the variable energy absorption rate and is calibrated for myocytes. In particular, we modify the transition phases using a piecewise linear Arrhenius model, that changes slope at 55°C, a decision supported by experimental evidence (Dewey et al 1977, Pearce 2013. The calibration of the model for cardiac myocytes is then twofold. On the one hand, the model returns a response of the cardiac tissue to thermal exposure in line with available experimental data (Haines 2011, Zaltieri et al 2021. On the other hand, when coupled with our previously validated computational model (Petras et al 2019), the model reproduces lesion sizes comparable to the in vitro experiments in Petras et al (2019) for standard ablation protocols and the reported lesions in Kotadia et al (2019) for High Power Short Duration (HPSD) protocols. As an illustration of its capabilities, the calibrated model is then applied to assess the thermal damage in the case of repeated HPSD ablations and in the case of the wall-induced movement of the catheter during the cardiac cycle.

Materials and methods
We consider here the three-state thermal cell death model initially proposed in Park et al (2016) as a starting point. The model consists of three cellular states: the native (N) state, representing a normally functional cell, the unfolded (U) state for cells that suffered some damage that might be reversible, and the denaturated (D) state, where the cell is irreversibly damaged. The interaction of the three states can be summarized as where k i (T), i = 1, 2, 3 are temperature-dependent transition rates among the different states, T denoting the temperature. The interactions in (1) are governed by the system of ordinary differential equations The system is initialised by considering all cells to be in native state, i.e. N = 1, U = 0 and D = 0. Due to the limited cardiomyocyte proliferation in mammalian hearts (He and Zhou 2017), no cell proliferation is considered in the model and a conservation equation holds at all times for the three states (N + U + D = 1).
In the original formulation of Park et al (2016), the transitions among the phases are modelled using an Arrhenius type model , where A i is the frequency factor (in s −1 ), ΔE i is the activation energy (in Jmol −1 ) and R is the universal gas constant. However, this approach does not account for the difference in energy absorption observed following a temperature increase (Qin et al 2014). In particular, a change in the rate of energy absorption of cells has been reported, typically occurring at around 43°C (Dewey et al 1977, Pearce 2013. Several studies have considered a piecewise linear energy absorption rate for thermal damage, with higher energy absorption rates for temperatures below some given thresholds, which vary from 43°C to 53°C (Pearce 2013). Other data reports much higher absorption rates at 50°C than at 70°C, which appear consistent with experiments using different cell types (Qin et al 2014). In particular, figure 7(b) from Qin et al (2014) shows a large decrease in the absorption rates between 50°C and 60°C and a smaller one as the temperature increases further to 70°C.
Following this experimental evidence, we consider a discontinuity in the energy absorption parameters of the transition k 1 between the native and the unfolded states at 55°C. The model parameters are calibrated as described in section 3.1, and summarized in figure 1 panel A. As shown in panel B of the same figure, the resulting transition rate k 1 between the N and U states is continuous across T = 55°C, and exhibits a slope change at T = 55°C, reflecting the regime change in absorption rates.
Following Park et al (2016), the overall frequency factor and activation energy are given as A = A 1 A 2 /A 3 and ΔE = ΔE 1 + ΔE 2 − ΔE 3 respectively. For T 55°C we obtain A = 2.97 × 10 70 s −1 and ΔE = 448.4 × 10 3 J mol −1 , which is consistent with the experimental data for 50°C in Qin et al (2014). In the case T > 55°C, the overall frequency factor and activation energy become A = 1.19 × 10 19 s −1 and ΔE = 125.5 × 10 3 J mol −1 respectively, which are within the values for 70°C in Qin et al (2014).
Different intrinsic mechanisms can contribute to the slow death of a damaged cell, including apoptosis, necroptosis and autophagy. Models of apoptosis and the caspase mechanisms involve complex protein dynamics and would be suitable for the estimation of delayed thermal cell death (Pearce 2013). In this work, we follow a simpler approach similar to O'Neill et al (2011), which uses a threshold value for the cells in native state N, below which cell death is expected to occur. Since no data is available for cardiac myocytes, we consider the irreversible damage threshold to be N 0.8 as in O'Neill et al (2011).

Model calibration I: physiological response
The parameters in figure 1 have been first calibrated to provide physiological response for myocardial cells. In particular, the following conditions should be satisfied.
(1) No damage should be observed at body temperatures of 37°C.
(3) The damage is reversible for any ablation treatment time at temperatures up to 48°C (Zaltieri et al 2021).
To test these cases, we simulated a commonly used experimental setup, in which the cells are placed in a thermal bath at a given temperature (Rylander et al 2005, O'Neill et al 2011. A constant temperature is applied to the tissue for a predefined period of time, followed by a thermal relaxation at body temperature of 37°C. Figure 2 shows the model response by applying three different temperatures for 60 s followed by a thermal relaxation period at a constant temperature of 37°C. We allow the system to reach a stable state by using a final time of 25 min, in which the variations of the different states are below 5 × 10 −6 . A 25 min exposure to a temperature at 37°C (left panel) returns a variation in the level of native cells of about 1%, predicting no damage at all at body temperature. At temperatures of both 43°C (central panel) and 48°C (right panel), some initial damage is observed, but it always remains reversible as expected for myocytes. Note that we chose a duration of 60 s for temperature administration, which corresponds to typical cardiac ablation applications. In our numerical experiment we do not consider a temperature ramp, which is what is typically obtained from ablation studies (Toupin et al 2017, Zaltieri et al 2022, thus overestimating the exposure at the given temperatures. Despite this, the model still predicts a reversible damage as expected.

Model calibration II: lesion size assessment
For the second calibration step, we coupled the three-state model with our previously developed computational model for RFA, that is outlined in appendix B. For a comprehensive description of the model, we refer the reader to Petras et al (2019). The simulation protocol depends on the applied power and the flow rate of the cooling saline. It is organised in three phases as shown in table 1.
The duration of the ablation phase depends on the applied protocol (30 s for standard, 4 s for HPSD), while the thermal relaxation continues until the rate of change of the three states drops below 5 × 10 −6 . Irreversible damage is assessed at the end of the thermal relaxation and it is identified by the threshold N 0.8. The results are compared against the in vitro experimental results provided in Petras et al (2019) for two standard ablation protocols and against experimental data from Kotadia et al (2019) for high-power short-duration (HPSD) (see table 2).
For standard protocols we simulated virtual ablations of 20 W for 30 s on a porcine tissue, with applied forces of 10 g, 20 g, using a 6-pore open-irrigated catheter with spherical electrode tip (see Petras et al 2019 for the electrode specifications). The resulting lesion dimensions are in good agreement with the experimental data from Petras et al (2019), particularly for the depth. The lesion width W is consistently underestimated, but such shortcoming follows from a limitation of the RFA model in Petras et al (2019) which does not account for anisotropy in the electrical and thermal conductivities due to the cardiac fiber orientation.
For high-power short-duration (HPSD), we performed a virtual ablation of 90 W for 4 s on a simulated human ventricular tissue using a force of 10 g and a 6-pore open-irrigated catheter with a cylindrical tip (see the supplementary material of Petras et al (2021) for the electrode specifications). The resulting depth compares well with the experimental data from Kotadia et al (2019), however the lesion width is underestimated by the model, for the same reasons as before. Figure 2. Response of the three-state tissue model to a 60 s exposure to a constant temperature (37°C, 43°C, 48°C) followed by a relaxation at body temperature (37°C). The native state N lies above the 0.8 threshold after the relaxation and the inflicted damage is always reversible.  Figure 3 takes a closer look at the results of the 90 W, 4 s, HPSD ablation protocol. As expected, the 50°C isotherm predicts a larger irreversible thermal damage than the three-state model (panel A). Panel B presents the temporal evolution of the temperature in three different locations, whose positions are shown in panel A: one at the edge of the 50°C isotherm (Node I), the other two closer to the electrode (Nodes II and III). A thermal latency effect is observed deep in the tissue, as shown by the temperature evolution for Node I (solid line): after the ablation stops, the temperature in Node I keeps increasing for a few seconds and remains above 50°C before starting to drop. This effect may be suggesting that the lesion size is underestimated when the 50°C isotherm is used. On the contrary, our simulations show that despite the thermal latency, the irreversible thermal damage is still smaller than the one predicted using the 50°C isotherm. Panel C collects the solution in time of the threestate model for the three selected nodes, showing that the damage in Node I, that lies at the edge of the 50°C isotherm, remains reversible, while Nodes II and III are irreversibly damaged.

Repeated HPSD ablations
Similarly to Petras et al (2021), we consider an ablation protocol with 2 subsequent HPSD applications of 90 W for 4 s with an in-between pause of 8 s, followed by 30 s of thermal relaxation. Using the model in Petras et al (2021) with a contact force of 5 g, the resulting evolution of the depth of the predicted damage is displayed in figure 4 using the 50°C isotherm contour and the three-state model. The depth of the damaged area predicted by the isotherm becomes zero before the second application, which is a non-physiological result. In comparison,  the thermal damage depth predicted from the three-state model decelerates, however the pause between the ablations is too small to observe a shrinkage. The final lesion size after the two applications is once again overestimated by the 50°C isotherm.

Ablation with catheter movement
We test the proposed method on the estimation of the irreversible damage region with a computational model that accounts for the catheter movement (as described in appendix B).

Discussion
In this paper, we propose a modification of the three-state thermal cell death model in Park et al (2016) calibrated for cardiac myocytes. The model features an Arrhenius transition with variable absorption rate from native to unfolded states accounting for the reduction in the energy absorption as the temperature increases (Qin et al 2014). A similar approach has been reported in the literature for other types of cells (Dewey et al 1977, Pearce 2013). The proposed method is capable of capturing the reversible damage that occurs due to intrinsic cellular protection mechanisms, which have been linked to the expression of heat shock proteins HSP70 (Rylander et al 2005). Irreversible damage due to direct protein denaturation is modeled by (1) The use of the long relaxation times until the rate of change of the three states is below 5 × 10 −6 is reasonable considering the lengthy process of lesion evaluation after the ablation process in the experimental setups. This technique can also be used in vivo, since the irreversible damage area should be predicted after the tissue healing. However, other factors may affect the cellular viability in vivo, including the extent of the edema formation after the ablation (Yano et al 2020).
The three-state model is capable of capturing the inflicted damage on the tissue in cases of ablation models with catheter movement as well as models that use ablation protocols with multiple applications. For such models, the 50°C isotherm may provide in some cases a crude estimation for the lesion size, however it can be very inaccurate and may lead to a serious overestimation of the irreversible damage area. The use of the threestate model allows for the introduction of alterations in the thermomechanical modelling of the tissue (Molinari et al 2022a) and its biophysical properties due to the inflicted thermal damage. These coupling aspects are part of our ongoing work.

Acknowledgments
This work has been partially supported by the State of Upper Austria. This research was funded in part by the Austrian Science Fund (FWF) P35374N. For the purpose of Open Access, the author has applied a CC BY public copyright licence to any Author Accepted Manuscript (AAM) version arising from this submission.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary information files). Data will be available from 23 February 2023. The irreversible thermal damage area given by the 50°C isotherm after the ablation and the three-state model after the relaxation time at the undeformed configuration. The temperature evolution is displayed on two nodes, one close to the surface and one deeper in the tissue. The evolution of the superficial thermal damage area as estimated by the 50°C isotherm and the three-state model. Table A1 collects a comparison with other commonly used lesion size assessment methods. The original threestate model (Park et al 2016) (calibrated for bovine chordae tendinae), with parameters A 1 = 3.68 × 10 30 s −1 , A 2 = 5.68 × 10 3 s −1 , A 3 = 2.58 × 10 5 s −1 , ΔE 1 = 210 × 10 3 J mol −1 , ΔE 2 = 38.6 × 10 3 J mol −1 and ΔE 3 = 47.2 × 10 3 J mol −1 was used in Molinari et al (2022aMolinari et al ( , 2022b. The study in Molinari et al (2022a) uses a physiological threshold of N thr = 0.8, but lacks validation, as stated in the limitations therein: it is a proof of concept for thermomechanical remodelling during ablation and produces quite small lesions. Hence, it will not be used here. In Molinari et al (2022b), the threshold cutoff for irreversible thermal damage is raised to N thr = 0.98 in order to fit their modelling results to the available experimental data.
Another popular approach uses only two states, the native and the denaturated state. As described in Pérez et al (2022), the transition follows the Arrhenius model with A = 7.39 × 10 39 s −1 and ΔE = 257.7 × 10 3 J mol. The irreversible damage is identified using the Ω = 1 isoline 90 s after the beginning of the ablation, corresponding to 63% dead cells, where This isoline is the one used by Pérez et al (2022), despite the lack of physiological background behind this choice. In fact, Ω = 1 was introduced to identify a trans-epidermal necrosis for skin burn (Henriques and Moritz 1947), but it has been used as a threshold for thermal injury to several other types of cells without any prior physiological assessment (Pearce 2013). All models predict a lesion depth that lies within the experimental data for both SP and HPSD protocols. Note that in all cases a consistent underestimation of the width is observed due to the lack of anisotropy in the thermal and electrical conductivities of the RFA model ( where ρ is the tissue density, Σ( · ) is the electrical conductivity tensor, Φ is the electrical potential, c( · ) is the specific heat, T is the temperature, v is the velocity, K( · ) is the thermal conductivity tensor, p is the pressure scaled by the density and T(v, p) = 2μρ −1 (∇v + ∇v T )/2 − pI is the stress tensor of the Navier-Stokes equation, I is the identity tensor and μ is the dynamic viscosity of the blood. The irrigated saline is delivered to Ω fluid through the inner electrode channels at a temperature of 28°C (Squara et al 2014). A dispersive electrode is placed at the bottom of the domain, and a voltage drop is imposed on the interface between the electrode and the catheter shaft. For the remaining boundary conditions of the system we refer the reader to Petras et al (2019).
Model (3) is solved with finite elements through a self-developed code implemented in FEniCSx (fenicsproject.org).

B.3. Catheter movement
The catheter movement is induced by the movement of the cardiac wall that, following the cardiac cycle, features a rapid contraction due to the ejection phase, followed by a slow relaxation and filling phase. Figure B1 shows the wall displacement, based on the left ventricular radial velocity data in Codreanu et al (2014) and using a reference of 85 heartbeats per minute corresponding to a normal heartbeat of a swine (Paslawska et al 2014). During the cardiac cycle, the cardiac wall is stiffening during the isovolumetric contraction (IVC) and remains stiff during the ejection phase (Couade et al 2009, Urban et al 2013. The wall relaxation occurs during the diastole phase (Couade et al 2009, Urban et al 2013, in which it reaches the lowest elasticity value as shown in figure B1. Thus, in opposition to the wall movement, the catheter is located deeper in the tissue during the diastole and the indentation becomes shallow during the ejection phase. Using a contact force that ranges from 5 to 25 g, during the diastole and systole respectively, following values reported in the literature (Aizer et al 2018), the tip displacement is shown in figure B1.
The deformations of the cardiac tissue induced by the catheter tip movement are described by where u is the displacement, σ is the stress tensor and Ω tissue is the tissue subdomain, which is modelled as a nearly incompressible, hyperelastic, and orthotropic material (Usyk et al 2000). The relationship between σ and u is, as usual for hyperelastic materials s = ¶Y ¶ J F F 1 , T with F the deformation gradient, = J F det and the strain energy density The fibers are considered parallel with the x-axis, and only the mechanical properties are considered anisotropic following the fiber distribution. The Fluid-Structure Interaction in our model is a one-way interaction, with the cardiac tissue moving freely as a consequence of the beating of the heart and the fluid domain following the imposed deformation. Algorithmically, we solve a time step of the tissue deformation, followed by a time step of the Navier-Stokes equations and finally a time step of the coupled potential and heat equations.
We have a fixed contact surface, so we simply impose the displacement due to contact as a Dirichlet boundary condition. This is a reasonable approximation for a cylindrical tip electrode (Petras et al 2022).
Equation (4) is paired with system (3) and the coupled problem is defined, at time t, on the deformed geometry Ω ≡ Ω t . In section 3.5 a porcine cardiac tissue is considered, and no blood inflow is imposed, due to the uncertainty in the blood flow velocity: the blood movement is thus driven by the wall movement during the cardiac cycle. In this experiment a cylindrical electrode is used to limit the computational complexity thanks to its fixed contact area with the tissue. The rest of the simulation protocol follows (Petras et al 2019). The coupled model (3)-(4) is solved with finite elements through a self-developed code implemented in FEniCSx (fenicsproject.org).

B.4. Parameters
The parameters of the models are based on Petras et al (2019) and Petras et al (2021) for standard and HPSD ablation protocols respectively. The thermal and electrical properties are isotropic and temperature dependent as detailed therein.