A particular model of cardiac electrophysiology

Mathematical models of electrical signal propagation in cardiac electrophysiology consist of suitable spatially dependent formulations of specified cell models reproducing the response of a single excitable cell to an applied electrical stimulus. All cell models describe the temporal evolution of the transmembrane potential v of the particular excitable cell considered and the changes in v caused by the opening and closing of the various ion channels present in the cell membrane, determining the movement of ions into and out of the cell. When an electrical stimulus is applied to a cell, v quickly increases (depolarizes) and moves away from its equilibrium value, known as the resting potential. If the applied stimulus is strong enough to push v above a certain threshold level, once the stimulus is removed the potential does not immediately return to rest but rather undergoes a large excursion—an action potential (AP)—before eventually reaching the resting state.

In order to account for pulse propagation, spatial dependence is typically introduced via the bidomain formulation of the considered cell model (for example, see [12]). The physiologically grounded bidomain model is often replaced in practice by a less computationally demanding formulation [13], known as the monodomain, based on some simplifying assumptions on the conductivity properties of the tissue in which electrical propagation occurs. In this work we focus on the monodomain formulation of the cell model originally developed by Beeler and Reuter [14] in the study of the excitability properties of ventricular myocardial fibres. The mathematical model considered can be written as the following coupled system: (1) (2) where C_m is the membrane capacitance per unit area, χ is the cell surface-to-volume ratio, z is a suitable vector of secondary variables used in the description of the dynamics of the ion channels in the cell membrane, I_ion is the sum of all transmembrane ionic currents, D is an effective conductivity tensor, and I_stim is the electrical stimulus.

The specification of the vector z, the vector-valued function f right-hand side of Eq (2), and the analytic expression of I_ion in Eq (1) depend on the particular type of cell model considered and the level of detail that one wants to incorporate into the cell model. The Beeler–Reuter (BR) model is probably one of the simplest cell models that captures the correct form of the cardiac AP of ventricular cells in a satisfactory and meaningful way, and its mathematical formulation consists of a system of eight ordinary differential equations (see S1 Appendix for the complete model description). Considering a higher level of detail in the cell model and the corresponding bidomain formulation (as in the case of state-of-the-art models of electrical propagation in cardiac tissue) comes at a considerable computational cost and is outside of the scope of the presented study. However, the numerical methodology proposed here can be easily generalised to more sophisticated electrophysiological models.

Fractional Laplacian on an insulated finite domain

Space-fractional differential equations are differential equations in which the classical integer-order differential operators in space are substituted by non-integer counterparts, often referred to as fractional operators. As previously mentioned, differential equations involving non-integer order derivatives have shown to be very powerful tools in the description of transport phenomena whose characteristics substantially deviate from the classical Gaussian and Markovian assumptions and have been adopted as mathematical models for a variety of practical applications.

Among fractional operators in space, the fractional Laplacian (−Δ)^α/2 [15] of order α plays an important role due to its probabilistic interpretation. This operator is naturally defined on the entire , n = 1, 2, 3, and only when α = 2 one recovers the local character of the definition. However, especially when dealing with applications of practical interest, relationships between the observed quantities often hold only within specified bounded regions , and hence suitable mathematical models for these applications must be defined on finite domains.

The restriction of fractional operators to finite domains is by no means trivial and a fundamental role in the definition and implementation of the solution of the problem is played by the particular boundary conditions to which the model is coupled. With fractional operators the mere specification of a local condition at the boundaries is no longer sufficient for the problem to be well-posed, and even though well-established methodologies have been developed for the case of homogeneous Dirichlet boundary conditions (see for example [16]), little is known about all other boundary constraints.

For the particular application considered in this paper, typical boundary conditions are defined under the assumption that the finite domain considered is insulated and hence, in the standard case, they are specified as homogeneous Neumann boundary conditions.

On an insulated one-dimensional domain [0, L], the continuous fractional Laplacian (−Δ)^α/2 can be defined via the spectral mapping theorem and the eigenvalue decomposition of the continuous standard Laplacian (−Δ), as originally proposed by Ilić et al. [17, 18]. This definition of the fractional operator correctly represents the insulation properties of the finite domain in both the standard and the fractional cases (as proved by Cusimano [19]) and will be used in this work to formulate the space-fractional modification of the BR monodomain model and to design an efficient solution strategy for the fractional problem considered.

