CC BY
67
Информационные технологии Вестник Нижегородского университета им. Н.И. Лобачевского, 2012, №9 5 (2), с. 438-447
УДК 621.391.2:519.72
ПРОГРАММНЫЙ КОМПЛЕКС «ВИРТУАЛЬНОЕ СЕРДЦЕ»
© 2012 г. В.С. Петров1, А.В. Вильдеманов, С.А. Григорьева1, Е.А. Козиной, М.А. Комаров1,
В.А. Костин1’2, А.К. Крюков1, Т.А. Леванова1, И.Б. Мееров1, Г.В. Осипов1
1 Нижегородский госуниверситет им. Н.И. Лобачевского
2
Институт прикладной физики РАН, Н. Новгород
levanova. tatiana@gmail .com
Поступила в редакцию 12.10.2012
Создан программный комплекс «Виртуальное сердце» для моделирования сердечной активности. Комплекс позволяет рассчитывать текущее состояние сердечной активности, его эволюцию, интегральные характеристики, в том числе виртуальную электрокардиограмму (ВЭКГ). В результате использования комплекса создана база данных, содержащая виртуальные электрокардиограммы, снимки пространственно-временной активности, временные реализации напряжений, полученные при нормальной работе сердца и при различных аритмиях.
Ключевые слова: моделирование активности сердца, виртуальная электрокардиограмма, КТ/МРТ, конечно-элементные сетки, электро-механические взаимодействия.
Введение
Изучение механизмов развития различного рода аритмий, разработка методов их диагностики и способов их предотвращения и лечения являются сейчас исключительно важными вследствие того, что в экономически развитых странах сердечно-сосудистые заболевания являются основной причиной смертности. На первом месте среди смертей от сердечно-сосудистых заболеваний стоит так называемая «внезапная сердечная смерть» (ВСС). Вот несколько фактов о ВСС [1]:
- каждый год в мире от ВСС погибает около
3 млн. человек;
- каждый день в России от ВСС погибает около 700 человек;
- не существует специфических предвестников внезапной остановки сердца.
Интенсивные исследования сердечных аритмий ведут медики и биологи, физики и специалисты в области математического моделирования. Так как сердце является динамической системой и происходящие в нем процессы могут быть описаны как эволюция некоторых переменных состояний: электрических мембранных потенциалов, проводимостей ионных каналов, ионных токов, то описание его работы можно получить анализируя соответствующие математические модели. Упрощенно сердечную ткань можно рассматривать как среду, состоящую из автоколебательных и возбудимых элементов-клеток. Тогда математическая модель сердца -
это система очень большого числа (до нескольких десятков миллиардов) обыкновенных дифференциальных уравнений. Эти уравнения описывают электрическую активность сердца. Но для адекватного анализа требуется учесть и механические движения сердечной ткани. А это еще десятки миллиардов уравнений. При этом надо проводить длительные вычисления. Поэтому необходимый компьютерный анализ невозможен без привлечения методов параллельных вычислений и использования современной вычислительной техники.
Сейчас наступает новый этап исследований сердечной деятельности, обусловленный наличием двух факторов: существованием достаточно реалистических моделей сердечной мышцы, базирующихся на новейших данных электрофизиологии, и доступностью высокопроизводительной вычислительной техники, позволяющей эффективно выполнять вычисления, необходимые при моделировании сложных пространственно-временных процессов в сердце. В мире существует несколько исследовательских групп, занимающихся исследованием сердечной деятельности [2-15] и создавших программные комплексы для моделирования сердца.
В настоящей работе описывается программный комплекс «Виртуальное сердце», предназначением которого является моделирование сердечной активности. Данный комплекс можно использовать в научных исследованиях для проведения крупномасштабных вычислительных экспериментов по анализу пространствен-
но-временной динамики в сердце. Комплекс позволяет рассчитывать как текущее состояние сердечной активности, его эволюцию, так и интегральные характеристики, например виртуальную электрокардиограмму (ВЭКГ). Другим важным аспектом использования комплекса является создание с его помощью базы данных, содержащих виртуальные электрокардиограммы, снимки пространственно-временной активности, временные реализации напряжений, полученные при нормальной работе сердца и при различных аритмиях. Такая база данных, оснащенная хорошим интерфейсом, должна помочь врачам при диагностике заболеваний.
1. Функциональное назначение и структура программного комплекса
Для использования программного комплекса требуются специально обработанные данные МРТ. На основе этих данных проводится выделение сердца, построение сетки и сегментация сердца на различные основные составные части: желудочки и предсердия, синоатриальный и атриовентрикулярный узлы, пучок Гиса и волокна Пуркинье, фибробласты и др.
Программный комплекс «Виртуальное сердце» реализует следующую основную функциональность:
- чтение геометрической конечно-элементной модели сердца, полученной на основе данных МРТ;
- ассемблирование конечно-элементных матриц масс и жесткости;
- решение конечно-элементных линейных уравнений;
- решение нелинейных обыкновенных дифференциальных уравнений, описывающих электрофизиологию;
- экспорт конечно-элементных матриц в заданном формате;
- экспорт результатов проведенных вычислений в формате УТК;
- визуализация данных и расчет ЭКГ.
Общую структуру программного комплекса
удобно схематично представлять в виде последовательного конвейера, где каждый следующий этап использует данные и результаты предыдущего шага. Схема данного конвейера представлена на рис. 1.
На первой стадии выполнения конвейера выполняется загрузка геометрии сердца, для которой будут выполняться последующие расчеты. Загружаемая геометрия имеет структуру конечно-элементной сетки. В программном комплексе поддерживаются тетраэдрические линейные конечные элементы. Для реализации данной функциональности создан специальный
Рис. 1. Принципиальная схема программного комплекса
класс MeshReader. Чтение производится из текстовых файлов node и ele. Файл node состоит из строк, которые имеют следующий формат:
Nn x y z,
где Nn - номер узла; x, y, z - координаты узла в трехмерном пространстве. Файл ele состоит из строк, которые имеют следующий формат:
Ne id\ id2 id3 id4 flag, где Ne - номер элемента, id\, id2, id3, id4 - номера узлов, образующих данный элемент. Опциональный параметр flag может описывать, например, принадлежность данного элемента к разным частям сердца.
Результатом работы MeshReader является инициализация внутренних структур, описывающих заданную сетку. Эти структуры реализованы в виде классов MeshNode и MeshElement. Далее загруженная сетка переходит на этап нормализации, которая упорядочивает нумерацию узлов и элементов сетки так, чтобы в этой нумерации не было пропусков. Описанная процедура значительно упрощает работу с данными сетки в дальнейшем. Нормализацию выполняет метод normalize класса Mesh, хранящего сетку во внутреннем формате. Необходимость в нормализации возникает, например, когда из текстовых файлов считываются элементы лишь с определенным значением параметра flag. В этом случае, очевидно, в нумерации считанных элементов могут возникать пропуски.
После нормализации сетка переходит на этап инициализации ассемблера метода конечных элементов (МКЭ ассемблера), реализованного в виде абстрактного класса FemAssembler. Абстрактность класса в данном случае удобна, т.к. конкретные реализации ассемблера могут быть сделаны различными способами. Например, на данный момент может быть использован МКЭ ассемблер, написанный на основе библиотеки GetFem++ [16]. Класс, реализующий эту функциональность, называется GetFemPpAssembler. Так как данная библиотека использует свой собственный формат хранения сетки, то стадия инициализации включает в себя также конвертацию считанной ранее сетки в формат, воспринимаемый GetFem++. Более того, нумерация узлов и элементов в данном формате отличается от исходной, вследствие чего становится необходимым также построение отображения одной нумерации в другую. Это выполняет метод constructGFIndex класса GetFemPpAssembler.
Также нами был разработан и реализован наш собственный алгоритм ассемблирования матриц жесткости и массы для указанной задачи. Для этого создан класс UnnAssembler.
Его функциональность аналогична классу GetFemPpAssembler. Преимущество данного метода состоит в отсутствии необходимости выполнять дополнительные отображения индексов узлов сетки.
После инициализации МКЭ ассемблера начинается непосредственно процедура ассемблирования конечно-элементных матриц: матрицы массы и матрицы жесткости. Матрицы, полученные в результате процедуры ассемблирования, являются сильно разреженными. Вследствие этого эффективным является их представление и хранение в специальном CSR-формате: compressed sparse row - разреженная матрица, сжатая по строкам. Конвертацию матриц в данный формат выполняют ряд методов шаблонного класса FemAssembler.
После того как МКЭ-матрицы прошли конвертацию в CSR-формат, они переходят на этап инициализации центрального модуля данного программного комплекса - модуля vHeart, реализованного в виде класса VirtualHeart. Кроме МКЭ-матриц модуль vHeart инициализирует электрофизиологическую компоненту комплекса, а именно выделяет память, создает объекты, описывающие клетки сердечной ткани, и вызывает методы инициализации данных объектов. Эти объекты являются классами, потомками базового абстрактного класса Cell, и описывают динамические переменные той или иной модели клетки, а также функцию, вычисляющую производные динамических переменных по времени. На данный момент программный комплекс поддерживает любую из моделей, представленных в открытом хранилище CellML, и, кроме этого, имеет ряд предустановленных моделей: LuoRudy1991 [17, 18], GrandiPasqualiniBers2010 [18], TenTusscherPanfilovNobleNoble2006 [18], MaltsevLakatta2009 [18, 19]. Кроме того, модуль vHeart также принимает на вход инициализированный решатель обыкновенных дифференциальных уравнений (класс DESolver; в текущей версии реализован ОДУ-решатель на основе прямого метода Эйлера с адаптивным дроблением временного шага). Также на этапе инициализации основного модуля vHeart выполняется чтение конфигурационного файла vheart.conf, в котором пользователем могут быть объявлены переменные, определяющие некоторые параметры расчетов (проводимости сердечной ткани и межклеточного пространства, максимальное время счета, временной шаг численного метода ОДУ-решателя и др.). Формат строк в файле конфигурации имеет вид: Variable Value,
где Variable - имя параметра, Value - значение.
Далее управление переходит непосредственно к выполнению численных расчетов. За это отвечает метод Solve. Этот метод принимает в качестве аргументов 3 входных параметра: интервал времени в миллисекундах, в течение которого будет выполняться расчет, шаг интегрирования ОДУ-решателя и его тип. Первые два параметра могут быть также заданы в файле конфигурации. В этом случае предпочтение будет отдано именно значениям, считываемым из файла. Третий параметр - тип решателя - может принимать три значения: VirtualHeart::SERIAL, VirtualHeart::TBB и VirtualHeart::MPI, которые соответствуют последовательному и параллельному режимам расчета. Распараллеливание для систем с общей памятью выполняется с использованием Intel Threading Building Blocks (TBB) [20], для систем с распределенной памятью используется интерфейс MPI.
Как видно из схемы на рис. 1, процесс численного расчета представляется циклическим процессом, состоящим из нескольких стадий. Каждая итерация данного цикла соответствует расчету на одном шаге интегрирования ОДУ-решателя. На первой стадии на основе заданного потенциала V в каждом узле сетки вычисляется лапласиан данного скалярного поля. Данная величина используется при подсчете суммарного внешнего тока Isum, входящего в домен кардиомиоцитов в узлах сетки. Затем Isum используется в качестве правой части при решении задачи Пуассона V((am + of + ов)Дфе) = = Isum, где фе - электрический потенциал в межклеточном пространстве. После этого выполняется численное интегрирование дифференциальных уравнений, описывающих электрофизиологию, правые части которых содержат дополнительное слагаемое, описывающие связь между различными клетками. Именно в это слагаемое входят вычисленные раннее величины: AV, ф.
Сохранение результатов вычислений, выполняемое в процессе расчета, производится в формате vtk, совместимом с программой визуализации Paraview.
2. Технологические решения
Основными специфическими технологическими решениями при создании программного комплекса является использование мультидо-менной модели сердца, высокореалистичных биологически релевантных моделей для различных сегментов сердца, а также учет электро-
механических взаимодействий. подробно на этих решениях.
Остановимся
2.1. Мультидоменная модель. Структура миокарда (сердечной ткани) представлена клетками различных типов, электрофизиологиче-ские параметры которых определяются их физиологической функцией. Клетки сердца соединены друг с другом проводящими и упругими связями. Связь кардиомиоцита с фибробластом осуществляется посредством электрического gap-контакта, через который протекает ионный ток. Собственная динамика клетки сердца зависит от обмена ионами с экстрацеллюлярным пространством. Описание и моделирование электрической динамики системы в целом является весьма сложной задачей, в том числе из-за сложной пространственной организации клеток различных типов в миокарде. Ещё один возможный подход состоит в представлении сердечной ткани как совокупности доменов различных типов, которые перекрываются в пространстве и взаимодействуют посредством ионных токов. Выделяются два больших домена: внутриклеточный и внеклеточный. Внутриклеточный домен разделяется на два домена: домен миоцитов и домен фибробластов. В нашем случае динамические процессы в системе обусловлены различными междоменными ионными токами. На рис. 2 представлена структура трехдоменной модели сердца.
s,fib
Рис. 2. Структура трехдоменной модели сердца
На рисунке Ітуо/іЬ - ионный ток между доменами миоцитов и фибробластов, Ітуое и І/іЬіЄ -ионный ток между доменами внеклеточного пространства и доменами миоцитов и фибробластов соответственно, І8,туо, Іц/іЬ, Іце - внешние стимулирующие токи в доменах миоцитов, фибробластов и внеклеточного пространства. Детали можно найти в [21]. Электрическое поле в домене миоцитов может быть представлено следующим образом [21]:
^(^т^Фт ) ^туо + в туоІтуо,я +
+ Р
myo, fib myo, fib
(i)
где cm - внутриклеточный тензор проводимости, ф,„ - внутриклеточный потенциал, -
число gap-контактов между фибробластами и миоцитами в единице объема, pmyo - доля миоцитов на единицу объема. Электрические поля в домене фибробластов и межклеточном пространстве описываются следующим образом:
V(CT f ^ f ) = - Is, fib +
P fib^fib,e Pmyo, fib^myo, fib,
V(°eVi ) = - Is,e P myoImyo,e
P fiblfib,e,
(2)
при этом значения параметров остаются теми же. Связь между внутри- и внеклеточными доменами осуществляется посредством токов через мембрану Ітуо,е и 1/ь,е:
дV
I = C m +1
myo,e myo ^ ion,m ’
дt
I fib,e = Cfib '
(3)
дt
- +1
ion, f
Здесь Vm = ф,„ - фe - трансмембранное напряжение домена миоцитов, Vf = фу - фe - трансмембранное напряжение домена фибробластов. Связь на уровне двух разных внутриклеточных доменов реализуется посредством тока ионов
фm -ф f
Rf
(4)
myo, fib
fib,myo
Здесь Rfib,myo - сопротивление gap-контакта между фибробластом и миоцитом. Из уравнений (1)-(4) получаем следующие выражения:
V(am VVm )+V(cf VVf )+V((Om +C f +Ce )x
f
xVфe ) = - Is,myo - Is, fib - Is,e .
(5)
Таким образом, интегрирование системы может быть выполнено за три шага:
1. Определение общего внеклеточного тока
I = V(c VV ) + V(c fVVf) +
sum,e \ m m' f
+ Is,myo +Is,fib +Is,e .
2. Нахождение внеклеточного потенциала фє (решение уравнения Пуассона)
V((°m + f + CTe )') = -Isum,e .
3. Нахождение напряжений на мембранах клеток в доменах миоцитов и фибробластов Vf,
Vm с учетом конкретных электрофизиологиче-ских моделей динамики этих клеток (решение дифференциальных уравнений в частных производных)
дVm
i
дt C„
V.=J_
дt Cf
I s,myo +V(CmVVm ) + V(<JmV^ )-P
2.2. Электродинамика сердца. Одним из основных компонентов математической модели сердца является модель одной клетки, воспроизводящая эволюцию мембранных потенциалов, ионных токов, концентраций различных ионов. Для различных животных разработаны специализированные модели, так или иначе базирующиеся на модели Ходжкина - Хаксли [22].
Данная модель была разработана в 1952 году с целью математического описания процессов генерации и передачи электрических импульсов в нейроне. Модель является системой обыкновенных дифференциальных уравнений, описывающих динамику изменения трансмембранного потенциала действия, а также так называемых воротных переменных. Данные переменные описывают способность различных ионных каналов клетки к проведению специфичных ионных токов. Уравнения модели имеют вид:
Ст^ = §Ыат^И(и - ЕМа ) + (и - ЕК ) +
т
+ §ь(и - Е1X т = а т (и)(1 - т) в т (и)т, п = а п (и)(1 - п)-в п (и)п,
И = аИ (и)(1 - И)-в И (и)И.
Приведем краткий перечень основных моделей клеток сердца человека.
2.2.1. Кардиомиоцит желудочков. Для моделирования электрической активности кардио-миоцита желудочка человека создано достаточно большое число биологически релевантных моделей. Мы отметим две из них, которые являются наиболее новыми, точными и подробными: модель Финка - Нобла и др. 2008 года [23] и модель Гранди - Берса 2010 года [6]. При создании модели Финка авторы отталкивались от модели Тен - Тушер - Панфилова [24]. В своей работе авторы дали существенно более точное, улучшенное описание двух К-токов, !К1 и ЫБЯО (или 1Кг). Результирующая математическая модель дает улучшенное представление временных реализаций калиевых токов, сильно влияющих на реполяризацию. Таким образом, данная модель, в комплексе с реальной геометрией сердца, потенциально способна воспроизводить Т-зубец ЭКГ с большей точностью.
Наиболее современной моделью кардиомио-цита желудочка человека является модель Гранди - Берса [6]. Ионные каналы и транс-
\
m e myo, fib myo, fib
ч Pmyo
4, fib + V(CT f VVf ) + V(CT f ^e ) + Pmyo, fibImyo, fib
- Iio
P
- Ii
fib
ion, f
портные молекулы моделируются здесь на основе наиболее новых экспериментальных данных миоцитов желудочков человека. Быстро и медленно инактивируемые компоненты тока 10 выражаются таким образом, что позволяют моделировать клетки как эндокарда, так и эпикарда. Также допускается моделирование трансмуральных градиентов протеинов, управляющих кальцием, и натриевых насосов. Данная модель была протестирована авторами в отношении соответствия экспериментальным данным длительности потенциала действия, кривой восстановления (реституции), частотно-зависимого увеличения изменения кальция.
2.2.2. Кардиомиоцит предсердий. В качестве моделей кардиомицита предсердия человека используются две модели, в достаточной мере схожие, однако основанные изначально на разных более ранних моделях. Первая модель - это модель Коуртеманша и соавторов 1998 года [25]. В качестве основы для модели авторы выбрали модель Луо - Руди [26] 1994 года. Однако исходная модель претерпела достаточно глубокие изменения. В частности, они касаются формулировок ионных токов: быстрый натриевый ток, кальциевый ток Ь-типа, выходной калиевый ток и др. Всего данная модель описывается 21 переменной: 12 воротными переменными для описания ионных токов, 3 переменными для моделирования выброса кальция из саркоплаз-матического ретикулума, переменными, описывающими концентрации кальция в цитоплазме, КБЯ и 1БЯ, а также внутриклеточные концентрации натрия и кальция.
Вторая модель кардиомицита предсердия человека, выбранная нами, - модель Найгрена и др. [27] также 1998 года. Данная модель, в отличие от модели Коуртеманша, основывается на модели Линдблада [28]. В ее описание включено 29 переменных: 12 воротных переменных, внутриклеточные концентрации натрия и кальция, концентрации кальция в ограниченном сабсаркомерном пространстве, в КБЯ и в 1БЯ, внеклеточные концентрации натрия, калия и кальция, 2 переменные описывают выброс кальция из БЯ и 5 переменных используются для описания буферизации кальция тропонином и кальмодулином в цитоплазме и кальсеквестри-ном в БЯ. Всего модель включает описание 12 ионных токов.
2.2.3. Клетки Пуркинье. Существует две наиболее актуальные модели клетки Пуркинье человека. Первая модель была опубликована Стюартом и соавторами в 2009 году [29]. Данная модель получена путем модификации моделей Тен - Тушер [24, 30] клеток эндокарда человека.
Авторы основывались на экспериментальных данных Хана и др. [31]. При создании модели было добавлено два дополнительных ионных тока: ток, активируемый гиперполяризацией, I/, поддерживаемый калиевый ток, Т^. Таким образом, всего в модель входят 14 ионных токов. Кроме введения новых токов, были даны новые формулы для ряда токов: 1К1, I\а. Кроме этого, были изменены значения максимальных проводимостей для токов: 1Кг, 1К!!, 1Ыа. Такие модификации обоснованы существенной разницей в кинетике ионных каналов клеток Пуркинье и клеток эндокарда.
2.2.4. Клетки синоатриального и атриовентрикулярного узлов. В качестве наиболее точной модели клетки синоатриального узла сердца человека можно выделить модель Хопер и соавторов 2003 года [32]. Данная модель во многом схожа с модель клетки САУ Жанга [33]. Для моделирования клетки атриовентрикулярного узла используется также модель Хопер, но с коэффициентами масштабирования времени.
2.3. Электро-механические взаимодействия. Необходимым условием адекватного моделирования динамики электрического потенциала в сердечной мышце является самосогласованный учет влияния механических деформаций ткани на текущие в ней токи. Для этого в модель вводится дополнительный ток 1$, отвечающий за электро-механический отклик [34]. Величина добавочного тока 1$ зависит от локальной деформации ткани и равна нулю, если деформации отсутствуют. Зависимость этого тока от деформаций среды может быть сведена к зависимости от локального относительного растяжения среды 1 + ё, причём обычно полагается, что ток пренебрежимо мал в случае локального сжатия среды [34]. В большинстве случаев используется выражение для дополнительного тока вида
15 = о£(т)ё(у-Е5), (6)
где 0$ и Е$ - проводимость и референтный потенциал для этого добавочного тока 1$, 0(ё) -единичная функция Хевисайда.
Величина т характеризует локальное относительное изменение линейных размеров ткани и в случае малых деформаций линейным образом связана с вектором смещения и [35]:
т=Уи. (7)
Сам вектор смещения находится из условий упругого равновесия ткани при наличии в ней активных напряжений (то есть напряжений, связанных с пространственно-временной динамикой электрического потенциала). В предпо-
ложении изотропных механических свойств ткани и изотропных активных напряжений вектор смещения удовлетворяет следующему уравнению [32]:
цУ2и + У[(Х + ц)Уи + Та] = 0. (8)
Здесь X и ц - коэффициенты Ламе изотропной деформируемой среды (коэффициент ц предполагается однородным, Уц = 0), Та - величина активного напряжения в ткани [34, 36]:
дТад = г(У)(ктУ - Та), (9)
где у - отклонение потенциала от равновесного значения; е(У) - пороговая функция, она равна 1 при значениях У, меньших некоторого порогового значения, и близка к 0 при больших значениях; кт - постоянный коэффициент, управляющий величиной возникающих активных напряжений.
Если тангенциальные напряжения пренебрежимо малы (ц = 0), то уравнение (8) может быть аналитически проинтегрировано:
ХУи + Та = -р, где р - внешнее давление на ткань, и в соответствии с (6) и (7)
1$ = -[1 - 0(Та + р)0Г‘(Та + р)( V - Е$). (10)
Таким образом, описанные выше модели для динамики электрического потенциала могут быть расширены уравнениями (6)-(9) или (9), (10) для анализа влияния механических деформаций сердечной мышцы на эту динамику.
3. Характеристики
Программный комплекс «Виртуальное сердце» обладает следующими характеристиками:
- кроссплатформенность (поддержка операционных систем МБ Windows и Ппих);
- поддержка параллельных вычислений в системах с общей памятью с использованием библиотеки 1П;е1 ТВВ;
- поддержка параллельных вычислений на кластерных системах с использованием интерфейса МР1;
- подключение решателей разреженных систем линейных алгебраических уравнений РЕТБс
[37], Intel MKL PARDISO [38], а также решателя, разработанного в [39];
- работа с наиболее детализированными и точными биологическими моделями электрофизиологии клеток. Поддержка любых моделей из архива CellML[18];
- вывод получаемых результатов в формате VTK, что позволяет использовать в качестве инструмента визуализации ПО с открытым исходным кодом - Paraview [40];
- возможность выполнения пред- и постобработки данных:
— чтение и очистка МРТ слайсов;
— сегментация различных тканей сердечной мышцы;
— вычисление ВЭКГ;
- построение конечно-элементной сетки из сегментированных слайсов;
- графический пользовательский интерфейс;
- кроссплатформенная система сборки scons;
- таймеры и замеры производительности;
- возможность сборки в виде библиотеки;
- проверка корректности решения.
4. Результаты применения
Результатом работы программного комплекса «Виртуальное сердце» являются значения интересующих биофизических величин, например потенциала действия мембраны в узлах конечно-элементной сетки для различных моментов времени. Другими словами, на выходе комплекса пользователь получает информацию о пространственно-временной динамике сердца. Например, комплекс позволяет моделировать процесс распространения электрического возбуждения в ткани сердца. В нормальных условиях возбуждение распространяется по стандартному пути: от САУ по ткани желудочков к АВУ, затем по проводящим путям в межжелу-дочковой перегородке в сторону желудочков, и после этого возбуждение охватывает ткань левого и правого желудочков. Данный процесс был успешно смоделирован при помощи нашего программного комплекса. Моментальные снимки данного процесса приведены на рис. 3.
Рис. 3. Распространение возбуждения в сердце с течением времени
Кроме этого, как уже отмечалось выше, программный комплекс позволяет вычислять ВЭКГ процесса. Детали данного вычисления приводятся далее на примере одномерной среды. Данные результаты легко обобщаются на трехмерный случай.
4.1. Вычисление ВЭКГ. Постановка задачи. В данном разделе будут представлены результаты построения виртуальной электрокардиограммы (ВЭКГ) по данным компьютерного моделирования. Целью данной работы является разработка базового метода вычисления ЭКГ по пространственно-временным реализациям элек-трофизиологических процессов, протекающих в сердце. На первом этапе нами были разработаны методы получения ВЭКГ на основе одномерной модели проводящей системы сердца.
Одномерная модель проводящей системы сердца. Как известно, распространение электрического возбуждения в сердце регулируется работой проводящей системы сердца. Проводящая система сердца представлена на рис. 4 (источник: Jaakko Malmivuo, Robert Plonsey. Bioelectromagnetism. Oxford University Press, 1995).
Рис. 4. Проводящая система сердца
Как видно из рисунка, проводящая система представляет собой, в сущности, одномерную среду, по которой распространяется волна возбуждения. Разные участки данной среды обладают различными электрофизиологическими характеристиками, что обусловливает различную форму потенциала действия в этих местах. В состав проводящей системы входят синоатриальный узел (САУ), атриовентрикулярный узел (АВУ), пучки Гиса и волокна Пуркинье (ГП). Кроме этого, возбуждение также распространяется в предсердиях (П) и желудочках (Ж). Таким образом, можно построить эквивалентную одномерную модель проводящей системы
сердца, представленную на рис. 5.
Моделирование. При построении описанной выше одномерной модели проводящей системы сердца использовалась одномерная цепочка связанных динамических элементов, описываемых системой Луо - Руди [17]. При различных конфигурациях контрольных параметров данная модель способна воспроизводить все необходимое разнообразие динамических режимов и форм потенциала действия клеток проводящей системы. В частности, значения параметра GK1 в диапазоне от 0 до 0.07 отвечают автоколебательному поведению системы, причем различные величины параметра соответствуют различным собственным частотам колебаний. Это используется для моделирования САУ и АВУ.
Результаты моделирования: пространственно-временная реализация. В данном разделе представлены результаты моделирования одномерной проводящей системы сердца. Моделирование выполнялось в течение 40 секунд модельного времени, из которых лишь в течение последних 3 секунд выполнялась запись данных. Первые 37 секунд выделялись под переходные процессы в системе. На рис. 6 (сверху) представлена пространственно-временная диаграмма процесса.
Как видно, полученный паттерн активации проводящей среды сердца соответствует биологическим данным (можно сравнить, например, с рис. 3). В частности, наблюдается замедление скорости распространения возбуждения в районах АВУ и САУ и, наоборот, ускорение в области ГП.
Вычисление ВЭКГ. Нами была предложена методика построения ВЭКГ по данным моделирования одномерной проводящей системы сердца. В настоящем разделе приводится описание данной методики. Исходными данными для вычисления ВЭКГ являются значения трансмембранного потенциала клеток ткани в каждый момент времени в каждой точке пространства У(х, і). Вычисление ВЭКГ Е(і) сводится к вычислению интеграла:
E (t) =
i
■ дV (X, t)
V(X, t)dX
max x {V (x, t )}J dt
Другими словами, ВЭКГ представляет собой среднее по пространству значение производной по времени от потенциала действия с весом, равным потенциалу действия. При этом функция E(t) нормируется, приобретая конечный вид. Результат применения данного преобразования к входным данным представлен на рис. 6 (снизу). Видно, что полученная ВЭКГ отображает все ключевые фазы возбуждения сердца,
2
САУ —► Предсердия
АВУ
Гис -Пуркинье
Желудочки
Рис. 5. Одномерная модель проводящей системы сердца
Рис. 6. Вычисленная ВЭКГ
а именно Р-волну, характеризующую возбуждение предсердий; ^ДО-комплекс, описывающий возбуждение желудочков; 7-волну, возникающую во время реполяризации желудочков.
Обобщение на трехмерный случай. Описанная в данном разделе методика построения ВЭКГ может быть легко обобщена на случай использования трехмерной модели сердца. В данном случае функция потенциала действия в каждый момент времени будет представлять собой трехмерное скалярное поле У(х,у, 2, ?), а одномерный интеграл в преобразовании должен быть заменен на трехмерный интеграл по объему.
5. Выводы
Созданный программный комплекс, построенный на основе параллельных алгоритмов и использующий в качестве моделей современные биологически релевантные системы обыкновенных дифференциальных уравнений, позволяет проводить крупномасштабные вычисления, с высокой точностью воспроизводящие динамические процессы в сердце. Высокая степень сегментации сердца позволяет получать достоверные пространственно-временные эффекты, интегральные характеристики работы сердца (например, виртуальную электрокардиограмму), дает возможность получения информации о влиянии тех или иных электрических, механических и медикаментозных воздействий на
сердце. Мультидоменная модель, лежащая в основе электродинамики, позволяет эффективно проводить исследования возможностей управления динамическими процессами с помощью глобальных и локализованных внешних электрических стимуляций. На основе этих исследований возможна разработка принципиально новых способов борьбы с аритмиями. Принципы построения комплекса позволяют эффективно подстраивать модели (меняя значения параметров) под конкретных пациентов и тестировать методики возможного лечения.
Таким образом, применение высокопроизводительной вычислительной техники с использованием параллельных алгоритмов обеспечивает принципиально новые возможности исследования динамики больших систем обыкновенных дифференциальных уравнений, моделирующих сложные процессы в живых системах, и в частности в сердце.
Работа выполнена при поддержке ФЦП «Исследования и разработки по приоритетным направлениям развития научно-технологического комплекса России на 2007-2013 годы» ГК № 11.519.11.2015 и ГК № 11.519.11.2022
Список литературы
1. http://www.ep-society.org/vss2/
2. Maltsev V.A., Lakatta E.G. // Am. J. Physiol. Heart Circ. Physiol. 2009. V. 296. P. H594.
3. Maltsev V.A., Lakatta E.G. // Am. J. Physiol. Heart Circ. Physiol. 2010. V. 298. P. H2010.
4. Maltsev V.A., Lakatta E.G., Vinogradova T.M. // Circ. Res. 2010. V. 106. P. 659.
5. Grandi E., Puglisi J.L., Wagner S. et al. // Biophys. J. 2007. V. 93(11). P. 3835.
6. Grandi E., Pasqualini F.S., Bers D.M. // Journal of Molecular and Cellular Cardiology. 2010. V. 48. P. 112-121.
7. Livshitz L.M., Rudy Y. // American Journal of Physiology: Heart and Circulatory Physiology. 2007. V. 292. P. H2854.
8. Moreno A.P., Abildskov J.A., Sachse F.B. // Annals of Biomedical Engineering. 2008. V. 36(1). P. 41.
9. Stewart P., Aslanidi O.V., Noble D. et al. // Philosophical Transactions of the Royal Society. 2009. V. 367. P. 2225.
10. Weiss J.N., Chen P.S., Qu Zh. et al. // Circ. Res. 2000. V. 87. P. 1103.
11. Hooke N., Henriquez C.S., Lanzkron P., Rose D. // Crit. Rev. Biomed. Eng. 1992. V. 120. P. 127.
12. Trayanova N.A. // Circ. Res. 2011. V. 108. P. 113-128.
13. Ten Tusscher K.H., Hren R., Panfilov A.V. // Circ. Res. 2007. V. 100. P. 87-103.
14. Jie X., Trayanova N.A. // Heart Rhythm. 2010. V. 7. P. 379.
15. Ashihara T., Constantino J., Trayanova N.A. // Circ. Res. 2008. V. 102. P. 737.
16. Библиотека GetFem++. http://download.gna.org/ getfem/html/homepage.html
17. Luo C.H., Rudy Y. // Circulation Research. 1991. V. 68. P. 1501.
18. http://models.cellml.org/electrophysiology
19. Maltsev V.A., Lakatta E.G. //American Journal of Physiology. 2009. V. 286. H594.
20. Библиотека Intel® Threading Building Blocks. http://software.intel.com/ru-ru/intel-tbb
21. Sachse F.B., Moreno A.P., Seemann G., Abildskov J.A. // Ann. Biomed. Eng. 2009. V. 37. P. 874.
22. Hodgkin A., Huxley A. // J. Physiol. 1952. V. 117. P. 500.
23. Fink M., Noble D., Virag L. et al. // Progress in Biophysics and Molecular Biology. 2008. V. 96 (1-3). P. 357. PubMed ID: 17919688.
24. Ten Tusscher K.H.W.J., Noble D., Noble P.J.,
Panfilov A.V. // American Journal of Physiology. 2004. V. 286. P. H1573. PubMed ID: 14656705.
25. Courtemanche M., Ramirez R.J., Nattel S. // American Journal of Physiology. 1998. V. 275. P. H301. PubMed ID: 9688927.
26. Luo C., Rudy Y. // Circ. Res. 1994. V. 74. P. 1071.
27. Nygren A., Fiset C., Firek L. et al. // Circulation Research. 1998. V. 82. P. 63.
28. Lindblad D.S., Murphey C.R., Clark J.W., Giles W.R. // Am. J. Physiol. Heart Circ. Physiol. 1996. V. 271(4). P. H1666.
29. Stewart P., Aslanidi O.V., Noble D. et al. // Philosophical Transactions of the Royal Society. 2009. V. 367. P. 2225. PubMedID: 19414454120.
30. Ten Tusscher K.H.W.J,. Panfilov A.V. // Am. J. Physiol. Heart Circ. Physiol. 2006. V. 291. P. 1088.
31. Han W., Zhang L., Schram G., Nattel S. // Am. J. Physiol. Heart Circ. Physiol. 2002. V. 283. P. 2495.
32. Hoper C. Modellierung der Anatomie und Elek-trophysiologie des menschlichen Vorhofs. Diploma Thesis, Universitaat Karlsruhe, Institut fur Biomedizinische Technik, 2003.
33. Zhang H., Holden A. V., Kodama I. et al. // Am. J. Physiol. Heart. Circ. Physiol. 2000. V. 279. P. H397.
34. Weise L.D., Panfilov A.V. // Phys. Rev. Lett. 2012. V. 108. P. 228104.
35. Ландау Л.Д., Лифшиц Е.М. Теория упругости. М.: Наука, 1987.
36. Panfilov A., Keldermann R., Nash M. // Phys. Rev. Lett. 2005. V. 95. P. 258104.
37. Библиотека PETSc. http://www.mcs.anl.gov/ petsc/index.html
38. Библиотека Intel Math Kernel Library. http://software.intel.com/ru-ru/intel-mkl
39. Козинов Е.А., Лебедев И.Г., Лебедев С.А. и др. Новый решатель для алгебраических систем разреженных линейных уравнений с симметричной положительно определенной матрицей // Вестник ННГУ им. Н.И. Лобачевского. 2012. №5 (2) (В печати).
40. Мультиплатформенное приложение с открытым исходным кодом для анализа и визуализации данных Paraview. http://www.paraview.org/
PROGRAM COMPLEX «VIRTUAL HEART»
V.S. Petrov, A.A. Vildemanov, S.A. Grigoryeva, E.A. Kozinov, M.A. Komarov, V.A. Kostin,
A.K. Kryukov, T.A. Levanova, I.B. Meyerov, G.V. Osipov
Program complex «Virtual heart» that might be used for heart activity modeling was created. The complex allows to calculate current state of the heart activity, its evolution, integral curves, including a virtual electrocardiogram (VECG). Database containing virtual electrocardiograms, spatio-temporal activity snapshots, voltage temporal realizations that were obtained from normal heart activity and during various arhytmia results from our complex.
Keywords: heart activity modeling, virtual electrocardiogram, CT/MRI, finite element mesh, electro-mechanical interactions.
