На основе предложенной ранее модели левого желудочка сердца построен численный метод для задачи расчета электрофизиологической активности миокарда. Эта задача сводится к решению системы дифференциальных уравнений реакции-диффузии. Численный метод позволяет найти Лапласиан и учесть классические краевые условия задачи. Метод основан на записи уравнений в специальной системе координат, связанной с желудочком. Строится равномерная по специальным координатам сетка, которая подвергается разрежению в декартовых координатах. Краевое условие - отсутствие потока потенциала через границу желудочка - также записывается в терминах специальных координат и учитывается с помощью метода фиктивных узлов. Приведен пример расчета на конкретной реализации модели левого желудочка сердца человека.
2.1
Модель формы ЛЖ и специальная система координат Для задания формы ЛЖ была использована специальная система координат ( ; ; ); где переменная 2 [ 0; 1] [0; 1] соответствует положению точки в толще стенки ЛЖ, = 0 эпикард, = 1 эндокард; 2 [0; =2] аналог географической широты, = 0 самая верхняя плоская часть модельного ЛЖ (зона клапанов и фиброзного кольца), = =2 верхушка ЛЖ; 2 [0; 2 ) аналог географической долготы.
Переход от специальных координат к цилиндрическим ( ; ; z) осуществляется согласно формулам (рис. 1): ( ; ; ) = layer( ; r0( ) d0( ); re( )
de( ); l0( ) + (l1( ) l0( )); p( )); z( ; ) = Z (Z layer( ; r0; re; l; p) = <8 re (re r0) 1
h l pi : re 1 =2 l ;
l h ) sin + z; p ; если иначе; < l; (1) где r0;e( ) радиус ЛЖ на эпикарде основания, экватора; d0;e( ) толщина стенки ЛЖ на эпикарде основания, экватора; l0;1( ) ¾широта¿ экватора на эпикарде, эндокарде; Z высота ЛЖ; h толщина стенки ЛЖ на верхушке; p > 1.
По аналогии с симметричной моделью [2], для подгонки углов наклона волокон с помощью величин 0; 1 мы пересчитываем входные параметры следующим образом: параметры, зависящие от координаты : радиусы и толщины стенки ЛЖ: широты экватора: параметры, не зависящие от меридиана: толщина стенки ЛЖ на верхушке: смещение ЛЖ по вертикали: высота ЛЖ: l0 = ( 1 1)l1i + (1 0)l0i
; r0 = r0i +
d0 = re = rei + de = Индекс i обозначает входное (input), или эффективное, значение параметра.
Рис. 1: меридиональное сечение модели (рисунок из [1]). 2.2
Приближение формы желудочка в сечениях Модель ЛЖ строится на основе данных медицинского обследования сердца: ЭхоКГ (УЗИ), КТ, МРТ или ДТ-МРТ. Вначале выделяется вертикальная ось ЛЖ, ось Oz. Она должна проходить через верхушку от эпикарда до эндокарда и не должна проходить через миокард вне верхушки. Полуплоскости, проходящие через ось, называются меридиональными. Данные о форме стенки ЛЖ задаются в нескольких меридианах. При обработке снимков УЗИ обычно рассматривают 4 меридиана; данные томографии позволяют рассмотреть большее количество меридианов (порядка 20-30). Подгонка формы модели к экспериментальным данным проводится с помощью подбора параметров модели r0;e, d0;e, l0;1, p в каждом меридиане по отдельности. После этого модель, построенную в меридианах, нужно достроить между ними. Достройка производится интерполяцией или аппроксимацией подобранных параметров по периодической переменной 2 [0; 2 ) – долготе. 2.3
От меридиональных сечений к цельной трёхмерной модели Пусть заданы значения долготы i 2 [0; 2 ), i = 0; 1; : : : ; n 1, сечений, в которых были проведены измерения. В результате измерений были получены значения yi некоторого параметра модели. Нужно построить приближающую функцию y( ); 2 [0; 2 ). Желательно, чтобы аппарат приближения в такой задаче обладал следующими свойствами.
1. Интерполяционность, если мы предполагаем, что обработка данных в меридианах осуществлена очень точно. Если такого предположения нет и известна погрешность измерения каждого отдельного параметра, можно ограничиться аппроксимацией с гарантией непревышения указанной погрешности. 2. Периодичность, так как независимая переменная периодическая. 3. 2-гладкость. 4. Неотрицательность, или равносильное ей требование: если оказалось, что зависимая переменная yi лежит на некотором отрезке [a; b], то и приближающая функция y( ) должна лежать на [a; b]. Это требование вызвано тем, что некоторые параметры по смыслу задачи неотрицательны или лежат на определенном отрезке. В ходе работы с моделью мы использовали три способа приближения.
1. Кусочно-линейные интерполяционные сплайны. Они обладают свойствами 1, 2 и 4, но не являются гладкими. В некоторых случаях трёхмерная модель ЛЖ, построенная по данным УЗИ, имела острые края на эндокарде, что не соответствует экспериментальным данным и плохо отражается на устойчивости счёта.
2. Периодические кубические интерполяционные сплайны. Эти сплайны обладают свойствами 1, 2 и 3, но могут выходить за отрезок, на котором лежат исходные данные.
3. Периодические тригонометрические сплайны из [3]. В этом случае сплайн удовлетворяет свойствам 2 и 4. Достаточно ли точно он приближает исходные данные, приходится проверять визуально или автоматически.
Результаты приближения для одной конкретной модели ЛЖ представлены на рис. 2. Приближаемая величина на графике А по своему смыслу неотрицательна, и мы видим, что кубический сплайн нарушает это требование. В то же время негладкий кусочно-линейный и гладкий тригонометрический сплайн остаются неотрицательными. Тригонометрический сплайн, будучи неинтерполяционным, существенно отличается от исходных данных в некоторых точках. Рис. 2: сравнение трёх способов приближения одного из параметров модели левого желудочка: кусочнолинейный сплайн (синяя сплошная линия), кубический сплайн (коричневая пунктирная линия), тригонометрический сплайн (красная штрихпунктирная линия). Квадраты – исходные данные. А – модель строится по 4-м меридианам, Б – по 24-м.
Пример модели ЛЖ – на рис. 3. 3 3.1 Моделирование электрофизиологической активности миокарда желудочка Система реакции-диффузии Электрофизиологические процессы в миокарде обычно описывают с помощью системы дифференциальных уравнений реакции-диффузии, имеющей следующий общий вид: = div(D grad u) + f (u; w); Рис. 3: модель левого желудочка, построенная по 4-м меридиональным сечениям. Сетка треугольников построена с помощью пакета ANI-3D [4], визуализация в Tetview [5]. где u = u(r; t) – трансмембранный потенциал клетки в точке r в момент времени t, D – матрица диффузии (3 3), w – вектор остальных фазовых переменных модели, f (u; w); g(u; w) – функции, зависящие от модели кардиомиоцита.
Задача ставится на некотором теле , имеющем границу .
Краевое условие – электрическая изолированность миокарда, то есть отсутствие потока потенциала u через границу : nD grad u = 0; (2) где n – вектор нормали к поверхности, D – матрица диффузии. Также заданы начальные условия – область начальной стимуляции.
Для решения начально-краевой задачи в случае трехмерного тела, имеющего простую форму (прямоугольного параллелепипеда, например), можно воспользоваться равномерной в Декартовых координатах расчетной сеткой. Если же тело имеет сложную форму, то обычно следуют одним из трёх путей. 1. Проводят разбиение тела на тетраэдры и решают задачу методом конечных элементов. 2. Строят равномерную в Декартовых координатах сетку в прямоугольном параллелепипеде (или кубе), целиком содержащем тело. В дальнейшем применяется метод конечных разностей. Достоинством этого пути является простота сетки, а недостатком – сложность учёта краевых условий, так как граница области обычно не проходит через узлы сетки и не является координатной поверхностью.
3. Связывают с телом специальную систему координат, в которой его граница является набором координатных поверхностей. В этом случае относительно легко учесть краевое условие, но требуется обратить особое внимание на качество сетки, потому что при применении явных конечно-разностных методов для довольно жёстких задач электрофизиологии миокарда расстояния между соседними узлами сетки должны находиться в достаточно узких пределах. Известно, что если узлы будут слишком близки, то для достижения устойчивого счёта потребуется уменьшать шаг по времени. Если слишком далеки, то волна возбуждения может остановиться. Если расстояния будут меняться в слишком широких пределах, то скорость волны будет слишком неравномерной, в то время как она должна быть постоянной. 3.2
Построение расчётной сетки Пусть равномерная по координатам ( ; ; ) сетка имеет N + 1 узел (с индексами i от 0 до N ) по , N + 1 узел (с индексами j от 0 до N ) по , N + 1 узел (с индексами k от 0 до N ) по .
Сетка, равномерная в специальных координатах, является чрезвычайно неравномерной в Декартовых координатах. Более того, все узлы с i = 0, j = N , k = 0; 1; : : : ; N совпадают; узлы с i = 1, j = N , k = 0; 1; : : : ; N совпадают и так далее. Эта особенность сетки вынуждает нас искать соседей узлов особым образом: если ближайший сосед в каком-то направлении расположен слишком близко, нужно рассмотреть более далёких соседей. Пусть заданы ограничения [dmin; dmax] на расстояния между соседними узлами сетки. Для каждого узла нужно найти 6 его соседей: 2 по , 2 по и 2 по – в положительном и отрицательном направлениях.
Если среди всех узлов на соответствующем координатном луче соседа на подходящем расстоянии не нашлось, то такой узел будем называть интерполируемым и пометим специальной меткой ¾Out¿. Она означает, что значения фазовых переменных в этом узле мы будем искать не методом Эйлера, а с помощью интерполяции на основе значений в узлах без этой метки. Мы воспользовались линейной интерполяцией по : u[i][j][k] = u[i][j][k0] + (u[i][j][k1] u[i][j][k0]) k k0 k1 k0 ; где k0 – индекс ближайшего неинтерполируемого соседа узла (i; j; k) в отрицательном направлении по , k1 – то же в положительном направлении. Очевидным образом вводится необходимая поправка на периодичность.
Если ограничениям удовлетворяет несколько узлов, то нужно выбрать тот, расстояние до которого ближе всего к заданной величине dopt. 3.2.1
Алгоритм построения сетки Следующий алгоритм для каждого узла пытается определить его соседей по каждой из трёх координат ( , , ) по двум направлениям (увеличение и уменьшение индекса). Кроме того, алгоритм выделяет среди узлов интерполируемые.
Исходные параметры: N , N , N – максимальные индексы узлов сетки по , , (нумерация с 0); dmin; dopt; dmax – минимальное, оптимальное, максимальное расстояние между соседями. Результат: массивы размера [N + 1][N + 1][N + 1]:
Out (логический) – является ли узел интерполируемым, Q , Q , Q и Q+, Q+, Q+ (целые неотрицательные) – разности индексов по координатам , , . Инициализация: все элементы массива Out Ложь; цикл для узлов сетки i = 0; 1; : : : ; N , j = 0; 1; : : : ; N 1, k = 0; 1; : : : ; N выполнять цикл для обоих направлений dir 2 f+; g по всем координатам var 2 f ; ; g выполнять Пусть рассматривается массив Qvdiarr: Qvdiarr[i][j][k] F ind(i; j; k; dir; var) если Qvdiarr[i][j][k] = 1, то
Out[i][j][k] Истина; выйти из циклов по dir, var конец конец конец цикл для узлов сетки i; j; k с меткой Out[i][j][k] =Истина выполнять Найти узлы с Out =Ложь, ближайшие по в обоих направлениях, и записать расстояния до них в Q+; [i][j][k] конец // Все узлы верхушки – интерполируемые. j N цикл для узлов сетки i = 0; 1; : : : ; N , k = 0; 1; : : : ; N выполнять
Out[i][j][k] Истина конец 3.2.2
Алгоритм F ind(i; j; k; dir; var) 3.3
Вычисление Лапласиана Лапласиан переменной u в анизотропной среде с матрицей диффузии D, как было сказано выше, вычисляется по формуле Здесь элементы матрицы D равны
div(D grad u): Dij = D2 ij + (D1
D2)vivj; где i; j = 0; 1; 2; D1 и D2 – коэффициенты диффузии соответственно вдоль и поперёк волокон; – символ Кронекера; v = (v0; v1; v2) – единичный вектор направления волокна.
Вслед за работой [6], для нахождения Лапласиана в криволинейной системе координат ( 0; 1; 2) рассмотрим матрицы
J = (Jij) = W = (Wij) = где v = (v0; v1; v2) – единичный вектор направления волокна. Эти матрицы связаны между собой следующими соотношениями:
T mkp =
W = SJ; X Jkl(JT HlJ)mp:
l Тогда Лапласиан выражается следующей формулой: где qkl – элементы матрицы Q: pk = D2 tr Tk + (D1
D2) (Jv)k tr(SJ) + (JSJv)k + vTTkv ;
Q = JDJT: (3) 3.3.1
Вычисление Лапласиана на сетке Применение вышеуказанных формул и построенной сетки приводит к следующему выражению для первой производной потенциала u по специальной координате (для , – аналогично): [i][j][k] = u[i + Q+][j][k] u[i Q ][j][k]
Q+ + Q 1 N 0 : Выражения для вторых производных: [i][j][k] = u[i Q ][j + Q+][k] + u[i + Q+][j + Q+][k] + u[i Q ][j Q ][k] u[i + Q+][j
Q+ u[i + Q+][j][k] u[i][j][k] u[i Q ][j][k] u[i][j][k]
1 Q+ + Q 1 N 0 2 ; другие вторые производные вычисляются аналогично.
(Q+ + Q ) (Q+ + Q ) N 1 +
N 3.4
Учёт краевых условий Напомним (см. (2)), что краевое условие записывается как где n – вектор нормали к поверхности Это уравнение преобразуется к виду
желудочка. где J ; ;' – вектор-столбцы производных специальных координат по Декартовым (см. (3)): nTDJ + nTDJ nD grad u = 0;
0 x1 (4) (5) и так далее. Краевые условия мы учитываем с помощью метода фиктивных узлов: вне модели ЛЖ строится дополнительный слой фиктивных узлов, и значения потенциала в них мы находим, исходя из производных первого порядка по соответствующим специальным координатам. Этот метод позволяет вычислять Лапласиан во всех нефиктивных узлах единообразно, будь то узлы на границе ЛЖ или в его толще.
Для эндокарда и эпикарда уравнение (5) нужно решить относительно @@u и найти значение потенциала в фиктивном узле за эндокардом (в полости ЛЖ) или за эпикардом соответственно.
Для базы ЛЖ уравнение (4) мы приводим к виду
X k 0
D20 x + D21 y + D22 z Требуемые производные специальных координат по декартовым можно найти численно. Фиктивный узел в этом случае будет расположен над базой ЛЖ. 4
Программная реализация
Рассматриваемый в настоящей статье численный метод реализован как часть программы для расчёта
электромеханической активности на трёхмерной модели ЛЖ сердца. Программа написана на языке Си.
Пользователь может выбрать требуемые модели клетки (Алиева–Панфилова (упрощённая и полная [7]),
TNNP [
Наиболее вычислительно затратные циклы распараллелены с помощью технологии OpenMP.
Исходные данные задаются в файле ini (работа с ним ведётся с помощью сторонней библиотеки iniparser
[
Пример расчета По данным ЭхоКГ сердца человека была построена модель ЛЖ в 4 меридианах, а затем она была достроена до полной с помощью тригонометрических сплайнов. Использована упрощенная модель Алиева– Панфилова = div(D grad u) ku(u a)(u 1) uv; где u = u(~r; t) – безразмерный трансмембранный потенциал клетки в точке ~r в момент времени t, v – безразмерная проводимость клетки для калиевого тока, k = 8, a = 0:03. Функция (u) = 1, если u < a, иначе (u) = 0:1.
Коэффициент диффузии вдоль волокон был равен D1 = 12, поперек – в 9 раз меньше, что соответствует физиологической анизотропии миокарда, в котором скорость распространения потенциала действия вдоль волокон в 3 раза выше, чем поперек. Параметры 0 = 0:1, 0 = 0:85 обеспечили угол вращения волокон в стенке ЛЖ, близкий к физиологическому, 130 .
Шаг по времени составлял dt = 0:0016666 ед. модельного времени, по пространству – 0.5 ед. модельного
расстояния. С помощью методики, описанной в [
На рис. 4 показан вид ЛЖ с раскраской по потенциалу. Мы видим, что спиральная волна начинается в области верхушки. Её передний фронт – верхняя зелёная кривая, а задний – нижняя зелёная кривая. Волна движется вверх и влево, то есть её направление вращения – против часовой стрелки. На рис. 4, справа, показаны две спиральные волны – на эндо- и эпикарде. Видно, что возбуждённая область связна и тянется от начала спирали, огибая желудочек, до его основания. Рис. 4: Спиральная волна на несимметричной модели ЛЖ сердца человека, построенной по данным ЭхоКГ. Слева – виден только эпикард, справа – полупрозрачный желудочек, виден эпикард и эндокард. Электрически возбужденные области получили цвета от красного до зеленого, состояние покоя показано синим цветом. 6 6.1 Обсуждение
Недостатки модели ЛЖ 1. Результаты кардиологического обследования пациентов и лабораторные эксперименты показывают, что стенка ЛЖ может иметь форму более сложную, чем предусмотрено моделью. В некоторых случаях подобрать параметры модели формы не удаётся и приходится сглаживать, возможно, существенные особенности стенки желудочка.
2. Базальная часть ЛЖ действительно может быть отсечена от предсердий по плоскому, как в модели, фиброзному кольцу. Но модель призвана приближать форму только миокарда без соединительной ткани, а предположение о плоском основании миокарда ЛЖ не вполне верно. Работы анатомов (например, классическая книга Ф. Торрента Гуаспа [12]) показывают, что миокардиальное основание желудочков сердца закруглено.
3. В модели волокна не касаются основания, но, по данным анатомов, около половины волокон проходит основание по касательной, не врастая в фиброзное клапанное кольцо. (Этот недостаток связан с предыдущим.)
4. Реальная архитектоника ЛЖ на верхушке, особенно геометрия слоёв миокарда, качественно отличается от модельной. 5. При малом числе меридианов, в которых заданы параметры сечений модели, между меридианами приближение формы ЛЖ оказывается довольно далеким от реального.
6. Вторые производные Декартовых координат x; y по специальной координате (элементы матриц H01;1 и H11;1) могут быть неограниченными в окрестности точки = l при p 2 (1; 2) в силу используемой формулы (1). Этот недостаток требует добавления в программу проверки при численном нахождении этих производных по трём узлам, не находятся ли узлы по по разные стороны от точки = l.
Тем не менее, пока сложно оценить, насколько все эти упрощения и недостатки модели влияют на результаты электрофизиологических расчётов. 6.2
Недостатки численного метода Недостатком изложенного метода является непостоянство коэффициента [1] S. Pravdin. Non-axisymmetric mathematical model of the cardiac left ventricle anatomy. Russian
Journal of Biomechanics, 17, 4(62): 75-94, 2013. [2] 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, 12(54), 2013. [3] K.V. Kostousov, V.T. Shevaldin. Approximation by local trigonometric splines.
77(3):354-363, 2005.
Mat. Zametki, [4] Advanced Numerical Instruments 3D – ANI-3D. https://sourceforge.net/projects/ani3d/ [5] TetView – A Tetrahedral Mesh and Piecewise Linear Complex Viewer. wias-berlin.de/software/ tetgen/tetview.html [6] 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):e93617, 2014. [7] R.R. Aliev, A.V. Panfilov. A simple two-variable model of cardiac excitation. Chaos, Solitons and
Fractals, 7(3):293-301, 1996. где ku – конечная разность u в формуле производной первого порядка по переменной с индексом k, k;lu – конечная разность u в формуле производной второго порядка по переменным с индексами k; l, k – приращение по специальной координате k. Представляется, что от этого недостатка можно избавиться путём использования корректирующих множителей, приводящих задачу с коэффициентом диффузии D и шагом по пространству dopt к задаче с изменённым коэффициентом диффузии D0 при фактическом расстоянии мужду узлами dfact 6= dopt. Благодарности
Автор выражает благодарность Т.В. Чумарной (ИИФ УрО РАН, г. Екатеринбург) за предоставленные данные ЭхоКГ; А.А. Данилову (ИВМ РАН, г. Москва) за консультации по использованию пакета ANI-3D; К.С. Ушенину (УрФУ) за предоставленную иллюстрацию (рис. 4).
Исследование выполнено за счет гранта Российского научного фонда (проект УрФУ 14-35-00005). При проведении работ были использованы суперкомпьютеры ИММ УрО РАН (¾Уран¿) и УрФУ. Список литературы при конечной разности переменной u в Лапласиане (dr – шаг по пространству). В трёхмерном случае Лапласиан имеет вид (см. (3)) div(D grad u)
X pk k ku + X qkl k k;l k;lu ; k
l D dr2
A method of solving reaction-diffusion problem on a non-symmetrical model of the cardiac left ventricle
Sergei F. Pravdin Krasovskii Institute of Mathematics and Mechanics (Yekaterinburg, Russia) Ural Federal University (Yekaterinburg, Russia)
Keywords: electrophysiology, computational cardiology, anisotropy, myocardium, left ventricle, mathematical biology, biophysics.