Nonlinear diffusion&thermo-electric coupling in a two-variable model of cardiac action potential

This work reports the results of the theoretical investigation of nonlinear dynamics and spiral wave breakup in a generalized two-variable model of cardiac action potential accounting for thermo-electric coupling and diffusion nonlinearities. As customary in excitable media, the common Q10 and Moore factors are used to describe thermo-electric feedback in a 10-degrees range. Motivated by the porous nature of the cardiac tissue, in this study we also propose a nonlinear Fickian flux formulated by Taylor expanding the voltage dependent diffusion coefficient up to quadratic terms. A fine tuning of the diffusive parameters is performed a priori to match the conduction velocity of the equivalent cable model. The resulting combined effects are then studied by numerically simulating different stimulation protocols on a one-dimensional cable. Model features are compared in terms of action potential morphology, restitution curves, frequency spectra and spatio-temporal phase differences. Two-dimensional long-run simulations are finally performed to characterize spiral breakup during sustained fibrillation at different thermal states. Temperature and nonlinear diffusion effects are found to impact the repolarization phase of the action potential wave with non-monotone patterns and to increase the propensity of arrhythmogenesis.

This work reports the results of the theoretical investigation of nonlinear dynamics and spiral wave breakup in a generalized two-variable model of cardiac action potential accounting for thermo-electric coupling and diffusion nonlinearities. As customary in excitable media, the common Q 10 and Moore factors are used to describe thermo-electric feedback in a 10-degrees range. Motivated by the porous nature of the cardiac tissue, in this study we also propose a nonlinear Fickian flux formulated by Taylor expanding the voltage dependent diffusion coefficient up to quadratic terms. A fine tuning of the diffusive parameters is performed a priori to match the conduction velocity of the equivalent cable model. The resulting combined effects are then studied by numerically simulating different stimulation protocols on a onedimensional cable. Model features are compared in terms of action potential morphology, restitution curves, frequency spectra and spatio-temporal phase differences. Two-dimensional long-run simulations are finally performed to characterize spiral breakup during sustained fibrillation at different thermal states. Temperature and nonlinear diffusion effects are found to impact the repolarization phase of the action potential wave with non-monotone patterns and to increase the propensity of arrhythmogenesis. Systemic temperature is kept approximately constant through many delicate physiological feedbacks. In the heart, temperature changes greatly affect the features of the excitation wave and occur because of pathological conditions, during surgery or because of various unfortunate events. It is therefore practically relevant to determine how thermo-electric feedbacks can alter normal cardiac pacing or sustain fibrillating scenarios. Recent studies have demonstrated that structural heterogeneity of cardiac tissue is a key ingredient for understanding dispersion of repolarization and tendency to arrhythmogenesis. Here we investigate whether a generalized two-variable reaction-diffusion model of cardiac electrophysiology can highlight the dual nonlinear effects induced by thermo-electricity and diffusive nonlinearities under several simulated tests. Our aim is to provide evidence of non negligible differences in the spatiotemporal behavior of the system when these contributions are considered and to stimulate the investigation of more reliable physiological cardiac models.

