Modeling and simulation of hypothermia effects on cardiac electrical dynamics

Previous experimental evidence has shown the effect of temperature on the action potential duration (APD). It has also been demonstrated that regional cooling of the heart can prolong the APD and promote the termination of ventricular tachycardia. The aim of this study is to demonstrate the effect of hypothermia in suppressing cardiac arrhythmias using numerical modeling. For this purpose, we developed a mathematical model that couples Pennes’ bioheat equation and the bidomain model to simulate the effect of heat on the cardiac action potential. The simplification of the proposed heat–bidomain model to the heat–monodomain model is provided. A suitable numerical scheme for this coupling, based on a time adaptive mesh finite element method, is also presented. First, we performed two-dimensional numerical simulations to study the effect of heat on a regular electrophysiological wave, with the comparison of the calculated and experimental values of Q10. Then, we demonstrated the effect of global hypothermia in suppressing single and multiple spiral waves.


Introduction
Several experimental studies have demonstrated the significant effect of induced hypothermia on cardiac and neurological outcomes for patients (see [1] for a review). Hypothermia is now recommended as a therapeutic treatment for cases of spinal cord and brain injuries (see [2] and [3]), and it is used as a standard treatment for cardiac arrest [4]. Numerical modeling can provide valuable contribution for the understanding of the role of temperature effects in the cardiac electrical dynamics, which is the main aim of this paper.
In literature, the modeling of the effects of temperature on the cardiac electrical wave has previously been performed mostly by modifying the ionic activity. For cells, an ionic model that considers the temperature dependence of electrical parameters was presented in [5]. It has been shown that the variations in the cellular responses due to changes in temperature can have profound effects on the behavior of the transmembrane potential, including action potential durations (APDs) (see [6] and [7]). In 2006, a modified FitzHugh-Nagumo monodomain model combined with the Pennes' equation was proposed [8] to include the influence of temperature on the behavior of a simulated nerve. In 2009, this work was extended [9]  including the effect of temperature in the Hodgkin-Huxley model to further increase the accuracy of the description of this behavior. Similarly, temperature dependence was added to ionic intestine models [10] to estimate its possible effects during surgery. Considering the cardiac dynamics, a cell model has recently been developed [11,12] that includes the influence of temperature on ventricular electrical activity. Regional cooling has been studied in different tissue sizes [12] and it has been shown that it can be suitable as an anti-arrhythmic therapy for small tissue sizes and pro-arrhythmic therapy for large tissue sizes. However, the effect of temperature has been included only in ionic cardiac activity. To the best of our knowledge, the effect of global hypothermia on the cardiac electrophysiological wave, which is the main contribution of this paper, has not been previously reported in literature. Therefore, based on the idea of Bini and coworkers [8,9,11,12], we have developed a mathematical model that combines Pennes' bioheat equation with the bidomain model. The simplification of the proposed heat-bidomain model to the heat-monodomain model is presented. In the proposed heat-bidomain model, heat source terms considering the Joule heating effects in the intra-and extracellular regions are added. These source terms induced by the Joule effect strongly depend on the gradient of the action potential, and in cardiac tissue, this potential has a steep slope during the de-and repolarization phases; thus, the simulation becomes challenging. Owing mostly to this challenge, the spatial temperature effects on cardiac tissue have not been previously studied. This difficulty, reported in [9], limits the generalization of the methodology to the case of cardiac tissue. Therefore, in this paper, we also present a time-dependent adaptive mesh algorithm to address this challenge. The main advantage of this method is that it can concentrate the spatial resolution along the areas with large gradients, enabling the simulation of the heat-bidomain coupling. In this paper, numerical simulations are presented to study the effect of temperature variations on a regular electrophysiological wave and to investigate the effect of hypothermia in suppressing cardiac arrhythmias. This paper is organized as follows. Section 2 discusses the proposed mathematical models for the coupling of the heat with the bi-and monodomain models. The finite element discretization and the time-dependent adaptive algorithm is presented in section 3. Finally, multiple simulations are provided in the last section, to study the effects of temperature variations on the action potential in two-dimensional tissues.

