=Paper=
{{Paper
|id=Vol-3132/Short_8.pdf
|storemode=property
|title=Construction of a Mathematical Model of the Heart Muscle
|pdfUrl=https://ceur-ws.org/Vol-3132/Short_8.pdf
|volume=Vol-3132
|authors=Oleksii Byckov,Evgen Gurko,Denys Khusainov,Andriy Shatyrko,Bedrich Puza,Veronika Novotna
|dblpUrl=https://dblp.org/rec/conf/iti2/BychkovGKSPN21
}}
==Construction of a Mathematical Model of the Heart Muscle==
https://ceur-ws.org/Vol-3132/Short_8.pdf
Construction of a Mathematical Model of the Heart Muscle
Oleksii Byckova, Evgen Gurkoa, Denys Khusainova, Andriy Shatyrkoa, Bedrich Puzab
and Veronika Novotnab
a
Taras Shevchenko University of Kyiv, 64, Volodymyrska str., Kyiv, 01033, Ukraine
b
Brno University of Technology, 4, Koleyni, Brno, Index, Czech Republic
Abstract
According to statistics, the problem of cardiovascular disease ranks first among diseases [1-3].
The purpose of this work is to build a mathematical model of the heart muscle, which will more
accurately reflect the behavior of the heart, for example in pathological contractions, in which
the reduction of arbitrary areas of arbitrary fibers or external mechanical impact. The main idea
of the model of the ventricle and the outer shell of the heart is the division of the surface into
longitudinal flexible fibers connected by elastic ligaments. The model is presented in the form
of an elastic system, which is described in terms of a system of differential equations [4-8].
They simulate the contraction forces of the heart muscle. A static mathematical model was also
used, which representing the left ventricle and the outer surface of the heart as a paraboloid of
rotation and an elliptical paraboloid, respectively, taking into account the spiral course of
muscle fibers. Euler's method was used to numerically solve the system of differential
equations. Based on the above models, software was created that reflects the dynamic behavior
of the heart muscle. In the course of the work it was possible to select such parameters at which
the dynamics of the model became close to the dynamics of the real heart muscle and other
models of heart work.
Keywords 1
Mathematical model, differential equation system, software, medicine, heart muscle
1. Introduction
In modern life, the problem of cardiovascular disease is quite acute - according to statistics, they
have the first place by cause of death. Therefore, there are constant attempts to improve the situation in
this direction. One method of improvement is to build a mathematical model of heart contraction.
Solving this problem will help to more accurately predict the work of the heart, which is important in
crisis situations (eg, pre-infarction), in post-crisis rehabilitation, with medication, will detect problems
in the early stages and more [9-12]. The main purpose of this work is to build a mathematical model of
the heart muscle, which will more accurately reflect the behaviour of the heart, for example, in
pathological contractions, in which the reduction of arbitrary areas of arbitrary fibers or external
mechanical impact. Of course, scientists has tried to solve this problem many times, but there are still
no exact mathematical models of heart muscle contraction [13-16]. The fact is that it has always been
believed that the heart contracts like an enema, that is, just pushes blood. In fact, recent studies have
shown that in this case, the aneurysm occurs rapidly, and the contraction is spiral, and the blood already
enters the vessels in a twisted form. Thus, the task of modeling the behavior of the heart muscle under
the condition of spiral contraction and the ability to simulate the disconnection of parts of the heart from
the contraction process (in the future).
Information Technology and Implementation (IT&I-2021), December 01-03, 2021, Kyiv, Ukraine
EMAIL: oleksiibychkov@knu.ua (A. 1); bos.knu@gmail.com (A. 2); d.y.khusainov@gmail.com (A. 3); shatyrko.a@knu.ua (A. 4);
puza@fbm.vutbr.cz (A. 5); novotna@fbm.vutbr.cz (A. 6)
ORCID: orcid.org/0000-0002-9378-9535 (A. 1); 0000-0001-5855-029X (A. 3); 0000-0002-5648-2999 (A. 4); 0000-0002-2949-4708 (A. 5);
0000-0001-9360-3035 (A. 6)
©️ 2022 Copyright for this paper by its authors.
Use permitted under Creative Commons License Attribution 4.0 International (CC BY 4.0).
CEUR Workshop Proceedings (CEUR-WS.org)
272
2. Formulation of the problem
Using the results of studies of the structure of the heart muscle and the basics of the geometry of
contraction of the left ventricular muscle and mathematical models based on them, a new mathematical
model is built to reproduce the physical behavior of the heart muscle in spiral contraction. The
possibility of exclusion from reduction of some sites of a myocardium is considered.
The paper uses the results of research by the Institute of Cardiology N.D. Strazhesko Academy of
Medical Sciences of Ukraine and scientific works of Taras Shevchenko National University of Kyiv.
2.1. Model of the ventricle and the outer shell of the heart as an elastic
surface
The main idea of the model of the ventricle and the outer shell of the heart is to break the surface
into longitudinal flexible fibers, which are connected by elastic ligaments. This most fully reflects the
physical nature of the muscles. Based on the results of previous research and mathematical models
based on them, it is possible to build a new model, which is based on the "fibrous" nature and more
fully reflects the behavior of the muscle in different situations.
2.2. Cardiac muscle fiber model as an elastic chain
The fiber model is as follows. Let the fiber be part of some curve in space. On this curve, points are
sequentially selected (Fig. 1). They are connected in series by segments (ribs), which will act as
connections between them. These connections themselves have spring-like properties. That is, when
you change the length of the rib, the force that tries to return to the original length of the connection
arose. Also, every two consecutive connections create an elastic system. When you change the position
of one rib relative to another, forces that try to return them to their original position arose.
Figure: 1 Geometric model of cardiac muscle fiber
Each point in this chain is a material point with mass m i . Suppose that at the moment of time t ,
the i -th point of the chain has a velocity vi (t ) , and it is acted upon by the forces of elasticity of the
bonds Fiu (t ) and Fi d (t ) co-directed with the corresponding edges connected with the given point. Also
at this point acts the force of preservation of the relative position Ri (t ) (Fig. 2). Let the equivalent of
these forces is
Fi (t ) Fiu (t ) Fid (t ) Ri (t ) (1).
Then you can write the following
dvi Fi (t )
(2)
dt mi
273
Figure: 2 Distribution of elastic and preservation forces acting on the point of the heart muscle fiber
The problem is that the value Fi (t ) becomes known only at a time t that corresponds to the
condition of "autonomy" of the system, the possibility of changing parameters for the operating system
over time, including the sudden exclusion (included) of muscle areas, and so on. Therefore, we find the
numerical solution of equation (2) for each i -th point by Euler's method.
We can construct the system of equations for all points:
dv1 F1 (t )
dt m
1
dv2 F2 (t )
dt m2 (3)
..................
dvn Fn (t )
dt m
n
We will show how to search Fi j (t ) and Ri (t )
Fi j (t ) k ij (t )(l ij (t ) l0ij (t )) , (4)
ij ij
where k (t ) - the coefficient of elasticity of this edge, l (t ) - the length between the i -th and j -th
ij
point at the time t ; l 0 (t ) - the length of the rib at rest between the i -th and j -th point.
For example, Fid (t ) Fii 1 (t ) , Fiu (t ) Fii 1 (t ) , if there exist, else Fi j (t ) 0
Ri (t ) k ri (t )(( pi pu , pi p d ) 0i (t ))) (5)
where k ri (t ) - the coefficient maintaining the relative position of the edges associated with the i -th
point; 0i (t ) - the output angle between the edges connected to the i -th point, the angle to which the
relative position of the edges goes.
Let is denote chain length as
n 1
L(t ) l i ,i 1 (t ) (6)
i 1
flexible chain as a triple G(t ) ( P(t ),V (t ), M ) , where P(t ) - position of the points of the chain
pi (t ) P(t ) , pi (t ) ( xi (t ), y i (t ), z i (t )) ; V (t ) - speed of chain points vi (t ) V (t ) ; M - mass of points.
Thus, we obtained a model of dynamic structure that tries to reach equilibrium and can change both
its own parameters (stiffness, length, etc.) and respond to external factors. We use it to build a more
complex structure that simulates the heart muscle.
2.3. The model of the heart muscle as a set of fibers
We represent a muscle as a set of longitudinal fibers connected by elastic ligaments. Let the surface
be s chains ( G 1 , G 2,, G s ). Then to each i -th point of the j -th chain, pij , you need to add "external"
connections connecting the points of adjacent chains. Let each point be connected to only one point of
the adjacent chain. Let two such edges be formed.
274
Then additional forces will appear Fijl (t ) and Fijr (t ) , even those that will hold the connected points
of adjacent chains, and will change the force Ri (t ) responsible for maintaining the position of the
edges relative to each other will change (Fig. 3).
Figure: 3 Distribution of forces acting on the point of the heart muscle fiber, taking into account the
set of longitudinal fibers
To this force must be added the force responsible for maintaining the relative position between the
edges connecting adjacent chains and belonging to the i -th point of the j -th chain (similarly to (5)).
Then, the equilibrium force will change
Fij (t ) Fiju (t ) Fijd (t ) Fijl (t ) Fijr (t ) Rij (t ) (7)
In general, the scheme will look like this
dv11 F11 (t )
dt m11
dv F (t )
12 12
dt m12 (8)
..................
dv1 F1 (t )
s s
dt m1s
2.4. Euler's method for solving the problem
We show how to find vij by the system (8). Let i -th point of the j -th chain have an initial velocity
vij (t 0 ) . Then we get the Cauchy problem, which can be solved by Euler's method.
h
vij (t h) vij (t ) Fij (t ) , (9)
mij
where h - step of method.
We modify the method, taking into account the distribution (scattering) of energy and, as a
consequence, the attenuation of velocity over time.
h
vij (t h) r vij (t ) Fij (t ) , (10)
mij
where r 1 - attenuation coefficient.
Also by Euler's method we find the position of the node at the time t
pij (t h) pij (t ) hvij (t ) (11)
275
2.5. The structure of a muscle as an elastic system of nodes
It is important to define the structure of the muscle, ie all the connections in the fibers and between
them. To do this, we use a static mathematical model that represents the left ventricle and the outer
surface of the heart as a paraboloid of rotation and an elliptical paraboloid, respectively, taking into
account the spiral course of muscle fibers. Based on the results of the study of the geometry of the
contraction of the left ventricle, the basis of his mathematical model was based on the paraboloid of
rotation y k ( x 2 z 2 ) or in the coordinate system ( , , y) , y k 2 , k 0 , which is bounded by
the plane y H on axis Oy . The fibers, according to the chosen model, lie completely on the surface
of the paraboloid and are spiral-shaped spatial curves, the beginning of which coincides with the center
O(0,0,0) , and the end belongs to the bounding plane y H .
A family of logarithmic spirals having equations in polar coordinates of the form
exp( l ( 0 )), l ln was taken as projections of fibers on a plane O . The fact that they
often occur in various natural phenomena speaks in favor of the choice of logarithmic spirals. From the
whole family we will choose a curve with a phase 0 0 . Based on these data, we can show that,
tg
l , or tg 4kl - the tangent of the angle is directly proportional to the current radius. The
4k
logarithmic spiral has only one drawback for this case - asymptoticity at a point O(0,0,0) . But this
problem is easy to solve, assuming that 0 -around the point O(0,0,0) of the fibers go into the vortex
and go to the opposite wall of the paraboloid.
Based on research, the outer shell of the heart can also be modeled as a paraboloid, but no longer
circular, but elliptical: y k x x 2 k z z 2 , ( k x 0 , k z 0 ), for simplicity, we believe that k x k z . The
fibers of the outer shell can no longer be projected by logarithmic spirals and at the same time rise
monotonically upwards. Therefore, the model is based on the calculation that before all calculations, an
affine transformation of space is performed - scaling along the axis Ox (multiplication by a factor
kz kx
), and after them - the inverse transformation of scaling (multiplication by a factor ). Thus,
kx kz
after scaling, the paraboloid becomes a body of rotation, and the fibers in the projection will have the
form of logarithmic spirals, after which it will be possible to apply the above modeling methods.
The right ventricle is modeled as a cavity between the outer surface of the heart and the left ventricle.
This is because it is much weaker than the left ventricle.
Based on the above model, it is possible to construct curves on the surface of paraboloids that will
correspond to the fibers, and on their basis to construct the structure of chains, and link them together.
On the surface of the paraboloid we choose the s curves that come out of the point O and give
logarithmic spirals in the projection. Let these curves divide the slice of the paraboloid y h
( 0 h H ) into equal s parts. We will consider these curves as fibers of a heart muscle. Draw n slices
y hi , ( 0 hi H , i 1, n ).The points pij of intersection of the j -th curve and this slice - hi , choose
as nodes on the j -th curve. Since the connections between nodes (points) in the chain do not differ
from the connections between nodes from different chains, so we represent all connections between
nodes as a connected graph, the vertices of which are points pij , and the edges are connections between
them. We also add a point O to the set of vertices of the graph and denote it p O0 . Add links ( pO0 , p1 j )
, j 1, s to the set of edges. For everyone j 1, s add the edges that connect the points in the j -th
fiber ( pi j , pi 1 j ) , i 1, (n 1) . Also, add the ribs that connect the nodes of adjacent fibers (in a circle)
– to i 1, n add the ribs ( pi j , pi j 1 ) , j 1, s and the rib ( pis , pi1 ) that "loops" the cross section. The
276
resulting graph defines the structure of the muscle, and formulas (6), (8), (9) define the behavior of this
structure over time
2.6. Muscle contraction control
One of the important aspects of the system is the control of muscle contraction. Namely, setting
rules that will allow this system of nodes to reduce or increase the volume of the internal cavity. This
ij
control is performed by setting the functions k ij (t ) - the stiffness coefficients of the rib and l 0 (t ) -
the length of the rib, which it tries to achieve (i, j {1..n} {1..s} {(0,0)}) . As you increase k ij (t ) and
ij
decrease l 0 (t ) , there is a force that tries to reduce the length of the rib. Thus, if this behavior is set for
all edges, it will reduce the linear size of the system, and thus reduce the volume of the cavity, which is
limited by this set of edges in space. The above steps can be used to simulate muscle contraction. To
ij
simulate its relaxation, it is necessary to increase l 0 (t ) and gradually decrease k ij (t ) at the right time,
which will correspond to a slower relaxation than reduction, as it actually happens.
We need to find this point in time when the heart needs to relax. Studies have shown that heart fibers
reduce their linear size by about 15%. That is, when, for example, the average fiber length reaches 85%,
s
L j (t )
j 1
the heart will begin to relax. That is 0.85 L j (t0 ) , where L j (t ) is the length of the j -th chain
s
(fiber). When the length of the fibers again reaches approximately the original size, there is a reduction.
In fact, this process is much more complicated, but for simulation so far this approach has been chosen
ij
for simplification. In the future, this process can be complicated by the choice of others l 0 (t ) and
k ij (t )
2.7. Realization
Based on the above-described mathematical model of the elastic system, software was created that
reflects this model and its behavior.
For the internal representation of the graph of connections between nodes, an incidence list was used
(ie a list in which for each vertex the vertices incident to it are stored), which allowed to calculate each
step of Euler's method for time O(wns) , where w is the number of edges leaving the node, and, since
in this case w not large (except for the node p 0 does not exceed 4 ), we can assume that the time is
proportional O(ns) . This is achieved by calculating the next step of the Euler method to go through all
ns 1 nodes, and calculate the forces arising on the edges and the force of maintaining the relative
position. The C# programming language and the Microsoft.NET platform were used to write this
software application, which significantly accelerated the design and development.
2.8. Selection of model parameters
One of the most important parts of modeling is the correct selection of parameters and functions, so
that the simulated process as closely as possible reflects the actual behavior of the muscle.
The parameters of the static model of the left ventricle and the outer shell of the heart are taken from
previous work on the model with paraboloids [1].The coefficients in the equation of the outer shell of
the heart are equal to k x 0.5, k z 1.0 . The coefficient in the equation of the left ventricle is k 1.5
. The center of the left ventricle is displaced relative to the center of the outer shell by a magnitude
(0.9;0.5;0.0) . The height of the outer shell H out 9.0 , the height of the left ventricle - H in 8.5
respectively. The initial coefficients for the equation of logarithmic spirals are: for the outer shell - 2.08,
for the left ventricle - 1.45. At such values, the fibers of the outer shell make one curl, and the left
ventricle - two (as it really is).
277
In the course of many experiments with system parameters, it was empirically found that the
following values and functions are the most suitable:
i i
n 20 , s 20 , y hi H ( ) 2 , mij 0.00012 , k ij (t ) 25 , k i r (t ) 5 , h 0.001 , r 0.6 ,
n n
0.7𝑙0 𝑖𝑗 (𝑡0 ),when decreasing
𝑙0 𝑖𝑗 (𝑡) = {
1.2𝑙0 𝑖𝑗 (𝑡0 ), when relaxing
Figure: 4. The surface of the heart is built by the program.
3. Results and comparisons
This model was compared with the model of the ventricle and the surface of the heart, which
represent them as paraboloids, in which the coefficients at the coordinates and angles of inclination of
the fibers change over time. Visually, the work of both models is very similar [17-21].
The following method was used to calculate the volume of the cavity bounding this system of rib.
We draw an axis through a point pO0 and a parallel axis Oy . Select two points on this axis A0 pO0
s
pnj
j 1
and An . Actions on points mean actions on the corresponding coordinates. The segment
s
A0 An is divided into n successive parts. Let the points be chosen Ai , i 1, n 1 . Then the volume of
this figure can be divided into n volumes bounded by the following sets:
Vgi (t ) V ({ pi j | j 1..s} Ai { pi1 j | j 1..s} Ai1 ) (12)
where pOj pO0 , j 1..s .
The total volume is as follows
n
V g (t ) V gi (t ) (13)
i 1
According to the given formulas and parameters found in the previous section, the volume dynamics is
found, which is presented in Fig.5 and Fig.6.
4. Conclusion
In the course of this work, a mathematical model of the heart muscle was built, which is quite close
in physical content to its actual structure. Muscle is represented as a set of elastic fibers that are
interconnected by elastic ligaments and together form an elastic surface, the behavior of which can be
controlled to some extent. The surface consists of nodes with a certain mass and elastic connections
between them. The system seeks to achieve a position of equilibrium relative to the internal forces that
arise when changing control parameters or external action on the surface.
278
Figure 5: Graph of the inner volume of the outer shell of the heart
Figure 6: Graph of the internal volume of the left ventricle of the heart
In the course of work it was possible to select such parameters that the dynamics of the model is
close to the dynamics of the real heart muscle and other heart models.
Building a mathematical model of the heart muscle will be widely used in cardiology. The results of
the work can be used to build a model of pathological contraction, in which arbitrary areas of arbitrary
fibers fall out of the contraction. Also, one of the ways to improve the model is to set the initial structure
based on the real geometric parameters of the heart, rather than paraboloids.
5. Acknowledgements
We are thanks’ to Erasmus project 2020-1-PL01-KA203-082197 (iBIGworld – «Innovations for Big
Data in a Real World» for collaboration and information about source data for problem analysis. Work
is conducted under the Agreement on scientific cooperation between Taras Shevchenko National
University of Kyiv and the Faculty of Business and Management, Brno University of Technology.
6. References
[1] V.N. Kovalenko “Functional morphology of the ventricles of the heart as the basis for changing
the geometry of contraction” - Institute of Cardiology named after N. D. Strazhesko AMS of
Ukraine,2002 (in Russian)
[2] Kukanova, B. and B. Mravec, Complex intracardiac nervous system. Bratisl Lek Listy, 2006.
107(3): p. 45-51.
[3] Fredrickson, B.L., et al., What good are positive emotions in crises? A prospective study of
resilience and emotions following the terrorist attacks on the United States on September 11th, 2001.
Journal of Personality and Social Psychology, 2003. 84(2): p. 365-376.
[4] Beckham, A.J., T.B. Greene, and S. Meltzer-Brody, A pilot study of heart rate variability
biofeedback therapy in the treatment of perinatal depression on a specialized perinatal psychiatry
inpatient unit. Arch Womens Ment Health, 2013. 16(1): p. 59-65.
279
[5] Teixeira, F.S., Neufeld, E., Kuster, N. et al. Modeling intracranial aneurysm stability and
growth: an integrative mechanobiological framework for clinical cases. Biomech Model Mechanobiol
19, 2413–2431 (2020).
[6] Guccione, J., Kassab, G. S., eds. (2020). Mathematical Modeling of Cardiovascular Systems:
From Physiology to the Clinic. Lausanne: Frontiers Media SA. doi: 10.3389/978-2-88963-323-4
[7] Pauli, J. Nam, M. Pasquali and M. Behr, Transient Stress- and Strain-Based Hemolysis
Estimationin a Simplified Blood Pump, accepted to International Journal for Numerical Methods in
Bio-medical Engineering (2013). DOI: 10.1002/cnm.2576
[8] Quarteroni, A., Dede', L., Manzoni, A., & Vergara, C. (2019). Modelling the heart. In
Mathematical Modelling of the Human Cardiovascular System: Data, Numerical Approximation,
Clinical Applications (Cambridge Monographs on Applied and Computational Mathematics, pp. 102-
152). Cambridge: Cambridge University Press. doi:10.1017/9781108616096.007
[9] K. C. Wong, L. Wang, H. Zhang, H. Liu, and P. Shi, “Integrating functional and structural
images for simultaneous cardiac segmentation and deformation recovery,” in Proceedings of the
Medical Image Computing and Computer-Assisted Intervention (MICCAI ’07), vol. 4791, pp. 270–
277, 2007.
[10] R. A. B. Burton, G. Plank, J. E. Schneider et al., “Three-dimensional models of individual
cardiac histoanatomy: tools and challenges,” Annals of the New York Academy of Sciences, vol. 1080,
pp. 301–319, 2006.
[11] R. Hren, J. Nenonen, and B. M. Hora´cek, “Simulated ep- ˇ icardial potential maps during paced
activation reflect myocardial fibrous structure,” Annals of Biomedical Engineering, vol. 26, no. 6, pp.
1022–1035, 1998.
[12] D. Rohmer, A. Sitek, and G. T. Gullberg, “Reconstruction and visualization of fiber and laminar
structure in the normal human heart from ex vivo diffusion tensor magnetic resonance imaging
(DTMRI) data,
[13] Q. Long, R. Merrifield, X. Y. Xu, P. Kilner, D. N. Firmin, and G. Z. Yang, “Subject-specific
computational simulation of left ventricular flow based on magnetic resonance imaging,” Proceedings
of the Institution of Mechanical Engineers H, vol. 222, no. 4, pp. 475–485, 2008.
[14] V. V. Fedorov, A. V. Glukhov, R. Chang et al., “Optical mapping of the isolated coronary-
perfused human sinus node,” Journal of the American College of Cardiology, vol. 56, no. 17, pp. 1386–
1394, 2010.
[15] R. Bordas, K. Gillow, Q. Lou et al., “Rabbit-specific ventricular model of cardiac
electrophysiological function including specialized conduction system,” June 2011 Progress in
Biophysics and Molecular Biology 107(1): pp. 90-100 DOI:10.1016/j.pbiomolbio.2011.05.002
[16] A. Atkinson, S. Inada, J. Li et al., “Anatomical and molecular mapping of the left and right
ventricular His-Purkinje conduction networks,” The Journal of Molecular and Cellular Cardiology, vol.
51, no. 5, pp. 689–701, 2011.
[17] F. Gaudenzi and A. P. Avolio, "Lumped parameter model of cardiovascular-respiratory
interaction," 2013 35th Annual International Conference of the IEEE Engineering in Medicine and
Biology Society (EMBC), 2013, pp. 473-476
[18] Daisuke Ogawa, Shinji Kobayashi, Kenji Yamazaki, Tadashi Motomura, Takashi Nishimura,
Junichi Shimamura, Tomonori Tsukiya, Toshihide Mizuno, Yoshiaki Takewa, Eisuke Tatsumi,
Tomohiro Nishinaka, Evaluation of cardiac beat synchronization control for a rotary blood pump on
valvular regurgitation with a mathematical model, Artificial Organs, 10.1111/aor.13795, 45, 2, (124-
134), (2020).
[19] Selim Bozkurt, Pressure, Flow Rate and Operating Speed Characteristics of a Continuous Flow
Left Ventricular Assist Device During Varying Speed Support, IEEE Access,
10.1109/ACCESS.2020.3009264, 8, (129830-129841), (2020).
[20] Roel Meiburg, Wouter Huberts, Marcel C. M. Rutten, Frans N. Vosse, Uncertainty in model‐
based treatment decision support: Applied to aortic valve stenosis, International Journal for Numerical
Methods in Biomedical Engineering, 10.1002/cnm.3388, 36, 10, (2020).
[21] Albanese A, Cheng L, Ursino M, Chbat NW. An integrated mathematical model of the human
cardiopulmonary system: model development. Am J Physiol Heart Circ Physiol. 2016;310: H899–921
280