I. INTRODUCTION
Understanding how transition from normal rhythm to fibrillation occurs within the heart has been for many years a primary focus of research in cardiac electrophysiology, theoretical physics, mathematical modeling as well as clinical practice [21][22][23][30][31][32][33][34]36,58 . Cardiac modeling, in fact, experienced a fast-paced growth of sophisticated mathematical models 20 that can accurately reproduce several electrophysiological features observed in isolated cardiomyocytes, small tissues as well as whole organs both in healthy and pathological conditions 21,23,37,52 . However, several open questions still need to be answered 14 , e.g.: what type of electrical activity induce spiral onset and provoke spiral breakup? Experimental evidences of the complex spatio-temporal alternans dynamics arising during fast pacing and their relation with tissue heterogeneity 27,40,41 , in particular, call for generalized theoretical approaches and mathematical modeling tools 11 able to provide a mechanistic description of these phenomena.
Virtually all mathematical models of excitable biological media derive from the seminal work of Hodgkin and Huxley on the squid giant axon and are typically formulated in terms of local kinetics and diffusive spatial features, either in their monodomain or bidomain version 1,2, 35,43,45 . The resulting set of equations constitute a nonlinear reaction-diffusion (RD) system in which spatial propagation of electrical potentials are coupled with evolution equations describing the ion-channel, exchanger and pump gate dynamics 39 . On these bases, spatio-temporal irregularities in excitable media have been studied from several points of view 4,55 and multiple approaches have been proposed to understand and control spiral waves dynamics in heterogeneous active media 5,6, 19,46,47,49,54 .
Cardiac tissue, as active biological medium, is inherently characterized by multiphysics couplings and a multiscale structure 48 . Cardiomyocytes, in fact, communicate via voltagesensitive gap junctions proteins 16 affecting the microscopic propagation timing 51 and inducing serious cardiac diseases when genetically modified 53 . In addition, local heterogeneities and anisotropies play a major role in action potential spread and propagation 44 . Temperature, in particular (as primary theme of the present work), affects the time constant of the local kinetic reactions thus inducing an enhanced level of heterogeneity to the tissue 12 . For these reasons, we focus here on a well established thermo-electric coupling known for greatly affecting spatio-temporal dynamics of excitable biological tissues [24][25][26] and for being of primary importance in cardiology 3,28,50 . Other approaches have been recently proposed in order to incorporate local spatial features within a computational framework, e.g. fractional diffusion operators 8-10,15 or statistical mechanics of cell-cell interaction 17,18 .
In healthy cardiac tissue dispersion usually has a weak role and the pulse speed is insensitive to the post-repolarization state due to a previous pulse. However, this condition does not hold when the tissue is subject to fast pacing and the diastolic interval (the time laps between the end of pulse and the beginning of the next pulse) is short.
As theoretical study, in this work we deliberately assume a minimal formulation comprehensive of generalized multiphysics couplings. In detail, we make use of the two-variable phenomenological Karma model 32 thus incorporating a Taylor expanded diffusion function up to quadratic terms and introducing the classical Q 10 and Moore factors describing the thermo-electric coupling 24,25 . The object of the numerical analysis is to investigate the emerging spatio-temporal features of linear versus nonlinear diffusion within a thermoelectric framework. We will demonstrate that the combined nonlinear effects modify the spatio-temporal behaviors of the system without impacting its local dynamics.
From the numerical point of view, we provide a fine tuning of the diffusive parameters with respect to a fixed conduction velocity and we show the ability of our formulation to avoid numerical oscillations as well as artificial discontinuities present in the classical porous medium approach 29 . We analyze a number of electrophysiological features associated with the action potential (AP) wave thus comparing linear versus nonlinear diffusion (LD, NLD) models at different thermal states: i) action potential duration (APD) and ii) conduction velocity (CV) restitution curves, iii) time and space phase differences (PD) and iv) dynamic averages of the membrane voltage ( V (t) ). One-dimensional (1D) cable, ring and twodimensional (2D) computational domains are considered in varying domain sizes, boundary conditions and local distribution of the temperature parameter T .

