User:Katje95/sandbox

Simulation of a realistic heart beat.

Schematic diagram of normal sinus rhythm for a human heart as seen on ECG.

The forward problem of electrocardiology is a computational and mathematical approach to study the electrical activity of the heart through the body surface. The principal aim of this study is to computationally reproduced an electrocardiogram (ECG), which has important clinical relevance to defined cardiac pathologies such as ischemia and infarction, or to test pharmaceutical intervation. Given the important functionalities and the relative small invasiveness of them, the electrocardiography techniques are used quite often as clinical diagnostic tests. Thus, it is natural to proceed to computationally reproduce an ECG, which means to mathematically model the cardiac behaviour inside the body.

The three main part of a forward model for the ECG are:

• a model for the cardiac electrical activity.
• a model for the diffusion of the electrical potential inside the torso, which represents the extracardiac region.
• some specific heart-torso coupling conditions.

Thus, to obtained an ECG, a mathematical electrical cardiac model must be considered, coupled with a diffusive model in a passive condutctor that describe the electrical propagation inside the torso.

The coupled model is usually a three-dimensional model made of partial differential equations (PDE). Such model is solved considering finite element method (FEM) for the solution's space evolution and semi-implicit numerical schemes for the solution's time evolution. However, the computational costs of such techniques, especially with three dimensional simulations, are quite high. Thus, a simplify model are often considered, solving for example the heart electrical activity separately from the problem on the torso. To provide realistic results, three dimensional anatomically-realistic models of the heart and the torso must be used.

Another possible simplification is a dynamical model made of three ordinary differential equations.

Heart tissue models

The electrical activity of the heart is caused by the ions flow across the cell membrane, between the inside and the outside of the cell, which determine a wave of excitation along the heart muscle that coordinates the cardiac contraction and, thus, the pump of the blood to the body. The modeling of cardiac electrical activity is thus related to the modelling of the ions flow on a microscopic level, and on the propagation of the exitation wave along the muscle fibers on a macroscopic level.

Between the mathematical model on the macroscopic level, Willem Einthoven and Augustus Waller defined the ECG through the conceptual model of a dipole rotating aroung a fixed point, whose projection on the lead axis determined the lead recordings. Then, a two-dimensional reconstruction of the heart activity in the frontal plane was possible using the Einthoven's limbs leads I, II and III as theoretical basis. Later on, the rotating cardiac dipole was considered inadequate and was substiduted by multipolar sources moving inside a bounded torso domain. Unfortunately, the short-comings of the methods used to quantified these sources are the lack of detailes, which are important to simulate cardiac phenomena.

On the other hand, microscopic models try to represent the behaviour of single cells and to connect them considering their electrical properties. Unfortunately, these techniques present some problems associated with the different scales that need to be connected and, expecially for large scale phenomenon such as re-entry or body surface potential, the collective behaviour of the cells is more importan than that of every single cell.

The third option to model the electrical activity of the heart is to consider a so called middle-out approach, where the model wants to incorporate both lower and higher level of details. This option considers the behaviour of a block of cells, called a continuum cell, thus avoiding scale and detail problems. The model obtained is called a bidomain model which is often replaced by its simplification, the monodomain model.

Bidomain model

Stylized representation of a human torso that describes the domain and the notation considered for the forward problem of electrocardiography. Two separated areas and their boundary are considered, which represent the heart and human torso around it.

The basic idea of the bidomain model is that the heart tissue can be dived in two ohmic conducting media, connected but separated throught the cell membrane. This media are called intracellular and extracellular regions, the first of which represents the cellular tissues, and the second one the space between them.

The standard formulation of the bidomain model, which considers also a dynamical model for the ionic current, is the following

${\begin{cases}\nabla \cdot (\mathbf {\sigma } _{i}\nabla V_{m})+\nabla \cdot (\mathbf {\sigma } _{i}\nabla u_{e})=A_{m}\left(C_{m}{\frac {\partial V_{m}}{\partial t}}+I_{ion}(V_{m},w)\right)+I_{app}&{\text{in }}\Omega _{H}\\\nabla \cdot \left((\mathbf {\sigma } _{i}+\mathbf {\sigma } _{e})\nabla u_{e}\right)+\nabla \cdot (\mathbf {\sigma } _{i}\nabla V_{m})=0&{\text{in }}\Omega _{H}\\{\frac {\partial w}{\partial {t}}}+g(V_{m},w)=0&{\text{in }}\Omega _{H}\end{cases}}$