Models
The bi and monodomain models are widely used in electrocardiology to simulate the spatial propagation of the transmembrane potential in the myocardium. These classical models do not consider the effect of temperature on the electrical wave; therefore, heat-bidomain and heat-monodomain couplings are required. conserved, the following elliptic equation can be written: At each point of the myocardial tissue, the transmembrane potential is considered as the difference between the intra-and extracellular potentials, V m = ϕ i − ϕ e . For any biological membrane, the transmembrane current is the sum of the ionic and capacitive currents given by where C m is the capacitance of the cell membrane, χ is the ratio of the membrane surface area to the volume, and according to the current conservation law, The total ionic current across the membranes, I ion , depends on the ionic models. In this paper, we consider two simplified two-variable ionic models. One is the Mitchell-Schaeffer model [13], introduced in 2003, which is derived from the Fenton-Karma ionic model [14]. This model quantitatively reproduces the behavior of the ventricular action potential, and it has similar de-and repolarization isochrons to those obtained by the more complex cardiac ionic models. In this work, the parameters for this model are adjusted to have a fast upstroke, which enables this model to be suitable for our study. The other is the Aliev-Panfilov model [15], introduced in 1996, which can reproduce more realistic shapes of the cardiac action potential as well as the APD restitution characteristic observed in experiments. This model has widely been used for simulating specific types of cardiac arrhythmia, such as those characterized by rotating waves, which are also explored in this study. In this paper, we scaled both the Mitchell-Schaeffer and the Aliev-Panfilov models to obtain physiologically interpretable values, based on the following equations. The Mitchell-Schaeffer model is given by The Aliev-Panfilov model is given by Assuming that the transmembrane potential V m is ϕ i − ϕ e , the bidomain model considered in this paper is given by wC m @V m @t À r � ðG i rV m Þ ¼ r � ðG i r� e Þ À wI ion ðV m ; WÞ; The bidomain model is not capable of studying any temperature effects. To address this problem, we first consider Pennes' bioheat equation [16]. This equation, first introduced in 1948, describes the transfer of heat in biological tissues, and it can be written as where k is the thermal conductivity, T is the temperature of the tissue, T � is the temperature of the supplied arterial blood, b c is the strength of the heat sink due to blood perfusion, ρ is the density of the tissue, and c p is the heat capacity of the tissue. However, this equation does not consider the effect of the heat due to the propagation of the electrical wave across the tissue. This response can be included simply by considering the Joule effect (see [17] and [8]). In our work, we consider the heat generation rate per unit volume from an electric field in the intraand extracellular regions p i and p e , respectively, as These equations are incorporated as internal source terms in Eq (5) as follows: Thus, using the expressions of J i and J e given in 1, the coupling between the heat transfer in the tissue and the transmembrane potential in the myocardium can be given by equation Nevertheless, the temperature also has a significant effect on the ionic behavior in the cells. The above ionic models still do not consider the effect of temperature; therefore, they need to be modified to accurately represent the effects of temperature variation on the APD. For this purpose, we refer to the original study by Hodgkin and Huxley [18], in which the effect of the temperature on the rate of change of the conductance variables is considered. Subsequent modifications were performed [6], and the linear changes of the ionic conductances with respect to temperature were also considered [7]. Based on these studies, the temperature properties can be added directly to the Mitchell-Schaeffer and Aliev-Panfilov models as

Heat-monodomain model
A widely applied method of reducing the computational time of the bidomain is to reduce the two-by-two set of partial differential Eq (4) to a scalar partial differential equation, resulting in the monodomain model. A similar technique can be used for the heat-bidomain model mentioned above. However, only the heat source due to Joule effect terms in Pennes' bioheat Eq (8) needs to be considered in this study. Thus, we assume that the wavefront of the transmembrane potential is not fully curved (see [19], page 110); thus the sum of the intra-and extracellular current densities is required to be equal to zero: Therefore, the heat source term in Eq (8) can be expressed as Under the assumption of equal anisotropy ratios, G i = λG e and from Eq (12) we can write, Thus, the coupling of the heat-bidomain model in Eq (11) can be simplified as where the conductivity tensor G can be considered as 1 1þl G i . Although the derivation of the heat-monodomain was presented, all our numerical simulations were performed using the heat-bidomain model. In addition, the assumption in Eq (12) is restricted to non-curved wavefronts, which is not universal. To overcome this limitation, we suggest, in the case of the bidomain model, the consideration of the Joule heat source generated by the membrane voltage, G i rV m � rV m , instead of that generated by the intra-and extracellular regions, G i rϕ i � rϕ i + G e rϕ e � rϕ e . Therefore, the Pennes' bioheat Eq (8) can be replaced by The numerical experiments [20] show that Eq (14) also provides appropriate results when combined with the bidomain model. For the sake of simplicity, in the following, the heat-bidomain model used for all numerical simulations is given by All parameters for both the heat-bidomain model in Eq (15) and its simplified heat-monodomain version in Eq (13) are determined in the next subsection, enabling these models to be more suitable for studying the effect of temperature on the cardiac transmembrane potential.