A. Generalized RD Model
We present the generalized formulation of a two-variable RD cardiac model where the variable V is a dimensionless representation of the transmembrane voltage and the variable n plays the role of a slow recovery current. Here ∇ and ∇· denote the gradient and divergence operators, respectively. According to Karma 32,33 , the nonlinear reaction functions f (V, n), g(V, n) identify a generalized FitzHugh-Nagumo model qualitatively that reproduce generic restitution and dispersion properties of cardiac tissue: where θ(x) denotes the standard Heaviside step function (i.e. θ(x) = 0 for x ≤ 0 and θ(x) = 1 for x > 0), the nondimensional membrane voltage V varies in [0 ÷ 4] and the rest state of the system corresponds to (V, n) = (0, 0). The restitution function, R(n), controls the relationship between two time intervals: the time-frame between the end of an action potential pulse and the beginning of the next one (diastolic interval, DI) and the duration of the next action potential pulse (APD). The dispersion function, D(n), prescribes the relation between the instantaneous speed of the AP front-end at a given spacial point and the time elapsed since the back-end of a previous pulse passed the same location. Provided the specific functions (2d), the restitution is controlled by the parameter Re while the dispersion is controlled by the parameter M . We account for a self diffusion quadratic law, D(V ), and thermal coupling, τ V (T ), τ n (T ), via the following constitutive relations: Such an approach generalizes the porous medium formulation 29 both in terms of RD properties and thermo-electric effects 24 and extends the restitution and dispersion properties of the original model. In detail, the functional form (3a) accounts for a basic constant diffusion D 0 enriched by linear, D 1 , and quadratic, D 2 , self diffusion terms that contribute to the diffusive flux according to the voltage field level V . When the voltage is low (close to the resting state) the additional contributions are negligible and the equation boils down to the standard cable equation with constant diffusion D 0 . However, when the voltage is at the plateau phase (depolarized state) then these two additional contributions affect the diffusive flux and consequentially speed and wavelength of the propagating wave. This particular feedback, well known in the context of animal dispersal and population growth 13,38 , is here introduced to mimic, in a phenomenological way, the nonlinear dependence of gap junctions conductance on the transjunctional voltage occurring at the sub-micron scale 51 . It is worth noticing here that we do not aim at reproducing the exact physical nature of gap junction couplings at the macroscopic level (whose description would require a more complex multiscale homogenization procedure), but we investigate the properties of our generalized formulation by analyzing its nonlinear behaviors via the classical tools used for dynamical systems.
On a similar phenomenological ground, well established is the nonlinear coupling adopted for the thermo-electric model (3b). In particular, an Arrhenius exponential law, Q 10 , is introduced into the dynamics of the gating variable τ n (T ) (3b) 2 whereas a linear deviation (referred to as Moore term) is considered for the time constant τ V (T ) of the membrane voltage (3b) 1 . These additional thermal factors contribute to adapt the restitution and dispersion time scales according to the selected temperatures. They greatly affect the duration of the AP wave as well as its conduction velocity 24 . Therefore, the concomitant effect of nonlinear diffusion and thermo-electric coupling will result in enhanced properties for the dynamical systems that we aim at exploring in the present contribution.

B. Parameter's Space Analysis
In order to avoid a wide and unstructured range of free parameters, we set the thermal coefficients according to the experimental fitting reported in Fenton et al. 24 and relative to endocardial optical mapping recordings 27 . Accordingly, we study the features of our thermo-electric dynamical system for a fixed set of model parameters such that extended comparisons are performed between linear and nonlinear diffusion operators only at different thermal states. To this aim, a preliminary fine tuning of the nonlinear diffusion parameters D 0 , D 1 , D 2 was conducted at physiological temperature, i.e. T = 37 • C, in order to reproduce the original maximum conduction velocity of 30 cm/s computed from the linear model with constant diffusion equal toD = 1.1 cm 2 /s. We remark that for T = 37 • C the thermo-electric coupling has no effects (both Q 10 and Moore terms are equal to 1).
As preliminary tuning, several combinations of the three parameters were tested by adopting as ultimate limiting condition the order of magnitude of the three coefficients, i.e. D 0 D 1 > D 2 . We remark that, for the nondimensional phenomenological model at hand, the physical dimensions of these coefficients are homogeneous, i.e. [cm 2 /s]. In the calculation we started using the value D 0 =D from the original model, then we introduced one coefficient at a time (maintaining the other at 0) and lowering the value of D 0 . Finally, we analyzed the three parameters together and we chose the optimum combination respecting the expected orders of magnitude. Such a fitting strategy, though not unique, is plausible since we aimed at introducing small deviations with respect to the linear model (cable equation) in a minimal simplified setting. The complete list of model parameters is provided in Tab. I.