where $V_{m}$ and $u_{e}$ are the transmembrane and extracellular potentials respectively, $I_{ion}$ is the ionic current, which depends also from a gating variable $w$, and $I_{app}$ is an external current applied to the domain. Moreover, $\mathbf {\sigma } _{i}$ and $\mathbf {\sigma } _{e}$ are the intracellular and extracellular conductivity tensors, $A_{m}$ is the surface to volume ratio of the cell membrane and $C_{m}$ is the membrane capacitance per unit area. Here $\Omega _{H}$ represents the heart domain.

The boundary conditions for this version of the bidomain model are obtained throught the assumption that there is no flow of intracellular potential outside of the heart, which means that

$\mathbf {\sigma } _{i}\nabla V_{m}\cdot \mathbf {n} +\mathbf {\sigma } _{i}\nabla u_{e}\cdot \mathbf {n} =0\quad \quad {\text{on }}\Sigma$

where $\mathbf {n}$ is the outward unit normal to $\Omega _{h}$ and $\Sigma$ represents the boundary of the heart domain.

Monodomain model

The monodomain model is a simplification of the bidomain model but netherless it is able to represent realistic electrophysiological phenomena in the sense of the transmembrane potential $V_{m}$.

The standard formulation is a partial differential equations with only one unkonwn, which is the transmembrane potential:

$\chi C_{m}{\frac {\partial V_{m}}{\partial t}}-\nabla \cdot \left(\mathbf {\sigma } _{i}{\frac {\lambda }{1+\lambda }}\nabla V_{m}\right)+\chi I_{ion}=I_{app}\quad \quad {\text{in }}\Omega _{H}$

where $\lambda$ is a parameter that related the intracelluar and extracellular conductivity tensors.

The boundary condition used for this model is

$\left(\mathbf {\sigma } _{i}\nabla V_{m}\right)\cdot \mathbf {n} =0\quad \quad {\text{on }}\partial \Omega _{H}.$

Torso tissue model

In the forward problem of electrocardiography, the torso is seen as a passive conductor and its model can be found starting from the Maxwell's equations under a quasi-static assumption.

The standard formulation considers a partial differential equation with one unknown, the torso potential $u_{T}$. Basically, the torso model is a generalized Lapalce equation

$\nabla \cdot (\mathbf {\sigma } _{T}\nabla u_{T})=0\quad \quad {\text{in }}\Omega _{T},$

where $\mathbf {\sigma } _{T}$ is a conductivity tensor and $\Omega _{T}$ is the boundary domain, i.e. the human torso.

Derivation

As for the bidomain model, the torso model can be derived from the Maxwell's equations and the continuity equation after some assumptions. First of all, since the electrical and magnetical activity inside the body is generated at low level, a quasi-static assumption can be considered. Thus, the body can be viewed as a passive conductor, which means that its capacitive, inductive and propagative effect can be ignored.

Under quasi-static assumption, the Maxwell's equations are

${\begin{aligned}\nabla \cdot \mathbf {E} &={\frac {\rho }{\epsilon _{0}}}\\\nabla \times \mathbf {E} &=0\\\nabla \cdot \mathbf {B} &=0\\\nabla \times \mathbf {B} &=\mu _{0}\mathbf {J} \end{aligned}}$

and the continuity equation is

$\nabla \cdot \mathbf {J} =0.$

Since its curl is zero, the electrical field can be represented by the gradient of a scalar potential field, the torso potential

$\mathbf {E} =-\nabla u_{T}$

(1)

where the negative sign means that the current flows from higher to lower potential regions.

Then, the total current density can be expressed in terms of the conduction current and other different applied currents so that, from continuity equation,

$\nabla \cdot \mathbf {J} =\nabla \cdot (\mathbf {J} _{app}+\mathbf {\sigma } _{T}\mathbf {E} )=0$

(2)

Then, substituting (1) in (2)

