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° 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.

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.

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, and mathematical modeling as well as clinical practice.21–23,31–35,37,59 Cardiac modeling, in fact, experienced a fast-paced growth of sophisticated mathematical models20 that can accurately reproduce several electrophysiological features observed in isolated cardiomyocytes and small tissues as well as whole organs both in healthy and pathological conditions.21,23,38,53 However, several open questions still need to be answered,14 e.g.,: what type of electrical activity induces spiral onset and provokes spiral breakup? Experimental evidence of the complex spatio-temporal alternans dynamics arising during fast pacing and their relationship with tissue heterogeneity,27,41,42 in particular, calls for generalized theoretical approaches and mathematical modeling tools11 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,36,44,46 The resulting set of equations constitute a nonlinear reaction-diffusion (RD) system in which spatial propagation of electrical potentials is coupled with evolution equations describing the ion-channel, exchanger, and pump gate dynamics.40 On these bases, spatio-temporal irregularities in excitable media have been studied from several points of view4,56 and multiple approaches have been proposed to understand and control spiral waves dynamics in heterogeneous active media.5,6,19,47,48,50,55

Cardiac tissue, as active biological medium, is inherently characterized by multiphysics couplings and a multiscale structure.49 Cardiomyocytes, in fact, communicate via voltage-sensitive gap junction proteins16 affecting the microscopic propagation timing52 and inducing serious cardiac diseases when genetically modified.54 In addition, local heterogeneities and anisotropies play a major role in action potential (AP) spread and propagation.45 Temperature, in particular (as a 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 tissues24–26 and for being of primary importance in cardiology.3,29,51 Other approaches have been recently proposed in order to incorporate local spatial features within a computational framework, e.g., fractional diffusion operators8–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 (DI) (the time laps between the end of pulse and the beginning of the next pulse) is short.

As a 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 model33 thus incorporating a Taylor expanded diffusion function up to quadratic terms and introducing the classical Q10 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 (LD, NLD) 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 (CV) and we show the ability of our formulation to avoid numerical oscillations as well as artificial discontinuities present in the classical porous medium approach.30 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 two-dimensional (2D) computational domains are considered in varying domain sizes, boundary conditions, and local distribution of the temperature parameter T.

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

{Vt=·D(V)V+f(V,n)τV(T)+Iext(1a)dndt=g(V,n)τn(T),(1b)
(1)

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,33,34 the nonlinear reaction functions f(V,n),g(V,n) identify a generalized FitzHugh-Nagumo model qualitatively that reproduces generic restitution and dispersion properties of cardiac tissue

{f(V,n)=V+[V*D(n)]h(V)(2a)g(V,n)=R(n)θ(VVn)[1θ(VVn)]n(2b)h(V)=[1tanh(VVh)]V22(2c)R(n)=1(1eRe)n1eRe,D(n)=nM,(2d)
(2)

where θ(x) denotes the standard Heaviside step function (i.e., θ(x)=0 for x0 and θ(x)=1 for x > 0), the nondimensional membrane voltage V varies in [04] 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), describes the relationship 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 relationships:

{D(V)=D0+D1V+D2V2(3a)τV(T)=τV01+β(TT0),τn(T)=τn0Q10(TT0)/10.(3b)
(3)

Such an approach generalizes the porous medium formulation30 both in terms of RD properties and thermo-electric effects24 and extends the restitution and dispersion properties of the original model.

In detail, the functional form (3a) accounts for a basic constant diffusion D0 enriched by linear, D1, and quadratic, D2, 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 D0. 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,39 is here introduced to mimic, in a phenomenological way, possible nonlinear functional dependencies occurring at the cellular and sub-cellular scale and affecting the excitation wave dynamics, e.g., gap junction voltage-conductance coupling.52 We refer, in particular, to recent contributions adopting a fractional diffusion approach,8,10,15 that clearly state the need for extended formulations encompassing the microstructure of the tissue. However, it is worth noting 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 nonlinear behavior of the present formulation, employing classical analysis tools for dynamical systems. In this perspective, we aim at reducing the complexity of the mathematical model while maintaining key novel aspects, which lead to consider a Taylor expansion of the diffusivity upon the membrane voltage.

On a similar phenomenological basis, well established is the nonlinear coupling adopted for the thermo-electric model (3b). In particular, an Arrhenius exponential law, Q10, is introduced into the dynamics of the gating variable τn(T) (3b)2 (inducing a change of 0.981.05 in the 10-degrees range of temperatures), whereas a linear deviation (referred to as Moore term) is considered for the time constant τV(T) of the membrane voltage (3 b)1 (leading to a marked change of 0.881.28 in the same range of temperatures). 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.

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 D0,D1,D2 was conducted at physiological temperature, i.e., T=37°C, in order to reproduce the original maximum conduction velocity of 30cm/s computed from the linear model with constant diffusion equal to D̃=1.1cm2/s. We remark that for T=37°C the thermo-electric coupling has no effects (both Q10 and Moore terms are equal to 1).

As preliminary tuning, several combinations of the three parameters were tested by adopting as the ultimate limiting condition the order of magnitude of the three coefficients, i.e., D0D1>D2. We remark that, for the nondimensional phenomenological model at hand, the physical dimensions of these coefficients are homogeneous, i.e., [cm2/s]. In the calculation we started using the value D0=D̃ from the original model, then we introduced one coefficient at a time (maintaining the other at 0) and lowering the value of D0. Finally, we analyzed the three parameters together and we chose the optimum combination respecting the expected orders of magnitude. Such a fitting strategy, although 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 Table I.

TABLE I.

Model parameters.

Vh   Vn  
V *   1.5415  Q10   1.5 
τ V 0   2.5 ms  β   0.008 
τ n 0   250 ms  D0   0.85 cm2/s 
Re   0.9  D1   0.09 cm2/s 
M   D2   0.01 cm2/s 
Vh   Vn  
V *   1.5415  Q10   1.5 
τ V 0   2.5 ms  β   0.008 
τ n 0   250 ms  D0   0.85 cm2/s 
Re   0.9  D1   0.09 cm2/s 
M   D2   0.01 cm2/s 

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 Iext at the left boundary of the cable with a constant amplitude equal to 2, a constant duration of 2ms but varying its duty cycle. For a fixed pacing cycle length (CL), ten 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 of the sequence: CL =800to2to100ms, a cable of length L=15cm 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.1ms, and fixed dx=0.025cm (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 CLS1=450ms followed by a single stimulation spaced with a shorter period CLS2< CLS1. The sequence of CLS2 was: 450 to −50 to 350 ms, 350 to −10 to 200 ms. Also for this case, since the spatial distribution of AP is necessary for the analysis, a cable of length L=10cm was considered. In this case, the minimum CL ensuring correct pacing was 280ms, higher than those reached during the full restitution. This protocol was implemented in finite elements with mesh size dx=0.025cm and maximum time step dt=0.1ms.

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 a 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.

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 Secs. III AIII C 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, and black, respectively. In addition, when not specified, LD and NLD quantities will be shown as dashed and solid line, respectively.

Figures 1(a) and 1(b) show the morphology of a single propagating AP wave generated at the left end of the 1D cable in time and space, respectively, comparing the LD and NLD models 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;34 (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.

FIG. 1.

Finite element numerical simulations of a 1D cable model of length L=25cm stimulated on the left boundary via an external squared current Iext (Δt=2ms, amplitude 2). Space and time snapshots are taken when the wave is fully developed and travels at constant velocity (mesh size dx=0.03cm and maximum time step dt=0.01ms 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 models 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°C). (d) Wave length calculation based on panel (c) for different threshold levels, Thr=3.25,3,2.25,0.25 of the AP amplitude. (e) Line integral of the normalized AP spatial distribution based on panel (c). The shadowed blue region highlights the range of temperatures in which the highest differences are observed between the LD and NLD models.

FIG. 1.

Finite element numerical simulations of a 1D cable model of length L=25cm stimulated on the left boundary via an external squared current Iext (Δt=2ms, amplitude 2). Space and time snapshots are taken when the wave is fully developed and travels at constant velocity (mesh size dx=0.03cm and maximum time step dt=0.01ms 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 models 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°C). (d) Wave length calculation based on panel (c) for different threshold levels, Thr=3.25,3,2.25,0.25 of the AP amplitude. (e) Line integral of the normalized AP spatial distribution based on panel (c). The shadowed blue region highlights the range of temperatures in which the highest differences are observed between the LD and NLD models.

Close modal

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[310.540]°C. Linear and nonlinear diffusion is 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 (linear spatial measure of the AP wave at a selected threshold) for any value of the intersecting threshold chosen [see Fig. 1(d)]. Four different threshold levels of the AP wave amplitude are compared, Thr=3.25,3,2.25,0.25. In addition, this information is confirmed by looking at 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. In this regard, the occurrence of numerical artifacts is indeed excluded, since a comparison of different numerical resolutions and schemes was carried out (including finite difference and finite element discretizations). Alternatively, the effect of interest could be a consequence of the additional nonlinearities appearing in that particular range of temperatures. Furthermore, richer information can be recovered when the system undergoes fast pacing, or when it sustains rotating spirals.

Figures 2(a) and 2(b) provide the normalized APD and CV restitution curves34 obtained via the full pacing down protocol conducted over a 1D cable (APD and CV values are taken from the last activation wave) for long run simulations in which more than 10 stimulations at constant CL were delivered in order to ensure steady state conditions. According to the discussion provided for the point-wise temporal analysis in Fig. 1(a), LD and NLD do not induce significant differences in the APD restitution curves that, in fact, usually fit experimental data with good accuracy.57 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,30 Since in our simplified case such a nonlinearity is moderate, the resulting APD restitution curves overlap.

FIG. 2.

Normalized restitution curves obtained via a full pacing down stimulation protocol comparing linear and nonlinear diffusion at different temperatures. (a) and (b) APD-DI and CV-CL normalized maps with respect to the maximum value of APD (APD-Max) and of CV (CV-Max) for the four selected temperatures using >10 stimulations at constant CL in order to ensure steady state conditions. In panel (a) the gray dashed line indicates slope 1. (c) Spatial average value of the excitation wave, V(t), calculated during the full restitution protocol (with 5 stimulations at constant CL) for the nonlinear model (the linear model behaves similarly). Each panel refers to a specific temperature value indicated in the legend.

FIG. 2.

Normalized restitution curves obtained via a full pacing down stimulation protocol comparing linear and nonlinear diffusion at different temperatures. (a) and (b) APD-DI and CV-CL normalized maps with respect to the maximum value of APD (APD-Max) and of CV (CV-Max) for the four selected temperatures using >10 stimulations at constant CL in order to ensure steady state conditions. In panel (a) the gray dashed line indicates slope 1. (c) Spatial average value of the excitation wave, V(t), calculated during the full restitution protocol (with 5 stimulations at constant CL) for the nonlinear model (the linear model behaves similarly). Each panel refers to a specific temperature value indicated in the legend.

Close modal

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 30s of the restitution protocol with 5 repetitive stimulations 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 oscillation amplitude is reached, V(t)1.8, for all the selected temperatures. 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 evidence27 and more sophisticated thermo-electric models,24 these fast pacing regimes are in fact characterized by a variety of alternating behaviors and bifurcation properties58 we do not capture in the present two-variable formulation.

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 (550 ms) and calculate the algebraic difference for the sole S2 activation (from 500 to 280 ms with smaller steps at shorter CLs), thus resulting in the traces shown in the top panel of Fig. 3. We consider then the sole maxima and minima values resulting from the phase differences (top panel of Fig. 3) 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 (middle and central panels of Fig. 3 with black dots indicating each S2 stimulation). 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.

FIG. 3.

Temporal (top left) and spatial (top right) synchronization of the LD and NLD AP signals and their algebraic difference for a representative S1S2 stimulation. The maximum and minimum phase difference is highlighted in blue and red symbol, respectively. Interpolated contour plots (spline Matlab routine) of the temporal (left) and spatial (right) maxima and minima voltage phase differences computed on the S2 wave of the S1S2 restitution protocol (see top panel). The color code refers to the amplitude of the difference at the four selected temperatures. Simulated data are provided as black circles.

FIG. 3.

Temporal (top left) and spatial (top right) synchronization of the LD and NLD AP signals and their algebraic difference for a representative S1S2 stimulation. The maximum and minimum phase difference is highlighted in blue and red symbol, respectively. Interpolated contour plots (spline Matlab routine) of the temporal (left) and spatial (right) maxima and minima voltage phase differences computed on the S2 wave of the S1S2 restitution protocol (see top panel). The color code refers to the amplitude of the difference at the four selected temperatures. Simulated data are provided as black circles.

Close modal

As usual in the context of cardiac alternans and arrhythmias,59 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 010Hz 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 intermediate peak at about 5Hz is present underlying an increased heterogeneity of the tissue affecting the characteristic time scales of the system during arrhythmias.34 

FIG. 4.

Comparison of linear and nonlinear diffusion (dashed, solid line) on a ring of L=6.55cm (the length is chosen according to the same analysis performed in the original Karma model34). (a) Transition from T=37°C to hyper- and hypothermia static states (left); V(t) time course of the system (color code as in previous plots) and corresponding power spectral density (right). (b) Static spatial distribution of temperatures for three different scenarios (left). Corresponding time course of V(t) and frequency spectrum (right).

FIG. 4.

Comparison of linear and nonlinear diffusion (dashed, solid line) on a ring of L=6.55cm (the length is chosen according to the same analysis performed in the original Karma model34). (a) Transition from T=37°C to hyper- and hypothermia static states (left); V(t) time course of the system (color code as in previous plots) and corresponding power spectral density (right). (b) Static spatial distribution of temperatures for three different scenarios (left). Corresponding time course of V(t) and frequency spectrum (right).

Close modal

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 3s; we keep the new thermal state for 3s (light blue region) and we finally recover the reference temperature in 3s (gray region). Although longer time profiles may be investigated, the selected 3s 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,43 but assume the thermal transition as instantaneous.

FIG. 5.

Comparison of V(t) (center) and corresponding spectra (right) for linear (LD) and nonlinear (NLD) diffusion models on a ring (L=6.55cm according to Ref. 34) for different spatial thermal profiles (left) (p1–p5 stand for the selected thermal profile–Gaussian warmer, Gaussian colder, linear, multimodal warmer, multimodal colder). A smooth dynamical transition from the reference physiological temperature, 37°C, is applied in 3s, kept for 3s (light blue region) and finally removed in 3s (gray region) such to recover the reference temperature. The central plots report the last 6s of V(t) in order to appreciate signal differences. The spectra are referred to the whole time course of the V(t) signal.

FIG. 5.

Comparison of V(t) (center) and corresponding spectra (right) for linear (LD) and nonlinear (NLD) diffusion models on a ring (L=6.55cm according to Ref. 34) for different spatial thermal profiles (left) (p1–p5 stand for the selected thermal profile–Gaussian warmer, Gaussian colder, linear, multimodal warmer, multimodal colder). A smooth dynamical transition from the reference physiological temperature, 37°C, is applied in 3s, kept for 3s (light blue region) and finally removed in 3s (gray region) such to recover the reference temperature. The central plots report the last 6s of V(t) in order to appreciate signal differences. The spectra are referred to the whole time course of the V(t) signal.

Close modal

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 Sec. III D.

We conclude the numerical investigation of the nonlinear dynamical system with long run simulations (60s of physical time) of sustained 2D spiral breakup at different temperatures. The re-entry onset protocol was the same for the four temperatures. We initiated a plane wave and then stimulated a small region of the wave back (during the repolarization phase) such to block the propagation in the direction of the previous wave and create a spiral. In order to compensate for the changing repolarization timing, it was necessary to fine tune the timing of the second stimulus for each temperature. In particular, a delay of the second stimulus is necessary at lower temperatures since the wave length widens.

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 (79Hz) and high (1316Hz) peaks are consistent with similar analyses25,34 although a shifted dominant frequency of about 1Hz is observed between linear and nonlinear cases. This result generalizes the analyses conducted in Secs. III AIII C 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.

FIG. 6.

Long run 2D spiral breakup simulations on a domain size 6.5×6.5cm2 implemented via Fortran routines with an explicit finite differences integration scheme (dt=0.1ms,dx=0.025cm) with Neumann zero flux boundary conditions. Comparison of LD (black line) and NLD (red line) models at four selected temperatures is shown. (Top) V(t) power spectrum density (PSD) in the frequency range 520Hz; (Center) numerosity of spiral cores identified in the 2D mapped field over the simulated time (060s); (Bottom) spatial distribution of the excitation waves with different spiral breakup at the final simulated time (60s). The color code refers to yellow–blue as maximum–minimum voltage V.

FIG. 6.

Long run 2D spiral breakup simulations on a domain size 6.5×6.5cm2 implemented via Fortran routines with an explicit finite differences integration scheme (dt=0.1ms,dx=0.025cm) with Neumann zero flux boundary conditions. Comparison of LD (black line) and NLD (red line) models at four selected temperatures is shown. (Top) V(t) power spectrum density (PSD) in the frequency range 520Hz; (Center) numerosity of spiral cores identified in the 2D mapped field over the simulated time (060s); (Bottom) spatial distribution of the excitation waves with different spiral breakup at the final simulated time (60s). The color code refers to yellow–blue as maximum–minimum voltage V.

Close modal

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.

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 us 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,34 the present contribution suggests 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 setting21,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.28 In this perspective, three- or four-variable phenomenological models are the best candidates to explore the key features of such a complex dynamical system.

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.

1.
M.
Bendahmane
and
K. H.
Karlsen
, “
Analysis of a class of degenerate reaction-diffusion systems and the bidomain model of cardiac tissue
,”
Networks Heterog. Media
1
,
185
218
(
2006
).
2.
M.
Bendahmane
,
R.
Bürger
, and
R.
Ruiz-Baier
, “
A multiresolution space-time adaptive scheme for the bidomain model in electrocardiology
,”
Numer. Methods Partial Differ. Equations
26
,
1377
1404
(
2010
).
3.
W. G.
Bigelow
,
J. C.
Callaghan
, and
J. A.
Hopps
, “
General hypothermia for experimental intracardiac surgery; the use of electrophrenic respirations, an artificial pacemaker for cardiac standstill and radio-frequency rewarming in general hypothermia
,”
Ann. Surg.
132
,
531
539
(
1950
).
4.
V. N.
Biktashev
,
I. V.
Biktasheva
,
A. V.
Holden
,
M. A.
Tsyganov
,
J.
Brindley
, and
N. A.
Hill
, “
Spatiotemporal irregularity in an excitable medium with shear flow
,”
Phys. Rev. E
60
,
1897
(
1999
).
5.
I. V.
Biktasheva
, “
Drift of spiral waves in the complex Ginzburg-Landau equation due to media inhomogeneities
,”
Phys. Rev. E
62
,
8800
(
2000
).
6.
I. V.
Biktasheva
,
H.
Dierckx
, and
V. N.
Biktashev
, “
Drift of scroll waves in thin layers caused by thickness features: Asymptotic theory and numerical simulations
,”
Phys. Rev. Lett.
114
,
068302
(
2015
).
7.
A.
Bueno-Orovio
,
E. M.
Cherry
, and
F. H.
Fenton
, “
Minimal model for human ventricular action potentials in tissue
,”
J. Theor. Biol.
253
,
544
560
(
2008
).
8.
A.
Bueno-Orovio
,
D.
Kay
, and
K.
Burrage
, “
Fourier spectral methods for fractional-in-space reaction-diffusion equations
,”
BIT Numer. Math.
54
,
937
954
(
2014
).
9.
A.
Bueno-Orovio
,
D.
Kay
,
V.
Grau
,
B.
Rodriquez
, and
K.
Burrage
, “
Fractional diffusion models of cardiac electrical propagation: Role of structural heterogeneity in dispersion of repolarization
,”
J. R. Soc. Interface
11
,
20140352
(
2014
).
10.
A.
Bueno-Orovio
,
I.
Teh
,
J. E.
Schneider
,
K.
Burrage
, and
V.
Grau
, “
Anomalous diffusion in cardiac tissue as an index of myocardial microstructure
,”
IEEE Trans. Med. Imaging
35
,
2200
2207
(
2016
).
11.
D. D.
Chen
,
R. A.
Gray
,
I.
Uzelac
,
C.
Herndon
, and
F. H.
Fenton
, “
Mechanism for amplitude alternans in electrocardiograms and the initiation of spatiotemporal chaos
,”
Phys. Rev. Lett.
118
,
168101
(
2017
).
12.
E. M.
Cherry
and
F. H.
Fenton
, “
A tale of two dogs: Analyzing two models of canine ventricular electrophysiology
,”
Am. J. Physiol.-Heart Circ. Physiol.
292
,
H43
H55
(
2007
).
13.
C.
Cherubini
,
S.
Filippi
, and
A.
Gizzi
, “
Electroelastic unpinning of rotating vortices in biological excitable media
,”
Phys. Rev. E
85
,
031915
(
2012
).
14.
R. H.
Clayton
,
O.
Bernus
,
E. M.
Cherry
,
H.
Dierckx
,
F. H.
Fenton
,
L.
Mirabella
,
A. V.
Panfilov
,
F. B.
Sachse
,
G.
Seemann
, and
H.
Zhang
, “
Models of cardiac tissue electrophysiology: Progress, challenges and open questions
,”
Prog. Biophys. Mol. Biol.
104
,
22
48
(
2011
).
15.
N.
Cusimano
,
A.
Bueno-Orovio
,
I.
Turner
, and
K.
Burrage
, “
On the order of the fractional laplacian in determining the spatio-temporal evolution of a space-fractional model of cardiac electrophysiology
,”
PLoS One
10
,
e0143938
(
2015
).
16.
S.
Dehin
, “
Cardiovascular gap junctions
,” in
Advances in Cardiology
(
Karger
,
2006
), Vol.
42
.
17.
V. S.
Deshpande
,
R. M.
McMeeking
, and
A. G.
Evans
, “
A bio-chemo-mechanical model for cell contractility
,”
PNAS
103
,
14015
14020
(
2006
).
18.
V. S.
Deshpande
,
M.
Mrksich
,
R. M.
McMeeking
, and
A. G.
Evans
, “
A bio-mechanical model for coupling cell contractility with focal adhesion formation
,”
J. Mech. Phys. Solids
56
,
1484
1510
(
2008
).
19.
M.
Dupraz
and
V.
Jacquemet
, “
Geometrical measurement of cardiac wavelength in reaction-diffusion models
,”
Chaos
24
,
033133
(
2014
).
20.
F. H.
Fenton
and
E. M.
Cherry
, “
Models of cardiac cell
,”
Scholarpedia
3
,
1868
(
2008
).
21.
F. H.
Fenton
and
A.
Karma
, “
Vortex dynamics in three-dimensional continuous myocardium with fiber rotation: Filament instability and fibrillation
,”
Chaos
8
,
20
(
1998
).
22.
F. H.
Fenton
,
E. M.
Cherry
,
H. M.
Hastings
, and
S. G.
Evans
, “
Multiple mechanisms of spiral wave breakup in a model of cardiac electrical activity
,”
Chaos
12
,
852
892
(
2002
).
23.
F. H.
Fenton
,
E. M.
Cherry
,
A.
Karma
, and
W. J.
Rappel
, “
Modeling wave propagation in realistic heart geometries using the phase-field method
,”
Chaos
15
,
013502
(
2005
).
24.
F. H.
Fenton
,
A.
Gizzi
,
C.
Cherubini
,
N.
Pomella
, and
S.
Filippi
, “
Role of temperature on nonlinear cardiac dynamics
,”
Phys. Rev. E
87
,
042717
(
2013
).
25.
S.
Filippi
,
A.
Gizzi
,
C.
Cherubini
,
S.
Luther
, and
F. H.
Fenton
, “
Mechanistic insights into hypothermic ventricular fibrillation: The role of temperature and tissue size
,”
Europace
16
,
424
434
(
2014
).
26.
A.
Gizzi
,
C.
Cherubini
,
S.
Migliori
,
R.
Alloni
,
R.
Portuesi
, and
S.
Filippi
, “
On the electrical intestine turbulence induced by temperature changes
,”
Phys. Biol.
7
,
016011
(
2010
).
27.
A.
Gizzi
,
E. M.
Cherry
,
R. F.
Gilmour
, Jr.
,
S.
Luther
,
S.
Filippi
, and
F. H.
Fenton
, “
Effects of pacing site and stimulation history on alternans dynamics and the development of complex spatiotemporal patterns in cardiac tissue
,”
Front. Physiol.
4
,
71
(
2013
).
28.
A.
Gizzi
,
A.
Loppini
,
E. M.
Cherry
,
C.
Cherubini
,
F. H.
Fenton
, and
S.
Filippi
, “
Multi-band decomposition analysis: Application to cardiac alternans as a function of temperature
,”
Physiol. Meas.
38
,
833
(
2017
).
29.
The Eurowinter Group
. “
Cold exposure and winter mortality from ischaemic heart disease, cerebrovascular disease, respiratory disease, and all causes in warm and cold regions of Europe
,”
Lancet
349
,
1341
1346
(
1997
).
30.
D. E.
Hurtado
,
S.
Castro
, and
A.
Gizzi
, “
Computational modeling of non-linear diffusion in cardiac electrophysiology: A novel porous-medium approach
,”
Comput. Methods Appl. Mech. Eng.
300
,
70
83
(
2016
).
31.
A.
Karma
, “
Universal limit of spiral wave propagation in excitable media
,”
Phys. Rev. Lett.
66
,
2274
2277
(
1991
).
32.
A.
Karma
, “
Scaling regime of spiral wave propagation in single-diffusive media
,”
Phys. Rev. Lett.
68
,
397
400
(
1992
).
33.
A.
Karma
, “
Spiral breakup in model equations of action potential propagation in cardiac tissue
,”
Phys. Rev. Lett.
71
,
1103
(
1993
).
34.
A.
Karma
, “
Electrical alternans and spiral wave breackup in cardiac tissue
,”
Chaos
4
,
461
472
(
1994
).
35.
A.
Karma
,
H.
Levine
, and
X.
Zou
, “
Theory of pulse instabilities in electrophysiological models of excitable tissues
,”
Physica D
73
,
113
127
(
1994
).
36.
J.
Keener
and
J.
Sneyd
,
Mathematical Physiology
(
Spinger–Verlag
,
2009
).
37.
V.
Krinsky
and
A.
Pumir
, “
Models of defibrillation of cardiac tissue
,”
Chaos
8
,
188
203
(
1998
).
38.
S.
Luther
,
F. H.
Fenton
,
B. G.
Kornreich
,
A.
Squires
,
P.
Bittihn
,
D.
Hornung
,
M.
Zabel
,
J.
Flanders
,
A.
Gladuli
,
L.
Campoy
,
E. M.
Cherry
,
G.
Luther
,
G.
Hasenfuss
,
V. I.
Krinsky
,
A.
Pumir
,
R. F.
Gilmour
, Jr.
, and
E.
Bodenschatz
, “
Low-energy control of electrical turbulence in the heart
,”
Nat. Lett.
475
,
235
239
(
2011
).
39.
J. D.
Murray
,
Mathematical Biology
(
Springer International Publishing
,
2002
).
40.
A. A.
Niederer
 et al, “
Verification of cardiac tissue electrophysiology simulators using an n-version benchmark
,”
Philos. Trans. R. Soc. A
369
,
4331
4351
(
2011
).
41.
J. M.
Pastore
and
D. S.
Rosenbaum
, “
Role of structural barriers in the mechanism of alternans-induced reentry
,”
Circ. Res.
87
,
1157
1163
(
2000
).
42.
J. M.
Pastore
,
S. D.
Girouard
,
K. R.
Lautira
,
F. G.
Akar
, and
D. S.
Rosenbaum
, “
Mechanism linking t-wave alternans to the genesis of cardiac fibrillation
,”
Circulation
99
,
1385
1394
(
1999
).
43.
H. H.
Pennes
, “
Analysis of tissue and arterial blood temperatures in the resting human forearm
,”
J. Appl. Physiol.
1
,
93
122
(
1948
).
44.
M.
Potse
,
B.
Dubé
,
J.
Richer
,
A.
Vinet
, and
R. M.
Gulrajani
, “
A comparison of monodomain and bidomain reaction-diffusion models for action potential propagation in the human heart
,”
IEEE Trans. Biomed. Eng.
53
,
2425
2435
(
2006
).
45.
M.
Potse
,
R.
Coronel
,
S.
Falcao
,
A.-R.
LeBlanc
, and
A.
Vinet
, “
The effect of lesion size and tissue remodeling on ST deviation in partial-thickness ischemia
,”
Heart Rhythm
4
,
200
206
(
2007
).
46.
A. J.
Pullan
,
M. L.
Buist
, and
L. K.
Cheng
,
Mathematically Modelling the Electrical Activity of the Heart: From Cell to Body Surface and Back Again
(
World Scientific
,
2015
).
47.
A.
Pumir
,
A.
Arutunyan
,
V.
Krinsky
, and
N.
Sarvazyan
, “
Genesis of ectopic waves: Role of coupling, automaticity, and heterogeneity
,”
Biophys. J.
89
,
2332
2349
(
2005
).
48.
A.
Pumir
,
S.
Sinha
,
S.
Sridhar
,
M.
Argentina
,
M.
Hörning
,
S.
Filippi
,
C.
Cherubini
,
S.
Luther
, and
V. I.
Krinsky
, “
Wave-train-induced termination of weakly anchored vortices in excitable media
,”
Phys. Rev. E
81
,
010901
(
2010
).
49.
Z.
Qu
,
G.
Hu
,
A.
Garfinkel
, and
J. N.
Weiss
, “
Nonlinear and stochastic dynamics in the heart
,”
Phys. Rep.
543
,
61
162
(
2014
).
50.
T.
Quail
,
A.
Shrier
, and
L.
Glass
, “
Spatial symmetry breaking determines spiral wave chirality
,”
Phys. Rev. Lett.
113
,
158101
(
2014
).
51.
J. B.
Reuler
, “
Hypothermia: Pathophysiology, clinical settings, and management
,”
Ann. Intern. Med.
89
,
519
527
(
1978
).
52.
S.
Rohr
, “
Role of gap junctions in the propagation of the cardiac action potential
,”
Cardiovasc. Res.
62
,
309
322
(
2004
).
53.
P.
Sáez
and
E.
Khul
, “
Computational modeling of acute myocardial infarction
,”
Comput. Methods Biomech. Biomed. Eng.
19
,
1107
1115
(
2016
).
54.
N. J.
Severs
,
A. F.
Bruce
,
E.
Dupont
, and
S.
Rothery
, “
Remodelling of gap junctions and connexin expression in diseased myocardium
,”
Cardiovasc. Res.
80
,
9
19
(
2008
).
55.
S.
Takagi
,
A.
Pumir
,
D.
Pazó
,
I.
Efimov
,
V.
Nikolski
, and
V.
Krinsky
, “
Unpinning and removal of a rotating wave in cardiac muscle
,”
Phys. Rev. Lett.
93
,
058101
(
2004
).
56.
A.
Vinet
,
D. R.
Chialvo
, and
J.
Jalife
, “
Irregular dynamics of excitation in biologic and mathematical models of cardiac cells
,”
Ann. N. Y. Acad. Sci.
601
,
281
298
(
1990
).
57.
A.
Vinet
,
D. R.
Chialvo
,
D. C.
Michaels
, and
J.
Jalife
, “
Nonlinear dynamics of rate-dependent activation in models of single cardiac cells
,”
Circ. Res.
67
,
1510
1524
(
1990
).
58.
M. A.
Watanabe
,
F. H.
Fenton
,
H. M.
Evans
,
S. J.
Hastings
, and
A.
Karma
, “
Mechanisms for discordant alternans
,”
J. Cardiovasc. Electrophysiol.
12
,
196
206
(
2001
).
59.
A. T.
Winfree
,
When Time Breaks down
(
Princeton University Press
,
1987
).