PDF
Electrical Signals in a Heart
This example was kindly provided by Prof. Simonetta Filippi and Dr Christian Cherubini from Università Campus Biomedico di Roma, Italy.
Introduction
Modeling the electrical activity in cardiac tissue is an important step in understanding the patterns of contractions and dilations in the heart. The heart produces rhythmic electrical pulses, initiated from a point known as the sinus node. The electrical pulses, in turn, trigger the mechanical contractions of the muscle. In a healthy heart these electrical pulses are damped out, but a number of heart conditions involve an elevated risk of re-entry of the signals. This means that the normal steady pulse is disturbed, a severe and acute condition often referred to as arrhythmia.
In this example, different aspects of electrical signal propagation in cardiac tissue are studied using the FitzHugh–Nagumo equations and the Ginzburg–Landau equations, both of which are solved on the same geometry. Interesting patterns emerging from these types of models are, for example, spiral waves, which, in the context of cardiac electrical signals, can produce effects similar to those observed in cardiac arrhythmia.
Excitable Media and the FitzHugh–Nagumo Equations
It has been shown that many important characteristics of electrical signal propagation in cardiac tissue can be reproduced by a class of equations which describe excitable media, that is, materials consisting of elementary segments or cells with the following basic characteristics:
•
Well-defined rest state
•
Threshold for excitation
•
A diffusive-type coupling to its nearest neighbors
Excitable media is a rather general concept, which is useful for modeling of a number of different phenomena, including nerve pulses, the spreading of forest fires, and certain types of chemical reactions, in addition to the electrical signals in cardiac tissue. One of the most important qualitative characteristics displayed by excitable media, and equally a common denominator between the diversity of phenomena mentioned above, is the almost immediate damping out of signals below a certain threshold. On the other hand, signals exceeding this threshold propagate without damping.
The heart works by passing ionic current inside the muscle, thus triggering the rhythmic contractions that pump blood in and out. The ions move through small pores or gates in the cellular membrane, which can be either open (excitation state) or closed (rest state).
In nerve cells and cardiac cells the three abstract characteristics of excitable media are manifested as:
•
Rest cell membrane potential
•
Threshold for opening the ionic gates in the cellular membrane
•
The diffusive spreading of the electrical signals
The state of the membrane gates is random on a microscopic scale, but the probability of a given state can be modeled as a continuous function of the voltage, thus allowing an averaged macroscopic continuum description of the current flow.
The FitzHugh–Nagumo equations for excitable media describe the simplest physiological model with two variables, an activator and an inhibitor. In these heart models the activator variable corresponds to the electric potential, and the inhibitor is a variable that describes the voltage-dependent probability of the pores in the membrane being open and ready to transmit ionic current.
Chaotic Dynamics and the Ginzburg–Landau Equations
The Ginzburg–Landau equations provide a relatively simple way of modeling some aspects of the transition, displayed by many dynamical systems under the influence of strong external stimulus, from periodic oscillatory behavior into a chaotic state with gradually increasing amplitude of oscillations and decreasing periodicity.
Although their first use was to describe the theory of superconductivity, the Ginzburg–Landau equations are also generic in their nature (as are the FitzHugh–Nagumo equations), and examples of dynamical systems that you can model successfully using these equations are:
•
The formation of vortices behind a slender obstacle in transversal fluid flow
•
Oscillating chemical reactions of the Belousov–Zhabotinsky type
In this example, the Ginzburg–Landau equations simulate the dynamics of the spiral waves in excitable media.
Model Definition
The geometry here is a simplified 3D model of a heart with two chambers, represented with semispherical cavities1.
Figure 1: Model geometry.
The FitzHugh–Nagumo Equations
The equations are the following:
Here u1 is an action potential (the activator variable), and u2 is a gate variable (the inhibitor variable). The parameter α represents the threshold for excitation, ε represents the excitability, and β, γ, and δ are parameters that affect the rest state and dynamics of the system.
The boundary conditions for u1 are insulating, using the assumption that no current is flowing into or out of the heart. The initial condition defines an initial potential distribution u1 where one quadrant of the heart is at a constant, elevated potential V0, while the rest remains at zero. The adjacent quadrant has instead an elevated value ν0 for the inhibitor u2. It is convenient to implement this initial distribution using the following logical expressions, where TRUE evaluates to 1 and FALSE to 0:
Here d is equal to 10−5, and it is included in the expressions to shift the elevated potential slightly off the main axes.
The Ginzburg–Landau Equations
The Ginzburg–Landau equations are:
The two variables v1 and v2 are the activator and inhibitor, respectively. The constants c1 and c3 are parameters reflecting the properties of the material. These constants also determine the existence and nature of the stable solutions.
As in the previous model, the boundary conditions are kept insulating. The initial condition, which gives a smooth transition step near z = 0, are the following:
Notes About the COMSOL Implementation
The simplified geometry is quite straightforward to create using the drawing tools in COMSOL Multiphysics. The FitzHugh–Nagumo and Ginzburg–Landau equations are also readily entered in one of the PDE interfaces2.
It is important to note that these equations are strongly nonlinear. It is therefore necessary (especially in full 3D models like these) to use a much finer mesh or use higher element order than in this example to get results with some degree of reliability for the time intervals of interest. This is particularly important in solving the Ginzburg–Landau equations, which describe inherently chaotic phenomena. They are highly sensitive to perturbations in the initial value and similarly to numerical errors during the course of the time-dependent solution. We recommend the use of the fourth-order Hermite element for the Ginzburg–Landau equation.
For the reasons above, the results presented here are only intended as a first rough estimate of the qualitative behavior that you can expect the system to show under a given stimulus. Consequently, higher-order elements, finer meshing, and smaller relative and absolute time dependent tolerances clearly give quantitatively more correct simulation results. The implementation of these improvements leads to longer computational time than that for the rough model described here which takes a few minutes to compute on a standard PC. When attempting these types of large models we strongly recommend the use of 64-bit platforms.
To avoid warnings about inconsistent units, turn off unit support.
Results and Discussion
The FitzHugh–Nagumo Equations
The plots in Figure 2 below show the action potential u1. To visualize the solution on the inside, a quarter of the outside shell of the heart and one of the chamber surfaces are suppressed in the plot.
The parameters used in the model along with the initial pulse lead to a reentrant wave which travels around the tissue without damping in a characteristic spiral pattern.
Figure 2: Solution to the FitzHugh–Nagumo equations at times t = 120 s (top) and t = 500 s (bottom).
The Ginzburg–Landau Equations
Figure 3 below shows the species v1 at different times.
Figure 3: Solution to the Ginzburg–Landau equations at times t = 75 s (top) and t = 45 s (bottom).
The equation parameters and initial condition used here lead the diffusing species (v1) to display characteristic spiral patterns with growing complexity over time.
References
1. F.H. Fenton, E.M. Cherry, H.M. Hastings, and S.J. Evans, “Real-time computer simulations of excitable media: JAVA as a scientific language and as a wrapper for C and FORTRAN programs,” BioSystems, vol. 64, pp. 73–96, 2002.
2. Y. Kuramoto, Chemical Oscillations, Waves and Turbulence, Dover Publications, 2003.
3. J. Keener and J. Sneyd, Mathematical Physiology, Springer, 1998.
Application Library path: COMSOL_Multiphysics/Equation_Based/heart_electrical
Modeling Instructions
From the File menu, choose New.
New
In the New window, click  Model Wizard.
Model Wizard
1
In the Model Wizard window, click  3D.
2
In the Select Physics tree, select Mathematics > PDE Interfaces > General Form PDE (g).
3
Click Add.
4
In the Number of dependent variables text field, type 2.
5
Click  Study.
6
In the Select Study tree, select General Studies > Time Dependent.
7
Click  Done.
Root
1
In the Model Builder window, click the root node.
2
In the root node’s Settings window, locate the Unit System section.
3
From the Unit system list, choose None.
The equations in this model are given in dimensionless form. Turning off unit support avoids warnings about inconsistent use of units.
Geometry 1
Follow the instructions below to recreate the geometry illustrated in Figure 1.
Sphere 1 (sph1)
1
In the Geometry toolbar, click  Sphere.
2
In the Settings window for Sphere, locate the Size section.
3
In the Radius text field, type 54.
4
Click  Build Selected.
Ellipsoid 1 (elp1)
1
In the Geometry toolbar, click  More Primitives and choose Ellipsoid.
2
In the Settings window for Ellipsoid, locate the Size and Shape section.
3
In the a-semiaxis text field, type 54.
4
In the b-semiaxis text field, type 54.
5
In the c-semiaxis text field, type 70.
6
Click  Build Selected.
Next create the egg-shaped solid by fusing the top half of the sphere with the bottom half of the ellipsoid. Use blocks to remove the bottom half of sphere and the top half of ellipsoid.
Block 1 (blk1)
1
In the Geometry toolbar, click  Block.
2
In the Settings window for Block, locate the Size and Shape section.
3
In the Width text field, type 110.
4
In the Depth text field, type 110.
5
In the Height text field, type 110.
6
Locate the Position section. From the Base list, choose Center.
7
In the z text field, type 55.
8
Click  Build Selected.
Difference 1 (dif1)
1
In the Geometry toolbar, click  Booleans and Partitions and choose Difference.
2
Select the object elp1 only.
3
In the Settings window for Difference, locate the Difference section.
4
Click to select the  Activate Selection toggle button for Objects to subtract.
5
Select the object blk1 only.
6
Click  Build Selected.
Block 2 (blk2)
1
In the Geometry toolbar, click  Block.
2
In the Settings window for Block, locate the Size and Shape section.
3
In the Width text field, type 110.
4
In the Depth text field, type 110.
5
In the Height text field, type 110.
6
Locate the Position section. From the Base list, choose Center.
7
In the z text field, type -55.
8
Click  Build Selected.
Difference 2 (dif2)
1
In the Geometry toolbar, click  Booleans and Partitions and choose Difference.
2
Select the object sph1 only.
3
In the Settings window for Difference, locate the Difference section.
4
Click to select the  Activate Selection toggle button for Objects to subtract.
5
Select the object blk2 only.
6
Click  Build Selected.
Union 1 (uni1)
1
In the Geometry toolbar, click  Booleans and Partitions and choose Union.
2
Click in the Graphics window and then press Ctrl+A to select both objects.
3
In the Settings window for Union, click  Build All Objects.
Now create the cavity inside the heart.
Scale 1 (sca1)
1
In the Geometry toolbar, click  Transforms and choose Scale.
2
In the Settings window for Scale, locate the Input section.
3
Select the Keep input objects checkbox.
4
Select the object uni1 only.
5
Locate the Scale Factor section. In the Factor text field, type 2/3.
6
Click  Build Selected.
7
Click the  Transparency button in the Graphics toolbar.
Difference 3 (dif3)
1
In the Geometry toolbar, click  Booleans and Partitions and choose Difference.
2
Select the object uni1 only.
3
In the Settings window for Difference, locate the Difference section.
4
Click to select the  Activate Selection toggle button for Objects to subtract.
5
Select the object sca1 only.
6
Click  Build Selected.
The next step is to create the wall separating the two chambers.
Cylinder 1 (cyl1)
1
In the Geometry toolbar, click  Cylinder.
2
In the Settings window for Cylinder, locate the Size and Shape section.
3
In the Radius text field, type 45.
4
In the Height text field, type 10.
5
Locate the Position section. In the x text field, type -5.
6
In the z text field, type -5.
7
Locate the Axis section. From the Axis type list, choose Cartesian.
8
In the x text field, type 1.
9
In the z text field, type 0.
10
Click  Build Selected.
Union 2 (uni2)
1
In the Geometry toolbar, click  Booleans and Partitions and choose Union.
2
Click in the Graphics window and then press Ctrl+A to select both objects.
3
In the Settings window for Union, locate the Union section.
4
Clear the Keep interior boundaries checkbox.
5
Click  Build All Objects.
This completes the geometry.
Global Definitions
Parameters 1
1
In the Model Builder window, under Global Definitions click Parameters 1.
2
In the Settings window for Parameters, locate the Parameters section.
3
In the table, enter the following settings:
Name
Expression
Value
Description
alpha
0.1
0.1
Excitation threshold
epsilon
0.01
0.01
Excitability
beta
0.5
0.5
System parameter
gamma
1
1
System parameter
delta
0
0
System parameter
V0
1
1
Elevated potential value
nu0
0.3
0.3
Elevated inhibitor value
d
1e-5
1E-5
Off-axis shift distance
Use the default Neumann boundary condition at all the boundaries.
General Form PDE (g)
General Form PDE 1
1
In the Model Builder window, under Component 1 (comp1) > General Form PDE (g) click General Form PDE 1.
2
In the Settings window for General Form PDE, locate the Conservative Flux section.
3
Specify the second Γ vector as
0
x
0
y
0
z
4
Locate the Source Term section. In the f text-field array, type (alpha-u1)*(u1-1)*u1-u2 on the first row.
5
In the f text-field array, type epsilon*(beta*u1-gamma*u2-delta) on the second row.
Initial Values 1
1
In the Model Builder window, click Initial Values 1.
2
In the Settings window for Initial Values, locate the Initial Values section.
3
In the u1 text field, type V0*((x+d)>0)*((z+d)>0).
4
In the u2 text field, type nu0*((-x+d)>0)*((z+d)>0).
Mesh 1
In the Model Builder window, under Component 1 (comp1) right-click Mesh 1 and choose Build All.
Study 1
Step 1: Time Dependent
1
In the Model Builder window, under Study 1 click Step 1: Time Dependent.
2
In the Settings window for Time Dependent, locate the Study Settings section.
3
In the Output times text field, type range(0,5,500).
4
In the Study toolbar, click  Compute.
Results
General Form PDE
The default plot shows a slice plot of the dependent variable. To create plots given in Figure 2, you need to suppress some of the boundaries.
Study 1/Solution 1 (sol1)
In the Model Builder window, expand the Results > Datasets node, then click Study 1/Solution 1 (sol1).
Selection
1
In the Results toolbar, click  Attributes and choose Selection.
2
In the Settings window for Selection, locate the Geometric Entity Selection section.
3
From the Geometric entity level list, choose Boundary.
4
Click  Paste Selection.
5
In the Paste Selection dialog, type 1 3 10-18 in the Selection text field.
6
Click OK.
7
In the Settings window for Selection, locate the Geometric Entity Selection section.
8
Click  Create Selection.
9
In the Create Selection dialog, click OK.
3D Plot Group 2
1
In the Results toolbar, click  3D Plot Group.
2
In the Settings window for 3D Plot Group, locate the Plot Settings section.
3
Clear the Plot dataset edges checkbox.
Surface 1
1
Right-click 3D Plot Group 2 and choose Surface.
2
In the 3D Plot Group 2 toolbar, click  Plot.
3
Click the  Transparency button in the Graphics toolbar.
3D Plot Group 2
1
In the Model Builder window, click 3D Plot Group 2.
2
In the Settings window for 3D Plot Group, locate the Data section.
3
From the Time (s) list, choose 120.
4
In the 3D Plot Group 2 toolbar, click  Plot.
This completes the modeling of the FitzHugh-Nagumo equations. To model the Complex Landau-Ginzburg equations use the following set of instructions.
Add Physics
1
In the Home toolbar, click  Add Physics to open the Add Physics window.
2
Go to the Add Physics window.
3
In the tree, select Mathematics > PDE Interfaces > General Form PDE (g).
4
Click to expand the Dependent Variables section. In the Field name (1) text field, type v.
5
In the Number of dependent variables text field, type 2.
6
In the Dependent variables (1) table, enter the following settings:
v1
7
Find the Physics interfaces in study subsection. In the table, clear the Solve checkbox for Study 1.
8
Click the Add to Component 1 button in the window toolbar.
9
In the Home toolbar, click  Add Physics to close the Add Physics window.
Add Study
1
In the Home toolbar, click  Add Study to open the Add Study window.
2
Go to the Add Study window.
3
Find the Physics interfaces in study subsection. In the table, clear the Solve checkbox for General Form PDE (g).
4
Find the Studies subsection. In the Select Study tree, select General Studies > Time Dependent.
5
Click the Add Study button in the window toolbar.
6
In the Home toolbar, click  Add Study to close the Add Study window.
Global Definitions
Parameters 1
1
In the Settings window for Parameters, locate the Parameters section.
2
In the table, enter the following settings:
Name
Expression
Value
Description
c1
2
2
PDE parameter
c3
-0.2
-0.2
PDE parameter
General Form PDE 2 (g2)
The accuracy of the solution can be improved by using cubic Hermite elements instead of the default quadratic Lagrange elements.
1
In the Model Builder window, under Component 1 (comp1) click General Form PDE 2 (g2).
2
In the Settings window for General Form PDE, click to expand the Discretization section.
3
From the Shape function type list, choose Hermite.
General Form PDE 1
Use the default Neumann boundary condition at all the boundaries also in this case.
1
In the Model Builder window, under Component 1 (comp1) > General Form PDE 2 (g2) click General Form PDE 1.
2
In the Settings window for General Form PDE, locate the Conservative Flux section.
3
Specify the first Γ vector as
-v1x+c1*v2x
x
-v1y+c1*v2y
y
-v1z+c1*v2z
z
4
Specify the second Γ vector as
-c1*v1x-v2x
x
-c1*v1y-v2y
y
-c1*v1z-v2z
z
5
Locate the Source Term section. In the f text-field array, type v1-(v1-c3*v2)*(v1^2+v2^2) on the first row.
6
In the f text-field array, type v2-(c3*v1+v2)*(v1^2+v2^2) on the second row.
Initial Values 1
1
In the Model Builder window, click Initial Values 1.
2
In the Settings window for Initial Values, locate the Initial Values section.
3
In the v1 text field, type tanh(z[1/m]).
4
In the v2 text field, type -tanh(z[1/m]).
Tighten the relative and absolute tolerances by a factor of 10.
Study 2
Step 1: Time Dependent
1
In the Model Builder window, under Study 2 click Step 1: Time Dependent.
2
In the Settings window for Time Dependent, locate the Study Settings section.
3
In the Output times text field, type range(0,5,75).
4
From the Tolerance list, choose User controlled.
5
In the Relative tolerance text field, type 0.001.
Solution 2 (sol2)
1
In the Study toolbar, click  Show Default Solver.
2
In the Model Builder window, expand the Solution 2 (sol2) node, then click Time-Dependent Solver 1.
3
In the Settings window for Time-Dependent Solver, click to expand the Absolute Tolerance section.
4
From the Tolerance method list, choose Manual.
5
In the Absolute tolerance text field, type 0.0001.
6
In the Study toolbar, click  Compute.
Results
Follow the instructions to create plots given in Figure 3.
Study 2/Solution 2 (sol2)
In the Model Builder window, under Results > Datasets click Study 2/Solution 2 (sol2).
Selection
1
In the Results toolbar, click  Attributes and choose Selection.
2
In the Settings window for Selection, locate the Geometric Entity Selection section.
3
From the Geometric entity level list, choose Boundary.
4
From the Selection list, choose Explicit 1.
3D Plot Group 4
1
In the Results toolbar, click  3D Plot Group.
2
In the Settings window for 3D Plot Group, locate the Data section.
3
From the Dataset list, choose Study 2/Solution 2 (sol2).
4
Locate the Plot Settings section. Clear the Plot dataset edges checkbox.
Surface 1
1
Right-click 3D Plot Group 4 and choose Surface.
2
In the Settings window for Surface, click Replace Expression in the upper-right corner of the Expression section. From the menu, choose Component 1 (comp1) > General Form PDE 2 > v1 - Dependent variable v1 - 1.
3
In the 3D Plot Group 4 toolbar, click  Plot.
3D Plot Group 4
1
In the Model Builder window, click 3D Plot Group 4.
2
In the Settings window for 3D Plot Group, locate the Data section.
3
From the Time (s) list, choose 45.
4
In the 3D Plot Group 4 toolbar, click  Plot.
1
Note that it is possible to import more realistic human/animal heart geometries, with anisotropies and inhomogeneities as well as proper dimensions, into COMSOL Multiphysics using the CAD import capabilities.

2
With this stage completed, it would also be straightforward to replace the rather simple equations with a physiologically more realistic mathematical model.