Model parameters setting
To study the effects of temperature on the cardiac action potential, first, we need to assign reasonable values to all the necessary constants in the heat-bidomain model in Eq (11). Specifically, these parameters need to be suitable to reproduce the realistic stiffness of the slope in the de-and repolarization phases of the cardiac transmembrane potential. This is crucial to this study, as the source term in the heat equation in the model in Eq (11) strongly depends on the gradient of the action potential. As shown in [21], this can be achieved using the Mitchell-Schaeffer model and, by tuning the parameters of this model, it is possible to obtain an action potential with nearly realistic amplitude, duration, and upstroke velocity. In this study, following the adjustment of the Mitchell-Schaeffer parameters, the suitable values presented in Table 1 provide an upstroke of *1 ms and a conduction velocity of *0.7 m/s at 37˚C. In addition, for both the Mitchell-Schaeffer and the Aliev-Panfilov models, we use a resting voltage of V rest = −85 mV, a peak value of V peak = 40 mV, and a total amplitude of the action potential as V amp = V peak − V rest , to obtain a dimensional version of the electrophysiological wave.
In the case of the bidomain tissue model, several values for the constants are available. In this study, we use the calibrated parameters obtained from [22] for the surface to volume ratio χ, the capacitance C m , and the conductances G i and G e (see Table 2). In [23], a full discussion of the experimentally measured values of the conductances can be found.
For the Pennes' bioheat Eq (8), the available experimental data for the heat properties of cardiac tissue is limited [9]; however, the values of thermal conductivity k, density ρ, and heat capacity c p , were experimentally determined for cardiac muscle [24]. We were unable to find an equivalent experimental value for the metabolic and blood perfusion term b c for cardiac tissue. Nevertheless, in [25] and [26], the heat due to the propagation of the electric wave in both the olfactory and the myelinated nerve fibers was rapidly reabsorbed by the medium (*30 ms). Therefore, despite the possible differences in the behavior of myocytes and nerve cells, we performed a parametric study to reproduce this re-absorption as performed in [9]. The values used in Pennes' bioheat equation are given in Table 3.
As the Mitchell-Schaeffer and Aliev-Panfilov models have not been used previously to study temperature effects, constants A, B, and Q in Eqs (9) and (10) are not available in literature in the context of these cell models. The values chosen for the simulations, shown in Table 4, are within the range of the experimental values found in [11].

Finite element discretization
A second-order mixed finite element formulation in both space and time was used for the bidomain model. For space discretization quadratic polynomials were used. A second-order fully implicit backward scheme (Gear method) was employed for the time derivative discretization. For instance, starting from T n−1 and T n , the Gear scheme gives @T @t ðt ðnþ1Þ Þ ' 3T ðnþ1Þ À 4T ðnÞ þ T ðnÀ 1Þ 2Dt : Using ψ v , ψ ϕ , ψ w , and ψ T as test functions, the overall algorithm for solving the proposed heat-bidomain model on a fixed mesh is as follows.
At each time step, the Newton's method was used to solve the non-linear system mentioned above. The linear systems resulting from the Newton's method were solved by iterative methods using incomplete LU decomposition (ILU) generalized minimal residual (GMRES) solver [27] from the PETSc library [28].

