Simulations of Magnetocardiographic Signals Using Realistic Geometry Models of the Heart and Torso

Although the first measurement of the cardiac magnetic field was reported almost half a century ago mag-netocardiography (MCG) is not yet widely used as a clinical diagnostic technique. With the development of a new generation of magnetoelectric sensors it is believed that MCG will become widely accepted in the clinical diagnosis. Our goal is to build a computer-based tool for medical diagnosis and to use it for the clarification of open electro-physiological questions. Here we present results from modelling of the cardiac electrical activity and computation of the generated magnetic field. For the simulations we use MRT-based anatomical models of the human atria and ventricles where the shape of the action potential is determined by ionic currents passing through the cardiac cell membranes. The monodomain reaction-diffusion equation is chosen for the description of the heart's electrical activity. This equation is solved for the transmembrane voltage which is in turn used to calculate current densities at discrete time instants. In subsequent simulations these current densities represent primary sources of magnetostatic fields arising from a volume conduction problem. In these simulations the heart is placed in a realistic torso model where the lungs are also considered. Both, the volume conduction problem as well as the reaction-diffusion problem are modelled using Finite-Element techniques.


Introduction
The human body can be seen as a volume conductor where the electrical activity of the heart gives rise to small electric currents in the body with associated potential distributions.The electrocardiography (ECG) records the variation in time of the electrical potential produced by the heart at the surface of the body.The graphical representation of the recorded potentials depicts the well known electrocardiogram.However, the same currents in the body give also rise to weak magnetic fields.The magnetocardiography (MCG) records the variation in time of the magnetic fields produced by the heart outside the body and the graphical representation of the recorded magnetic signals is called magnetocardiogram.Whereas the ECG is a widely used technique in the clinical diagnosis, the MCG is not yet.The reasoning for this may lay in the fact that currently used magnetic sensors (Superconducting QUantum Interference Devices -SQUIDs) require extremely low temperatures as well as magnetically shielded environments.Therefore, electrocardiography (ECG) is considerably more attractive for clinicians, though MCG and ECG supposedly supply disjunctive information and MCG offers several advantages (e.g.no electrodes needed, magnetic transparency of tissue).Efforts are nowadays directed to the overcoming of the above mentioned disadvantages of MCG by using magnetoelectric composite materials (Jahns et al., 2011) with giant magnetoelectric effect.In the same time computational tools for modelling the electrical activity of the heart as well as for solving the inverse problem of magnetocardiography are developed as a framework for studies on open electrophysiological questions such as the non-response of some patients to bi-ventricular stimulation.
In this paper we present a modelling of the electrical activity of the heart at both, cellular and organ level and a method to compute the associated magnetic field produced by a healthy heart.

Biological and electro-physiological background
The heart is a hollow muscular organ being structured on four chambers i.e. two atria and two ventricles.The function of the heart is to pump blood.This function is achieved by contractions of the muscles composing the atria and ventricles (working myocardium).The contraction of the cardiac muscles is electrically controlled by the excitation-conduction system of the heart.This system is made of cells that are able to fire impulses -a property called autorhythmicity.
Published by Copernicus Publications on behalf of the URSI Landesausschuss in der Bundesrepublik Deutschland e.V.  gram and the repolarization of the ventricles that determines the T-wave.

linkenbusch: Simulations of Magnetocardiographic Signals
The atrial repolarization takes place during the QRScomplex but it is masked by the much stronger process of ventricular depolarization.

