Parallel Left Ventricle Simulation Using the FEniCS Framework Timofei Epanchintsev1,2,3 and Vladimir Zverev3,1 1 Krasovskii Institute of Mathematics and Mechanics, Yekaterinburg, Russia 2 Institute of Immunology and Physiology of the Ural Branch of the Russian Academy of Sciences, Yekaterinburg, Russia 3 Ural Federal University, Yekaterinburg, Russia Abstract. Heart simulation is complex task that requires multiscale modeling on cell, tissue and organ levels. Such structure makes difficult high performance code development and its maintenance. In this paper, we evaluate how scientific software could be used for heart simulation. An overview of existing frameworks for automated scientific computing is presented. The FEniCS framework was chosen since it supports auto- mated solution of differential equations by the finite 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 propaga- tion of membrane potential alternation over a cardiac left ventricle using the electrophysiological model of a left ventricle and the Ekaterinburg- Oxford cell model. The FEniCS framework showed good performance and near-linear scalability up to 240 CPU cores. Keywords: FEniCS · heart simulation · parallel computing · finite ele- ment method 1 Introduction 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 acceler- ators such as GPU and Xeon Phi. However, porting a complicated multilevel simulation code to a new computational architecture requires a long time (of- ten 3–5 years) during which the architecture may become obsolete. In addition, adaptation for parallel computing architectures often leads to significant changes of code. Consequently, it is very difficult to determine which mathematical mod- els and numerical methods are used in the optimized code. Hence, reflecting 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 scientific computing frameworks. Such frameworks allow the development of simulation software using 30 Timofei Epanchintsev and Vladimir Zverev programming languages with near-mathematical notation. Traditionally, a sig- nificant disadvantage of such frameworks was their low simulation performance. However, some modern implementations use advanced tools such as highly ef- ficient mathematical libraries, just-in-time compilers, parallel execution, and so on, which improve their performance. Still, it is not clear if the modern auto- mated scientific computing frameworks are efficient enough for real multiscale simulation tasks such as heart simulation. In this paper, we evaluate how the automated scientific computing frame- works could be used for heart simulation on parallel computing systems. Heart simulation is an attractive task for evaluation of the performance of automated scientific 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 Scientific Computing Frameworks Nowadays, finite element method is the most popular method for numerical in- vestigation of systems with complex geometry, such as human heart. Several frameworks for automated scientific computing using finite element method ex- ist. 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 finite volume framework with solvers that are known to be efficient 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 applica- tion of finite element analysis techniques to a variety of complex bioengineer- ing problems. For distributed and parallel computing, OpenCMISS uses MPI, OpenMP, and ParMETIS.However, the OpenCMISS project is still under devel- opment, 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 aris- ing in biology and physiology [5, 6]. Current functionality includes the tissue and cell level electrophysiology, the discrete tissue modeling, and the soft tis- sue 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. Parallel Left Ventricle Simulation Using the FEniCS Framework 31 FEniCS [7] is a collaborative project for the development of the tools for automated scientific computing with a particular focus on the automated solu- tion of differential equations by the finite element method. Implementation of a finite 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 pro- vides the automated solution of differential equations by finite element methods, supports automatic parallel execution using MPI, and has a good documenta- tion. 3 Description of the Heart Model We simulated the electrical activity of a human heart, which is a result of spa- tial and temporal propagation of the electrical signal from each cardiac cell. We used the Ekaterinburg-Oxford cell model [8] with a detailed description of electri- cal, 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. Fig. 1. The schemes of the electric potential and some of the underlying ionic current 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 concentra- tions of calcium, sodium and potassium, extracellular potassium, the kinetics of calcium complexes, and calcium kinetics in the organelles. The mechanical block 32 Timofei Epanchintsev and Vladimir Zverev of the model includes equations describing the voltage of the cell that depends on its length. The differential equations of the electrical, chemical, and mechanical processes can be presented in the following simplified form: ∂S = g(V, S), (1) ∂t 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 defined 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-diffusion equation” [9]: ∂V = D∇2 V + Iions (V, S), (2) ∂t where D is the diffusion coefficient, ∇2 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 differential 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 first order (Marchuk–Yanenko method) [11]. Let us consider time domain t ∈ [0, T ] and the uniform grid tn = ht n, 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 first step, we solve the following partial differential equation: Vn∗ − Vn−1 ∗ = D∇2 V ∗ , V ∗ (t = tn−1 ) = V (tn−1 ), t ∈ [tn−1 , tn ]. (3) ht At the second step, we should solve the following system of ordinary equa- tions: Parallel Left Ventricle Simulation Using the FEniCS Framework 33 Fig. 2. An example of 3D mesh of left ventricle (asymmetric model) Vn∗∗ − Vn−1 ∗∗ = Iions , V ∗∗ (t = tn−1 ) = V ∗ (tn ), ht (4) Sn∗∗ − Sn−1 ∗∗ = g(V ∗ , S ∗∗ ), S ∗∗ (t = tn−1 ) = S(tn−1 ). ht 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 34 Timofei Epanchintsev and Vladimir Zverev configuration parameters of the computational nodes are presented in Table 1. The FEniCS version 1.6.0 was used. Table 1. Configuration of computational nodes Configuration parameter Value CPU 2 x Intel(R) Xeon(R) CPU X5675 @ 3.07GHz RAM 192 GB Interconnect Infiniband DDR (20 Gbit) Operating System CentOS Linux 7.2 200 180 160 The simulation time, min 140 120 100 80 60 40 20 0 16 24 32 40 48 56 64 72 80 88 96 8 104 112 120 132 156 168 180 192 204 216 228 240 144 Number of CPU cores Fig. 3. Simulation time depending on the number of CPU cores We estimated the scalability by conducting the simulation on varying num- bers of CPU cores, from 1 to 240. Fig. 3 shows the simulation time using different numbers of CPU cores and Fig. 4 presents the achieved speedup. 6 Discussion The FEniCS framework demonstrated good performance and near-linear scal- ability due to the data parallelism. The mesh was distributed among the com- putational 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. Parallel Left Ventricle Simulation Using the FEniCS Framework 35 100 90 80 The simulation speedup 70 60 50 40 30 20 10 0 16 24 32 40 48 56 64 72 80 88 96 1 8 104 112 120 132 156 168 180 192 204 216 228 240 144 Number of CPU cores Fig. 4. Simulation speedup depending on the number of CPU cores Our previous manual implementation of the same model in the LeVen sys- tem [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 modified 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 Phys- iome 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 files on each stage due to the com- plexity 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 efficient tool for heart simulation on parallel com- puting systems because it provides convenient near-mathematical notation, high simulation performance, and scales well. In comparison to our previous manual 36 Timofei Epanchintsev and Vladimir Zverev 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 scientific 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 Founda- tion (grant no. 14-35-00005). Our study was performed using the Uran super- computer of the Krasovskii Institute of Mathematics and Mechanics. References 1. Kerckhoffs, 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., Bollensdorff, C., Mason, F., et al.: Generation of histo- anatomically representative models of the individual heart: tools and application. Philosophical Transactions of the Royal Society of London A: Mathematical, Phys- ical 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: effects 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 fiber. Russian Journal of Numerical Analysis and Mathematical Modelling 29(5) (2014) 275–284 Parallel Left Ventricle Simulation Using the FEniCS Framework 37 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 efficient split-operator scheme for 2-D advection-diffusion sim- ulations using finite elements and characteristics. Applied Mathematical Modelling 13(4) (1989) 248–253 12. Geuzaine, C., Remacle, J.F.: Gmsh: a three-dimensional finite element mesh gen- erator 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. SIMU- LATION 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)