Time-dependent adaptive algorithm
Although a second-order mixed finite element formulation method for both space and time was used, the above-mentioned algorithm presents complex computational challenges. In addition to computational difficulties related to the bidomain model, the source term in the Pennes' bioheat Eq (8) strongly depends on the gradients of the action potential, and as the transmembrane potential is a traveling wave with a very sharp depolarization front in cardiac tissues, the numerical simulation is more intensive computationally and requires extremely fine meshes. Therefore, an accurate numerical method is necessary for suitable simulations of the temperature effect on the electrical waves in the human heart.
Mesh adaptative methods can be more suitable to address these difficulties. The main advantage of these methods is that finer mesh cells can be located near the front, while a coarser mesh can be used away from the front. Therefore, the accuracy of the prediction of the electrical wavefronts is enhanced, and the total number of mesh elements can be greatly reduced, along with the computational time. Several adaptive strategies have been introduced in the context of simulating the cardiac electrical activity (see [29], [30], [31], [32] and the references therein). An adaptive method, based on a hierarchical error estimator, was developed for the two-dimensional simulation of electrical waves in the heart [33]. A three-dimensional adaptive method, based on the definition of edge lengths using a solution-dependent metric, was employed in [34]. In this work, we adopt an adaptive mesh technique based on the definition of edge lengths using a solution-dependent metric. A complete description of this technique is presented in [34,35] and not repeated here. The computational efficiency of our adaptive algorithm was previously assessed on single and complex cardiac wave dynamics [36]. However, the heat-bidomain model presented in this paper requires a different timedependent algorithm which is described in the following.
The overall adaptive algorithm has the following steps: 8. Solve the system in Eq (17) on the mesh M ðnþ1Þ for T n+1 . 9. Next time step: go to step 2.
According to step 5, adapting the mesh using all different numerical solutions at each time step is more favorable, which depends on the time discretization scheme. In this work, a second-order fully implicit backward scheme was employed for the time stepping. Thus, the mesh is required to represent the solutions at times t (n−1) , t (n) , and t (n+1) . We adopted a common metric on which the mesh is adapted by considering the errors on all these solutions. Moreover step 5 indicates that 12 variables need to be considered. This is costly and therefore to reduce the computational time, a linear combination of the different solutions can be adapted. In our experience, the following four variables provide satisfactory results:

Results
In this section, we discuss the performance of the model and the numerical method introduced in the previous sections. First, we present the effect of heat on regular electrophysiological wave. The cardiac transmembrane potential is demonstrated at different temperatures and the results are validated by comparing the calculated and the experimental values of Q 10 , which is a commonly used quantification in biology to describe the rate of change of a system subjected to a temperature increase of 10˚C. Then, the effect of hypothermia in suppressing cardiac arrhythmias is investigated.

Effect of heat on action potential duration and conduction velocity
In this section, a two-dimensional regular wave is presented. The simulation was performed for 750 ms on a sheet of size of 10 cm × 10 cm. The heat-bidomain model combined with the Mitchell-Schaeffer model was used for the ionic representation. All physical parameters are summarized in Tables 1, 2  Six different values were considered for T � : 28, 31, 34, 37, 40, and 43˚C. The initial temperature was set to be equal to the value of T � . A time step of t = 0.3 ms was employed in all calculations. This time step was relatively large because a fully implicit method combined with an adaptive mesh technique was used. However, a smaller time step is be required when using regular meshes or at T � = 43˚C, where the wavefront is sharper. The effects of these changes in temperature on the APD and the rise time is shown in Fig 1. As can be seen in Fig 1a), the increase in temperature significantly decreases the APD, while a decrease in temperature increases it. This is the expected non-linear behavior, compared to the experimental results presented in [37]. Fig 1b) shows that the rise time decreases when the temperature is increased. Several studies were performed using the Q 10 temperature coefficient. The Q 10 value can be used to describe the temperature sensitivity of properties in cardiac tissue, such as the temperature dependence of the APD and the conduction velocity. The experimental value of Q 10 (CV) = 2.3 was determined in rat cardiomyocytes [38] and a value of Q 10 (APD) = 2.5 [39] was determined through measurements of guinea pig ventricular myocytes. Using the results of our simulations, it is possible to measure these values based on the following equations: The APDs were calculated as the duration at which the voltage is maintained above −80 mV at point (6,5). Additionally, the conduction velocity was calculated using the time difference between the first front crossing −80 mV at locations near (6,5). Table 5 shows the obtained calculated values, which are in good agreement with the experimentally measured values in [38] and [39].
Furthermore, owing to the inclusion of Pennes' bioheat equation to the bidomain model, it is possible to study the effect of heat on the electrical wave propagation. Fig 2 shows the transmembrane potentials at T � = 31˚C and T � = 37˚C obtained at a single time instant. As can be seen in the figure, the area of the excited potential is larger for the warmer tissue. This is due to the higher propagation speed at higher temperatures. For these simulations, the Mitchell-Schaeffer model was employed, and at 37˚C, the applied parameters result in a conduction velocity of *0.7 ms and a rise time of 1.07 ms. Higher upstroke and conduction velocities were obtained when the temperature was increased; thus, this model is suitable to demonstrate the performance of the presented adaptive strategy. Fig 3 shows these adapted meshes at different time instants for both T � = 31˚C and T � = 37˚C. This demonstrates that elongated elements are obtained at the appropriate position to capture the transmembrane potential wave. The adapted mesh evolves with time, and at each time step, the adapted meshes are obtained to accurately determine the stiffness of the depolarization wavefront, which enables the feasibility of the technique for this study.