Since the theoretical results obtained in Cusimano [19] have been proved in the one-dimensional case and for α ∈ (1, 2], in this study we focus only on this particular range of the fractional order and on one spatial dimension, leaving the generalisation of the proposed study to higher dimensions as future extension of this work.

Space-fractional Beeler–Reuter model

Let us consider the BR monodomain formulation on an insulated finite interval [0, L], that is, system Eqs (1)–(2) in which I_ion, z and f are given by the BR cell model as defined in S1 Appendix. We introduce the following space-fractional modification of system Eqs (1)–(2): (3) (4) where D is now simply a constant and (−Δ)^α/2 is the one-dimensional fractional Laplacian of order α ∈ (1, 2] on [0, L], defined via its spectral expansion and embedding in its definition the correct representation of insulating boundary conditions at both ends of [0, L]. A physiological justification for the use of fractional diffusion models in the description of structurally heterogeneous excitable media, as a representation of the modulation of the electrical field of a homogeneous conductor by the secondary electrical sources associated with tissue inhomogeneities, has been presented in [11].

Specifically, (−Δ)^α/2 is the continuous operator with eigenfunctions , defined for x ∈ [0, L] as (5) with , for j ≥ 1, and corresponding eigenvalues , that is, the fractional powers of the eigenvalues of the standard continuous Laplacian (−Δ), coupled to homogeneous Neumann boundary conditions on [0, L], defined as for all j ≥ 0.

Throughout the rest of the paper we assume that the cell model description in system Eqs (3)–(4) is given by the BR model. However, most of the considerations made in this section do not depend on the particular specification of I_ion, z, f, and hence could be easily extended to the monodomain formulation of any other cell model of interest.

Solution methodology

In order to compute the numerical solution of system Eqs (3)–(4) on a given temporal interval [0, t_f], with t_f > 0, we adapt the strategy proposed by Whiteley [20] for the standard monodomain formulation of an excitable media model and combine it with a spectral approach for the solution of the space-fractional reaction-diffusion Eq (3). The main idea behind Whiteley’s approach [20] is to introduce a temporal grid for [0, t_f] and starting from a given initial condition, to compute the numerical solution at each time iteration in two steps. First, the value of the transmembrane potential v is updated by solving Eq (3), while keeping the value of z fixed at the previous time iteration. Second, the vector z is updated by solving the system of differential Eq (4) in which the vector-valued function f has been evaluated by using the newly computed value of v.

The only difference between the strategy presented here and the one proposed by Whiteley [20] for the standard case is the way in which the value of the transmembrane potential v is updated at each time step. Here, to compute the solution of Eq (3) we exploit the spectral definition of the fractional Laplacian and extend the spectral approach presented in [19] to the case of a space-fractional reaction-diffusion equation (similarly to that done by Bueno-Orovio et al. [21]).

In particular, we look for the solution v(x, t) as a linear combination of the orthonormal eigenfunctions of the one-dimensional fractional Laplacian (−Δ)^α/2 on the insulated domain, that is, in the form , where each is the integral coefficient (these eigenfunctions can be shown to form a basis for a suitable Hilbert space in which the solution to the problem is sought—see for example Hutson et al. [22] for more details). Let us introduce a grid of time points t_k with uniform time step Δt. Let v^k: = v(x, t_k), and z^k: = z(x, t_k). The solution of Eq (3) is computed at each time step via the semi-implicit Euler scheme in which the spatial derivative of v is evaluated at t_k+1, while the ionic term is considered fixed at t_k, that is, by solving (6) By using the spectral definition of the fractional Laplacian, that is, the fact that , and due to the orthonormality of the eigenfunctions φ_j, the computation of v^k+1 can be fully diagonalised and each integral coefficient can be updated independently at each time step as follows: (7) where once again is the set of eigenvalues of the fractional operator (−Δ)^α/2 and is the j-th integral coefficient of the spectral expansion of the source term evaluated at v^k and z^k.

