High performance computing (HPC) allows one to run in parallel large amount of independent numerical experiments for computationally intensive simulations of a complex system. Results of such experiments can be used to derive dependencies between functional characteristics of simulated system and parameters of the computational model. In this paper, we implemented this HPC approach with using a computational model of the electrical activity in the left ventricle of human heart. To illustrate possibilities of the approach, we analyzed dependencies of electrophysiological characteristics of the left ventricle on the parameters of its geometry. Particularly, we identi ed a dependence of the dynamics of activated myocardium part during excitation on the model parameters of the myocardial ber orientation in the ventricular wall.
heart Modern progress in systems biology and high performance computing (HPC) have made it possible to perform computer simulations of the physiological functions of complex living systems. Such simulation can be utilized to derive the dependencies between hidden, unobservable in clinical data parameters and functional characteristics of the systems. These dependencies can be used for solving the inverse problem and make the clinical diagnostics better.
However, deriving such dependencies for computational models of complex living system is complicated process because simulation requires long time and contemporary models have large number of parameters that in uence functional characteristics. To solve this problem, we developed an HPC-based approach to study the parameter dependency of living system models.
According to suggested approach, large amount of independent numerical experiments for computationally intensive simulations are executed on HPC cluster in parallel. The simulation results are used to derive dependencies between functional characteristics of simulated system and parameters of the computational model.
The approach could be used even for the simulation software that is not support parallel computing, because large amount of independent simulation tasks with di erent parameter values are executed on the HPC cluster.
To illustrate possibilities of the approach, we analyzed dependencies of electrophysiological characteristics of the left ventricle (LV) of human heart on the parameters of its geometry. 2
Our approach is based on the execution of large amount of independent simulation tasks on computational cluster (each task with di erent parameter values) and statistical analysis of simulation results. We assume that simulation software for the investigated computational model is already exists. In addition, we assume that simulation takes a few hours or days. Hence, the parallel computing support for one simulation task is not necessary.
The approach consists of the following steps: 1. Choosing a set of model parameters that we plan to change and determining the functional characteristics that are dependent on those parameters. 2. Determining the domain of parameters values, which will be used during the experiments, and preparing the input data for the simulation software. 3. Performing simulation on the HPC cluster. Large amount of independent computational tasks using existing simulation software with di erent parameter values are executed in parallel (embarrassing parallelism). 4. Performing statistical analysis of experiment results. The results are represented in the form yj = f (xj1; xj2; ::; xjN ), where f is the dependency between parameters and functional characteristic. To determine the form and coe cients of the dependencies we use the regression analysis method.
This research approach gives us some advantages in comparison with qualitative descriptions of experiments. First, the obtained result is a function. A type of this function presents a hypothesis regarding the process behavior, and its coe cients can be used for distinguishing criterion that separate normal and pathological cases. In addition, through the principal component analysis, we can nd physiological parameters the mostly a ecting on processes. Second, the approach provides the possibility to restore hidden parameters of the physiological system from clinical data. 3
To be able to use the suggested approach, we implemented the software system that executes tasks on the computational cluster, collects the results, and provides tools for statistical analysis and visualization. The architecture of the system is shown in Fig. 1.
Two types of parallel execution of the simulation tasks on the cluster are implemented: based on MPI and SLURM resource manager tools.
The MPI implementation is based on the master-slave concept with load balancing over task queue. This implementation shows itself to be well, but when applying in more complicated experiments, cause several problem. When some combination of parameters leads to insoluble problem formulations the simulation software quits unexpectedly. This, in turn, can lead to an unexpected end of all MPI tasks. We do not describe this complex experiment in the paper.
The SLURM implementation is the most convenient, as an error in one task has no e ect on other tasks. Furthermore, the SLURM's task distribution is more exible in comparison with the MPI implementation.
For reduction of the output data, the simulation software must be compiled as a library and wrapped in Python with module ctypes. Result arrays are saved in a hard drive in compressed binary le through standard numpy command savez compressed. The wrapper structure is shown in Fig. 2.
After this, a data is moved into the machine with GUI and is processed by modules numpy and scipy.
This approach allows us to divide numerical experiments and results processing, for maximal exibility of our research. 4
In order to illustrate the possibilities of our approach, we analyzed dependences
of electrophysiological characteristics of the LV on the parameters of its
geometry. Particularly, we identi ed a dependence of the dynamics of activated
myocardium part during excitation on the model parameters of the myocardial
ber orientation in the ventricular wall.
The research is based on the continuum description of the LV tissue. The
simulation includes an anatomical model for description of the shape and the ber
direction, and an electrophysiological models for description of the excitation
properties of the LV tissue. The Euler method is used for calculation.
Anatomical model: We used a simpli ed anatomical model, which was
developed in our group earlier [
Electrophysiological model: We use the Aliev-Pan lov [
uv + r (Dru) v_ = (u)(8u v) where (u) = 0:1 if u > a and (u) = 1 if u 1; k = 8; a = 0:03. For introducing the material anisotropy, we used a symmetrical di usion tensor: D11 = 12; D22 = 1:33; D12 = D21 = 0:0. Thus, the excitation wave propagation speed along the bers is 3 times faster than across them.
Initial activation The model does not account the conduction system of the
heart. We used the initial activation on endocardium apex applying the model
coordinates: 0:7 2 ; 0:6. This solution is based on the previous work [
In our experiments, we used an existing implementation of these models [
During the experiment, we measured a part of the excited myocardium of the LV depending on the time for every model. For every point of the electrical mesh, we calculated the corresponding myocardium volume and the arrival time of the excitation wave. The part of the excited myocardium Vact(t) is a percent of the activated LV before the instant t. The choice of this parameter is connected with the fact that the transmural angle of bers has in uence on this parameter. In practice, it can be calculated through the activation of the ventricle's map during the physiological experiment or clinical electrophysiological study of the heart. Moreover, Vact(t) may be interpreted as the index of myocardium activation coherence.
We selected parameters 0 and 1 for changing in the series of numerical experiment. These parameters determine the ber direction in the LV. Parameters of in uence on the ventricle's shape stayed xed.
For calculation of Vact(t), the time interval of the wave spread [0; Tmax] is divided into 10,000 parts. For every sub-interval [tk; tk + 1), we calculated the activated volume as sum of activated point volumes. The activated myocardium volume at the instant t is calculated as the sum of volumes in groups t1; ::; tk + 1, if t 2 [tk; tk + 1]. Thus, the arrival time of waves into activated points group is equal to average arrival time of wave in these points. So, we get two sequences of data 10,000 elements in each. Xtime is an average time of point groups' activation. Yvol is the volume of activated myocardium in time moment Xtime. 4.3
After computation of the Xtime and Yvol sequences, we use the ordinary least squares (OLS) method. The Levenberg{Marquardt algorithm was used for nding the minimum of the OLS discrete functional. We tried to use polynomials, logarithmic, sedate, and exponential regression for function approximation. Finally, we have selected the 4-th degree polynomial (1). The value p0 is xed on rst activated points group's level: in this case it matches the initial activation. The results of approximation are shown in Fig. 3.
Vact(t) = p4t4 + p3t3 + p2t2 + p1t1 + p0 (1)
For coe cients p0; p1; p2; p3, we searched dependence on the initial changeable parameters 0 and 1. By analogy with previous variants, we use the least square method and the Levenberg { Marquardt algorithm.
pi = B3i 1 0 + B2i 1 + B1i 0 + B0 i (2)
Approximation by the polynomial (2) is the most simple and accurate. Results are shown on Fig.(4). The nal formulas are represented below in the Results section.
We insert coe cients of equation (2) into (3) to get the nal equation (3), that describes the activated myocardium volume dependence. This equation is based on 55 experiments, and has the following form:
There are two ways to verify the obtained results. For the rst method, the researchers can measure the real shape and the activation map of a mammalian heart. After that, they can run the simulation with the parameters of the measured shape and obtain the activation map. These results can be used for the veri cation of the dependence we obtained. The second method of the veri cation is similar to the rst, but the researchers can use data from the electrophysiological study.
Inverse functions cannot be obtained for the deduced formulas in strictly
mathematical terms. Nevertheless, with the activation map data, the result still
can be used for creating hypotheses about the ber direction. Thus, we can
construct the anatomical model [
The usage of statistical processing of the patients clinical data is widespread
in medicine [
The parameters similar to Vact(t) is used in work [
Large number of experiments for dependencies research is used in paper [
We presented and approach to study the dependency of computational model on parameters that is based on an execution of large amount of simulations with di erent parameter values on the HPC cluster. To derive dependencies between functional characteristics of simulated system and parameters, the statistical analysis of simulation results is used.
To provide an example of application of our approach to living system modeling, we studied the dependence of the dynamics of activated myocardium part during excitation on the model parameters of the myocardial ber orientation in the ventricular wall.
Future works include a study of the dependencies between geometric and electrophysiological parameters of the LV to the dynamics of spiral waves. Acknowledgments. The research was supported by the grant under the Program of the Presidium of RAS II.4P \Fundamental program of mathematical modeling" no. 22. Our work was performed using the \Uran" supercomputer from Institute of Mathematics and Mechanics UrB RAS and computational cluster of Ural Federal University.