Dynamical modelling of cardiac electrical activity using bidomain approach : The effects of variation of ionic model parameters

This work presents the dynamical modelling of cardiac electrical activity using bidomain approach. It focuses on the effects of variation of the ionic model parameters on cardiac wave propagation. Cardiac electrical activity is governed by partial differential equations coupled to a system of ordinary differential equations. Numerical simulation of these equations is computationally expensive due to their non-linearity and stiffness. Nevertheless, we adopted the bidomain model due to its ability to reflect the actual cardiac wave propagation. The derived bidomain equations coupled with FitzHugh-Nagumo’s ionic equations were time-discretized using explicit forward Euler method and space-discretized using 2-D network modelling to obtain linearized equations for transmembrane potential Vm, extracellular potential φe and gating variable w. We implemented the discretized model and performed simulation experiments to study the effects of variation of ionic model parameters on the propagation of electrical wave across the cardiac tissue. Time characteristic of transmembrane potential, Vm, in the normal cardiac tissue was obtained by setting the values of ionic model parameters to 0.2, 0.2, 0.7 and 0.8 for excitation rate constant ε1, recovery rate constant ε2, recovery decay constant γ and excitation decay constant β respectively. Changing the values of ε1, ε2 to 0.04 and 0.28 respectively, the obtained Vm showed a time dilation at 0.04 indicating cardiac arrhythmia but no significant change to Vm was observed at 0.28. Also, changing β to 0.3 and 1.1 and γ to 0.4 and 1.2 sequentially, there was no significant change to the time characteristic of Vm. The obtained results revealed that only decrease in є1, є2 impacted significantly on the cardiac wave propagation.


INTRODUCTION
Cardiovascular disease is the number one leading cause of death and disability globally.According to Kwunife and Aguwa [1], in Nigeria, cardiovascular diseases alone accounted for 9.2% of the total deaths (WHO, 2001), killing even more than malaria (WHO, 2002).The global status report on non-communicable diseases by the World Health Organisation (WHO) [2] revealed that 80% of the 17 million deaths due to cardiovascular disease occur in low and middle income countries.
Cardiac arrhythmia is a name for a large family of cardiac behaviour that shows abnormalities in the electrical behaviour of the heart [3] and has been considered the main source of sudden cardiac death, preventing blood flow to different parts of the body.Prompt diagnosis and appropriate treatment can remedy the devastating physiological consequences due to cardiac arrhythmia.
Mathematical modelling and simulation of cardiac electrical activities can play a vital role in diagnosing cardiac electrical abnormalities by providing valuable information about the functional status of the heart.Mathematical models of the heart can be used to simulate heart conditions and the effects of certain drugs designed to treat them [4].The development of drugs for treatment of cardiovascular and other diseases is often very expensive.Computer simulation can reduce this cost by reducing the number of physical experiments needed in designing a drug [5] as well as offer an easy way of assessing the medical well-being of patients and detecting any potential side effects of drugs on cardiac rhythms.It is a cost effective method for understanding the functioning of not only heart, but also other human diseases and biological systems.Usually, models of cardiac electrical activities are governed by differential equations consisting of partial differential equations (PDEs) coupled to a system of ordinary differential equations (ODEs) [6,7].While the PDEs describe the propagation of the electrical signal through the cardiac tissue, the ODEs describe the electrochemical reaction in the cell.
Cardiac wave propagation within the cardiovascular system is a multiscale problem with the electrical activity through ion channels of the individual cells driving a wavefront over the whole cardiac tissue, the heart.The model that gives complete description of such a complex process is the anisotropic bidomain model [6].It consists of a system of two non-linear partial differential equations coupled to a system of ordinary differential equations.
Several researches have been conducted to describe the electrical characteristics of human heart using an efficient mathematical model such as the bidomain model [8][9][10][11][12].Despite the versatility of the model to give realistic simulation of the electrical wave propagation across the cardiac tissue, its implementation is still considered a demanding scientific computing problem due to the requirement of very fine computational grids size.One of the remedies to these computational challenges is the use of the monodomain model.This model consists of a single non-linear partial differential equation coupled with the same system of ordinary differential equations for the ionic currents.Although, it has been reported that the central processing unit requirements are reduced when simplifying the bidomain model to a monodomain model, but both models still encounter computational difficulties because of the need for fine meshes and small time-steps [13].To lessen the computational requirements of the bidomain model, the present work takes advantage of discrete simplification of the model where the cardiac tissue is represented by interconnected network of cells, each individually described by a given system of ionic model.
This work presents dynamical modelling of cardiac electrical activity using bidomain approach.The main focus is on studying the effects of variation of the ionic model parameters on the propagation of electrical wave across the cardiac tissue.The rest of this article is organised as follows.The mathematical model of the cardiac electrical activity is derived and explained in Section 2. In Section 3, the numerical simulation of the model is presented.This is followed by the simulation results and discussion in Section 4. The article ends with a conclusion in Section 5.