Once the updated value of v has been computed, for the vector of secondary variables z we adopt Whiteley’s method without the need of any modification and hence consider a fully implicit Euler scheme to solve the ODE system Eq (4). Due to the particular form of the right-hand side of the differential equation of six out of seven of the secondary variables in the BR model (linear in the considered component of the vector z and independent from all other secondary variables—see S1 Appendix), the backward differentiation formula allows us to obtain an unconditionally stable explicit expression for the updated value. When the right-hand side of the equation is not linear, as in the case of the calcium concentration for the BR cell model, an inexact Newton method (in which the Jacobian is “frozen” at the first iteration) is used instead.

Note that for practical reasons the eigenfunction expansion of the solution and the source term must be truncated after a finite number of terms, and the integral coefficients and must be somehow approximated. Let us consider a uniform spatial grid of N + 1 nodes x_i = ih, for the discretisation of the insulated finite domain [0, L], with , and i = 0, 1, …, N. In all our numerical simulations we will consider the first N + 1 eigenfunctions of the spectral expansion and approximate all integral coefficients by using the trapezoidal rule on the considered nodes. This approach can be shown to be equivalent to computing the solution on [0, L] via the discrete cosine transform, assuming symmetry with respect to both ends of the spatial interval and periodicity with period 2L for the extended solution defined on .

Convergence

In order to assess convergence of the numerical solution produced by our solution method and to identify acceptable values of the key parameters N and Δt for our numerical simulations, we follow a similar approach to the one proposed by Pathmanathan et al. [12]. We consider successive refinements of a key parameter q (either q = N or q = Δt) and compare the resulting solutions via four relative metrics against a reference solution (computed with the finest spatial mesh, smallest uniform time step). Specifically, we consider the relative Euclidean norm of the temporal profile of the solution computed at a selected node P_i ∈ [0, L], the relative error in the AP duration at 50% and 90% repolarization at the node P_i (APD₅₀ and APD₉₀, respectively), and the relative change in the conduction velocity (CV) measured as the ratio of the distance between two nodes of [0, L] and the difference in the corresponding activation times. Activation time of a given node was defined as the time when the voltage reached 10% of its full depolarization. APD₅₀ and APD₉₀ were defined as the time intervals between activation and the time when the voltage (at the considered point) reached 50% and 90% of repolarization to its resting value, respectively. In both cases, linear interpolation was used to obtain better resolved time values.

In order to compute the quantities of interest for our metrics, we focus on three equally spaced nodes, namely P₁(x = 0.25), P₂(x = 0.5), and P₃(x = 0.75). In Fig 1, we display the above mentioned metrics as a function of the key parameter N for four different values of α. For each α we use as reference solution the solution computed with the considered α on a grid of 3201 spatial nodes.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. The four relative metrics as a function of the key parameter N (the number of eigenfunctions), for four selected values of α.

The solution is computed for all values of α on the time interval [0,500 ms] with a uniform time step Δt = 0.01 ms. All relative errors are evaluated by comparing the solution for a specific combination of α and N with the solution obtained for the same α when N = 3200. In all simulations, the stimulus I_stim = 40 μA⋅cm⁻³ was applied on [0,0.05 cm] at t_stim = 10 ms for two consecutive milliseconds and then removed.

https://doi.org/10.1371/journal.pone.0143938.g001

Fig 1(A), 1(B) and 1(C) display the value of the relative Euclidean norm and the relative APD₅₀ and APD₉₀ only at the mid-point P₂ but very similar qualitative results were obtained for these metrics also at P₁ and P₃ (results not shown). In order to compute the relative velocity of Fig 1(D) we use the activation times of and the distance between P₁ and P₃.

As is clearly visible in Fig 1(A), as N increases, for all values of α the approximation of the solution over the entire temporal interval improves. However, as we reduce α, a larger value of N is necessary in order to achieve a given accuracy. The other three metrics exhibit a similar overall trend, even though the error decay is not always monotonic and in some instances for a fixed N a lower value of the relative error is obtained when using a smaller α.

In all simulations performed to compute the errors displayed in Fig 1 we considered Δt = 0.01 ms and only varied N. There is a certain amount of a posteriori analysis in this choice and the particular value chosen is the result of a series of simultaneous refinements of both key parameters, N and Δt.

