В данном исследовании при помощи численных экспериментов мы сравниваем процессы активации, деполяризации и реполяризации, происходящие в левом желудочке (ЛЖ) сердца и в тонком срезе ЛЖ. В компьютерных моделях геометрия ЛЖ и направление волокон являются реалистичными, персонифицированными по данным диффузионно-тензорной магнито-резонансной томографии. Было показано общее сходство электрофизиологических процессов в тонком срезе и в трехмерном теле. Наблюдаемые различия во времени активации было более явным, чем во времени реполяризации, причем данные различия были наиболее выражены в отдельных регионах. Значительное расхождение было получено в значениях потенциалов внеклеточного домена во время распространения волны возбуждения. Основываясь на полученных результатах, мы можем предположить, что различия в геометрии и ходе волокон влияют на процессы деполяризации и реполяризации в миокарде.
Тонкие срезы миокардиальной ткани (100-200 мкм) в экспериментальных исследованиях электрофизиологических свойств сердца начали использовать в конце 20 века. Так, Бурнашевым и соавт. был применен метод локальной фиксации потенциала (patch-clamp) для записи ионных токов на одиночных калиевых каналах на тонком срезе сердца новорожденных крыс без ферментативной обработки миокарда [2].
Позднее было показано, что тонкие срезы сердечной ткани у новорожденных крыс обладают более репрезентативными электрофизиологическими характеристиками, чем монослойные культуры [3]. В настоящее время тонкие срезы применяют для исследования вклада различных типов клеток в функцию сердечной ткани, а также для разработки биологических кардиостимуляторов и методов регенерации сердечной ткани с использованием стволовых клеток [4,5,6,11]. Было показано, что такие срезы могут являться адекватной экспериментальной моделью для изучения сложных процессов, происходящих в миокарде, представляя собой компромисс между моделью изолированных клеток и исследованием сердца in vivo [13,14,24,25,27,28].
Тем не менее, при исследовании на тонких срезах возникает ряд новых эффектов. Потенциал покоя клеток в тонком срезе выше, чем в других препаратах [2]. Также для анализа процессов имеет значение толщина среза [1]. У исследователей часто возникает вопрос, насколько записи на тонком, практически плоском срезе соответствуют процессам, протекающим в трехмерном сердце.
При помощи математических моделей сердца, построенных на основании экспериментальных данных, поведение сердца можно проанализировать на различных уровнях организации -от клеточного до органного. При этом, в экспериментах in silico появляется возможность контролировать все внешние условия и регистрировать все интересующие исследователя параметры, что невозможно в исследованиях ex vivo и in vivo. Таким образом, математическое моделирование выступает в роли метода для оценки отличий электрофизиологических характеристик в тонком срезе и в интактном органе. Также при помощи математических моделей можно построить новые протоколы физиологических экспериментов для коррекции внешних условий (температуры, частоты стимулирующего тока, концентрации ионов в растворе в перфузируемой ванне и др.) таким образом, чтобы поведение среза можно было интерпретировать на уровне целого органа с меньшей ошибкой.
Современные компьютерные трехмерные (3D) модели электрической функции левого желудочка (ЛЖ) сердца могут учитывать такие анатомические и физиологические особенности, как реалистичную геометрию, направление волокон, пространственную неоднородность свойств клеток, а также включать свойства среды, в которой находится орган. Дальнейшее развитие таких 3D моделей выявило недостатки существующих численных методов для решения дифференциальных уравнений, на основании этих исследований были выработаны тестовые примеры для верификации таких моделей [17]. Большинство исследований, проведенных на двумерных (2D) моделях, выполнены в методологии редукционизма и не учитывают реалистичную геометрию, направление волокон и т.д. Такие 2D модели, учитывающие неоднородность электрофизиологических свойств клеток, могут найти применение при разработке новых численных методов.
В настоящей статье мы предлагаем компьютерную модель 2D миокардиального среза ЛЖ для сравнения процессов деполяризации и реполяризации, происходящих в таком срезе и трехмерном ЛЖ. Наша модель включает анизотропию среды, соответствующую реалистичному ходу волокон в миокарде, полученному на основании данных диффузионно-тензорной магнито-резонансной томографии (ДТ МРТ).
Мы воспользовались данными ДТ МРТ сердца собаки, полученными ex vivo и предоставленными в открытый доступ университетом Джонса Хопкинса (Johns Hopkins University) [26].
Данные ДТ МРТ представляют собой дискретное трехмерное тензорное поле, полученное на основе измерения величины и направления диффузии молекул воды в миокарде. Каждый тензор является симметричным диадическим тензором диффузии, описывающим диффузию молекул воды в точке измерения. Расстояние между точками измерения по осям X и Y составляет 0.3 мм. Расстояние между точками измерения по оси Z составляет 0.8 мм. Общий размер трехмерного массива в элементах составляет 256×256×130.
Направление оси сердца в данных приблизительно совпадает с осью Z. Предполагается, что преимущественный ход волокон в миокарде в точке совпадает с направлением главной компоненты тензора диффузии в ней.
Рис. 1: (а) Поверхностная сетка, полученная применением алгоритма шагающих кубов, (б) поверхностная сетка, сглаженная с помощью пакета Blender, (в) часть тетраэдральной сетки, построенной в пакете GMSH, (г) совмещение исходных данных с построенной конечноэлементной сеткой, (д) распределение значений двугранных углов в градусах в трехмерной сетке, (е) разметка субэндокарда для дальнейшей активации током стимуляции в 3D модели (выделена фиолетовым цветом).
Для построения 3D расчетной сетки было сделано следующее. Была проведена очистка данных ДТ МРТ от шумов, из данных были удалены предсердия и правый желудочек, затем была получена воксельная геометрия ЛЖ сердца. У каждого тензора диффузии в точках, оставшихся после очистки от шумов, была взята длина главной компоненты, которая позволяет судить о принадлежности точки миокарду или окружающему пространству. Таким образом, 3D тензорное поле было преобразовано в 3D скалярное поле с сохранением местоположения точек измерения. Предсердия и правый желудочек были удалены из данных отсечением по плоскостям. Получившийся трехмерный массив обрабатывался по срезам вдоль оси Z. Таким образом, 3D воксельное изображение представлялось как серия 2D изображений, к которым можно применять классические методы обработки изображений без их модификации.
Значения тех пикселей, которые на изображении соответствовали сердечной ткани, отделялись от пикселей, соответствующих внешней среде, с помощью порогового метода Отсу [10]. Значение порога составило 3.737e-06. К получившемуся бинарному изображению были применены последовательно методы математической морфологии: удаление связных областей с площадью менее 300 пикселей, заливка замкнутых областей с площадью менее 6 пикселей, эрозия с квадратным окном 9 × 9, заливка замкнутых областей с площадью менее 500 пикселей, эрозия с квадратным окном 5 × 5, дилатация с квадратным окном 5 × 5 пикселей.
Первая поверхностная сетка строилась алгоритмом шагающих кубов (marching cube), рис. 1а. Затем поверхностная сетка, полученная на предыдущем этапе и отвечающая за геометрию ЛЖ, была улучшена применением фильтров Blender, рис. 1б. Конечноэлементная тетраэдральная сетка была построена в программе GMSH [8], рис. 1в. В полученной таким образом тетраэдральной сетке якобианы всех элементов положительны, двугранные углы тетраэдров лежат в диапазоне 17°-82°. Максимальная длина ребра -0.07 см, количество тетраэдров в сетке -9 630 862.
Для построения 2D расчетной сетки в качестве тонкого среза миокарда был взят продольный срез ЛЖ. Для этого поверхностная сетка ЛЖ была рассечена плоскостью, проходящей через верхушку ЛЖ и параллельной плоскости OY Z. На полученном 2D срезе была построена треугольная конечноэлементная сетка в программе GMSH. В полученной таким образом 2D модели тонкого среза максимальная длина ребра составила -0.07 см, а количество треугольников -22406. Геометрическое соответствие исходных данных и построенных сеток представлено на рис. 1г.
Для каждого элемента тетраэдральной и треугольной конечноэлементной сетки было задано направление волокна в нем путем покомпонентной интерполяции главных компонент ближайших тензоров диффузии на середину элемента. Для плоской треугольной сетки векторы, задающие направления волокон, были спроецированы на плоскость.
В качестве электрофизиологической модели кардиомиоцита использовалась модель генерации потенциала действия (ПД) клетками ЛЖ собаки Luo-Rudy I [15], описывающая 6 трансмембранных токов: I N aбыстрый натриевый ток, I si -медленный входящий кальциевый ток, I K -исходящий калиевый ток, I K1ток аномального выпрямления, I Kp -исходящий калиевый ток во время фазы плато, I b -ток утечки. Вместе со стимулирующим током I stim эти токи составляют общий ток
Миокардиальная ткань представлена как непрерывная среда, описываемая уравнениями бидоменной модели без бассейна [21]:
где емкость мембраны C m = 0.7 мкФ см 2 [23]; V -трансмембранный потенциал; w -вектор фазовых переменных; φ -внеклеточный потенциал; отношение площади поверхности мембраны к объему β m = 700 [22,20]. При данных параметрах скорость проведения вдоль волокна 0.76 м с и поперек волокна 0.44 м с , что соответствует экспериментальным данным по здоровому миокарду [7]. На границе ∂Ω расчетной области Ω устанавливали граничные условия Неймана с нулевым потоком через границу ∂V ∂n = 0, ∂φ ∂n = 0. Активация длительностью 2 мс проводилась в начальный момент времени. В физиологических экспериментах и при помощи математического моделирования было показано, что при активации возбудимых тканей, если деполяризующий потенциал больше порогового значения, развивается ПД, амплитуда и длительность которого не зависит от амплитуды стимулирующего тока [12]. Поэтому в нашей работе плотность стимулирующего тока была подобрана таким образом, чтобы довести деполяризующий потенциал до порогового значения, и составила −13.5 мA см 2 для каждого узла, попадающего в область стимуляции, в 2D и 3D моделях. Процедура подбора выполнялась вручную отдельно на 2D и 3D моделях с использованием метода деления пополам на отрезке от −70 мA см 2 до −3 мA см 2 . В качестве зоны активации для 3D модели ЛЖ и 2D модели тонкого среза были взяты точки, лежащие на глубине R = 0.1 см от размеченной ранее поверхности эндокарда (рис. 1е).
Для проведения вычислительных экспериментов был использован пакет Chaste [16]. Метод дискретизации уравнения и техническая реализация, используемые в этом пакете, подробно описаны в статье [18].
Для сравнения результатов моделирования в 2D и 3D моделях ЛЖ значения трансмембранного потенциала с тетраэдральной сетки ЛЖ были интерполированы на треугольную сетку среза ЛЖ методом линейной интерполяции. Таким образом, мы сравнивали два скалярных поля V 2D t,i и V 3D t,i , где t -момент времени, а i -номер узла, V 2D t,i -значение трансмембранного потенциала в i-м узле в момент времени t в 2D модели, V 3D t,i -интерполированное значение трансмембранного потенциала в i-м узле в момент времени t в срезе 3D модели ЛЖ.
Для оценки различий в значениях трансмембранных потенциалов ∆V t и потенциалов во внеклеточной среде ∆φ t в 3D и 2D моделях в момент времени t, мы сравнивали значения потенциалов в узлах конечно-элементной сетки, используя нормированное на максимум Евклидово расстояние:
где N -количество узлов в двумерной модели, t -момент времени.
Различие в значениях трансмембранного потенциала ∆V t и потенциала во внеклеточном домене ∆φ t между 2D и 3D моделями в каждый момент времени представлено на рис. 2. Максимальная разница наблюдалась в момент времени 13 мс, а минимальная -в момент времени 50 мс. Помимо этого, наблюдалось значительное возрастание ошибки на промежутке от 300 до 345 мс. Карты различий в значениях трансмембранных потенциалов в отдельные моменты времени (где они наиболее значительны) представлены на рис. 3. Отдельного внимания заслуживает карта потенциалов во внеклеточном домене, представленная на рис. 5. Максимальное различие между моделями составило 57 мВ в момент времени t = 6 мс.
В данной работе мы сравнили электрофизиологические характеристики компьютерных моделей ЛЖ и тонкого среза ЛЖ, основанных на данных ДТ МРТ ex vivo. Были построены карты трансмембранного потенциала, потенциала во внеклеточном домене, карты активации, длительности ПД и карты реполяризации, которые показывают разницу в поведении системы в определенные моменты времени и в определенных регионах модели.
Общее время активации, длительности ПД и времени реполяризации отличаются не более чем на 1-2 мс между 2D моделью миокардиального среза и 3D моделью ЛЖ (см. таблицу) и не более чем на 3-5 мс для отдельных точек в 2D модели и срезе 3D модели ЛЖ (рис. 4). Полученные отличия в 1-2 мс являются не столь значимыми для физиологических экспериментов: так, экспериментально измеряемые время активации, длительность ПД и др. имеют дисперсию более 1-2 мс [9,19]. Тем не менее, по рис. 2 и рис. 4 видно, что отличия могут быть более выраженными в отдельных регионах на 2D срезе и 3D модели. В целом, плоские срезы достаточно близко повторяют процессы в трехмерном теле, однако объемная геометрия и ход волокон приводят к количественным различиям в электрофизиологических свойствах. На картах активации можно заметить зоны, в которых разница максимальна, они локализованы небольшими группами, примерно по 1 см в длину (рис. 4). Данные различия, скорее всего, вызваны влиянием хода волокон. Заметим, что на карте реполяризации таких зон с максимальным отклонением значений не было обнаружено.
Так как в работе использовалась бидоменная модель без бассейна, то не учитывалась утечка заряда из внеклеточного домена. Этим можно объяснить существенную разницу в потенциалах внеклеточного домена (рис. 5). Данный результат должен быть учтен в постановке задач для будущих исследований.
Отметим, однако, что в физиологическом эксперименте можно воспроизвести ситуацию, когда поток
В данном исследовании мы сравнили процессы деполяризации и реполяризации в компьютерных моделях ЛЖ сердца и тонкого среза ЛЖ, в которых реалистичный ход волокон задавался на основании данных ДТ МРТ.
Мы показали, что электрофизиологические процессы в трехмерном теле и срезе схожи, однако существуют промежутки времени и регионы миокарда с существенными отличиями. Процесс активации проявляет больше различий, чем процесс реполяризации. Нами была обнаружена значительная разница в потенциалах во внеклеточном домене в момент распространения волны возбуждения.
Основываясь на полученных результатах, мы полагаем, что различия в геометрии и ходе волокон влияют на процессы деполяризации и реполяризации в миокарде.
Работа поддержана грантом Президиума РАН I.33П. и Постановлением Правительства РФ № 211 от 16.03.2013. Расчеты выполнялись на суперкомпьютере Уральского федерального университета и суперкомпьютере Уран Института математики и механики Уральского отделения Российской академии наук.
Comparison of depolarization and repolarization in mathematical models of the left ventricle and a thin myocardial slice Аnastasia E. Bazhutina 1 , Konstantin S. Ushenin 1,2 , Alexey A. Golovin 1 , Аnastasia D. Khokhlova 1,2 1 -Ural Federal University (Yekaterinburg, Russia) 2 -Institute of Immunology and Physiology (Yekaterinburg, Russia) Keywords: heart simulation, fiber orientation, heart slice, conduction velocity, repolarization.
In this simulation study we compare activation, depolarization and repolarization in the left ventricle (LV) of the heart and in a thin slice of the LV. In these models LV geometry and fiber orientations are realistic and personalized by diffusion tensor magnetic resonance imaging data. We show the general similarity of electrophysiological processes between the LV slice and the three-dimensional LV. The difference in the activation time was more pronounced than the difference in repolarization time, moreover, these differences were most pronounced in individual regions. A significant discrepancy between the models was found in values of the extracellular domain potentials during the excitation wave propagation. Based on the results obtained, we can assume that the differences in the heart geometry and fiber orientation affect the processes of depolarization and repolarization in the myocardium.