Программный комплекс для сравнения электрофизиологических свойств волокон миокарда на одномерных и трехмерных моделях Незлобинский Т.В.2 Правдин С.Ф.1,2 nezlobinsky@yandex.ru 1 – ИММ УрО РАН (Екатеринбург) 2 – УрФУ (Екатеринбург) Аннотация Актуальной задачей на сегодняшний день является изучение осо- бенностей проведения электрического возбуждения по сердечной мышце. Нами было проведено сравнение скорости распростране- ния волны возбуждения по трехмерному волокну модели левого желудочка сердца и по одномерному волокну такой же длины. Раз- работано программное обеспечение, позволяющее проводить расче- ты для волокон различной длины и различного пространственного расположения в толще стенки модели левого желудочка. Получен- ные результаты могут быть полезны в исследованиях физиологии сердечной мышцы. 1 Введение В данной статье описано разработанное программное обеспечение для постановки вычислительного экс- перимента по сравнению скорости прохождения волны возбуждения по трехмерному волокну, являюще- муся частью миокарда, и одномерному (изолированному) волокну. Подобная задача является достаточно сложной для реализации в натурном эксперименте, поэтому возникает необходимость в использовании математических моделей, позволяющих ставить вычислительный эксперимент без затраты значительных ресурсов. Существует ряд моделей, описывающих строение и работу сердца (их обзор см. в [1]). Среди них можно выделить класс анатомических моделей, описывающих геометрию сердца или отдельных его камер, а так- же систему мышечных волокон, и класс электрофизиологических моделей, описывающих электрические явления, происходящие в миокарде. Совместное использование анатомических и электрофизиологических моделей в рамках одного вычислительного эксперимента дает возможность получить результаты, близкие к физиологическим. 2 Анатомическая модель левого желудочка сердца Для моделирования распространения волны возбуждения по трехмерного волокну, мы воспользовались анатомической моделью левого желудочка сердца из работы [2]. Модель основана на идее представления левого желудочка в виде семейства вложенных спиральных поверхностей, заполненных кривыми – во- локнами. Ранее, для данной модели левого желудочка был разработан программный комплекс, который Copyright c by the paper’s authors. Copying permitted for private and academic purposes. In: A.A. Makhnev, S.F. Pravdin (eds.): Proceedings of the 47th International Youth School-conference “Modern Problems in Mathematics and its Applications”, Yekaterinburg, Russia, 02-Feb-2016, published at http://ceur-ws.org 268 позволяет использовать различные электрофизиологические модели кардиомиоцитов для расчета охвата левого желудочка волной возбуждения. Для моделирования электрических явлений, происходящих в серд- це, мы использовали электрофизиологическую модель кардиомиоцита Алиева – Панфилова, качественно описывающую потенциал действия кардиомиоцита [3]. Миокард как проводящая среда имеет анизотропные свойства, которые заключаются в разнице в скоро- стях распространения волны возбуждения вдоль и поперек волокон. Кроме того, имеет место вращательная анизотропия миокарда, вызванная трансмуральным вращением волокон в стенке левого желудочка, что также оказывает влияние на скорость проведения электрического возбуждения [4] (рис. 1). Рис. 1: Осесимметричная модель левого желудочка (слева) и трансмуральное вращение волокон в толще его стенки (справа) [5] 3 Описание программного комплекса Вычислительный эксперимент включал следующие этапы: 1. Выбор зоны начальной стимуляции (область миокарда, откуда происходит запуск волны возбужде- ния). 2. Проведение расчета на модели левого желудочка сердца со встроенной электрофизиологической мо- делью Алиева – Панфилова для выбранной зоны стимуляции. 3. Выбор волокна на спиральной поверхности модели левого желудочка и нахождение времени прихода волны возбуждения в каждый из узлов волокна. 4. Решение одномерной задачи для волокна той же длины, что и выбранное волокно на спиральной поверхности, с аналогичными параметрами модели Алиева – Панфилова. Результаты проведенного вычислительного эксперимента представлялись в виде графика скоростей про- хождения волны возбуждения по одномерному и трехмерному волокнам. Кроме того, для каждого экспе- римента строилась спиральная поверхность с наложенным на нее волокном, для оценки его расположения в толще стенки левого желудочка. Программное обеспечение (табл. 1) разработано на языке python 2.7, с использованием библиотек matplotlib (построение двумерных графиков) [6], mayavi (построение и работа с трехмерной графикой) [7] и pyside (создание графического интерфейса) [8]. Для проведения одного вычислительного эксперимента требуется в среднем 1 – 3 минуты в зависимости от производительности процессора и установленного мо- дельного времени эксперимента (при фиксированном числе узлов разбиения волокна и шаге по времени). 269 Таблица 1: Структура программного комплекса Модуль Назначение Решатель системы дифференциальных уравнений one_dim_solver.py в частных производных модели Алиева–Панфилова для случая одномерных волокон Обработка результатов расчета на модели левого же- лудочка. Построение волокна на спиральной поверх- three_dim_processing.py ности и подготовка данных для трехмерной визуали- зации interpolation.py Реализация метода трилинейной интерполяции gui.py Построение графического интерфейса пользователя Модуль запуска программы, осуществляющий вызов main.py остальных модулей и обеспечивающий обмен данны- ми между ними 3.1 Решение задачи распространения волны возбуждения по одномерному волокну Модель Алиева – Панфилова представляет собой систему из двух дифференциальных уравнений в частных производных: ∂u ∂2u = D1 2 − ku(u − a)(u − 1) − uv, ∂t ∂x ∂v = (u, v)(−v − ku(u − a − 1)), (1) ∂t µ1 v (u, v) = 0 + , u + µ2 где фазовая переменная u описывает трансмембранный потенциал (разность потенциалов между внешней и внутренней сторонами клеточной мембраны), фазовая переменная v – мембранный ток восстановления; D1 – коэффициент диффузии вдоль направления волокна; k, a, 0 , µ1 , µ2 – параметры модели клетки. Мышечное волокно миокарда может быть представлено в виде последовательности узлов x0 , x1 , . . . , xn . Будем считать, что в начальный момент времени стимулируется только первый из n узлов волокна: u(0, 0) = 1, u(xi , 0) = 0, i = 1, 2, 3, . . . , n. (2) xi = i∆x, i = 0, 1, 2, . . . , n. Граничные условия второго рода моделируют электрическую изолированность концов волокна: ∂u = 0, x = x0 , x = xn . (3) ∂x Решая систему (1) – (3), мы получим картину распространения волны возбуждения по одномерному волокну с постоянной скоростью (из-за постоянства коэффициента диффузии). Первое уравнение пред- ставляет собой уравнение теплопроводности с функцией внутренних источников. Для численного решения уравнения теплопроводности была выбрана безусловно устойчивая неявная разностная схема. Для числен- ного решения второго уравнения системы (ОДУ) мы воспользовались явным методом Эйлера. По ходу решения мы также фиксировали момент времени, когда потенциал превышал заранее установ- ленный порог (находили время прихода волны возбуждения). В поставленных вычислительных экспери- ментах этот порог был равен 0.5 (в модели Алиева – Панфилова безразмерный потенциал u ∈ [0, 1]). 3.2 Обработка результатов расчета на трехмерной модели левого желудочка Построение модели левого желудочка осуществляется в специальной системе координат (γ, ψ, φ), где γ – положение точки в толще стенки левого желудочка (γ = γ0 – эндокард, γ = γ1 – эпикард), ψ – аналог 270 географической широты (ψ = 0 – самая верхняя плоская часть модели левого желудочка, π/2 – верхушка, самая узкая часть модели левого желудочка), φ – аналог географической долготы (от 0 до 2π вокруг оси вращения левого желудочка). Точки в специальной системе координат достаточно просто переводятся в цилиндрическую и декартову системы координат [2]. Модель левого желудочка имеет следующие геомет- рические параметры (на рис. 2 представлено меридиональное сечение модели): высота левого желудочка Z, толщина стенки на основании d0 , толщина стенки на экваторе de , «широта» ψ экватора на эпикарде l0 , «широта» ψ экватора на эндокарде l1 , толщина стенки на верхушке h, радиус желудочка на эпикарде основания r0 , радиус желудочка на эпикарде экватора re . Рис. 2: Меридиональное сечение модели левого желудочка [2] Система дифференциальных уравнений в частных производных модели Алиева–Панфилова для трех- мерного случая имеет вид: ∂u = div(D grad u) − ku(u − a)(u − 1) − uv, ∂t ∂v = (u, v)(−v − ku(u − a − 1)), (4) ∂t µ1 v (u, v) = 0 + , u + µ2 где D – матрица диффузии (с элементами Di,j = D2 δi,j + (D1 − D2 )wi wj , где D1 – коэффициент диффузии вдоль волокон, D2 – коэффициент диффузии поперек волокон, δi,j – символ Кронекера, w ~ = (w0 , w1 , w2 ) – единичный вектор направления волокна). Для расчета на трехмерной модели левого желудочка используется ранее разработанный программный комплекс [9], обеспечивающий построение геометрической модели и расчет волны возбуждения на трех- мерной сетке. В результате данного расчета мы получаем данные о времени прихода волны возбуждения в каждый из узлов сетки левого желудочка в форме: координаты узла (γi , ψj , φk ) и время прихода волны в узел Tijk . Задача заключается в построении волокна на спиральной поверхности и нахождении времени прихода волны в точки волокна, исходя из полученных данных. Согласно используемой модели левого же- лудочка [2], прообразом волокна на спиральной поверхности является хорда Y = const полукруга радиуса P = 1 (рис. 3). Двигаясь равномерно по полярному углу от одного из концов хорды к другому, можно разбить хорду на конечное число узлов (отрезки разбиения равны не будут): Φ1 − Φ0 ∆Φ = , (5) n 271 Рис. 3: Горизонтальные хорды на полукруге. Φ0 и Φ1 – полярные углы соответственно правого и левого концов хорды [2] Φi = Φ0 + i∆Φ, i = 0, 1, 2, . . . , n, (6) Y Pi = , i = 0, 1, 2, . . . , n, (7) sin(Φi ) где ∆Φ – шаг разбиения по полярному углу, n – число узлов разбиения. Узлы хорды однозначно отобра- жаются в узлы волокна на спиральной поверхности: Φ γ(Φ) = , (8) π π ψ(P ) = (1 − P ) . (9) 2 По массиву координат узлов волокна мы вычисляем длину волокна, а также время прихода волны возбуждения в эти узлы, воспользовавшись известными значениями времени прихода волны Tijk в узлы сетки модели левого желудочка (методом трилинейной интерполяции). 3.3 Графический интерфейс пользователя На рис. 4 представлен графический интерфейс программного комплекса. В левой половине – панель управления, состоящая из блоков: «Fiber parameters» – задание числа узлов разбиения волокон и выбор волокна на спиральной поверхности; «Conditions» – загрузка файла, хранящего значения времен прихо- да волны возбуждения в узлы модели левого желудочка (Time of arrival), и конфигурационного файла, хранящего геометрические параметры модели левого желудочка, параметры модели Алиева – Панфилова и другие параметры вычислительного эксперимента (Initial parameters); «Control bar» – запуск или оста- новка вычислительного эксперимента; «Status screen» – окно, отображающее ход процесса вычисления. В правой половине – окно с графиком сравнения (для вкладки «Comparison graph»). Вкладка «Spiral surface» переключает на трехмерное изображение спиральной поверхности с волокном (рис. 5). 4 Пример расчета Для вычислительного эксперимента были выбраны следующие параметры геометрической модели ле- вого желудочка: высота левого желудочка Z = 100 мм, толщина стенки на основании d0 = 10 мм, толщина стенки на экваторе de = 10 мм, «широта» ψ экватора на эпикарде l0 = 0.2, «широта» ψ экватора на эн- докарде l1 = 0.2, толщина стенки на верхушке h = 8 мм, радиус желудочка на эпикарде основания r0 = 30 мм, радиус желудочка на эпикарде экватора re = 33 мм. Границы расчетной области: γ ∈ [0.1, 0.9]; ψ ∈ [0, π/2]; φ ∈ [0, 2π]. Параметры электрофизиологической модели Алиева – Панфилова были равны: k = 8, a = 0.03, µ1 = 0.12, µ2 = 0.3,  = 0.01. Коэффициент диффузии вдоль волокон D1 = 0.05, а поперек волокон – в девять раз меньше: D2 = D1 /9. Такое соотношение связано с различием скорости распространения вол- ны возбуждения вдоль и поперек волокон в миокарде по экспериментальным данным (поперек волокон в три раза меньше, чем вдоль волокон) [10]. 272 Рис. 4: Графический интерфейс программы с графиком времени прихода волны как функции расстояния по волокну от первого стимулированного узла Рис. 5: Графический интерфейс программы со спиральной поверхностью и волокном на ней. Цвет поверх- ности соответствует времени прихода волны возбуждения Рис. 6: Левый желудочек в разрезе. Красным выделена область стимуляции – субэндокард 273 Чтобы продемонстрировать основные наблюдения, выберем в качестве области начальной стимуляции зону субэндокарда – часть миокарда, примыкающую к полости желудочка (рис. 6). На рис. 4 представлен график времен прихода волны возбуждения в узлы волокна на одномерных и трехмерных моделях. Из графика следует, что волна возбуждения в трехмерном случае распространяется значительно быстрее, чем в одномерном, несмотря на равенство коэффициентов диффузии D1 вдоль во- локон. Данный эффект вызван тем, что фронт волны возбуждения способен достигать узлов трехмерного волокна через соседние узлы, то есть через тканевое окружение волокна. Для выбранной зоны стимуляции эффект отчетливо наблюдается из-за превращения анизотропной среды в псевдоизотропную [4]. 5 Заключение Разработанная методика проведения вычислительного эксперимента и соответствующее программное обеспечение позволяют отслеживать скорость распространения волны возбуждения по выбранному во- локну модели левого желудочка сердца, а также сравнивать полученные результаты с результатами на изолированном волокне такой же длины, что довольно сложно осуществить в условиях реального экспе- римента. В частности, нами было показано наличие существенного опережения скорости распространения волны возбуждения по трехмерному волокну по сравнению с одномерным, что можно объяснить наличием окружающей трехмерное волокно тканью миокарда, оказывающей значительное влияние на скорость рас- пространения возбуждения по волокну. Полученные результаты могут представлять интерес в вопросах, связанных с изучением электрических явлений, происходящих в миокарде во время сердечных сокращений. Благодарности Авторы выражают благодарность своим коллегам из Института иммунологии и физиологии УрО РАН: Л.Б. Кацнельсону, Н.А. Викуловой, А.Д. Хохловой, О.Э. Соловьевой. Работа поддержана Российским на- учным фондом №14-35-00005. Список литературы [1] N.A. Trayanova. Whole Heart Modeling: Applications to Cardiac Electrophysiology and Electromechanics. Circulation research, 108(1):113-128, 2011. [2] S.F. Pravdin. Nonaxisymmetric mathematical model of the cardiac left ventricle anatomy. Russian journal of biomechanics, 17, (62):84–105, 2013. [3] R.R. Aliev and A.V. Panfilov. A simple two-variable model of cardiac excitation. Chaos, Solitons and Fractals, 7(3):293–301, 1996. [4] S.F. Pravdin, H. Dierckx, L.B. Katsnelson, O. Solovyova, V.S. Markhasin, A.V. Panfilov. Electrical wave propagation in an anisotropic model of the left ventricle based on analytical description of cardiac architecture. PLOS ONE, 9(5), 2014. [5] S.F. Pravdin, V.I. Berdyshev, A.V. Panfilov, L.B. Katsnelson, O. Solovyova, V.S. Markhasin. Mathematical model of the anatomy and fibre orientation field of the left ventricle of the heart. BioMedical Engineering OnLine, 54(12):1–21, 2013. [6] Matplotlib. http://matplotlib.org/ [7] The MayaVi Data Visualizer. http://mayavi.sourceforge.net/ [8] PySide. https://wiki.qt.io/PySide [9] S.F. Pravdin. A method of solving reaction-diffusion problem on a non-symmetrical model of the cardiac left ventricle, Proceedings of the 47th International Youth School-conference «Modern Problems in Mathematics and its Applications»: 284–296, Yekaterinburg, Russia, 2016. [10] Kanai A., Salama G. Optical mapping reveals that repolarization spreads anisotropically and is guided by fiber orientation in guinea pig hearts. Circ Res, 77(4):784-802, 1995. 274 Software complex for the comparison of electrophysiological properties of myocardial fibers in 1D and 3D models Timur V. Nezlobinsky2 , Sergei F. Pravdin1,2 1 – Krasovskii Institute of Mathematics and Mechanics (Yekaterinburg, Russia) 2 – Ural Federal University (Yekaterinburg, Russia) Keywords: left ventricle, cardiac model, myocardial fibers, rotational anisotropy. Nowadays, investigations of the heart fibers and their role in myocardium conduction in relation to the electrical stimulation are actual in physiology. We compared the propagation speed of the excitation wave in three- dimensional fibers in mathematical model of the left ventricle with the propagation speed in one-dimensional fiber. We developed software that allows us to conduct experiments for fibers with different lengths and spatial locations. Results can be useful in investigations of the cardiac physiology. 275