Fig 2 shows the behaviour of the four metrics as a function of the key parameter Δt, for the same four values of α. Clearly, as Δt decreases the error in all relative metrics becomes smaller. However, similarly to that observed in Fig 1, for a fixed value of the key parameter Δt, the corresponding value of rel_norm_full_P₂ increases as α is reduced, indicating the need for refining the temporal mesh as α is reduced to achieve a given level of accuracy in the numerical approximation. The effect of α on the other three metrics is less significant. As the time step is varied, the overall trend exhibited by these three metrics is indeed quite similar and for a given Δt the magnitude of the error is not affected by the use of different values of α. Once again, we computed the value of the first three metrics also at P₁ and P₃ obtaining very similar results (not shown).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. The four relative metrics as a function of the key parameter Δt (the time step), for four selected values of α.

The solution is computed for all values of α on the time interval [0,500 ms] with N = 400. All relative errors are evaluated by comparing the solution for a specific combination of α and Δt with the solution obtained for the same α when Δt = 6.25⋅10⁻⁴ ms. In all simulations, the stimulus I_stim = 40 μA⋅cm⁻³ was applied on [0,0.05 cm] at t_stim = 10 ms for two consecutive milliseconds and then removed.

https://doi.org/10.1371/journal.pone.0143938.g002

In light of the proposed convergence results, we decided to set N = 400 (corresponding to a mesh size equal to h = 0.0025 cm) and Δt = 0.01 ms. With this choice of key parameters we observe that rel_norm_full ≤10⁻² for all the considered values of α in Fig 1, meaning that the percentage error in the solution approximation is smaller than 1% over the entire temporal interval [0, t_f]. Similarly, from Fig 2 we see that the error computed in the calculation of the CV with this choice of key parameters is smaller than 1% as well. Finally, Fig 2(B) and 2(C), show that when N = 400 the percentage error in the computation of both the APD₅₀ and APD₉₀ is smaller than 0.01% and hence, ensures that the value of these two biomarkers is computed with an even higher precision.

Effects of a reduced α on AP shape and propagation across the domain

In this section we investigate the effect of reducing the value of the fractional order α on the spatio-temporal propagation of the excitation wave triggered by the applied stimulus.

Temporal effects.

We begin by studying how changes in α modify the temporal profile of the solution and focus on the differences observed in the AP recorded at a given set of points in the spatial domain. Fig 3 shows the AP at P₁, P₂, and P₃ generated by the same applied stimulus I_stim = 40 μA⋅cm⁻³ for α = 2 and α = 1.5.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Temporal profile of the solution for t ∈ [0,500 ms] at three equally spaced nodes P₁ (x = 0.25), P₂ (x = 0.5), P₃ (x = 0.75) of the spatial domain [0, 1 cm].

(A) Solution obtained when α = 2. (B) Solution obtained when α = 1.5. In both cases the stimulus I_stim = 40 μA⋅cm⁻³ was applied on [0,0.05 cm] at t_stim = 10 ms for two consecutive milliseconds and then removed.

https://doi.org/10.1371/journal.pone.0143938.g003

In the standard case depicted in Fig 3(A), we see how the excitation wave triggered by the applied stimulus travels across the spatial domain and produces at each node the AP shape typical of the BR cell model, that is, characterised by the presence of a depolarization peak followed by a small concavity known as an early repolarization, and finally a large dome and a slow recovery to the transmembrane potential resting state. When a fractional value for the parameter α is considered (here set equal to α = 1.5), the electrical stimulus is still able to trigger an excitation wave travelling across the spatial domain and each considered node, once depolarized, goes through what seems to be the exact same AP (differences discussed below). However, it can be observed that there are interesting differences between the two plots, such as the larger times required for each of the three nodes to be depolarized in the fractional case and the more pronounced early repolarization phase when α < 2.

In order to better visualise these differences we focus on a single node in space and compare the AP shape generated at this node by I_stim = 40 μA⋅cm⁻³ for ten different values of α ∈ (1, 2]. The results obtained at the mid-point P₂ are plotted in Fig 4. Very similar qualitative results were also found at P₁ and P₃ (results not shown).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Temporal profile for t ∈ [0,500 ms] of the solution obtained with ten different values of α at the mid-point of the spatial domain [0, 1 cm].