Mathematical modelling at the cellular level
As already mentioned in Sect. 2 the electrical activity of the cardiac myocyte is mainly the result of the ionic currents that are crossing the cell membrane through ion channels and transporters.These currents determine the shape of the action potential of a cell.The available cellular ionic models at present can be divided into three classes (Clayton and Pan-  (Malmivuo and Plonsey, 1995).
Electrical impulses result from ionic activity at the cellular level.At rest, the cell membrane of the cardiac muscle cells (myocytes) composing the atria and ventricles, is electrically polarized in the sense that the inner side of the membrane is more electronegative than the outer side of the membrane.This produces the so called transmembrane voltage.When a strong enough stimulus is applied to the membrane, certain protein channels will open letting particular ionic species to cross the membrane along their electro-chemical gradients.This depolarization of the membrane is followed by a repolarization which brings the myocyte in the resting state.In contrast to the atrial and ventricular myocytes (of the working myocardium), the cells belonging to the excitationconduction system of the heart do not need an external stimulus for the depolarization because they are able to depolarize themselves initiating thus electrical impulses (pacemaker cells).The transmembrane voltage propagates from a myocyte to the neighboring myocytes in the form of an action potential.The spread of the action potential is supported by some rapid intercellular channels called gap junctions which couple the interiors of adjacent cells.Thus, at the organ level the excitation propagates in the form of isochronal surfaces also called wave-fronts.Electrically active cells of the heart are located within these wave-fronts and they are usually associated to current dipoles.It is these moving current dipoles with their accompanying conduction currents that gives rise to the magnetic field of the heart.
In contrast to the atrial and ventricular myocytes (of the working myocardium), the cells belonging to the excitationconduction system of the heart do not need an external stimulus for the depolarization because they are able to depolarize themselves initiating thus electrical impulses (pacemaker cells).The transmembrane voltage propagates from a myocyte to the neighboring myocytes in the form of an action potential.The spread of the action potential is supported by some rapid intercellular channels called gap junctions which couple the interiors of adjacent cells.Thus, at the organ level the excitation propagates in the form of isochronal surfaces also called wave-fronts.Electrically active cells of the heart are located within these wave-fronts and they are usually associated to current dipoles.It is these moving current dipoles with their accompanying conduction currents that give rise to the magnetic field of the heart.

The normal electrocardiogram
The excitation-conduction system of the heart (see Fig. 1) consists of the sinoatrial (SA) node, the atrioventricular (AV) node, the common bundle also called the bundle of His, the left and right bundle branches and the Purkinje fibres.
Under normal conditions, the SA node generates action potentials at the rate of about 70 pulses/min (Malmivuo and Plonsey , 1995) which spread through the right and left atria depolarizing them and determining the occurrence of the so called P-wave in the recorded electrocardiogram (see Fig. 2).The only conduction path between the atria and the ventricles is provided by the AV node which has a low conduction velocity and consequently introduces a delay in the propagation of the excitation that corresponds to the PQ-segment of the electrocardiogram.Once the AV node has been passed, the excitation propagates through a high conduction-velocity system (the His-Purkinje system) to the inner sides of the ventricular wall.The ventricular depolarization determines the QRS complex on the electrocardiogram.After the complete ventricular depolarization, it follows a plateau phase which corresponds to the ST-segment on the electrocardio-Fig.1.Schematic representation of a section through the human heart showing the excitation-conduction system.Picture adapted from (Malmivuo and Plonsey , 1995).The atrial repolarization takes place during the QRScomplex but it is masked by the much stronger process of ventricular depolarization.

Mathematical modelling at the cellular level
As already mentioned in Sect. 2 the electrical activity of the 125 cardiac myocyte is mainly the result of the ionic currents that are crossing the cell membrane through ion channels and transporters.These currents determine the shape of the action potential of a cell.The available cellular ionic models at present can be divided into three classes (Clayton and Pan- simplified two-variable models that reduce the behavior of all the ionic currents to two variables that describe excitation and recovery; biophysically detailed models that use experimental in-135 formation about the voltage and time dependence of ion channel conductances (first generation) as well as intracellular ionic concentrations (second generation) in order to reproduce the action potential;

The normal electrocardiogram
The excitation-conduction system of the heart (see Fig. 1) consists of the sinoatrial (SA) node, the atrioventricular (AV) node, the common bundle also called the bundle of His, the left and right bundle branches and the Purkinje fibres.
Under normal conditions, the SA node generates action potentials at the rate of about 70 pulses/min (Malmivuo and Plonsey, 1995) which spread through the right and left atria depolarizing them and determining the occurrence of the so called P-wave in the recorded electrocardiogram (see Fig. 2).The only conduction path between the atria and the ventricles is provided by the AV node which has a low conduction velocity and consequently introduces a delay in the propagation of the excitation that corresponds to the PQ-segment of the electrocardiogram.Once the AV node has been passed, the excitation propagates through a high conduction-velocity system (the His-Purkinje system) to the inner sides of the ventricular wall.The ventricular depolarization determines the QRS complex on the electrocardiogram.After the complete ventricular depolarization, it follows a plateau phase which corresponds to the ST-segment on the electrocardiogram and the repolarization of the ventricles that determines the T-wave.
The atrial repolarization takes place during the QRScomplex but it is masked by the much stronger process of ventricular depolarization.

Mathematical modelling at the cellular level
As already mentioned in Sect. 2 the electrical activity of the cardiac myocyte is mainly the result of the ionic currents that are crossing the cell membrane through ion channels and transporters.These currents determine the shape of the action potential of a cell.The available cellular ionic models at present can be divided into three classes (Clayton and Panfilov, 2008): -simplified two-variable models that reduce the behavior of all the ionic currents to two variables that describe excitation and recovery; s of one or two phases.The bi-domain model conthe cardiac muscle to be a two-phase medium conof intracellular and extracellular spaces separated at oint in space by the cellular membrane.Although this ry detailed model, it is also computationally expenhe simpler monodomain model is based on a single description of the cardiac muscle but it is still able oduce many of the experimentally observed phenomth a much higher numerical efficiency.The simulation oblem with the monodomain model can be by a facten or more faster than the same problem simulated e bi-domain model (Clayton and Panfilov , 2008).In rk we consider that the action potential propagates usion in the monodomain model of the cardiac tisth local excitation of the membrane voltage.The so monodomain reaction-diffusion model consists of a lic partial differential equation (PDE) coupled to a of ODEs (Whiteley , 2006): [ ∂V ] channel conductances (first generation) as well as intracellular ionic concentrations (second generation) in order to reproduce the action potential; -reduced cardiac models that simplify the underlying electrophysiology to a minimum number of currents for which the activation and inactivation is controlled by some empirical variables.
Choosing the appropriate cellular model for a certain problem is a matter of balancing its benefits and limitations.For this work we have chosen the Luo-Rudy phase 1 model of mammalian myocytes (Luo and Rudy, 1991).This cellular model belongs to the first generation of the biophysically detailed models and it is considered that this model offers a good balance between electrophysiological details and computational efficiency.Within this model the Hodgkin-Huxley type of formalism is used to describe the action potential.
The rate of change of the membrane voltage V m is given by: where C m is the membrane capacitance, I st is an applied stimulus current and I ion is the sum of six ionic currents: a fast sodium current (I Na ), a slow inward current (I si ), a timedependent potassium current (I k ), a time-independent potassium current (I k1 ), a plateau potassium current (I kp ) and a time-independent background current (I b ).In general, each current I S carried by an ionic species "S" through an ion channel is expressed as: with G S being the maximal conductance of the ion channel,

Method
The procedure adopted by us for computing the magneti field produced by the heart consists of three subsequent step where the output of one step is used as an input for the nex

225
step.The steps are also differentiated by the involved tim and space scales.
In the first step we simulate the propagation of electrica excitation in the isolated heart shown in the right-hand sid of Fig. 4. As a result we obtain the distribution of the trans 230 membrane voltage in the heart at discrete time instants during a complete heart beat.
In the second step a volume conduction problem is solved where the driving force is represented by impressed curren densities in the heart which are computed using the trans 235 membrane voltage resulted at each time instant considered in the simulations of the previous (first) step.The computa tional domain consists of the heart and lungs inside the torso as shown on the left-hand side of Fig. 4. The output of thes simulations is represented by the total (impressed and con 240 duction) current density in the entire torso.
In the third step the magnetic field outside the body pro duced by the total current density in the torso, is computed as a magnetostatic problem at each time instant considered in the previous (second) step.[0,1] and describing the activation, inactivation and recovery of the ion channel and V S being the membrane voltage of the species "S" alone.The gating variables are obtained by solving a system of eight nonlinear ordinary differential equations (ODEs) of the general form: where x S∞ is the steady-state value of a single gating variable x S and τ S is the time constant for the return of the variable x to its steady-state value after a voltage perturbation.The voltage dependencies of x S∞ and τ S are experimentally determined.
The shape of the action potential as given by the Luo-Rudy phase 1 cellular model and used in this work is shown in Fig. 3.

Mathematical modelling at the tissue level
At the tissue level the granular structure may be neglected and the cardiac muscle may be modelled as a functional syncytium in which the membrane voltage propagates smoothly.Under this assumption the action potential can be considered to propagate by diffusion in an excitable medium that consists of one or two phases.The bi-domain model considers the cardiac muscle to be a two-phase medium consisting of intracellular and extracellular spaces separated at each point in space by the cellular membrane.Although this is a very detailed model, it is also computationally expensive.The simpler monodomain model is based on a single phase description of the cardiac muscle but it is still able to reproduce many of the experimentally observed phenomena with a much higher numerical efficiency.The simulation of a problem with the monodomain model can be by a factor of ten or more faster than the same problem simulated  meshes by using the "Tetgen" -a freely available software code for research and non-commercial purposes (Si, 2006).

255
The atria are discretized with 1.23 millions mesh nodes with an average edge length of 0.8 mm and the ventricles contain 1.43 millions mesh nodes with an average edge length of 0.95 mm.
The simulations are based on the monodomain reaction- with the bi-domain model (Clayton and Panfilov, 2008).In this work we consider that the action potential propagates by diffusion in the monodomain model of the cardiac tissue with local excitation of the membrane voltage.The so called monodomain reaction-diffusion model consists of a parabolic partial differential equation (PDE) coupled to a system of ODEs (Whiteley, 2006): where σ (Sm −1 ) is the electrical conductivity tensor, V m (V) is the transmembrane voltage, β (m −1 ) is the surface to volume ratio of the cardiac cells, C m (Fm −2 ) is the specific membrane capacitance, I ion (Am −2 ) is the current flow through unit area of the cell membrane and I st (Am −2 ) is an applied stimulation current.The net membrane current is given by the cellular ionic models such as the one described in Sect. 4.

Anatomical models
Anatomical models of the heart, lungs and torso used in this work are those provided with the freely available software code for simulations of electrocardiograms called "ECGsim" (van Oosterom, 2004).The models are based on magnetic resonance tomography (MRT) data of a healthy young male of 22 yr.They are available in the form of triangulated surfaces as shown in Fig. 4.

Method
The procedure adopted by us for computing the magnetic field produced by the heart consists of three subsequent steps where the output of one step is used as an input for the next C.V. Motrescu and L. Klinkenbusch: Simula repolarization of the atria and ventricles are a the form of short films.

Computation of currents in the body
The electrical excitation in the heart gives step.The steps are also differentiated by the involved time and space scales.
In the first step we simulate the propagation of electrical excitation in the isolated heart shown in the right-hand side of Fig. 4. As a result we obtain the distribution of the transmembrane voltage in the heart at discrete time instants during a complete heart beat.
In the second step a volume conduction problem is solved where the driving force is represented by impressed current densities in the heart which are computed using the transmembrane voltage resulted at each time instant considered in the simulations of the previous (first) step.The computational domain consists of the heart and lungs inside the torso as shown on the left-hand side of Fig. 4. The output of these simulations is represented by the total (impressed and conduction) current density in the entire torso.
In the third step the magnetic field outside the body produced by the total current density in the torso, is computed as a magnetostatic problem at each time instant considered in the previous (second) step.
Thus slowly varying magnetocardiographic signals with time are obtained from a sequence of stationary computations.

Simulation of the excitation spreading in the heart
The spread of excitation in the heart is simulated with Finite-Element methods.For this purpose the surface meshes of the atria and ventricles composing the heart shown on the righthand side of the Fig. 4 are transformed to tetrahedral volume meshes by using the "Tetgen" -a freely available software code for research and non-commercial purposes (Si, 2006).The atria are discretized with 1.23 millions mesh nodes with an average edge length of 0.8 mm and the ventricles contain 1.43 millions mesh nodes with an average edge length of 0.95 mm.The simulations are based on the monodomain reactiondiffusion model described by Eq. ( 4).Each node of the computational domain contains the Luo-Rudy ionic cellular model described in Sect. 4. This cellular model with constant action potential duration (APD) is used for both, atria and ventricles.The cardiac muscles are assumed to be homogeneous and isotropic with the scalar electrical conductivity σ = 0.25 Sm −1 .Other parameter values used in the simulation are: C m = 1 µFcm −2 and β = 1400 cm −1 .
The simulations are performed within "CHASTE" -the "Cancer, Heart And Soft-Tissue Environment" -a multipurpose software library for simulations in the area of biology (Pitt-Francis et al., 2009) provided under the LGPL license by a team from the Computer Science Department of the University of Oxford, the UK.
Examples showing the membrane voltage distribution at three different time instants during the atrial and ventricular depolarization are provided in Fig. 5. Due to the constant action potential durations considered in the simulations, the repolarizations follow a similar time course.Rectangular current pulses of 2 ms duration are used for the initial stimulations.The areas of initial stimulations are taken from the ECGsim (van Oosterom, 2004) software where depolarization isochrones on the surfaces enveloping the atria and ventricles are obtained from the potentials on the surface of the body by following an inverse procedure.The atria are stimulated in an area corresponding to the position of the SAnode.In the case of the ventricles, three areas inside the left ventricle and one area inside the right ventricle are initially stimulated with some small time delays between them.Animated representations of the spread of electrical excitation during the depolarization and repolarization of the atria and ventricles are also provided in the form of short films.

Computation of currents in the body
The electrical excitation in the heart gives rise to currents in the entire body.Impressed current densities in the heart can be computed as J e = −σ ∇V m (Aoki et al., 1987)   ure, the so called butterfly plot is obtained as shown in Fig. 9.

Conclusions
The method used in this work proved to be able to provide realistic values for the computed magnetic fields of the heart.

335
We coupled different state-of-the-art existent software codes and realistic anatomical models in order to compute the magnetic field produced by the heart.We used the magneto-static approximation to compute slowly varying magnetic field with time at successive time instants.The Finite-Element 340 method (FEM) with tetrahedral elements was used for com- solving a stationary volume conduction problem described by the equation ∇J t = 0 where J t = J e + J c is the total current density.We compute the total current density in the entire torso shown on the left-hand side of Fig. 4 by using the commercially available software for Finite-Element analysis called COMSOL Multiphysics (COMSOL Multiphysics v. 4.2, 2011).
The electrical conductivity values used for the torso and the lungs are σ torso = 0.2 Sm −1 and σ lung = 0.04 Sm −1 , respectively.

Computation of magnetic fields (MCG signals)
The currents in the torso produce magnetic fields inside and outside the body.We compute these magnetic fields using Ampere's law ∇ ×H = J t at each time instant previously considered in the computation of the total current density in the torso J t .These simulations are coupled in COMSOL Multiphysics (COMSOL Multiphysics v. 4.2, 2011) with those for the solution of the volume conduction problems.All biological tissues are considered to be non-magnetic with the relative magnetic permeability µ r = 1.
For the results we consider a four by four magnetic sensor array positioned just in front of the torso as shown in Fig. 6.The nodes of the magnetic sensor array are spaced 8.2 cm and 7.2 cm apart from each other on the horizontal and vertical directions, respectively.
The variation in time of the magnetic field component perpendicular to the chest computed in each node of the magnetic sensor array is shown in  puting conduction currents in the inhomogeneous and irregularly shaped torso.Future work will focus on the reduction of the number of unknowns in order to reduce the computation time as well as for more conveniently solving the corre-345 sponding inverse problem.

Fig. 1 .
Fig.1.Schematic representation of a section through the human heart showing the excitation-conduction system.Picture adapted from(Malmivuo and Plonsey , 1995).

Fig. 1 .
Fig. 1.Schematic representation of a section through the human heart showing the excitation-conduction system.Picture adapted from(Malmivuo and Plonsey, 1995).

Fig. 3 .
Fig. 3. Action potential shape of the Luo-Rudy phase cellular ionic model.
x s being gating variables with values varying in the interval of 22 years.They are available in the form of triangulated 220 surfaces as shown in Fig. 4.

Fig. 4 .
Fig. 4. MRT-based anatomical mesh models provided with th freely available software code ECGsim (van Oosterom, 2004).Left heart and lungs inside the torso.Right: detailed view of the heart. 245

Fig. 4 .
Fig. 4. MRT-based anatomical mesh models provided with the freely available software code ECGsim (van Oosterom, 2004).Left: heart and lungs inside the torso.Right: detailed view of the heart.

Fig. 5 .
Fig. 5.The spread of electrical excitation in the heart during the atrial depolarization (top row), and ventricular depolarization (bottom row).

260Fig. 5 .
Fig. 5.The spread of electrical excitation in the heart during the atrial depolarization (top row), and ventricular depolarization (bottom row).

Fig. 6 .
Fig. 6.Position of the magnetic sensor array (red the body.

Fig. 6 .
Fig. 6.Position of the magnetic sensor array (red color) relative to the body.

Fig.Fig. 8 .
Fig. Magnetic flux density component perpendicular on the chest recorded as a function of time at the position of each node of the sensor array shown in Fig. 6.

Fig. 7 .
Fig. 7. Magnetic flux density component perpendicular on the chest recorded as a function of time at the position of each node of the sensor array shown in Fig. 6.

Fig. 7 .Fig. 8 .
Fig. 7. Magnetic flux density component perpendicular on the chest recorded as a function of time at the position of each node of the sensor array shown in Fig. 6.

Fig. 7 .
Fig. 7.In this figure all the magnetocardiographic signals are represented on the same scale.Detailed views of the signals recorded by the sensors 6 and 12 are given in Fig. 8.When plotting the signals recorded by all the sensors in a single fig-330

Fig. 8 .
Fig. 8. Detailed view of the magnetocardiographic signals recorded at the positions of the sensors denoted by 6 (left) and 12 (right) in Fig. 7.

Fig. 7 .
In this figure all the magnetocardiographic signals are represented on the same scale.Detailed views of the signals recorded by the sensors 6 and 12 are given in Fig. 8.When plotting the signals recorded by all the sensors in a single figure, the so called butterfly plot is obtained as shown in Fig. 9.

Fig. 9 .
Fig. 9. Butterfly plot of the magnetic flux density component perpendicular on the chest recorded at the position of each point sensor.