C. Restitution Protocols
The numerical analysis of the 1D model was conducted via two alternative activation protocols, i.e.: pacing down (or full restitution) and S1S2. The two protocols were used for computing the APD and CV restitution curves and for extracting phase differences between linear and nonlinear diffusion cases.
The full restitution protocol consisted of periodic electrical activations induced through the square wave current I ext at the left boundary of the cable with a constant amplitude equal to 2, a constant duration of 2 ms but varying its duty cycle. For a fixed pacing cycle length (CL), n = 10 stimulations were delivered before CL was decreased. In particular, the protocol was applied for each selected temperature and for LD and NLD models. In order to compute the fine transition of the CV-CL restitution, the applied protocol consisted in the sequence: CL= 800 ÷ −2 ÷ 100 ms, a cable of length L = 15 cm was adopted in order to ensure the complete evolution of the propagating wave. The protocol ended when capture of activation was lost and was implemented via finite differences Fortran routines with an explicit time integration scheme of first order, a fixed dt = 0.1 ms, and fixed dx = 0.025 cm (equivalent results were obtained with finite element implementations using piecewise linear approximations for every field and a semi-implicit discretization in time).
The S1S2 protocol consisted of the series of a constant stimulation at CL S1 = 450 ms followed by a single stimulation spaced with a shorter period CL S2 <CL S1 . The sequence of CL S1 was: 450 ÷ −50 ÷ 350 ms, 350 ÷ −10 ÷ 200 ms. Also for this case, since the spatial distribution of AP is necessary for the analysis, a cable of length L = 10 cm was considered. In this case, the minimum CL ensuring correct pacing was 280 ms, higher than those reached during the full restitution. This protocol was implemented in finite elements with mesh size dx = 0.025 cm and maximum time step dt = 0.1 ms.

D. Averaged Action Potential and Gradient Operator
In the following analysis we make use of two additional quantities related to the membrane voltage V and to the recovery variable n. We considered the spatial average of the membrane voltage, V (t) , computed as function of time on both 1D and 2D domains. Accordingly, the time course of the average signal was investigated by means of the FFT transform such that its spectrum was studied. In addition, the sharp gradients of the recovery variable n were computed for 2D simulations in order to identify the numerosity of localized regions with high gradients (spiral cores) during long run simulations thus to compare LD and NLD model spatio-temporal behaviors.

III. RESULTS
We start by comparing the one-dimensional spatio-temporal features of the AP wave arising from the nonlinear thermo-electric model in cable (Neumann zero flux) and ring (periodic boundary conditions) domains. We conclude with 2D long run simulations of spiral breakup dynamics. An extended analysis is conducted comparing linear (LD) and nonlinear (NLD) diffusion at selected temperatures based on the indicators described in the previous section and enriched by the temporal and spatial analysis of AP wave phase differences. In the description of the results we will make use of a unique color code associated with the temperature parameter: T = 40, 37, 34, 31 • C will be shown in red, green, blue, black, respectively. In addition, when not specified, LD and NLD quantities will be shown as dashed and solid line, respectively.
A. 1D Cable Analysis   Figures 1(a) and 1(b) show the morphology of the simulated AP waves in time and space, respectively, comparing LD and NLD model at the four selected temperatures. As expected from the model formulation: i) the shape of the AP wave does not change neither due to temperature nor to the nonlinear diffusion operator; ii) temperature acts on both the temporal duration and the spatial distribution of the AP wave, in particular widening the AP duration at lower temperatures; iii) NLD acts on the spatial properties of the AP wave and in particular on the repolarization phase of the wave 33 ; iv) the combined effect of NLD and thermo-electric coupling entails a non-unique variation in space of the AP wave and produces higher differences at lower temperatures. In order to emphasize the last observation, Fig. 1(c) shows a zoomed view of the back-end of the AP wave in space for a finer gradient of the temperatures: T ∈ [31÷0.5÷40] • C. Linear and nonlinear diffusion are provided (dashed, solid) and two temperatures 33, 34 • C are highlighted. In particular, we observe that the linear model does not differentiate 33, 33.5 • C and 34, 34.5 • C thus resulting in the same AP wave length for any value of the intersecting threshold chosen (see Fig. 1(d)). In addition, this information is confirmed by the line integral of the normalized AP wave in space as shown in Fig. 1(e). It is worth noticing that the combined effect of nonlinear diffusion and thermo-electric coupling seems to regularize the non-monotone behavior observed for the linear model when a single propagating wave is analyzed. Richer information can be recovered when the system undergoes fast pacing or sustains rotating spirals.
Figures 2(a) and 2(b) provide the normalized APD and CV restitution curves 33 obtained via the full pacing down protocol conducted over a 1D cable (APD and CV values are taken from the last activation wave). As expected from the point-wise temporal analysis shown in Fig. 1(a), LD and NLD do not induce significant differences in the APD restitution curves that usually fit experimental data 56 . Such a quantity, in fact, mainly relies on the front- end of the excitation wave which is known to be affected by nonlinear diffusion operators primarily on the foot of the AP wave spatial distribution 15,29 . Since in our simplified case such a nonlinearity is moderate, the resulting APD restitution curves overlap.
On the contrary, small deviations are observed on the CV restitution curves thus entailing the eventual role of diffusion in the spatio-temporal features of the system in accordance with Fig. 1(b). Temperature shifts the APD and CV curves: longer APDs and smaller CVs are in fact obtained at lower temperatures. A reduced minimum pacing CL and a faster decay of the CV curve are also observed at lower temperatures. It is important to note here that these small deviations may result in greatly different behavior for more complex models accounting for heterogeneities and anisotropies or when long run simulations are conducted.
The analysis of the restitution curves is here enriched with the dynamical calculation of the average action potential V (t) on the simulated domain during periodic pacing. Figure 2(c) highlights the last 30 s of the restitution protocol in terms of V (t) (the sole nonlinear case is shown since similar to the linear case). It is worth noting that a stable oscillating signal is present before a reduction of the oscillations amplitude is reached, V (t) 1.8, for all the selected temperature. This value corresponds to the physical condition for which the new excitation wave entering the domain from one boundary is balanced by the previous one leaving the domain from the opposite boundary. However, according to the restitution maps described before, the transition towards such a regime occurs at earlier times (e.g. at larger pacing CL) for lower temperatures. In addition, the last simulated times (shortest CLs) are characterized by an eventual increase of the oscillations just before loss of activation. At comparison with experimental evidences 27 and more sophisticated thermo-electric models 24 , these fast pacing regimes are in fact characterized by a variety of alternating behaviors and bifurcation properties 57 we do not capture in the present two-dimensional formulation.