Effect of hypothermia in suppressing cardiac arrhythmias
One of the main advantages of the proposed model and the methodology used in this study is that it enables the investigation of the effect of hypothermia and regional cooling on cardiac tissue. For this purpose, we first study a single spiral wave generated by a similar technique described in [36] and with an initial temperature of 37˚C. We consider the heat-bidomain model combined with the Aliev-Panfilov ionic model. All physical parameters are given in Tables 1, 2 and 3 except for a, which has a value of 0.1. The regional cooling is generated by the source term in the model in Eq (11). In particular, we set T � to 30˚C in a circle of radius of 2 cm centered in the computational domain, while 37˚C is maintained in the rest of the domain. The case of hypothermia is considered by setting T � to 30˚C in the entire computational domain. The results are shown in Fig 4. The first column in this figure shows the transmembrane potential when there is no heat effect (T � = 37˚C), the effect of regional cooling is presented in the second column, and the last column shows the case of hypothermia. As can be seen, the regional cooling prolongs the transmembrane potential, but it does not suppress the spiral wave. However, the hypothermia at T � = 30˚C allows the termination of the spiral during the first 1.68 s interval of the simulation. These results were investigated in more complex spiral waves as well. The study begun with three spiral waves and then the parameters of the Aliev-Panfilov model were set to the values After t = 600 ms we considered two cases of regional cooling, where T � was set to either 30˚C or 28˚C in a circle of radius of 2 cm centered in the computational domain, while 37˚C was maintained in the rest of the domain, and two cases of global hypothermia, where T � was set to either 30˚C or T � = 28˚C in the entire computational domain. The transmembrane potential in the case of regional cooling is shown in Fig 5, where the first column shows the transmembrane potential in the case of no temperature variation (T � = 37˚C), the effect of regional cooling using 30˚C is presented in the second column, and the last column shows the effect of regional cooling at 28˚C. Similarly, the transmembrane potential in the case of hypothermia is shown in Fig 6, where the first column shows the transmembrane potential in the case of no heat effect (T � = 37˚C), the first case of hypothermia (T � = 30˚C) is presented in the second column, while the last column shows the second case of hypothermia (T � = 28˚C). As can be seen, in the case of regional cooling (Fig 5) the spiral waves are not terminated and they are breaking up into multiple spirals, which can be considered pro-arrhythmic. However, hypothermia at T � = 30˚C promotes the prolongation of waves and reduces the number of spirals, while hypothermia at T � = 28˚C terminates all spirals during the first 1.35 s interval of the simulations.
To investigate if the spiral wave was terminated due to the physical boundary, we performed simulations in a region of 13 cm × 13 cm. As shown in Fig 7, the hypothermia can still suppress the spiral waves in this larger computational domain. Based on our numerical simulations, the spiral wave can remain in the tissue depending on the value of the cooling temperature and independent of the size of the tissue. In the case of one spiral wave, a temperature of 30˚C was sufficient to suppress the wave; however, for a more complex case, a temperature of 28˚C was required to terminate the waves. The heat induced by the Joule effect associated with the action potential is in the order of magnitude of μ˚C. Although this change is small, it is of similar magnitude as that previously found in nerves and intestinal tissues [8][9][10]. It was also found, that the heat due to the Joule effect accumulated in these tissues at the spiral tip; therefore, it could be potentially used as an alternative method for detecting the location of the spiral waves tips. Our simulations showed similar results for the cardiac tissue. Fig 8 shows the temperature field, T − T � . As can be seen, that in the case of no heat (T � = 30˚C), as the spiral wave evolves, an accumulation of heat develops near the tip of the spiral wave (shown in the first column in Fig 8). However, this behavior is not present in the case of hypothermia, (T � = 28˚C), where the temperature in the tissue is nearly constant in the entire computational domain (shown in the second column in Fig 8). The accumulation of heat was not observed in the case of hypothermia, which provides a potential explanation for the inhibition of spiral waves formation and for the promotion of termination of these waves by hypothermia.
All numerical results were obtained using the adaptive method described in section 3.2. The obtained meshes in the case of regional cooling and hypothermia are shown in Fig 9. As can be seen, using the adaptive mesh, accurate numerical results can be obtained. In the case of regional cooling (first column), the mesh is well adapted not only around the de-and repolarization fronts but also around the circle where T is set as 30˚C. Similarly, in the case of hypothermia (second column), the adapted mesh provides noticeably better results, as only few elements are obtained after the termination of the spirals waves.
As mentioned in the introduction, certain results in the literature include heat in the monodomain model via the ionic term only. Finally, we performed a comparison between the approach in the literature and the approach developed in our work. We simply replaced T by T � in Eqs (9) and (10). Then, the Pennes' bioheat equation does not affect the bidomain model and the temperature is included via the ionic term only. We repeated the simulation for the case of hypothermia for one spiral wave, as shown in Fig 4, and we compared it to the case where T was replaced by T � in the ionic model. The simulations are illustrated in Fig 10. As can be seen, there are considerable differences between the two methods as the approach in the literature failed to terminate the spiral wave. Although our approach can be relatively more expensive computationally, it can still be considered more logical, as it considers the heat over  the entire cardiac tissue rather than just the ionic terms. However, experimental results are needed to further validate the numerical results.