In all cases the stimulus I_stim = 40 μA⋅cm⁻³ was applied on [0,0.05 cm] at t_stim = 10 ms for two consecutive milliseconds and then removed.

https://doi.org/10.1371/journal.pone.0143938.g004

As α is reduced we observe:

1. a decrease in both the depolarization peak height and the dome peak;
2. a visibly more pronounced early repolarization phase;
3. an increase in the AP foot width (the AP foot is defined as the first portion of the quick AP depolarization);
4. a shift towards the right of the AP, corresponding to a later activation time for the considered node and indicating a reduction in the CV of the propagation pulse.

In order to quantify the effects observed qualitatively in Fig 4 and to establish specific relationships between the considered quantities and the fractional order α, we then fit the data obtained from our numerical simulations with suitable functions of the parameter α. Specifically, in Fig 5(A) we consider the data stored for the AP depolarization peak height and clearly see that a reduction in α produces a linear decrease of this quantity. On the other hand, we find that the relationship between α and the dome peak (Fig 5(B)) and between α and the early repolarization minimum (Fig 5(C)) is nonlinear, and in both cases it can be well described by a power law of the form f(α) = aα^b + c. The value of the parameters a, b, c, and the goodness of fit statistics for the presented fitting curves is reported in the S2 Appendix.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Curve fitting of three data sets obtained from the AP solution profile at the node P₂ for ten different values of α ∈ (1, 2].

(A) AP depolarization peak height, (B) dome peak height, (C) early repolarization minimum.

https://doi.org/10.1371/journal.pone.0143938.g005

A better visualisation of the effects produced by varying the fractional order α on the foot of the AP profile at P₂ is given in Fig 6. By aligning the AP obtained with different α so that the activation time of each solution profile coincides, we clearly see the change in AP foot width and decrease in curvature as α is reduced.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. AP foot of the solution profile obtained at the node P₂ for ten different values of α ∈ (1, 2].

To aid the visualisation of the differences produced by varying the fractional order, we align the AP foot of the ten solution profiles considered so that the activation time (AT) of the node P₂ coincides for all values of α.

https://doi.org/10.1371/journal.pone.0143938.g006

Spatial effects.

The differences in the activation time of P₂ observed for different values of α are a clear indication that the introduction of a fractional order does not only modify the AP profile but also has a sensible effect on the spatial spread of the excitation wave.

Fig 3 suggests that for a fixed α the depolarization front of the excitation wave triggered by the applied stimulus moves across the spatial interval (away from the stimulus site) with constant speed. In Fig 7(A) we plot the activation time of forty-eight equally spaced nodes in [0, L] (we exclude a small part of the domain surrounding the origin to avoid artificial effects due to the applied stimulus), for ten different values of α ∈ (1, 2]. We see that the excitation pulse travels linearly in space for all values of α (with very slight variations from linearity only in proximity of the boundary x = L, further away from the stimulus site). The CV of the triggered excitation wave, corresponding to the reciprocal of the gradient of each of these straight lines, is hence constant across the spatial domain. However, its value depends on α and decreases as α is reduced.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7. Activation time of forty-eight equally spaced nodes in the spatial domain [0.06 cm, 1 cm] and conduction velocity for ten different values of α.

The CV corresponding to the reciprocal of the gradient of each of the ten straight lines in (A) was considered as data in (B) and a quadratic dependence of these data points from the fractional order α was established.

https://doi.org/10.1371/journal.pone.0143938.g007

In particular, we find that the relationship between the CV and α is not arbitrary, but rather it can be expressed as a quadratic function of α as shown in Fig 7(B). Once again, see S2 Appendix for specific values of the parameters involved in the quadratic fit and the corresponding goodness of fit statistics.

Dispersion of action potential duration.