B. Phase Analysis
We consider now the S1S2 restitution protocol and analyze temporal and spatial phase differences (PD) of the AP wave. In detail, we synchronize the time course and spatial distribution of the AP signal on the S1 activation wave and calculate the algebraic difference for the sole S2 activation thus resulting in the plots shown in the top panel of Fig. 3. We consider then the sole maxima and minima values resulting from these differences and plot them in the phase space (CL,T) by running the analysis over the four selected temperatures and pacing periods adopted during the stimulation protocol. The resulting diagrams indicate regions (basins) with high variations between LD and NLD models. In particular, the temporal PD only indicates maximum differences at high pacing CLs and normal temperature, while spatial PD shows a strong correlation with the major differences obtained at 37, 31 • C for short CLs. This second observation is in agreement with the restitution analysis discussed before thus implying a wide spectrum of novel properties to arise from the combination of nonlinear diffusion and nonlinear thermo-electricity.  T=40°C -LD  T=40°C -NLD  T=37°C -LD  T=37°C -NLD  T=34°C -LD  T=34°C -NLD  T=31°C -LD  T=31°C -

C. 1D Ring Analysis
As usual in the context of cardiac alternans and arrhythmias 58 , in this section we investigate the dynamical properties of our thermo-electric system by analyzing a 1D cable with periodic boundary conditions (ring). In view of multi-dimensional analyses, we also consider different static and dynamic spatial distributions of the temperature. Figure 4 provides the time course of V (t) and the corresponding frequency spectrum for steady state uniform and distributed temperatures (a,b). The analysis shows that the rotating wave stabilizes on different average values according to the temperature parameter T . Persisting oscillations are present for the NLD case (solid) with respect to the LD one (dashed) and in particular at 40 • C. However, the dominant frequencies of the system in the range 0 ÷ 10 Hz are captured in a similar manner and curves overlap (right panel). However, if a spatial distribution of temperatures is considered in the ring (warmer and colder regions in panel b), then a clear shift between LD and NLD cases arises as well as the appearance of novel intermedia peaks in the corresponding frequency spectra. Of note an intermedia peak at about 5 Hz is present underlying an increased heterogeneity of the tissue affecting the characteristic time scales of the system during arrhythmias 33 .
An extended analysis of rotating waves in the 1D ring setting is further provided in Fig. 5. We compare different spatial distributions of the parameter T with the additional investigation of thermal dynamical variation. We impose a smooth dynamical transition from the reference physiological temperature, 37 • C, in 3 s; we keep the new thermal state for 3 s (light blue region) and we finally recover the reference temperature in 3 s (gray region). Though longer time profiles may be investigated, the selected 3 s intervals allow the system to stabilize in the new rotating spiral configuration. It is worth noting here that we do not consider the effect of thermal diffusion and tissue perfusion which would require the additional coupling of a Pennes bio-equation 26,42 , but assume the thermal transition as instantaneous.  Comparison of V (t) (center) and corresponding spectra (right) for linear (LD) and nonlinear (NLD) diffusion models on a ring (L = 6.55 cm) for different spatial thermal profiles (left) (p1-p5 stand for the selected thermal profile). A smooth dynamical transition from the reference physiological temperature, 37 • C, is applied in 3 s, kept for 3 s (light blue region) and finally removed in 3 s (gray region) such to recover the reference temperature. The plot reports the last 6 s of V (t) in order to appreciate signal differences.
The results confirm the shift between LD and NLD models and indicate that the LD case is able to recover the original average value in shorter time than the NLD one. Such a behavior, with additional oscillations observed for the NLD V (t) signal, anticipates the arrhythmogenic propensity of the nonlinear diffusion model that we will highlight via two-dimensional simulations in the next section.