$\nabla \cdot (\mathbf {\sigma } _{T}\nabla u_{T})=\nabla \cdot \mathbf {J} _{app}=I_{v}$

in which $I_{v}$ is the current per unit volume.

Finally, since aside from the heart there is not current source inside the torso, the current per unit volume can be set to zero, giving the generalized Laplace equation which represents the standard formulation of the diffusive problem inside the torso

$\nabla \cdot (\mathbf {\sigma } _{T}\nabla u_{T})=0\quad \quad {\text{in }}\Omega _{T}.$

Boundary condition

The boundary condition represents the property of the media surrounding the torso. Generally, the air have zero conductivity which means that the current can not flow outside of the torso. This is translated in the following equation

$\mathbf {\sigma } _{T}\nabla u_{T}\cdot \mathbf {n} _{T}=0\quad \quad {\text{on }}\Gamma _{T}$

where $\mathbf {n} _{T}$ is the unit outward normal to the torso and $\Gamma _{T}$ is the torso boundary, which means the torso surface.

Torso conductivity

Usually, the torso is considered to have isotropic conductivity, which means that the current can flow in the same way in all directios. However, the torso is not an empty or homogeneous enveloped, but containes different organes characterised by different conductivity coefficients which can be experimentally obtained. A simple example of conductivity parameters in a torso that consideres the bones and the lungs are reported in the following table.

Values of the torso conductivity.

$\sigma _{T}^{\text{lungs}}$	$\sigma _{T}^{\text{bones}}$	$\sigma _{T}^{\text{remaining regions}}$
(S/cm)	(S/cm)	(S/cm)
$2.4\times 10^{-4}$	$4\times 10^{-5}$	$6\times 10^{-4}$

Heart-torso models

The coupling between the electrical activity model and the torso model is done applying some boundary condition at the epicardium, which means at the inteface surface between the heart and the torso.

The heart-torso model can be fully coupled, which means that a perfect electrical transmission between the two domains are considered, or can be uncoupled, which means that the heart electrical model and the torso model can be solved separately but with an exchanged of informations between them.

Fully coupled heart-torso models

The complete coupling between the heart and the torso is obtained imposing a perfect electrical transmission between the heart and the torso. This is done considering two equations that establish a relationship between the extracellular potential and the torso portensial, which are

${\begin{aligned}u_{e}&=u_{T}\quad \quad &{\text{on }}\Sigma \\(\mathbf {\sigma } _{e}\nabla u_{e})\cdot \mathbf {n} _{e}&=-(\mathbf {\sigma } _{T}\nabla u_{T})\cdot \mathbf {n} _{T}\quad \quad &{\text{on }}\Sigma .\end{aligned}}$

This equations ensure the continuity of both the potential and the current across the epicardium.

Using these boundary conditions, it is possible to obtained two dirrent fully complete heart-torso models, one that considers the bidomain model for the heart electrical activity, and the other that considers the monodomain model for the same purpouse. Unfortunately, numerically, the two model are computationally high expensive with similar computation costs.

Alternative boundary conditions

Boundary conditions that represent a perfect electrical coupling between the heart and the torso are the most used and the classical ones. However, between the heart and the torso there is the pericardium, a sac with a double-walled that containes a serous fluid wich have a specific effect on the electrical transmission. Considering the capacitance $C_{p}$ and the resistive$R_{p}$ effect that the pericardium has, alternative boundary conditions that take into account this effect can be formulated

${\begin{aligned}R_{p}\mathbf {\sigma } _{e}\nabla u_{e}\cdot \mathbf {n} &=R_{p}C_{p}{\frac {\partial (u_{T}-u_{e})}{\partial t}}+(u_{T}-u_{e})\quad \quad &{\text{on }}\Sigma \\\mathbf {\sigma } _{e}\nabla u_{e}\cdot \mathbf {n} &=\mathbf {\sigma } _{T}\nabla u_{T}\cdot \mathbf {n} \quad \quad &{\text{on }}\Sigma \end{aligned}}$

Formulation with the bidomain model

The fully coupled heart-torso model which considered the bidomain model for the heart electrical activity in its complete form is