In this section we make some considerations on the action potential duration (APD), and as a measure of APD we consider the APD₉₀ previously introduced in the text. One important characteristic observed in standard simulations of excitable media models is the dependence of the repolarization front advance on the domain geometry and on the proximity to the boundaries [23]. These spatial variations result in differences of APD across the domain already in the standard case. Specifically, an initial reduction in the APD due to the proximity to the stimulus site is followed by the adjustment of the APD about an almost constant level before sensibly reducing in proximity of the insulated boundary x = L. In order to compare the effect of different values of the parameter α on the variations produced in APD across the domain, we measure the APD dispersion (computed here for each α as the difference between the value of APD at a given spatial node and the maximum APD value recorded for the same value of α over the entire spatial domain).

Differences in the APD across the spatial interval are present for all values of α (Fig 8). However, as α decreases we find a more pronounced effect in the dispersion across an increasing portion of the domain and the reduction in the APD in proximity of the boundary x = L becomes significantly larger.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 8. APD dispersion at forty-eight equally spaced nodes in [0.06 cm, 1 cm] for ten different values of α.

For a fixed α, the APD is computed from the AP solution profile recorded at each considered node in the spatial interval. Dispersion of APD is then defined as the difference between the computed value of APD at that node and the maximum APD value recorded for the same value of α over the entire spatial domain.

https://doi.org/10.1371/journal.pone.0143938.g008

The shortening in APD in proximity to the boundary is essentially linked to the finiteness of the spatial domain and to the increase in the repolarization electrotonic current due to the loss of the more depolarized neighbouring cells, as observed by Cherry and Fenton [23]. When α < 2, we add to the reflective nature of the insulated boundary the fractional character of the solution. Hence, as the repolarization front moves closer to the end of the spatial interval, the influence on the solution of the increasing number of more repolarized nodes (closer to the stimulus site) becomes stronger as α decreases. Therefore, the acceleration produced on the repolarization becomes greater and the shortening of APD more marked.

We stress that the presented results on APD dispersion are not a simple effect due to the small domain size chosen for our simulations. In fact, when longer cables are considered dispersion of APD is still present and becomes increasingly pronounced as the fractional parameter α is reduced (results not shown).

Reduced α versus reduced conductivity.

The introduction of a fractional operator in space in the BR monodomain formulation resulted in the solution of the model being characterised by a reduction of the CV and various changes in a number of other features of the AP shape. Clearly, a reduction in the CV could be obtained also by considering the standard BR monodomain formulation (α = 2) and simply reducing the conductivity parameter D in Eq (3). We conclude this section by providing a comparison between the solution obtained from the space-fractional modification of the BR monodomain model for a given fractional index and the one found by reducing the conductivity coefficient in the standard case in order to match the value of the CV obtained in the fractional case. Specifically, in Fig 9 we consider the case α = 1.5, D = 1 mS⋅cm⁻¹, generating an excitation wave with CV = 12.2903 cm⋅s⁻¹, and compare the AP generated at the usual three equidistant nodes P₁, P₂, P₃, to the solution obtained at the same nodes when α = 2 and D = 0.093 mS⋅cm⁻¹.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 9. Temporal evolution for t ∈ [0,400 ms] of the standard solution with reduced conductivity and the fractional solution with unit conductivity parameter at three equally spaced nodes P₁, P₂, P₃, in the spatial domain [0, 1 cm].

The continuous blue line corresponds to the standard solution obtained with α = 2, D = 0.093 mS⋅cm⁻¹. The dashed red line represents the fractional solution obtained by setting α = 1.5, D = 1 mS⋅cm⁻¹. In both cases the stimulus I_stim = 40 μA⋅cm⁻³ was applied on [0,0.05 cm] at t_stim = 10 ms for two consecutive milliseconds and then removed.

https://doi.org/10.1371/journal.pone.0143938.g009

Even though the choice of these two combinations of parameters α and D was made so that the CV of the corresponding solutions coincided, once again we see that the fractional solution (red) and the standard solution (blue) are quite different in the AP shape. In fact, the fractional solution is characterised by a wider AP foot, a lower depolarization peak, an accentuated early repolarization phase, and a slightly lower dome peak, highlighting the fundamentally different nature of the two cases considered. We also see that at P₃ the distinction between the repolarization phase of the AP generated in the standard and the fractional cases becomes visible, identifying also the more pronounced shortening of APD as another distinctive feature of the fractional modification of the model. Achieving a larger CV in the fractional case can be simply obtained by increasing the conductivity parameter D.