Depolarization waves in the heart are an important electrophysiological phenomenon which is described by the basic model for the cardiac reaction-diffusion system: The bidomain equations for electrical potentials in and between cells and the spread of the depolarization waves. The cardiac electrophysiological system triggers contractions of the heart muscle, which is modeled as a time-varying elastic body. The electrical cellular transmembrane voltages are given as solutions to nonlinear reaction-diffusion equations on a manifold. Depending on the accuracy of the model, this is coupled with an ODE model in every nodal point of every cell to determine the transmembrane currents.

video

## Mathematical modeling

The numerical simulation of the electrical behavior of the heart is best described by a mathematical model called the bidomain model. However, it has been shown that the reduced monodomain model is sufficient to describe the propagation of depolarization waves and electrical signals on the myocardium in many cases. The derivation of the monodomain model is based on the assumption of equal anisotropy ratios in the intracellular and extracellular space and the equations for the cardiac electrophysiology are given by: \[\begin{equation*} \beta \Big( C_\text{m} \frac{\partial v}{\partial t} + \big(I_\text{ion}(v,w)+I_\text{stim}\big)\Big) = \nabla \cdot (\boldsymbol{\sigma} \nabla v)\quad \text{in}~ \Omega. \end{equation*} \] Here \(\beta\) is the surface-to-volume ratio of the cells, \(I_\text{ion}(v,w)\) is the transmembrane current density, \(v\) is the membrane voltage, \(\boldsymbol{\sigma}\) is the conductivity and \(C_\text{m}\) is the capacitance of the membrane. A biophysically detailed model of human ventricular myocytes is described in ten Tusscher and Pantilof [9], and concretizes \(I_\text{ion}(v,w)\).

## Numerical treatment

To solve the monodomain equations we use an operator splitting method. As starting point we have \[\begin{equation*} \partial_t v = \frac{1}{\beta C_\text{m}} \nabla\cdot (\boldsymbol{\sigma} \nabla v) - \frac{1}{C_m}\Big( I_\text{ion}(v,w)+I_\text{stim}\Big), \end{equation*} \] and split it as follows: \[ \begin{align*} \partial_t v &= - \frac{1}{C_\text{m}} \Big(I_\text{ion} v,w)+I_\text{stim}\Big),\\ \partial_t v &= \frac{1}{\beta C_\text{m}} \nabla\cdot (\boldsymbol{\sigma} \nabla v). \end{align*} \] Then we have to solve a nonlinear system of ODEs and a linear PDE. We use the second order method by Strang [7].

## Verification of the monodomain setup

To test our implementation we use the benchmark proposed by Niederer et al. [3] which is realized by several international groups. We reproduce their setup regarding the model equations and their simplifications. Their geometry is selected to be very simple: a cuboid \(\Omega = (0,3)\times(0, 7)\times(0,20)\). The stimulus current \(I_\text{stim}(t,x)\) is delivered to a volume of \(\Omega_\text{stim}=(0,1.5)\times(0, 1.5)\times(0,1.5)\) located at a corner of the cuboid and is given by: \[ \begin{equation*} I_\text{stim}(t,x)=\begin{cases} 50 & 0\leq t \leq 0.002 ~\text{and}~ x\in \Omega_\text{stim},\\ 0 & \text{else.}\end{cases} \end{equation*} \]

video: Geometry of the cuboid with stimulation area

## Realistic simulation pattern

In the numerical examples above we use a short, local, external impulse \(I_\text{stim}\) to initiate depolarization waves. Here we have a detailed geometry of the ventricles equipped with labeled cells distributed in the endocardium to mimick the excitation by the Purkinje network. The Purkinje network is simulated by different excitation times at these cells to apply \(I_\text{stim}\). Nevertheless, the duration and the amplitude of the external stimulus is the same.

video?

## A four chamber model with circulation

We develop a lumped element model of the closed loop human vascular system together with embedding the heart model into a surrounding elastic tissue by forcing permanent contact at the pericardial surface. This leads to more realistic temporal behavior of atrial volumes and atrial pressure-volume curves. By estimating the cardiac pressures that lead to continuous flows across the model interfaces, we were able to accomplish a strong coupling between the heart and the vascular system.

Furthermore we have made first steps towards a coupled *More*{electro-mechanical model}.

It is determined by

- the electric potential \(V_\text{m}\);
- a vector of ion concentrations and gating variables \(w\);
- the deformation vector \(\boldsymbol{\varphi}: \Omega\times[0,T]\rightarrow\mathbb{R}^3\);