${\begin{cases}\nabla \cdot (\mathbf {\sigma } _{i}\nabla V_{m})+\nabla \cdot (\mathbf {\sigma } _{i}\nabla u_{e})=A_{m}\left(C_{m}{\frac {\partial V_{m}}{\partial t}}+I_{ion}(V_{m},w)\right)+I_{app}&{\text{in }}\Omega _{H}\\\nabla \cdot \left((\mathbf {\sigma } _{i}+\mathbf {\sigma } _{e})\nabla u_{e}\right)+\nabla \cdot (\mathbf {\sigma } _{i}\nabla V_{m})=0&{\text{in }}\Omega _{H}\\{\frac {\partial w}{\partial {t}}}+g(V_{m},w)=0&{\text{in }}\Omega _{H}\\\nabla \cdot (\mathbf {\sigma } _{T}\nabla u_{T})=0&{\text{in }}\Omega _{T}\\(\mathbf {\sigma } _{i}\nabla V_{m})\cdot \mathbf {n} +(\mathbf {\sigma } _{i}\nabla u_{e})\cdot \mathbf {n} =0&{\text{on }}\Sigma \\u_{e}=u_{T}&{\text{on }}\Sigma \\(\mathbf {\sigma } _{e}\nabla u_{e})\cdot \mathbf {n} =(\mathbf {\sigma } _{T}\nabla u_{T})\cdot \mathbf {n} &{\text{on }}\Sigma \\(\mathbf {\sigma } _{T}\nabla u_{T})\cdot \mathbf {n} =0&{\text{on }}\Gamma _{T}\end{cases}}$

where the first four equations are the partial differential equations representing the bidomain and the torso model, while the remaining ones represents the boundary conditions for the bidomain and torso model and the coupling conditions between them.

Formulation with the monodomain model

The fully coupled heart-torso model which considered the monodomain model for the elctrical activity of the heart is more complicated that the bidomain problem. The coupling conditions use, infact, the extracellular potential which is not computed by the monodomain model. Thus, it is necessary to used also the second equation of the bidomain model which gives

$\nabla \cdot \left((\mathbf {\sigma } _{i}+\mathbf {\sigma } _{e})\nabla u_{e}\right)+\nabla \cdot (\mathbf {\sigma } _{i}\nabla V_{m})=0\quad \quad {\text{in }}\Omega _{H}.$

In this way, the coupling conditions do not need to changed, and the complete heart-torso model is made of two different blocks:

• First the monodomain model with its usual boundary condition must be solve, which is

${\begin{cases}\chi C_{m}{\frac {\partial V_{m}}{\partial t}}-\nabla \cdot \left(\mathbf {\sigma } _{i}{\frac {\lambda }{1+\lambda }}\nabla V_{m}\right)+\chi I_{ion}=I_{app}&{\text{in }}\Omega _{H}\\\left(\mathbf {\sigma } _{i}\nabla V_{m}\right)\cdot \mathbf {n} =0&{\text{on }}\partial \Omega _{H}.\end{cases}}$

• Then, the coupled model that considers the recover of the extracellular potential and the torso model with the coupling conditions must be solved, which is

${\begin{cases}\nabla \cdot \left((\mathbf {\sigma } _{i}+\mathbf {\sigma } _{e})\nabla u_{e}\right)+\nabla \cdot (\mathbf {\sigma } _{i}\nabla V_{m})=0&{\text{in }}\Omega _{H}\\\nabla \cdot (\mathbf {\sigma } _{T}\nabla u_{T})=0&{\text{in }}\Omega _{T}\\\mathbf {\sigma } _{T}\nabla u_{T}\cdot \mathbf {n} =0&{\text{on }}\Gamma _{T}\\u_{e}=u_{T}&{\text{on }}\Sigma \\(\mathbf {\sigma } _{e}\nabla u_{e})\cdot \mathbf {n} =(\mathbf {\sigma } _{T}\nabla u_{T})\cdot \mathbf {n} &{\text{on }}\Sigma \end{cases}}$

Uncoupled heart-torso models

The fully coupled heart-torso models are very detailed models but they are also computationally expensive to solve. A possible simplification is the so called uncoupled assumption in which the heart is considered completely electrically isoleted from the heart. Mathematically, this is done imposing that the current can not flow across the epicardium, from the heart to the torso, namely

