CC BY
42
Изв. вузов «ПНД», т. 19, № 3, 2011
УДК 621.391.01
МОДЕЛИРОВАНИЕ СЕРДЕЧНОЙ АКТИВНОСТИ НА ОСНОВЕ ОТОБРАЖЕНИЙ*
Часть II. Ансамбль связанных элементов
Е.А. Павлов, Г.В. Осипов
На основе отображения, предложенного в работе [1], исследуется динамика ансамблей связанных элементов в приложении к моделированию пространственно-временных процессов в сердечной мышце. В частности, рассмотрена динамика двух связанных отображений, моделирующих взаимодействие пейсмекера (осцилляторной клетки) и миоцита (возбудимой клетки), а также двух пейсмекеров. Установлено возникновение синхронного режима через совпадение характерных временных характеристик - частот колебаний - при увеличении силы связи. Исследуется динамика цепочек и решеток связанных осцилляторных элементов со случайным распределением индивидуальных частот. Обнаружены эффекты кластерной синхронизации и переход к глобальной синхронизации при увеличении силы связи. Проанализировано распространение импульсов по цепочке, а также концентрических и спиральных волн в двумерной решетке связанных отображений. Изучены характеристики спиральной волны в зависимости от изменения индивидуальных параметров и связи. Представлены результаты сравнения вычислительной эффективности отображения и исходной модели в форме ОДУ
Ключевые слова: Пространственно-временная динамика, синхронизация, образование структур, концентрическая волна, спиральная волна.
Введение
Исследование коллективной динамики крупномасштабных ансамблей связанных элементов является сегодня одной из наиболее важных задач, имеющих как теоретическое, так и прикладное значение для понимания процессов, протекающих в различных сложных радиофизических, механических, биологических, химических системах. В таких системах могут наблюдаться всевозможные динамические эффекты, например, синхронизация колебаний, распространение концентрических и спи* Статья написана по материалам доклада на IX Международной школе «Хаотические автоколебания и образование структур», Саратов, Россия, 4-9 октября, 2010.
116
ральных волн, а также формирование различных пространственно-временных структур. К подобным сложным системам, демонстрирующим практически весь спектр описанных явлений, относится, в частности, сердце.
Безусловно, с развитием современных методов компьютерного моделирования и высокопроизводительной техники изучение таких систем стало возможным, однако даже сейчас исследование больших ансамблей (тем более с физиологически реалистичной топологией) представляет собой процесс, требующий значительных вычислительных ресурсов. По этой причине для моделирования сердечной динамики в [1] была предложена новая вычислительно-эффективная модель в форме точечного отображения, позволяющая воспроизводить основные динамические механизмы сердечной активности и допускающая описание на основе ионных токов.
Настоящая работа представляет результаты численных экспериментов по моделированию пространственно-временной динамики с использованием указанной модели, в том числе по анализу особенностей синхронизации и выявлению некоторых волновых эффектов в ансамблях связанных элементов. Кроме того, при воспроизведении тех или иных динамических эффектов анализируется связь между полученными численными результатами и процессами, происходящими в реальной сердечной ткани.
1. Моделирование электрической связи
Одной из важных компонент в моделировании электрической активности сердечной ткани является выбор типа связи между клетками. Как известно, проведение электрического импульса по миокарду характеризуется внеклеточным электрическим взаимодействием между близлежащими клетками, осуществляемым через гэп-контакты (gap junctions) [2]. Такой способ соединения клеток в миокарде осуществляется с помощью специальных каналов (коннексонов), образуемых шестью белками-коннексинами. Через гэп-контакты могут непосредственно передаваться от клетки к клетке электрические сигналы (потенциалы действия), а также некоторые молекулы.
Моделирование гэп-контактов на различных уровнях сложности предложено в работах [3-5] в приложении к моделям в форме ОДУ. В качестве таковой в большинстве случаев использовался диффузионный тип связи. Таким образом, как и для моделей ОДУ, для учета влияние гэп-контактов в случае модели в форме отображения перепишем уравнения для мембранного потенциала в виде
Vj i = F V) + dAd(Vn), (1)
где функция F(Vn1) описывает индивидуальную динамику ij-й клетки, а параметр d -силу электрической связи между близлежащими клетками. Здесь Ad - дискретный Лапласиан.
Для иллюстрации некоторых эффектов и особенностей сердечной динамики рассмотрим несколько достаточно простых одномерных и двумерных решеток связанных отображений с указанным типом связи (1), моделирующих сердечную ткань в случае нормального и патологического функционирования миокарда.
117
2. Взаимодействие возбудимых и осцилляторных клеток
Вначале рассмотрим динамику малого ансамбля, а именно, двух связанных отображений, моделирующих взаимодействие пейсмекера (осцилляторной клетки) и миоцита (возбудимой клетки). Необходимость детального анализа связанных сердечных клеток различных типов диктуется тем фактом, что данное взаимодействие имеет принципиальное значение при проведении электрического импульса от источника автоматической активности (пейсмекера) к рабочим клеткам миокарда, проводящим электрический импульс (кардиомиоцитам). При проведении моделирования мы продемонстрируем как ситуацию нормального проведения импульсов (1 : 1), так и ситуацию, когда в силу недостаточной силы связи будут различные типы синхронных откликов и, как следствие, проведения импульсов по сердечному волокну (цепочке отображений).
Синхронный отклик возбудимой клетки в отношении S : R оценивался через совпадение соответствующих временных характеристик (частот) для обеих клеток. Для вычисления частоты определялись плоскости пересечения как V^+i > УП,
Г
5-10'4
410'4 1Г) О II
з-ю-4 j 3:2 —Id = 2.5'
210'4 2:1
МО'4
0
1.175 1.180 1.185 d
Угп+1 = —30 и регистрировалось каждое пересечение траектории указанных плоскостей. Частота колебаний каждого элемента вычислялась по формуле
fi = (1г — 1)/Ьиг
(2)
Рис. 1. Взаимодействие пейсмекера (Id=2.5) и возбудимой клетки (Id =0.5)
Рис. 2. Временные реализации двух связанных клеток (пейсмекера (Id=2.5) и возбудимой клетки (Id=0.5)) при различных типах синхронного отклика: a - отклик 2:1; б - отклик 1 : 1
где 1г есть число пересений указанной плоскости, а Ьиг - число итераций между первым и последним пересечением. Установление синхронного режима в отношении S : R можно расценивать как R откликов возбудимой клетки в ответ на S импульсов, поступающих от пей-смекерной клетки. На рис. 1 приведены результаты одного из экспериментов по моделированию двух связанных клеток. Показано изменение частоты колебаний для обеих клеток и уменьшение отношения S : R при увеличении значения связи. На рис. 2 представлены временные реализации для обеих клеток при различных значениях S : R синхронного отклика, отмеченных на рис. 1.
Далее рассматривается цепочка связанных элементов, сконструированная следующим образом: к уже имеющимся двум пейсмекерной и возбудимой клеткам справа добавляется еще 48 возбудимых клеток. Это дает возможность изучить эффект распространения импульса от пейсмекера по возбудимой
118
Рис. 3. Пространственно-временная динамика цепочки связанных отображений при различных значениях силы связи d: a - 0.011807; б - 0.01182
среде. Пространственно-временная динамика такой цепочки при различных значениях d изображена на рис. 3. При достаточно малом d нет проведения колебаний по возбудимой части цепочки. При d=0.011807 (рис. 3, а) только каждый второй импульс от пейсмекера проходит по возбудимой среде. Наконец, при d>0.01182 (рис. 3, б) есть полное проведение импульсов по цепочке (1:1).
3. Цепочки и решетки осцилляторных элементов
Рассмотрим результаты анализа ансамблей, состоящих только из осцилляторных элементов. Исследование таких ансамблей связано, в первую очередь, с моделированием функционирования SA-узла. Известно, что SA-узел является основным водителем ритма в сердечной мышце человека. Он представляет собой особый участок на своде правого предсердия, расположенный в месте впадения верхней полой вены, и состоит из небольшого числа сердечных мышечных волокон, иннервированных окончаниями нейронов из вегетативной нервной системы. В SA-узле зарождается каждая волна возбуждения, которая приводит к сокращению сердечной мышцы и служит стимулом для возникновения следующей волны. В этом смысле синхронная активность клеток SA-узла является одним из ключевых моментов в процессе формирования единого сердечного ритма. И наоборот, рассогласованность (некогерентная активность) участков SA-узла влечет дисфункцию последнего, что может вызвать различного рода аритмии. Все это показывает, насколько важно изучение особенностей синхронизации в ансамблях связанных осцилляторных клеток в данном контексте приложения к моделированию электрофизиологии сердца.
В настоящей работе анализ ансамблей осцилляторных элементов проводился в несколько этапов. Сначала рассматривалась динамика двух связанных отображений, моделирующих взаимодействие двух пейсмекеров (осцилляторных клеток). Каждый из пейсмекеров обладает своей собственной индивидуальной частотой, задаваемой с помощью различных значений параметра Id. Так же как и в случае с возбудимой клеткой частота колебаний оценивалась по формуле (2).
Установление синхронных режимов через совпадение частот при увеличении силы связи d представлено на рис. 4.
Рис. 4. Взаимодействие двух пейсмекеров
119
На втором этапе анализа исследовалась цепочка осцилляторных клеток со свободными граничными условиями. Для этого была проведена серия численных экспериментов по моделированию цепочек из 30-ти локально диффузионно связанных осцилляторных элементов со случайным распределением индивидуальных частот (параметр Id был случайным образом распределен на интервале [2.8;3.6]). Ниже приводятся результаты одной из таких реализаций. На рис. 5 показаны пространственновременные диаграммы эволюции мембранных потенциалов клеток, взятые после временного интервала (ttrans = 106) переходных процессов в цепочке. Распределения соответствующих средних частот fi для каждого из 30-и элементов при различных значениях силы связи d представлены на рис. 6. Можно видеть, что с ростом связи коллективное поведение элементов в цепочке становится все более когерентным, то есть все большее количество элементов становится взаимно синхронизованным.
Рис. 5. Пространственно-временная динамика цепочки связанных отображений при различных режимах синхронизации. При связи d = 0.002 (а) существуют кластеры синхронизации для элементов вблизи концов цепочки и полностью несинхронизованные движения для элементов в середине цепочки. При d = 0.003 (б) имеет место пять кластеров синхронизованных элементов. При d = 0.005 (в) реализуется режим глобальной синхронизации
Рис. 6. Распределение средних частот колебаний для режимов синхронизации, представленных на рис. 5
120
Так, для коэффициента связи d = 0.003 (рис. 5, б) наблюдается пять кластеров синхронизованных элементов. При d = 0.004 - два кластера синхронизации. На границе двух соседних кластеров синхронизации имеют место проскоки (см. рис. 5, а, б), частота появления которых определяется разностью средних частот синхронизации в соседних кластерах. При d = 0.005 (рис. 5, в) существует режим глобальной синхронизации. Таким образом, при увеличении силы связи d наблюдается переход от полностью некогерентного состояния элементов в цепочке к установлению глобальной синхронизации через режим кластерной синхронизации.
Переходим теперь к анализу динамики двумерной решетки 100 х 100 связанных осцилляторных элементов. Так же как и для цепочки случайное распределение индивидуальных частот элементов обусловлено распределением параметра Id на указанном интервале. Граничные условия - свободные концы. На рис. 7 представлено распределение средних частот элементов и соответствующее мгновенное распределение мембранного потенциала при одной из реализаций. Видно, что при слабой силе связи d формируются кластеры синхронизации, однако они состоят из достаточно малого числа элементов и в целом динамика решетки выглядит некогерентной (рис. 7, а и б). При дальнейшем увеличении силы связи d динамика становится более регулярной, размеры кластеров увеличиваются (рис. 7, в и г) Наконец, при d = 0.005 почти во всей решетке имеет место единственный кластер синхронизации, за исключением малых «дефектов», характеризующихся отличными от усредненного значения частотами (рис. 7, д). Пространственная динамика решетки демонстрирует почти регулярные волны, распространяющиеся от лидера - пейсмекера, задающего частоту колебаний всем элементам - по всей решетке с препятствиями на указанных «дефектах» (рис. 7, е). Дальнейший рост значения силы связи d приводит к режиму глобальной синхронизации элементов решетки.
Рис. 7. Распределение средних частот элементов (а, в, д) и мгновенное распределение мембранного потенциала (б, г, е) в решетке при различных значениях силы связи d: 0.001 (а, б), 0.002 (в, г),
0.005 (д, е)
121
4. Решетки возбудимых элементов
В этом разделе рассмотрим ансамбли, состоящие из связанных возбудимых клеток, которые, как известно, играют фундаментальную роль в проведении электрического возбуждения по миокарду. Действительно, исследование таких ансамблей достаточно типично при моделировании проводящих участков сердечной ткани, а анализ волновых процессов, протекающих в таких ансамблях, позволяет адекватно описывать различные режимы сердечной активности.
Рассмотрим двумерную решетку 250 х 250 связанных элементов со свободными концами, моделирующую возбудимую среду. На рис. 8 приведен результат распространения концентрических волн в решетке возбудимых клеток. Источниками концентрических волн в решетке являются элементы с номерами (75,175) и (175, 50). Видно, что два движущихся навстречу друг другу фронта волн, исходящих из указанных источников, сталкиваются и гасят друг друга. Это происходит вследствие характерной для возбудимых сред сильной релаксационности, то есть превышения времени рефрактерности над временем возбуждения. В результате распространяющаяся в сердце волна имеет резкий передний фронт и длинный рефрактерный хвост заднего фронта, где возбуждение невозможно. Две волны поэтому не могут пройти одна сквозь другую и гаснут. Отметим, что концентрические волны, распространяющиеся от пейсмекерной зоны (SA-узла), являются результатом пространственно-временной организации нормальной деятельности миокарда, что обеспечивает нормальные условия сокращения сердечной мышцы.
При определенных условиях возбудимости взаимодействие фронтов волн с препятствиями может привести к разрушению фронта концентрической волны и образованию спиральных волн, частота вращения которых выше частоты нормального следования импульсов возбуждения от SA-узла. Активность спиральных волн может вызывать повторную самостоятельную активацию различных участков миокарда, что связано с появлением различного рода сердечных аритмий [6], таких как, например, пароксизмальная тахикардия, которая затем может перерасти в фибрилляцию.
Моделирование эффекта спиральных волн выполнено на решетке 500 х 500 связанных элементов. Влияние изменения параметров gca и gKd на характеристики спиральной волны и ее проведения в двумерной решетке представлено на рис. 9, а-в
~\V II 0 -20
Й'40
1-60
I'80
50 100 150 200 г
Рис. 8. Распространение концентрических волн в двумерной решетке
и г-е, соответственно, здесь d = 0.2. Рис. 9, ж-и представляет результаты моделирования волн при различных значениях параметра связи d для возбудимых клеток с gca = 0.09, gKd =
0.282. Приведенные результаты по моделированию спиральных волн в двумерной решетке получены с использованием протокола Si — S2 [7]. Вышеуказанный метод достаточно часто используется при моделировании и заключается в следующем. Сначала воздействуют на элементы с одной стороны решетки Si , запуская тем самым плоскую волну,
122
Рис. 9. Спиральные волны в двумерной решетке - пространственная динамика при вариации параметров gca: 0.11 (а), 0.13 (б), 0.15 (в); gKd: 0.27 (г), 0.24 (д), 0.21 (е) и связи d: 0.20 (ж), 0.24 (з),
0.28 (и)
распространяющуюся в одном направлении (слева направо). После того, как волна достигает середины решетки (в нашем случае i = 250), прикладывается дополнительный импульс S2 параллельно распространению фронта волны возбуждения. Этот второй стимул S2 используется с целью «задавить» половину распространяющейся плоской волны, свободный конец которой в этих условиях как раз и дает формирование спиральной волны в решетке.
5. Вычислительная эффективность
Как было указано в [1], важной отличительной особенностью предлагаемой модели является высокая вычислительная эффективность, что позволяет быстро производить численные эксперименты для крупномасштабных ансамблей связанных элементов. Напомним, что модель представляет собой точечное отображение и по этой причине не требует необходимости выбора как самого метода, так и шага численного интегрирования, как в случае моделей в форме ОДУ. В настоящей работе проведено сравнение вычислительной эффективности модели в форме отображения [1] и оригинальной модели Луо-Руди [8], на основе которой она была построена.
Для одного элемента сравнение скоростей счета были проведены для моделирования возбудимых клеток с приложенной стимуляцией. Модельное время со-
123
Сравнение вычислительного времени для различных размеров двумерных решеток
Таблица
Размер решетки (N х N) 100 х 100 250 х 250 500 х 500
Отображение 0.15 мин 0.9 мин 3.65 мин
Модель Луо-Руди 19 мин 115 мин 457 мин
ставляло 50.00 мс. Сравнительные вычисления были выполнены на процессоре Intel Pentium IV 2.80-GHz CPU с 1.00 Gb RAM с использованием шага численного интегрирования At = 0.01 мс для модели Луо-Руди. Для исходной модели (Луо-Руди) время вычисления составляло 52.1 с машинного времени, в то время как моделирование того же эффекта с использованием модели в форме отображения заняло всего 0.2 с. Безусловно, скорость интегрирования модели ОДУ (модель Луо-Руди) может быть увеличена за счет выбора более сложной схемы численного интегрирования или использования адаптивных сеток [9,10]. Однако даже тогда модель в форме отображения демонстрирует более высокую вычислительную эффективность в сравнении с моделями в форме ОДУ. Также было проведено сравнение вычислительной эффективности для двумерного случая, то есть в зависимости от различных размеров решеток, где каждый элемент описывался либо моделью Луо-Руди, либо моделью в форме отображения. Результаты сравнения времен счета для отображения и модели Луо-Руди приведены в таблице.
Заключение
Проведено исследование динамики ансамблей связанных отображений в контексте описания пространственно-временных процессов в сердечной мышце, в том числе ситуации нормального и патологического функционирования миокарда.
Проанализирована система из двух связанных отображений, моделирующих взаимодействие пейсмекера (осцилляторной клетки) и миоцита (возбудимой клетки). Представлены различные типы синхронных откликов S : R возбудимой клетки в ответ на импульсы, поступающие от пейсмекерной клетки. Показано уменьшение отношения S : R при увеличении значения связи d. Проанализировано распространение импульсов от пейсмекера по цепочке возбудимых элементов.
Исследована динамика цепочек и решеток связанных осцилляторных элементов со случайным распределением индивидуальных частот. В ходе численных экспериментов установлено, что в таких ансамблях возможно существование следующих режимов коллективной динамики:
• некогерентная активность элементов ансамбля;
• режимы кластерной синхронизации;
• глобальная синхронизация всех элементов в ансамбле.
Представлен эффект перехода от полностью некогерентного поведения связанных элементов к глобальной сихронизации через режим кластерной синхронизации при увеличении силы связи d.
124
Рассмотрена динамика двумерных решеток связанных возбудимых элементов. Продемонстрировано распространение концентрических и спиральных волн в решетке связанных отображений. Изучены характеристики спиральной волны (форма и величина фронта возбуждения) в зависимости от изменения индивидуальных параметров отображения: gca и дка, а также от силы связи d между элементами решетки.
Представлены результаты сравнения времен вычисления отображения и модели Луо-Руди, на основе которой оно было построено. Сравнение проведено как для одного элемента, так и для различных размеров двумерных ансамблей связанных элементов. Показана высокая вычислительная эффективность модели в форме отображения.
Данная работа выполнена при финансовой поддержке РФФИ (гранты № 08-02-92004, 08-02-97049, 10-02-00940) и при поддержке Федеральной целевой программы «Научные и научно-педагогические кадры инновационной России» на 2009-2013 гг. (контракты № П2308, 02.740.11.5138, П2018).
Библиографический список
1. Павлов Е.А., Осипов Г.В. Моделирование сердечной активности на основе отображений: Динамика одного элемента// Изв. вузов. Прикладная нелинейная динамика. 2011. Т. 19, № 3. С. 104.
2. Rohr S. Role of gap junctions in the propagation of the cardiac action potential // Cardiovasc. Res. 2004. Vol. 62. P. 309.
3. Joyner R.W., van Capelle ///.Propagation through electrically coupled cells. How a small SA node drives a large atrium // Biophys. J. 1986. Vol. 50. P. 1157.
4. Jamaleddine R.S., Vinet A., Roberge F.A. New frequency entrainment scenario in a pair of cardiac cells coupled through a variable resistance (Proceedings of the 18th Annual International Conference of the IEEE Engineering) // Medicine and Biology Society. 1996. Vol. 3. P. 1270.
5. Henriquez A.P., Vogel R., Muller-Borer B.J., Henriquez C.S., Weinggart R., Cascio W.E. Influence of dynamic gap junction resistance on impulse propagation in ventricular myocardium: A computer simulation study // Biophys. J. 2001. Vol. 81. P. 2112.
6. Kleber A.G., Rudy Y. Basic mechanisms of cardiac impulse propagation and associated arrhythmias // Physiol. Rev. 2004. Vol. 84. P. 431.
7. Bernus O., Wilders R., Zemlin C.W., Verschelde H., Panfilov A.V. A computationally efficient electrophysiological model of human ventricular cells // Am. J. Physiol. 2002. Vol. 282. H2296.
8. Luo C.H., Rudy Y A model of the ventricular cardiac action potential, depolarization, repolarization and their interaction // Circ. Res. 1991. Vol. 68. P. 1501.
9. Cherry E.M., Greenside H.S., Henriquez C.S. A space-time adaptive method for simulating complex cardiac dynamics // Phys. Rev. Lett. 2000. Vol. 84. P. 1343.
10. Qu Z., Garfinkel A. An advanced algorithm for solving partial differential equation in cardiac conduction // IEEE Trans. Biomed. Eng. 1999. Vol. 46. P. 1166.
Нижегородский государственный Поступила в редакцию 4.02.2011
университет им. Н.И. Лобачевского После доработки 12.04.2011
125
MODELING OF CARDIAC ACTIVITY ON THE BASIS OF MAPS: ENSEMBLES OF COUPLED ELEMENTS
E.A. Pavlov, G.V. Osipov
The dynamics of coupled maps’ ensembles is investigated in the context of description of spatio-temporal processes in the myocardium. Particular, the dynamics of two coupled maps is explored as well as modeling the interaction of pacemaker (oscillatory) cell and myocyte (excitable cell), and the interation of two pacemakers. Setting of synchronous regime by increasing of coupling strength is considered through a coincidence of their characteristic time scales (characteristic frequencies). Effects of cluster and global synchronization in 1D and 2D lattices of oscillatory cells with a random distribution of individual frequencies are discussed. Impulse propagation in the chain of excitable cells has been observed. Analysis of 2D lattice of excitable elements with target and spiral waves has been made. The characteristics of the spiral wave have been analyzed in depending on the individual parameters of the map and coupling strength between elements of the lattice. Comparative results of computational efficiency with the map-based model and original ODE model are presented.
Keywords: Spatio-temporal dynamics, synchronization, pattern formation, target wave, spiral wave.
126