D. Spiral Breakup
We conclude the numerical investigation of the nonlinear dynamical system with long run simulations (60 s of physical time) of sustained 2D spiral breakup at different temperatures. The results are shown in Fig. 6 in terms of V (t) spectrum, spiral core numerosity and spatial distribution of the excitation waves. The LD and NLD models exhibit similar trends of V (t) (not shown) but the resulting spectra are characterized by higher frequencies in the NLD case. In particular, both low (∼ 7÷9 Hz) and high (∼ 13÷16 Hz) peaks are consistent with similar analyses 25,33 though a shifted dominant frequency of about 1 Hz is observed between linear and nonlinear case. This result generalizes the analyses conducted in the previous sections on the ring and is corroborated by the core numerosity plots highlighting a higher propensity of the system to produce spiral breakup at higher temperatures since the wave lengths are shorter and the repolarization times are faster. Accordingly, the spatial distribution of spiral waves shows different fibrillating scenarios between LD and NLD cases.
It is important to observe here that the sole computed restitution curves are not able to predict these differences in the spatio-temporal dynamics. In fact, via the thermo-electric coupling we introduced an additional nonlinear effect in the system's temporal properties whereas a supplemental nonlinearity is enforced in the spatial flux via a dispersal model. The resulting dynamical system is therefore characterized by more complex emerging behaviors and additional nonlinearities and future studies are demanded in this direction on more physiological models.

IV. CONCLUSIONS
In summary, we have introduced a generalized two-variable model of action potential propagation accounting for thermo-electric coupling and nonlinear diffusion which reproduces characteristic restitution and dispersion properties of cardiac tissues. The simplified structure adopted allows to limit the parameter's space to the sole new elements ruling thermo-electricity within the framework of non-Fickian fluxes. However, several missing physiological factors are not captured and extended computational studies are envisioned with more physiological models 7,20 .
On this ground, we have investigated the dynamical origin of spiral wave breakup and compared long run simulations of linear and nonlinear diffusion under different thermal regimes. The main conclusion is that the onset of spiral wave breakup is enhanced at the highest and lowest temperatures tested by the nonlinear model. In addition, the resulting fibrillating scenario is affected by the non-Fickian flux in a non predictive way. By analogy with similar theoretical works 8-10,15,33 , the present contribution suggest a direct connection between the emerging behavior at the tissue level with the dynamical response of single cells.
The second main conclusion of the work is that the linear and nonlinear dynamical systems adjust with slightly different behaviors when the temperature assumes heterogeneous distributions or is modulated in time. It is worth noting that spatial heterogeneities and anisotropies may potentially play an important role in the present setting 21,22 and we expect to address this particular aspect in forthcoming contributions.
Of particular importance for understanding and control of cardiac arrhythmias, the study of cardiac alternans and their multiple transitions towards fibrillation would represent an additional important element to investigate within the context of thermo-electricity and nonlinear diffusion. In this perspective, three-or four-variable phenomenological models are the best candidates to explore the key features of such a complex dynamical system.