MODELLING OF CARDIAC ELECTRICAL ACTIVITY
Three fundamental electrical laws govern the physics of the electrical field model that describes the wave propagation in the thoracic volume [6,7,14]:  The electrical charge conservation law;  The electrical conduction law (Ohm's Law);  The consequence to the electromagnetic induction law.
Using these laws together with the divergence theorem and by treating thoracic volume as the volume of good conductors, the following fundamental equations emerge [15][16][17]: where j is the current density in A/m, σ is the conductivity in S/m, E is the electric field in V/m and ϕ is the electric potential in V.
Equations ( 1) to ( 4) are respectively called electrical charge conservation law, Ohm's law, consequence to electromagnetic induction law and modified Ohm's law.In order to fully adapt Equations ( 1) to (4) to this work, the bidomain model assumes the cardiac tissue as a homogenized two-phase ohmic conducting medium with one phase representing the intracellular space and the other, extracellular space.The phases are linked by a network of resistors and capacitors representing the ion channels and the capacitive current driven across the cell membrane due to a difference in potential respectively as shown in Figure 1 [18].
Considering a post homogenization process, the intracellular and extracellular domains can be assumed to be superimposed to occupy the whole heart volume Ω H [6,[19][20][21] and this also applies to the cell membrane.Hence, the average intracellular and extracellular current densities, i and e , conductivity tensors i j j  and e  and electric potentials i  and e  are defined in Ω H . Application of Equation ( 1) to the heart volume gives Equations ( 5) and ( 6) respectively: where j i is the intracellular current density in A/m 2 , j e is the extracellular current density in A/m 2 , m  is the surface-to-volume ratio of the cell membrane per meter (m -1 ), and m I is the cell membrane current in ampere (A).
The use of Equation ( 4) in Equation ( 5) yields: The transmembrane potential, m V , defined as difference in potential between intracellular and extracellular spaces is represented by: where i  is the intracellular electric potential in volt (V), and e  is the extracellular electric potential in volt (V).
Expressing Equation ( 7) in terms of the transmembrane potential V m , we obtain: This becomes: By extending the cell model formulated by Hodgkin and Huxley in 1952 as reported in Matthias [18] with its electric circuit equivalence diagram as shown in Figure 2 to the context of bidomain, we obtain Equation ( 11) which together with Equations ( 5) and ( 8) gives another fundamental equation of bidomain model represented by Equation (12).

 
, where m is the membrane capacitance in per area unit, I m is the membrane current in ampere (A), I app is the excitation current in ampere (A) and I ion is the ionic current in ampere (A).
The ionic variable w satisfies a system of ODE of the type given by: where g is a vector-valued function.
Hence, the bidomain equations describing the cardiac electrical activity within the cardiovascular system are summarised as: The bidomain model described by Equations ( 10), ( 12) and ( 13) depicts a non-linear elliptic equation for the extracellular potential e  coupled with the parabolic differential equation for the transmembrane potential V m as well as an ordinary differential equation representing the ionic current w.While Equations ( 10) and (12) give the description of the electrical propagation through the cardiac tissue, Equation ( 13) describes the electrochemical reaction in the cell.For complete definition of bidomain model, it must be coupled with an ionic model and appropriate initial and boundary conditions.

Ionic Model
Various ionic models exist for use in cardiac modelling problems.However, of importance to this work is Fitz-Hugh-Nagumo's ionic model and is of the form given by Equations ( 14) and ( 15) [22]: where 1 2 , , ,     are positive constant parameters respectively called excitation rate constant, recovery rate constant, recovery decay constant and excitation decay constant.It was considered principally in this work because of these essential parameters it contains making our dynamical analysis of cardiac electrical activity using bidomain approach a comprehensive study apart from its straightforwardness and wider theoretical and computational applications.

Initial and Boundary Conditions
The bidomain equations described by Equations ( 10), (12), and ( 13) are subjected to the initial conditions given by Equation ( 16): The boundary conditions imposed on this system of Equations ( 10), ( 12) and ( 13) is that of a sealed boundary where no current flows across the boundary between the intracellular and extracellular domains, that is: where n is the normal vector to the domain boundary.

NUMERICAL SIMULATION
With the complete system of differential equations describing the cardiac electrical activity fully derived, we describe in this section the numerical implementation for solving ϕe, Vm and w (the main variables of bidomain equations) using Equations ( 10), ( 12), ( 13), ( 14) and ( 15), and the initial and boundary conditions in Equations ( 16) and ( 17).

Discretization
Discretization of boundary value problems often lead to solving system of linear equations with matrices having bound, block or sparse forms [23].The bidomain equations are both time and space dependent and the process of discretization has to be dealt with separately.Various explicit and implicit time discretization techniques have been explored for use in many modelling works involving differential equations, though, implicit method offers greater stability than explicit method but the latter is a very simple and straightforward method and problem of instability can be reduced by making the time step size very small.On this account, we employed explicit forward Euler time discretization scheme to linearize Equations ( 12) and ( 13), which contain time derivatives and for the space-discretization of Equations ( 10) and ( 12), 2-D discrete (network) modelling was adopted making it possible for us to replace Del operators with intracellular and extracellular admittances G i and G e .The final discretized equations are given by Equations ( 18) to (20).(18) with m  and m assumed unity; G i , the intracellular admittance matrix equivalent to   where G e , the intracellular admittance matrix equivalent to e     .G i and G e were constructed by considering node arrays N x -by-N y defined in the 2-D network domain to be linked by network of resistors arranged along x-and y-direction with η ex , η ey, η ix, and η iy representing the extracellular and intracellular resistance values along these directions.These resistors arrays were then transformed into matrices in the implementation code.This procedure is illustrated here using 2 by 3 nodes in Figure 3 as an example and the result was generalized to the N x by N y nodes considered in this work.
Each resistor was represented by five indices; the x i and y i indices of one the nodes to which the resistor was connected, the x j and y j indices of the other nodes to which the resistor was connected and resistor value η.These two indices were then converted into a single index (encircled numbers in Figure 3) which now represented the first node   i x to which the resistor was connected and the second node   j y to which the resistor was connected.It is these two single-indexed numbers that were stored in the matrices of G i and G e to represent the positions p ij where the resistors are to be placed in the admittance matrices.
Scanning through the network of resistors in Figure 3 from left to right and right to left, it was observed that the resistor values were equal in both directions, that is, η 12 is equal to η 21 and so on.Based on this, the admitance matrices G i and G e (6 by 6 matrices) are as shown in the matrix given by Equation (21).where g x is 1 x  and g y is 1 y  The intracellular and extracellular admittance matrices G i and G e finally obtained represent homogeneous but anisotropic system since the resistors appeared the same everywhere but current flows in the two different directions x and y due to difference in the resistances in the x and y directions were different.
because of its enriched mathematical library and well designed Graphical User Interface (GUI) for displaying graphical representation of results.Another benefit of Java is its portability across various operating systems.The portability however comes with price of slow performance in comparison to C/C++.A flow chart for the implementation algorithms is shown in Figure 4 below.Simulation experiments were carried out on a 4GB RAM, Intel (R) Core (TM) i7 CPU M620 @ 2.67GHz and 32-bit operating system computer.Running time of one simulation experiment ranged between 3 and 5 minutes for unchanged and changed parameters.

Implementation
The discretized equations were implemented in Java programming language (Java 6.0 version).Java is an object-oriented programming language.We adopted Java

SIMULATION RESULTS AND DISCUSSION
The simulations carried out in this work were grouped into two categories, namely;  Simulation using the standard values of the basic parameters.
 Simulation using the varying values of the basic parameters.

Simulation Using the Standard Values of the Basic Parameters
We performed simulation experiment using the developed 2-D Java programme based on the linearized bidomain Equations ( 18), ( 19) and ( 20) and the parameters in Table 1.
The selected cells 4, 6, 8, 10 of the 50-by-50 nodes (cells) specified in 2-D network domain produced the propagated electrical waves in the normal cardiac tissue as shown in Figures 5(a

Simulation Using the Varying Values of the Basic Parameters
Here we present the results of simulations we performed to study the effects of the variation of the parameters є 1 , Table 1.Values of basic parameters [22].

Parameter Value
Excitation rate constant (є 1 ) 0. є 2 , β and γ that characterizes the adopted FitzHugh-Nagumo's ionic model on the electrical wave propagation in the normal cardiac tissue.The standard values considered for these parameters as presented in Table 1 are 0.2, 0.2, 0.7 and 0.8 respectively.

Effects of Variation of є 1 and є 2 on the Cardiac Electrical Wave Propagation
The parameters є 1 and є 2 basically control the dynamics between the transmembrane potential V m and gating variable w.In this work, the two parameters had been considered identical.With the value of є 1 , є 2 as 0.2, the electrical wave propagation in the normal cardiac tissue presented in Figure 5 above was obtained.However, when the value of є 1 , є 2 was changed to 0.04 and 0.28 respectively, some differences were observed.At є 1 , є 2 equals 0.04, the electrical wave obtained for cells 4, 6, 8 and 10 have near-infinite slope values when the electrical potential peaks and relaxes.Time dilation effect where the cardiac excitation cycle extends beyond the normal time was also observed.Figures 6(a)-( d) show this effect of decrease value of є 1 and є 2 on the normal cardiac electrical wave.Comparing Figure 5 and Figure 6, it was observed that each signal in Figure 6 had an extended time for relaxation (resting state) over the corresponding signal in Figure 5.The implication of this time dilation effect is that the cardiac tissue has a delayed repolarization, which is an indication of high risk of arrhythmia (abnormal cardiac electrical behaviour).If this delay persists, it may result in sudden cardiac death.
With the value of є 1 , є 2 changed to 0.28, no significant differences were observed in terms of excitation time before resting state and shape of signal obtained as compared to those of Figure 5.However, it was observed that the signal dwelled more into negative potential before resting.This situation can be regarded as normal since the cell may still be excitable and conduction through the cardiac tissue may not be delayed.Figures 7(a)-(d

Effects of Variation of β on the Cardiac
Electrical Wave Propagation In order to study the effect of variation of β on the normal cardiac electrical wave presented in Figure 5, the value of β was changed to 0.3 and 1.1 respectively from the considered value 0.7.It was observed during the simulation that variation in the value of β did not impact any significant effect on the signal at 0.3 and 1.1 respectively.The cardiac electrical signal for β equals 0.3 and 1.1 respectively are shown in Figures 8 and 9

Effects of Variation of γ on the Cardiac
Electrical Wave Propagation The effect of variation in the value of γ is similar to that of variation in the value of β.When the value of γ was changed from 0.8 to 0.4 and 1.2 respectively, it was observed that variation in the value of γ did not impact any significant effect on the electrical signal as presented in Figures 10 and 11.
Summarily, the results of these simulations showed that of the parameters whose effects of variation on the cardiac electrical wave propagation were studied, it was only decrease in є 1 , є 2 that impacted significantly on the cardiac wave propagation while increase in є 1 , є 2 , decrease or increase in β and γ merely had subtle effect on normal cardiac electrical wave propagation.

CONCLUSION
In this work, dynamical modelling of cardiac electrical activity using bidomain approach was presented.Apart from the fact that this work has been able to provide some insights into the electrical behaviour of human heart, revealing the nature of the electrical wave propagation pattern in the normal cardiac tissue, it also showed that cardiac electrical signal could be significantly affected by variation in the values of some of the basic parameters characterizing the adopted FitzHugh-Nagumo's ionic model coupling the bidomain model.The simulation results showed the excitation pattern in 2-D.Time dilation effect of the cardiac electrical wave resulting from decrease in value of є 1 and є 2 characterizing the ionic model is a sign of potential cardiac arrhythmias and the overall resulting effect is a delayed repolarization, which if persists for a long time may result in sudden cardiac death.The obtained results in this work are very useful in studying the characteristic properties of action potential (the time characteristics of transmembrane potential) as it propagates through the cardiac tissue and in effect detect any electrical wave abnormalities in the cardiac tissue.Despite the computational difficulty of numerical simulation of bidomain equations, we were still able to achieve the major objectives of this work by restricting ourselves to 2-D and FitzHugh-Nagumo's ionic model.In our future work, we will investigate potential multiscale and other numerical methods [26,27] for efficient simulation of bidomain equations in 3-D using other ionic models.

Figure 1 .
Figure 1.Schematic model of the bidomain space; the intracellular and extracellular domains are separated by cell membrane [18].
, and t  , time step size.

Figure 3 .
Figure 3. Network of resistors connecting the nodes.
Input parameters Sum matrices G i and G e Invert the sum of matrices G i and G e Construct the right hand side of Equaion 20 Solve for e  at timestep n in Equation 20 Solve for V m at timestep n+1 in Equaion 18 Stop Next timestep Impose initial conditions on V m and w Define arrays Construct admittance matrices G i and G e Keep matrix G i Plot results V m against Time Solve for w at timestep n+1in Equaion 19

Figure 4 .
Figure 4. Flow chart for the bidomain implementation code.
) to (d) respectively with depolarization, partial repolarization, plateau, repolarization and resting membrane potential phases identified using the values 0, 1, 2, 3, and 4 in Figure 5(d) for clarity.This electrical signal produced is called the action potential, with the highest period observed in this work around 600 ms.Figures 5(a) to (d) are typically of the same wave pattern, consistent with the theoretical standard and the experimental findings from other researchers [3,24, 25].