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

TL;DR

This study investigates how nonlinear diffusion and thermo-electric coupling influence cardiac action potential dynamics, spiral wave breakup, and arrhythmogenesis through numerical simulations of a generalized two-variable model.

Contribution

It introduces a nonlinear Fickian diffusion formulation and analyzes their combined effects on cardiac wave behavior and fibrillation susceptibility.

Findings

Temperature and nonlinear diffusion affect repolarization patterns.

Thermo-electric coupling increases arrhythmogenesis risk.

Spiral wave breakup is influenced by thermal and diffusion nonlinearities.

Abstract

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.

Tables

Table 1. Table 1: Model parameters.

$V_h$	3			$V_n$	1
$V^{\ast }$	1.5415			$Q_{10}$	1.5
${\tau }_V$	2.5	[ms]		$\beta$	0.008
${\tau }_n$	250	[ms]		$D_0$	0.85	[cm2/s]
$Re$	0.9			$D_1$	0.09	[cm2/s]
$M$	9			$D_2$	0.01	[cm2/s]

Equations

$$[left=\text{empheqlbrace}]$$

$$\frac{dn}{dt}$$

$$[left=\text{empheqlbrace}]f(V,n)$$

$$g(V,n)$$

$$h(V)$$

$$\mathcal{R}(n)$$

$$[left=\text{empheqlbrace}]D(V)$$

$${\tau }_V(T)$$

Full text

††preprint: AIP/123-QED

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 spatio-temporal 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 Winfree [1987], Karma [1991, 1992, 1993, 1994], Karma et al. [1994], Fenton and Karma [1998], Krinsky and Pumir [1998], Fenton et al. [2002, 2005]. Cardiac modeling, in fact, experienced a fast-paced growth of sophisticated mathematical models Fenton and Cherry [2008] that can accurately reproduce several electrophysiological features observed in isolated cardiomyocytes, small tissues as well as whole organs both in healthy and pathological conditions Fenton and Karma [1998], Fenton et al. [2005], Luther et al. [2011], Sáez and Khul [2016]. However, several open questions still need to be answered Clayton et al. [2011], 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 Pastore et al. [1999], Pastore and Rosenbaum [2000], Gizzi et al. [2013], in particular, call for generalized theoretical approaches and mathematical modeling tools Chen et al. [In Press] 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 Bendahmane and Karlsen [2006], Potse et al. [2006], Keener and Sneyd [2009], Pullan et al. [20015], Bendahmane et al. [2010]. 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 Niederer and et al. [2011]. On these bases, spatio-temporal irregularities in excitable media have been studied from several points of view Vinet et al. [1990a], Biktashev et al. [1999] and multiple approaches have been proposed to understand and control spiral waves dynamics in heterogeneous active media Takagi et al. [2004], Pumir et al. [2005, 2010], Biktasheva [2000], Dupraz and Jacquemet [2014], Quail et al. [2014], Biktasheva et al. [2015].

Cardiac tissue, as active biological medium, is inherently characterized by multiphysics couplings and a multiscale structure Qu et al. [2014]. Cardiomyocytes, in fact, communicate via voltage-sensitive gap junctions proteins Dehin [2006] affecting the microscopic propagation timing Rohr [2004] and inducing serious cardiac diseases when genetically modified Severs et al. [2008]. In addition, local heterogeneities and anisotropies play a major role in action potential spread and propagation Potse et al. [2007]. 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 Cherry and Fenton [2007]. For these reasons, we focus here on a well established thermo-electric coupling known for greatly affecting spatio-temporal dynamics of excitable biological tissues Gizzi et al. [2010], Fenton et al. [2013], Filippi et al. [2014] and for being of primary importance in cardiology Bigelow et al. [1950], Reuler [1978], Group [1997]. Other approaches have been recently proposed in order to incorporate local spatial features within a computational framework, e.g. fractional diffusion operators Bueno-Orovio et al. [2014b, a], Cusimano et al. [2015], Bueno-Orovio et al. [2016] or statistical mechanics of cell-cell interaction Deshpande et al. [2006, 2008].

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 Karma [1993] 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 Fenton et al. [2013], Filippi et al. [2014]. The object of the numerical analysis is to investigate the emerging spatio-temporal features of linear versus nonlinear diffusion within a thermo-electric 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 Hurtado et al. [2016]. 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 (${\langle V(t)\rangle}$). One-dimensional (1D) cable, ring and two-dimensional (2D) computational domains are considered in varying domain sizes, boundary conditions and local distribution of the temperature parameter $T$.

II Methods

II.1 Generalized RD Model

We present the generalized formulation of a two-variable RD cardiac model

[TABLE]

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 $\nabla$ and $\nabla\cdot$ denote the gradient and divergence operators, respectively. According to Karma Karma [1993, 1994], 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:

[TABLE]