Conclusions and discussion
A mathematical model obtained by coupling Pennes' bioheat equation with the bidomain model is presented to simulate the effects of temperature variations on the cardiac APD. The simplification of the proposed heat-bidomain model to the heat-monodomain model is also presented. Several numerical experiments were performed. We demonstrated the effect of temperature on the APD and the conduction velocity, and the obtained numerical results were compared to the experimental results in the literature. We also investigated the effects of hypothermia and regional cooling in suppressing cardiac arrhythmias. We clearly demonstrated that hypothermia of short duration could terminate the spiral wave breakup. To our knowledge, this is the first study that presents numerical modeling of the effect of hypothermia in cardiac arrhythmias. As mentioned in the introduction, the effect of temperature has been included in studies in the literature via the ionic terms. However, in our approach the effect of spatial temperature on the cardiac electrical activity has been introduced directly. The comparison between the two approaches presented in this paper demonstrates the advantages of the proposed methodology.
The heat generated by the Joule effect was shown to accumulate at the tips of the spiral waves under non-hypothermic conditions. Although this induced variation in temperature is very small (* μ˚C), such small temperature changes have previously been detected using thin-film synthetic pyroelectric material, polyvinylidene fluoride (PVDF), when the heat produced by electrical propagation in myelinated nerve fibers was measured [26]. Therefore, there is a potential to apply this heat-accumulation property as an alternative method for detecting the location of spiral tips in cardiac tissues.
The coupling between the heat and the bidomain model depends on the gradient of the transmembrane potential, which results in challenges in the simulation. Therefore, an appropriate numerical method is necessary for simulating the developed heat-bidomain model. In our simulations we used a time adaptive mesh algorithm, which is suitable for this coupling. However, any other appropriate numerical method can be used for simulating the developed heat-bidomain model, including parallel computing on fixed meshes, high-order methods, etc.
The overall methodology presented in this paper can be extended to further consider important temperature effects on cardiac tissues, such as the heat generated by ionic pumps and mechanical contraction. In this paper, a Mitchell-Schaeffer ionic model is used, which provides a fast upstroke and reproduces similar behavior to the ventricular action potential. However, more complex cardiac ionic models can be included to the proposed heat-bidomain model to study the heat generated by ionic pumps that maintain ionic gradients crucial for the action potential. Preliminary results were obtained, where the heat-monodomain model was coupled to the Luo-Rudy model [20]. However, this type of coupled model requires different numerical algorithms than the one used in our study, which are in the focus of a future work. Heat generation by mechanical contraction is also important. However, it requires more complex coupling where the heat-bidomain model needs to be coupled with a large-deformation mechanical model. This is the subject of a future work, where our aim is to study the effects of temperature and contraction on spiral wave breakup as well.