=Paper=
{{Paper
|id=Vol-1513/paper-08
|storemode=property
|title=An HPC-Based Approach to Study Living System Computational Model Parameter Dependency
|pdfUrl=https://ceur-ws.org/Vol-1513/paper-08.pdf
|volume=Vol-1513
|authors=Konstantin Ushenin,Dmitry Byordov
}}
==An HPC-Based Approach to Study Living System Computational Model Parameter Dependency==
https://ceur-ws.org/Vol-1513/paper-08.pdf
An HPC-Based Approach to Study Living
System Computational Model Parameter
Dependency
Konstantin Ushenin1,2,3 and Dmitry Byordov1
1
Ural Federal University, Yekaterinburg, Russia,
2
IMM UB RAS, Yekaterinburg, Russia,
3
IIP UB RAS, Yekaterinburg, Russia
kostaNew@gmail.com,berd gow@mail.ru
Abstract. High performance computing (HPC) allows one to run in
parallel large amount of independent numerical experiments for compu-
tationally intensive simulations of a complex system. Results of such ex-
periments can be used to derive dependencies between functional charac-
teristics of simulated system and parameters of the computational model.
In this paper, we implemented this HPC approach with using a computa-
tional 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 identified a dependence of the dynamics
of activated myocardium part during excitation on the model parameters
of the myocardial fiber orientation in the ventricular wall.
Keywords: high performance computing · parallel computing · heart
simulation
1 Introduction
Modern progress in systems biology and high performance computing (HPC)
have made it possible to perform computer simulations of the physiological func-
tions of complex living systems. Such simulation can be utilized to derive the
dependencies between hidden, unobservable in clinical data parameters and func-
tional 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 influence 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 ex-
periments for computationally intensive simulations are executed on HPC cluster
68 Konstantin Ushenin and Dmitry Byordov
in parallel. The simulation results are used to derive dependencies between func-
tional 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 different parameter values are executed on the HPC cluster.
To illustrate possibilities of the approach, we analyzed dependencies of elec-
trophysiological characteristics of the left ventricle (LV) of human heart on the
parameters of its geometry.
2 An HPC-based Approach to Study Computational
Model Parameter Dependency
Our approach is based on the execution of large amount of independent simula-
tion tasks on computational cluster (each task with different parameter values)
and statistical analysis of simulation results. We assume that simulation soft-
ware 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 different pa-
rameter values are executed in parallel (embarrassing parallelism).
4. Performing statistical analysis of experiment results. The results are repre-
sented in the form y j = f (xj1 , xj2 , .., xjN ), where f is the dependency between
parameters and functional characteristic. To determine the form and coeffi-
cients of the dependencies we use the regression analysis method.
This research approach gives us some advantages in comparison with qual-
itative descriptions of experiments. First, the obtained result is a function. A
type of this function presents a hypothesis regarding the process behavior, and
its coefficients can be used for distinguishing criterion that separate normal and
pathological cases. In addition, through the principal component analysis, we can
find physiological parameters the mostly affecting on processes. Second, the ap-
proach provides the possibility to restore hidden parameters of the physiological
system from clinical data.
3 Implementation
To be able to use the suggested approach, we implemented the software sys-
tem that executes tasks on the computational cluster, collects the results, and
Living System Computational Model Parameter Dependency 69
Fig. 1: Architecture of system with MPI or SLURM
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 effect on other tasks. Furthermore, the SLURM’s task distribution is
more flexible 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 file 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 process-
ing, for maximal flexibility of our research.
4 Computational Experiments
In order to illustrate the possibilities of our approach, we analyzed dependences
of electrophysiological characteristics of the LV on the parameters of its ge-
ometry. Particularly, we identified a dependence of the dynamics of activated
myocardium part during excitation on the model parameters of the myocardial
fiber orientation in the ventricular wall.
70 Konstantin Ushenin and Dmitry Byordov
Fig. 2: SLURM and Python wrapper above C (C99) code.
4.1 LV Models
The research is based on the continuum description of the LV tissue. The sim-
ulation includes an anatomical model for description of the shape and the fiber
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 simplified anatomical model, which was devel-
oped in our group earlier [3]. The anatomical shape parameters are fixed: LV
height Zb = 60 mm, wall thickness on apex h = 8 mm, radius from the axis to
epicardium on base Rb = 33 mm, wall thickness on base Lb = 12 mm, and cur-
vature ε = 0.85. Only parameters γ0 , γ1 are changed. This parameters influence
on fiber direction in the LV. In addition, the anatomical model provide a method
for generate a structured grid. This grid is used for calculate electrophysiological
differential equation by Euler’s method.
Electrophysiological model: We use the Aliev-Panfilov [1] model for simulation
of the electrophysiological activity of the heart myocardium at the cellular level.
The model equations are represented as follows:
u̇ = −ku(u − a)(u − 1) − uv + ∇ · (D∇u)
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 diffusion tensor:
D11 = 12, D22 = 1.33, D12 = D21 = 0.0. Thus, the excitation wave propagation
speed along the fibers 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 [4].
In our experiments, we used an existing implementation of these models [6]
in C language with OpenMP support.
Living System Computational Model Parameter Dependency 71
4.2 Experiment Protocol
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 fibers has influence 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 ex-
periment. These parameters determine the fiber direction in the LV. Parameters
of influence on the ventricle’s shape stayed fixed.
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 ∈ [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 Results
After computation of the Xtime and Yvol sequences, we use the ordinary least
squares (OLS) method. The Levenberg–Marquardt algorithm was used for find-
ing the minimum of the OLS discrete functional. We tried to use polynomials,
logarithmic, sedate, and exponential regression for function approximation. Fi-
nally, we have selected the 4-th degree polynomial (1). The value p0 is fixed on
first activated points group’s level: in this case it matches the initial activation.
The results of approximation are shown in Fig. 3.
Vact (t) = p4 t4 + p3 t3 + p2 t2 + p1 t1 + p0 (1)
For coefficients p0 , p1 , p2 , p3 , we searched dependence on the initial change-
able 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 + B0i (2)
Approximation by the polynomial (2) is the most simple and accurate. Re-
sults are shown on Fig.(4). The final formulas are represented below in the
Results section.
We insert coefficients of equation (2) into (3) to get the final equation (3),
that describes the activated myocardium volume dependence. This equation is
based on 55 experiments, and has the following form:
72 Konstantin Ushenin and Dmitry Byordov
(a) (b)
Fig. 3: Example of approximation of dependence the LV activated part from time
using the 4-th degree polynomial; (a): γ0 = 0.0, γ1 = 1.0 (b): γ0 = 5.0, γ1 = 7.0
Vact (t) = t4 (−0.0492γ0 γ1 + 0.0013γ0 + 0.0442γ1 − 0.0056)+
t3 (1.1161γ0 γ1 − 0.0487γ0 − 1.0192γ1 + 0.1559)+
t2 (−5.1988γ0 γ1 + 0.5706γ0 + 4.7070γ1 − 1.0035)+
t(12.9829γ0 γ1 − 0.8279γ0 − 11.8881γ1 + 4.7832)+
1.6822
(3)
There are two ways to verify the obtained results. For the first method,
the researchers can measure the real shape and the activation map of a mam-
malian 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 verification of the dependence we obtained. The second method of the
verification is similar to the first, 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 fiber direction. Thus, we can
construct the anatomical model [3].
5 Related Work
The usage of statistical processing of the patients clinical data is widespread
in medicine [5]. Our research differs from them by using numerical experiments
with mathematical models instead of clinical data.
The parameters similar to Vact (t) is used in work [7]. Authors of this work
use volume of activated myocardium as intermediate equation in mathematical
transformation.
Living System Computational Model Parameter Dependency 73
(a) (b)
(c) (d)
Fig. 4: Approximation of the coefficients p0 , p1 , p2 , p3 with polynomial B3i γ1 γ0 +
B2i γ1 + B1i γ0 + B0i . Images (a), (b), (c), and (d) correspond to coefficients B0 ,
B1 , B2 , B3 , respectively.
Large number of experiments for dependencies research is used in paper [2].
The authors of this paper studied structural differences in the myocardium of
different mammal species. In our work we are focused on another physiological
problem and suggest technical implementation of its decision.
6 Conclusion
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
different 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 mod-
eling, we studied the dependence of the dynamics of activated myocardium part
during excitation on the model parameters of the myocardial fiber 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.
74 Konstantin Ushenin and Dmitry Byordov
Acknowledgments. The research was supported by the grant under the Pro-
gram 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 clus-
ter of Ural Federal University.
References
1. Aliev, R.R., Panfilov, A.V.: A simple two-variable model of cardiac excitation.
Chaos, Solitons and Fractals 7(3), 293–301 (1996), www.scopus.com
2. Bueno-Orovio, A., Kay, D., Grau, V., Rodriguez, B., Burrage, K.: Fractional dif-
fusion models of cardiac electrical propagation: role of structural heterogeneity in
dispersion of repolarization. Journal of The Royal Society Interface 11(97) (2014)
3. Pravdin, S.F., Berdyshev, V.I., Panfilov, A.V., Katsnelson, L.B., Solovyova, O.,
Markhasin, V.S.: Mathematical model of the anatomy and fibre orientation field of
the left ventricle of the heart. BioMedical Engineering Online 12(1) (2013), www.
scopus.com
4. Pravdin, S.F., Dierckx, H., Katsnelson, L.B., Solovyova, O., Markhasin, V.S., Pan-
filov, A.V.: Electrical wave propagation in an anisotropic model of the left ventricle
based on analytical description of cardiac architecture. PLoS ONE 9(5) (2014),
www.scopus.com
5. Riffenburgh, R.: Statistics in Medicine. Elsevier Science (2012), https://books.
google.ru/books?id=Pd4KCgJeXeEC
6. Sozykin, A., Pravdin, S., Koshelev, A., Zverev, V., Ushenin, K., Solovyova, O.:
Leven – a parallel system for simulation of the heart left ventricle. In: Application of
Information and Communication Technologies (AICT), 2015 IEEE 9th International
Conference on. pp. 249–252 (2015)
7. Young, R.J., Panfilov, A.V.: Anisotropy of wave propagation in the heart can be
modeled by a Riemannian electrophysiological metric. Proceedings of the National
Academy of Sciences of the United States of America 107(34), 15063–15068 (2010)