where $\theta(x)$ denotes the standard Heaviside step function (i.e. $\theta(x)=0$ for $x\leq 0$ and $\theta(x)=1$ for $x>0$), the nondimensional membrane voltage $V$ varies in $[0\div 4]$ and the rest state of the system corresponds to $(V,n)=(0,0)$. The restitution function, ${\mathcal{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, ${\mathcal{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 (1bd), 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, $\tau_{V}(T),\tau_{n}(T)$, via the following constitutive relations:

[TABLE]

Such an approach generalizes the porous medium formulation Hurtado et al. [2016] both in terms of RD properties and thermo-electric effects Fenton et al. [2013] and extends the restitution and dispersion properties of the original model.

In detail, the functional form (1ca) 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 Murray [2002], Cherubini et al. [2012], 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 Rohr [2004]. 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 (1cb). In particular, an Arrhenius exponential law, $Q_{10}$, is introduced into the dynamics of the gating variable $\tau_{n}(T)$ (1cb)2 whereas a linear deviation (referred to as Moore term) is considered for the time constant $\tau_{V}(T)$ of the membrane voltage (1cb)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 Fenton et al. [2013]. 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.

II.2 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. Fenton et al. [2013] and relative to endocardial optical mapping recordings Gizzi et al. [2013]. 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^{\circ}$ C, in order to reproduce the original maximum conduction velocity of $30~{}cm/s$ computed from the linear model with constant diffusion equal to $\tilde{D}=1.1~{}cm^{2}/s$. We remark that for $T=37^{\circ}$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}\gg 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}=\tilde{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. 1.

II.3 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_{\rm 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\div-2\div 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${}_{\rm S1}=450~{}ms$ followed by a single stimulation spaced with a shorter period CL${}_{\rm S2}<$CLS1. The sequence of CLS1 was: $450\div-50\div 350~{}ms$, $350\div-10\div 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$.

II.4 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, $\langle V(t)\rangle$, 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.

Figure 1: Finite element numerical simulations of a 1D cable model of length $L=25~{}cm$ stimulated on the left boundary via an external squared current $I_{\rm ext}$ ($\Delta t=2~{}ms$, amplitude 2). Space and time snapshots are taken when the wave is fully developed and travels at constant velocity (mesh size $dx=0.03~{}cm$ and maximum time step $dt=0.01~{}ms$ ensure stable conduction velocity). Comparison of the AP wave morphology in time (a) and space (b) at selected temperatures. Linear (dashed) and nonlinear (solid) diffusion model are compared. Wave overlapping is obtained on the basis of the wave front threshold AP$=1$ (axes rescale is provided for comparative purposes). (c) Zoomed view of the back-end space morphology of the excitation wave (LD-dashed, NLD-solid) for a fine gradient of temperatures (step of $0.5^{\circ}$C). (d) Wave length calculation based on panel (c) for different threshold (Thr) levels. (e) Line integral of the normalized AP spatial distribution based on panel (c). Blue shadowed region highlight the range of temperatures in which the highest differences are observed between LD and NLD model.

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^{\circ}$ 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.

III.1 1D Cable Analysis

Figures II.4 and II.4 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 Karma [1994]; 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. II.4 shows a zoomed view of the back-end of the AP wave in space for a finer gradient of the temperatures: $T\in[31\div 0.5\div 40]^{\circ}$C. Linear and nonlinear diffusion are provided (dashed, solid) and two temperatures $33,34^{\circ}$C are highlighted. In particular, we observe that the linear model does not differentiate $33,33.5^{\circ}$C and $34,34.5^{\circ}$C thus resulting in the same AP wave length for any value of the intersecting threshold chosen (see Fig. II.4). In addition, this information is confirmed by the line integral of the normalized AP wave in space as shown in Fig. II.4. 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 and 2 provide the normalized APD and CV restitution curves Karma [1994] 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. II.4, LD and NLD do not induce significant differences in the APD restitution curves that usually fit experimental data Vinet et al. [1990b]. 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 Hurtado et al. [2016], Cusimano et al. [2015]. 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. II.4. 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 ${\langle V(t)\rangle}$ on the simulated domain during periodic pacing. Figure 2 highlights the last $30~{}s$ of the restitution protocol in terms of ${\langle V(t)\rangle}$ (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, ${\langle V(t)\rangle}\simeq 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 Gizzi et al. [2013] and more sophisticated thermo-electric models Fenton et al. [2013], these fast pacing regimes are in fact characterized by a variety of alternating behaviors and bifurcation properties Watanabe et al. [2001] we do not capture in the present two-dimensional formulation.

III.2 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^{\circ}$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.

III.3 1D Ring Analysis

As usual in the context of cardiac alternans and arrhythmias Winfree [1987], 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 ${\langle V(t)\rangle}$ 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^{\circ}$C. However, the dominant frequencies of the system in the range $0\div 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 Karma [1994].

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^{\circ}$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 Pennes [1948], Gizzi et al. [2010], but assume the thermal transition as instantaneous.

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 ${\langle V(t)\rangle}$ signal, anticipates the arrhythmogenic propensity of the nonlinear diffusion model that we will highlight via two-dimensional simulations in the next section.

III.4 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 ${\langle V(t)\rangle}$ spectrum, spiral core numerosity and spatial distribution of the excitation waves. The LD and NLD models exhibit similar trends of ${\langle V(t)\rangle}$ (not shown) but the resulting spectra are characterized by higher frequencies in the NLD case. In particular, both low ($\sim 7\div 9~{}Hz$) and high ($\sim 13\div 16~{}Hz$) peaks are consistent with similar analyses Karma [1994], Filippi et al. [2014] 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 Bueno-Orovio et al. [2008], Fenton and Cherry [2008].

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 Karma [1994], Bueno-Orovio et al. [2014b, a], Cusimano et al. [2015], Bueno-Orovio et al. [2016], 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 Fenton and Karma [1998], Fenton et al. [2002] 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.

Acknowledgements.

The authors thank Flavio H. Fenton (Georgia Institute of Technology, USA) and Elizabeth M. Cherry (Rochester Institute of Technology, USA) for helpful discussions and comments. This work was supported by funding from the Italian National Group of Mathematical Physics (GNFM-INdAM), the International Center for Relativistic Astrophysics Network (ICRANet) and the London Mathematical Society through its Grant Scheme 4 - Research in Pairs.