$\mathbf {\sigma } _{e}\nabla u_{e}\cdot \mathbf {n} =0\quad \quad {\text{on }}\Sigma .$

Applying this equation to the boundary conditions of the fully coupled models, it is possible to obtained two uncoupled heart-torso models, in which the electrical models can be solved separately from the torso model reducing the computational cost.

Uncoupled heart-torso model with the bidomain model

The uncoupled version of the fully coupled heart-torso model that uses the bidomain to represents the elctrical activity of the heart is made of two separated part:

• The bidomain model in its isolated form

${\begin{cases}\nabla \cdot (\mathbf {\sigma } _{i}\nabla V_{m})+\nabla \cdot (\mathbf {\sigma } _{i}\nabla u_{e})=A_{m}\left(C_{m}{\frac {\partial V_{m}}{\partial t}}+I_{ion}(V_{m},w)\right)+I_{app}&{\text{in }}\Omega _{H}\\\nabla \cdot \left((\mathbf {\sigma } _{i}+\mathbf {\sigma } _{e})\nabla u_{e}\right)+\nabla \cdot (\mathbf {\sigma } _{i}\nabla V_{m})=0&{\text{in }}\Omega _{H}\\{\frac {\partial w}{\partial {t}}}+g(V_{m},w)=0&{\text{in }}\Omega _{H}\\(\mathbf {\sigma } _{i}\nabla V_{m})\cdot \mathbf {n} +(\mathbf {\sigma } _{i}\nabla u_{e})\cdot \mathbf {n} =0&{\text{on }}\Sigma \\\mathbf {\sigma } _{e}\nabla u_{e}\cdot \mathbf {n} =0&{\text{on }}\Sigma .\end{cases}}$

• The torso diffusive model in its standard formulation plus the potential continuity condition

${\begin{cases}\nabla \cdot (\mathbf {\sigma } _{T}\nabla u_{T})=0&{\text{in }}\Omega _{T}\\\mathbf {\sigma } _{T}\nabla u_{T}\cdot \mathbf {n} =0&{\text{on }}\Gamma _{T}\\u_{T}=u_{e}&{\text{on }}\Sigma \\\end{cases}}$

Uncoupled heart-torso model with the modomain model

As in the case of the fully coupled heart-torso model which uses the monodomain model, the corrisponding uncoupled model needs to recover the extracellular potential which is not computed by the monodomain. In this case, three different problem must be solved:

• The monodomain model with its usual boundary condition

${\begin{cases}\chi C_{m}{\frac {\partial V_{m}}{\partial t}}-\nabla \cdot \left(\mathbf {\sigma } _{i}{\frac {\lambda }{1+\lambda }}\nabla V_{m}\right)+\chi I_{ion}=I_{app}&{\text{in }}\Omega _{H}\\\left(\mathbf {\sigma } _{i}\nabla V_{m}\right)\cdot \mathbf {n} =0&{\text{on }}\partial \Omega _{H}.\end{cases}}$

• The problem to recover the extracellular potential with the no intracellular current across the epicardium boundary condition

${\begin{cases}\nabla \cdot \left((\mathbf {\sigma } _{i}+\mathbf {\sigma } _{e})\nabla u_{e}\right)+\nabla \cdot (\mathbf {\sigma } _{i}\nabla V_{m})=0&{\text{in }}\Omega _{H}\\(\mathbf {\sigma } _{i}+\mathbf {\sigma } _{e})\nabla u_{e}\cdot \mathbf {n} =-(\mathbf {\sigma } _{i}\nabla V_{m})\cdot \mathbf {n} &{\text{on }}\Sigma \end{cases}}$

• The torso diffusive model with the potential continuity boundary condition to be imposed at the epicardium

${\begin{cases}\nabla \cdot (\mathbf {\sigma } _{T}\nabla u_{T})=0&{\text{in }}\Omega _{T}\\\mathbf {\sigma } _{T}\nabla u_{T}\cdot \mathbf {n} =0&{\text{on }}\Gamma _{T}\\u_{T}=u_{e}&{\text{on }}\Sigma \\\end{cases}}$

Electrocardiogram computation

Precordial leads in ECG

