Parallel Left Ventricle Simulation Using the FEniCS Framework

Abstract

Heart simulation is complex task that requires multiscale modeling on cell, tissue and organ levels. Such structure makes di cult high performance code development and its maintenance. In this paper, we evaluate how scienti c software could be used for heart simulation. An overview of existing frameworks for automated scienti c computing is presented. The FEniCS framework was chosen since it supports automated solution of di erential equations by the nite element method of parallel computing systems, provides near-mathematical notation, uses high performance backend and has comprehensive documentation. FEniCS performance was evaluated by simulation the space propagation of membrane potential alternation over a cardiac left ventricle using the electrophysiological model of a left ventricle and the EkaterinburgOxford cell model. The FEniCS framework showed good performance and near-linear scalability up to 240 CPU cores.

Keywords

• FEniCS
• heart simulation
• parallel computing ment method

-

Heart simulation is a complex task that requires multiscale modeling on cell, tissue, and organ levels [ 1 ]. Such problems are computationally intensive and require parallel computing, including the use of modern computational accelerators such as GPU and Xeon Phi. However, porting a complicated multilevel simulation code to a new computational architecture requires a long time (often 3{5 years) during which the architecture may become obsolete. In addition, adaptation for parallel computing architectures often leads to signi cant changes of code. Consequently, it is very di cult to determine which mathematical models and numerical methods are used in the optimized code. Hence, re ecting the changes of a mathematical model in the optimized code can be very complicated. As a result, complex multiscale simulation software is rarely adapted to modern parallel computing architectures.

An alternative approach is based on using automated scienti c computing frameworks. Such frameworks allow the development of simulation software using programming languages with near-mathematical notation. Traditionally, a signi cant disadvantage of such frameworks was their low simulation performance. However, some modern implementations use advanced tools such as highly efcient mathematical libraries, just-in-time compilers, parallel execution, and so on, which improve their performance. Still, it is not clear if the modern automated scienti c computing frameworks are e cient enough for real multiscale simulation tasks such as heart simulation.

In this paper, we evaluate how the automated scienti c computing frameworks could be used for heart simulation on parallel computing systems. Heart simulation is an attractive task for evaluation of the performance of automated scienti c computing frameworks because numerical simulations can be used as a virtual environment for testing and predicting tissue behavior in the cases where experimental techniques can not be applied [ 1, 2 ]. As a benchmark problem, we chose the investigation of the space propagation of the membrane potential alternation over the left ventricle of human heart. 2

Automated Scienti c Computing Frameworks

Nowadays, nite element method is the most popular method for numerical investigation of systems with complex geometry, such as human heart. Several frameworks for automated scienti c computing using nite element method exist. The most popular among them are OpenFOAM, OpenCMISS, Chaste, and FEniCS.

OpenFOAM (Open Source Field Operation and Manipulation) [ 3 ] is a free Computational Fluid Dynamics (CFD) software designed for solving problems in continuum mechanics. OpenFOAM is a C++ library that provides numerical schemes implemented in the traditional nite volume framework with solvers that are known to be e cient for continuum mechanics problems. OpenFOAM uses domain decomposition in order to provide parallel execution based on MPI.

OpenCMISS is a part of the global international project Physiome [ 4 ]. OpenCMISS is a mathematical modeling environment that enables the application of nite element analysis techniques to a variety of complex bioengineering problems. For distributed and parallel computing, OpenCMISS uses MPI, OpenMP, and ParMETIS.However, the OpenCMISS project is still under development, its documentation is incomplete, and it lacks examples.

Chaste (Cancer, Heart, and Soft Tissue Environment) is a general purpose simulation package aimed at multiscale, computationally intensive problems arising in biology and physiology [ 5, 6 ]. Current functionality includes the tissue and cell level electrophysiology, the discrete tissue modeling, and the soft tissue modeling. Chaste uses solvers from PETSc, mesh distribution algorithms from ParMETIS, and provides the ability to run parallel simulations using MPI. Chaste has already implemented models and methods, but to modify them a developer has to deal with sophisticated internal structure. Hence their further extension is complicated.

FEniCS [ 7 ] is a collaborative project for the development of the tools for automated scienti c computing with a particular focus on the automated solution of di erential equations by the nite element method. Implementation of a nite element method consists of several stages (obtaining equations' weak form, discretizing the equations in the weak form, assembling the values calculated on each element, and solving of the system of algebraic equations). Each stage is covered by a separate component of the FEniCS framework. It uses third-party high-performance libraries such as PETSc, Trilinos, uBLAS, or Intel MKL. In addition, FEniCS allows parallel simulation using MPI.

For our heart simulation task, we chose the FEniCS framework because it provides the automated solution of di erential equations by nite element methods, supports automatic parallel execution using MPI, and has a good documentation. 3

Description of the Heart Model

We simulated the electrical activity of a human heart, which is a result of spatial and temporal propagation of the electrical signal from each cardiac cell. We used the Ekaterinburg-Oxford cell model [ 8 ] with a detailed description of electrical, chemical, and mechanical processes. On the intracellular level, the electrical potential arises from a very complicated interaction among ionic currents and cell organelles (organized structures in cells), as presented in Fig. 1 (left). The trans-membrane potential is shown in Fig. 1 (right). Scaled windows presents fast processes in period between 0 and 0.015 seconds.

The electrical block of the model contains the equations of the membrane potential and dynamic parameters describing the opening and closing of ion channels. The chemical block describes the kinetics of intracellular concentrations of calcium, sodium and potassium, extracellular potassium, the kinetics of calcium complexes, and calcium kinetics in the organelles. The mechanical block of the model includes equations describing the voltage of the cell that depends on its length. The di erential equations of the electrical, chemical, and mechanical processes can be presented in the following simpli ed form: (1) (2) where V is the electrical potential, S is the vector of model variables that govern the ion currents, and g is the vector-valued function that describes the time evolution of each variable. The dimension of the vectors S and g is 30. System (1) is de ned at each point of the heart tissue, and, consequently, we should solve it for each node of the computational mesh. The space and time propagation of the electrical potential is governed by the \reaction-di usion equation" [ 9 ]: where D is the di usion coe cient, r2 is the Laplace operator, and Iions is the sum of the ionic currents related to the capacitance of cell membrane. Boundary conditions correspond to the condition of electrical isolation.

This model is a nonlinear system of di erential equations that can not be solved analytically and is a very computationally intensive task due to the large amount of variables in the 3D domain. 4

Benchmark Problem

During the experiments, we simulated the electrical activity of the human heart left ventricle (LV). We used the asymmetric model of LV that was previously developed in our group [ 10 ]. The important feature of the model is the ability to vary the size of the mesh elements. The personalized model parameters were captured with the help of ultrasound imaging. An example of 3D mesh for LV is presented in Fig. 2.

In order to solve model (1){(2), we use the operator splitting scheme of rst order (Marchuk{Yanenko method) [ 11 ]. Let us consider time domain t 2 [0; T ] and the uniform grid tn = htn, where ht = T =N and n is an integer that counts time level, 0 <= n <= N . We denote a quantity at time tn as Vn. The scheme of computing Vn and Sn consists of two steps. Let us assume that we have already calculated the values of V (t) and S(t) for t < tn. At the rst step, we solve the following partial di erential equation:

Vn

V ht n 1 = Dr V ; V (t = tn 1) = V (tn 1); t 2 [tn 1; tn]: 2 (3)

At the second step, we should solve the following system of ordinary equations:

Vn Sn

V ht

S ht n 1 = Iions; V

(t = tn 1) = V (tn); n 1 = g(V ; S ); S (t = tn 1) = S(tn 1): (4)

The solution of (4) gives us the values of V (tn) and S(tn) according to the rules V (tn) = V (tn) and S(tn) = S (tn). This method allows us to tackle task (3) using the implicit method and use the explicit time-scheme for task (4). In addition, we avoid the Newton-like iterations. The disadvantage of such approach is that we have to use a very small integration time step, in order to capture the fast electrochemical processes (Fig. 1).

For testing purposes, we chose activation of an entire LV (the potential is greater than 40 millivolt) as our initial condition. The simulation duration period was 0.3 seconds of physical time because the electrical activity tends to the equilibrium state without an external stimulus approximately after this period.

We use the tetrahedral mesh that was generated by the GMSH software [ 12 ]. The minimal length of the tetrahedrons was set to 2mm and maximal to 4mm. As a result, the mesh contained 7178 points and 26156 tetrahedrons. 5

Performance Evaluation

In order to estimate the performance and scalability of the LV simulation using the FEniCS framework, a series of experiments was performed. We used the Uran supercomputer of the Krasovskii Institute of Mathematics and Mechanics. The con guration parameters of the computational nodes are presented in Table 1. The FEniCS version 1.6.0 was used. The FEniCS framework demonstrated good performance and near-linear scalability due to the data parallelism. The mesh was distributed among the computational nodes before the launch of the simulation. Almost all computations were performed independently except for the transfer of boundary values of each mesh fragment between the nodes.

Number of CPU cores

Our previous manual implementation of the same model in the LeVen system [ 13 ], which uses the C language and OpenMP for parallelization, provides approximately the same performance. However, the FEniCS-based solution has better scalability: it scales up to 240 CPU cores while LeVen scales only up to 8 cores. In addition, the FEniCS implementation can be easily modi ed since it has near{mathematical notation.

Another problem we faced was the import of the model's description from CellML [ 14 ] into FEniCS. CellML is a language created within the Physiome Project to describe mathematical models of cellular biological functions in order to aid distribution and reuse of the models. We use the CellML description of the Ekaterinburg-Oxford model as a basis for our code. However, the standard tool for converting CellML descriptions to the simulation program from the Physiome Project generates non human-readable code and, therefore, is unsuitable for further use. We found a workaround by using the tools from the Gotran Project [ 15 ] and converting the CellML model description to the UFL language. However, we had to manually edit output les on each stage due to the complexity of models and the limitations of the just-in-time compilation algorithm of FEniCS. For example, we had to replace expressions such as x by e(ln( )x) because FEniCS does not support raising to a fractional power.

Conclusion and Future Work

The FEniCS framework is an e cient tool for heart simulation on parallel computing systems because it provides convenient near-mathematical notation, high simulation performance, and scales well. In comparison to our previous manual implementation FEniCS provides better scalability and can be easily utilized by biologists, chemists or physicists.

Possible directions of future work include: { Testing the scalability of FEniCS on thousands of CPU cores. { Applying FEniCS for real tasks such as simulation of scroll wave dynamics. { Evaluating the ability of FEniCS and other automated scienti c computing frameworks to use modern computational accelerators such as GPU and Xeon Phi.

{ Developing tools for automatic import of the CellML models to FEniCS. Acknowledgments. This work was supported by the Russian Science Foundation (grant no. 14-35-00005). Our study was performed using the Uran supercomputer of the Krasovskii Institute of Mathematics and Mechanics.

References

1. Kerckho s , R.C.P. , Healy , S.N. , Usyk , T.P. , McCulloch , A.D. : Computational methods for cardiac electromechanics . Proceedings of the IEEE 94 ( 2006 ) 769 { 783

2. Plank , G. , Burton , R.A. , Hales , P. , Bishop , M. , Mansoori , T. , Bernabeu , M.O. , Garny , A. , Prassl , A.J. , Bollensdor , C. , Mason , F. , et al.: Generation of histoanatomically representative models of the individual heart: tools and application . Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 367 ( 1896 ) ( 2009 ) 2257 { 2292

3. Jasak , H. , Jemcov , A. , Tukovic , Z. : OpenFOAM: a C++ library for complex physics simulations . In: International workshop on coupled methods in numerical dynamics. Volume 1000 . ( 2007 ) 1 { 20

4. Crampin , E.J. , Halstead , M. , Hunter , P. , Nielsen , P. , Noble , D. , Smith , N. , Tawhai , M. : Computational physiology and the physiome project . Experimental Physiology 89 ( 1 ) ( 2004 ) 1 { 26

5. Pitt-Francis , J. , Pathmanathan , P. , Bernabeu , M.O. , Bordas , R. , Cooper , J. , Fletcher , A.G. , Mirams , G.R. , Murray , P. , Osborne , J.M. , Walter , A. , et al.: Chaste: a test-driven approach to software development for biological modelling . Computer Physics Communications 180 ( 12 ) ( 2009 ) 2452 { 2471

6. Mirams , G.R. , Arthurs , C.J. , Bernabeu , M.O. , Bordas , R. , Cooper , J. , Corrias , A. , Davit , Y. , Dunn , S.J. , Fletcher , A.G. , Harvey , D.G. , et al.: Chaste: an open source C++ library for computational physiology and biology . PLoS Comput Biol 9 ( 3 ) ( 2013 ) e1002970

7. Aln s, M.S. , Blechta , J. , Hake , J. , Johansson , A. , Kehlet , B. , Logg , A. , Richardson , C. , Ring , J. , Rognes , M.E. , Wells , G.N. : The FEniCS project version 1.5. Archive of Numerical Software 3 ( 100 ) ( 2015 )

8. Solovyova , O. , Vikulova , N. , Katsnelson , L.B. , Markhasin , V.S. , Noble , P. , Garny , A. , Kohl , P. , Noble , D. : Mechanical interaction of heterogeneous cardiac muscle segments in silico: e ects on Ca 2+ handling and action potential . International Journal of Bifurcation and Chaos 13 ( 12 ) ( 2003 ) 3757 { 3782

9. Katsnelson , L.B. , Vikulova , N.A. , Kursanov , A.G. , Solovyova , O.E. , Markhasin , V.S. : Electro-mechanical coupling in a one-dimensional model of heart muscle ber . Russian Journal of Numerical Analysis and Mathematical Modelling 29 ( 5 ) ( 2014 ) 275 { 284

10. Pravdin , S. : Non-axisymmetric mathematical model of the cardiac left ventricle anatomy . Russian Journal of Biomechanics 17 ( 62 ) ( 2013 ) 75 { 94

11. Li , Y. , Chen , C. : An e cient split-operator scheme for 2-D advection-di usion simulations using nite elements and characteristics . Applied Mathematical Modelling 13 ( 4 ) ( 1989 ) 248 { 253

12. Geuzaine , C. , Remacle , J.F. : Gmsh: a three-dimensional nite element mesh generator with built-in pre- and post-processing facilities . International Journal for Numerical Methods in Engineering 79 ( 11 ) ( 2009 ) 1309 { 1331

13. Sozykin , A. , Pravdin , S. , Koshelev , A. , Zverev , V. , Ushenin , K. , Solovyova , O. : LeVen - a parallel system for simulation of the heart left ventricle . 9th International Conference on Application of Information and Communication Technologies, AICT 2015 Proceedings ( 2015 ) 249 { 252

14. Cuellar , A.A. , Lloyd , C.M. , Nielsen , P.F. , Bullivant , D.P. , Nickerson , D.P. , Hunter , P.J.: An overview of CellML 1.1, a biological model description language . SIMULATION 79 ( 12 ) ( 2003 ) 740 { 747

15. Hake , J.E. : A general ODE translator (gotran): Towards a versatile tool for general ODEs . CBC and CCI Workshop on Advancing Numerical Technologies in the Cardiac Domain, May 15 (May 2013 )