in the reference configuration \(\Omega\subset\mathbb{R}^3\) for the heart tissue and pericardium, such that \[ \begin{equation*} \beta C_\text{m} \partial_t^{\boldsymbol{\varphi}} V_\text{m} = \nabla^{\boldsymbol{\varphi}}\cdot\sigma \nabla^{\boldsymbol{\varphi}} V_\text{m} + \beta I_\text{ion}(V_\text{m},w;\boldsymbol{d}^{\boldsymbol{\varphi}}) + I_\text{ext} \end{equation*} \] with total derivative in time \( \partial_t^{\boldsymbol{\varphi}} w = \partial_t w+\nabla^{\boldsymbol{\varphi}}\!\cdot\!(w \partial_t \boldsymbol{\varphi}) \) and the covariant gradient \(\nabla^{\boldsymbol{\varphi}}\). The depolarization wave is initiated be an external impulse \(I_\text{ext}\) located at predefined nodes. The model is complemented by initial values and boundary conditions. For a detailed explanation see [8]. A further coupling to the circulatory system and the embedding in the pericardium is discussed in [2]. *end: more*

We do not use a strongly coupled, monolithic approach. Instead, we initialize both, the electrophysiology model and the elastomechanic model and exchange the tension development as well as the rotation, deformation and stretch between the models. To achieve this, we first update the electrophysiology model and continue by updating the elastomechanic model.

Both models operate on different spatial and temporal resolutions that differ by multiple magnitudes. To accommodate for these differences, we either subdivide the coarse mechanics mesh or use a new, finer mesh for the electrophysiology model.

Due to the high computational demand of simulating coupled multiphysics problems, we work on reducing the complexity of the mathematical models that describe the myocardium cells. We use the manifold boundary approximation method to systematically reduce highly detailed, physiological cell models to build a hierarchy of cell models that allow an efficient simulation of our coupled model. So far, we have managed to reduce one particular cell model, consisting of 41 ODEs, to a model with only 23 ODEs. In consequence, we reduced the computation time in a single cell to about 50%.

Additionally, we made various studies on different relevant aspects of cardiac electrophysiology; see for example [6, 5, 4].

## Outlook

Currently we work on the coupling of the mechanics and the electropysiology including the new implementations in M++:

- the verified mododomain model with the interface to the cell model library;
- more realistic material laws of Guccione or Ogden type;
- improving the full model with more realistic properties like fiber orientations (see [1] and a realistic excitation profile;
- characterization of cellular pacemaking and its synchronization in the human sinus node.

## References

[1] Lukas Baron, Thomas Fritz, Gunnar Seemann, and Olaf Dössel. Sensitivity study of fiber orientation on stroke volume in the human left ventricle. In *Computing in Cardiology Conference (CinC), 2014*, pages 681–684. IEEE, 2014.

[2] Thomas Fritz, Christian Wieners, Gunnar Seemann, Henning Steen, and Olaf Dössel. Simulation of the contraction of the ventricles in a human heart model including atria and pericardium. *Biomechanics and modeling in mechanobiology*, 13(3):627–641, 2014.

[3] Sander Land, Viatcheslav Gurev, Sander Arens, Christoph M. Augustin, Lukas Baron, Robert Blake, Chris Bradley, Sebastian Castro, Andrew Crozier, Marco Favino, et~al. Verification of cardiac mechanics software: benchmark problems and solutions for testing active and passive material behaviour. * Proc. R. Soc. A*, 471(2184):20150641, 2015.

[4] Heiko Lehrmann, Amir~S Jadidi, Jan Minners, Juan Chen, Björn Müller-Edenborn, Reinhold Weber, Olaf Dössel, Thomas Arentz, and Axel Loewe. Novel electrocardiographic criteria for real-time assessment of anterior mitral line block:“v1 jump” and “v1 delay”. *JACC: Clinical Electrophysiology*, page 655, 2018.

[5] Axel Loewe, Martin W. Krueger, Fredrik Holmqvist, Olaf Dössel, Gunnar Seemann, and Pyotr~G Platonov. Influence of the earliest right atrial activation site and its proximity to interatrial connections on p-wave morphology. *EP Europace*, 18(suppl_4):iv35–iv43, 2016.

[6] Axel Loewe, Mathias Wilhelms, Fathima Fischer, Eberhard P. Scholz, Olaf Dössel, and Gunnar Seemann. Arrhythmic potency of human ether-à-go-go-related gene mutations l532p and n588k in a computational model of human atrial myocytes. *Europace*, 16(3):435–443, 2014.

[7] Shev MacNamara and Gilbert Strang. Operator splitting. In *Splitting Methods in Communication, Imaging, Science, and Engineering*, pages 95–114. Springer, 2016.

[8] Alfio Quarteroni, Toni Lassila, Simone Rossi, and Ricardo Ruiz-Baier. Integrated heart—coupling multiscale and multiphysics models for the simulation of the cardiac function. *Computer Methods in Applied Mechanics and Engineering*, 314:345–407, 2017.

[9] Kirsten H. W. J. ten Tusscher and Alexander V. Panfilov. Alternans and spiral breakup in a human ventricular tissue model. *American Journal of Physiology-Heart and Circulatory Physiology*, 2006.