Solving the fully coupled or the uncoupled heart-torso models allows to obtained the electrical potential generated by the heart in all the human torso, and especially on all the surface of the torso. Definint the electrodes positions on the torso, it is possible to find the time evolution of the potential on such points. Then, the electrocardiograms can be computed, for example according to the 12 standard leads, considering the following formulas

${\begin{aligned}I&=u_{T}(L)-u_{T}(R)\\II&=u_{T}(F)-u_{T}(R)\\III&=u_{T}(F)-u_{T}(L)\\aVR&={\frac {3}{2}}(u_{T}(R)-u_{w})\\aVL&={\frac {3}{2}}(u_{T}(L)-u_{w})\\aVF&={\frac {3}{2}}(u_{T}(F)-u_{w})\\V_{i}&=u_{T}(V_{i})-u_{w}\quad i=1,\dots ,6\end{aligned}}$

where $u_{W}=(u_{T}(L)+u_{T}(R)+u_{T}(F))/3$ and $L,R,F,\{V_{i}\}_{i=1,\dots ,6}$ are the standard locations of the electrodes.

Numerical methods

The heart-torso models are made of partial differential equations in which the unknowns change both in space and in time. Moreover, in the case of the bidomain model, an ordinary differential equation is present. Thus, multiple numerical schemes can be used. Usually, finite element method is applied for the space discretization and semi-implicit schemes are used for the time discretization.

Uncoupled heart-torso model are the simplest to treat numerically because the electrical model can be solved separately from the torso one, so that classic numerical methods to solve each of them can be applied. This means that the bidomain and monodomain models can be solved for example whith a backward differentiation formula (BDF) for the time discretizion, while the recover of the extracellular potential and torso problem can be easily solved applying only the finite element method because they are time independent.

The fully coupled heart-torso models, instead, are more complexes and need more sofisticated numerical models. For example, the fully heart-torso model that uses the bidomain model for the electrical simulation of the cardiac behaviour can be solved considering domain decomposition techniques, such as a Dirichlet-Neumann domain decomposition.

Geometric torso model

Three-dimensional torso model including the most organs.

To simulate and electrocardiogram using the fully coupled or uncoupled models, a three-dimensional reconstruction of the human torso is needed. Today, diagnostic imaging techniques such as MRI and CT can provide a sufficient background to recostruct also detailed anatomycal human parts and, thus, appropriate torso geometry. For example the Visible Human Data is a useful dataset to create a threee dimensional torso model detailed with internal organs including the skeletal structure and muscle.

Dynamical model for the electrocardiogram

Even if the results are quite detailed, solving three-dimensional model is usually quite expensive. A possible simplification is a dynamical model based on three coupled ordinary differential equations.

The quasi-periodicity of the heart beat is reproduced by a three-dimensional trajectory around an attracting limit cycle in the $(x,y)$ plane. The principal peaks of the ECG, which are the P,Q,R,S and T, are described at fixed angles $\theta _{P},\theta _{Q},\theta _{R},\theta _{S}{\text{ and }}\theta _{T}$, which gives the following three ODEs

${\begin{aligned}x'&=\alpha z-\omega y\\y'&=\alpha y+\omega x\\z'&=-\sum _{i\in \{P,Q,R,S,T\}}a_{i}\Delta \theta _{i}{\text{exp}}\left(-{\frac {\Delta \theta _{i}^{2}}{2b_{i}^{2}}}\right)-(z-z_{0})\end{aligned}}$

with $\alpha =1-{\sqrt {x^{2}+y^{2}}}$, $\Delta \theta _{i}=(\theta -\theta _{i}){\text{mod}}(2\pi )$, $\theta ={\text{atan}}2(y,x)$

The equations can be easily solved with classical numerical algorithms like Runge-Kutta methods for ODEs.

See also

• Monodomain model
• Bidomain model
• Electrocardiography

References

[[Category:Cardiac electrophysiology]] [[Category:Cardiac procedures]] [[Category:Electrodiagnosis]] [[Category:Electrophysiology]] [[Category:Mathematics in medicine]] [[Category:Medical tests]] [[Category:Differential equations]] [[Category:Partial differential equations]] [[Category:Mathematical modeling]] [[Category:Numerical analysis]]