Индуцирование дрейфа спиральных волн в двумерной изотропной модели миокарда с помощью внешней стимуляции С.Ф. Правдин Т.В. Незлобинский ИММ УрО РАН (Екатеринбург) ИММ УрО РАН (Екатеринбург) УрФУ (Екатеринбург) УрФУ (Екатеринбург) А.В. Панфилов Гентский университет (Гент, Бельгия) Аннотация Спиральные волны электрического возбуждения возникают во мно- гих физических, химических и биологических активных средах, но особенно важно изучать спиральные волны в миокарде в силу того, что они возникают только при опасных аритмиях и поэтому долж- ны быть устранены. Известно, что при движении к границе среды спиральные волны исчезают. Индукция дрейфа спиральных волн в направлении границы является одной из перспективных стратегий лечения аритмий сердца. В данной работе мы моделируем динами- ку спиральных волн на изотропном квадрате с помощью феноме- нологических и биофизических моделей клеток сердечной мышцы. Внешняя стимуляция реализована через подачу стимулирущего то- ка с точечных и линейных электродов. Найдены траектории дрей- фа спиральной волны. Наблюдались два механизма устранения са- моподдерживающихся волн: дрейф в направлении границы среды и аннигиляция с новой спиральной волной. 1 Введение Механическая работа сердца как насоса, перекачивающего кровь из вен в артерии, неразрывно связана с электрическими процессами, протекающими в клетках миокарда. Каждая отдельная клетка сердечной мышцы, кардиомиоцит, может находиться в возбужденном, рефрактерном (невозбудимом) и покоящемся состоянии. При нормальной работе сердца по всему рабочему миокарду проходят волны электрического возбуждения, которые обеспечивают почти синхронные переходы из одного состояния в другое всех кардио- миоцитов предсердий и, с определенной задержкой, желудочков. Таким образом, рабочий миокард может быть рассмотрен как активная возбудимая среда, ритм волн возбуждения в которой задается внешним во- дителем ритма, синусным узлом в правом предсердии. Кроме вышеупомянутых плоских волн, в миокарде, Copyright c by the paper’s authors. Copying permitted for private and academic purposes. In: A.A. Makhnev, S.F. Pravdin (eds.): Proceedings of the International Youth School-conference «SoProMat-2017», Yekaterinburg, Russia, 06-Feb-2017, published at http://ceur-ws.org 268 как и в других возбудимых средах, могут появляться спиральные волны. Они являются самоподдержи- вающимися источниками патологического возбуждения в миокарде, подавляют сигнал от синусного узла и соответствуют опасным аритмиям сердца: пароксизмальной тахикардии и фибрилляции предсердий и желудочков. В некоторых случах приступ аритмии проходит самопроизвольно без лечения, но если этого не происходит в течение нескольких секунд или минут, то необходимо медицинское вмешательство для купирования приступа. Методы лечения аритмий, связанных со спиральными волнами, могут быть разделены на три крупные ветви: медикаментозные, хирургические и электрические. Лекарственная терапия аритмий направлена на общее изменение свойств кардиомиоцитов с тем, чтобы спиральные волны в такой среде дрейфовали в сторону электрической границы среды и исчезали либо чтобы спиральная волна в такой среде не мог- ла возникнуть. Если известен очаг в миокарде, где возникают разрывы плоских волн и формируются спиральные волны, то его в ряде случаев можно хирургически удалить или коагулировать, превратив в электрически невозбудимую соединительную ткань. Серьезными недостатками лекарств-антиаритмиков являются их побочные эффекты, а иногда – недостаточная эффективность. Основные недостатки опера- тивного лечения – это риски самой операции на сердце и наркоза. Альтернативой этим методам служит электротерапия аритмий. Традиционные методы электролечения аритмий – это дефибрилляция и кардиоверсия, то есть пропус- кание короткого электрического импульса очень высокого напряжения (сотни-тысячи вольт) и большой силы тока. Виды дефибрилляторов: наружные (их электроды накладываются на кожу грудной клетки и, иногда, спины), операционные (во время операции на открытом сердце электроды касаются внешней по- верхности сердца, эпикарда), имплантируемые (такие аппараты миниатюрны, находятся под кожей, а их электроды внедряются в миокард). Недостатками наружной дефибрилляции являются чрезвычайно высо- кие напряжение и сила тока, подаваемого аппаратом, что приводит к болезненности процедуры и повре- ждению миокарда, сходного с электроожогом. Операционная дефибрилляция, хоть при ней и используют ток меньшей силы и напряжения, требует вскрытия грудной клетки. Таким образом, самым щадящим ме- тодом дефибрилляции является такая, которая проводится с имплантированного прибора, что позволяет примерно на порядок сократить электрическое напряжение. Тем не менее, даже в этом случае напряжение составляет сотни вольт, а лечение пациенты сравнивают с очень сильным ударом в грудь. Низковольтная дефибрилляция-кардиоверсия (НВКД) – это современный активно изучаемый метод электротерапии аритмий. Его цель – вытеснить спиральную волну, используя низкое электрическое на- пряжение, подаваемое на электроды (порядка 10 В), таким образом сведя к минимуму повреждающее действие терапии на миокард. Важное применение НВКД находит в лечении предсердных аритмий: па- роксизмальной тахикардии, трепетания и фибрилляции предсердий. Особенность лечения предсердных аритмий в том, что пациент в сознании, в отличие от желудочковых аритмий, когда пациент без сознания. Вследствие этого снижение энергии импульса при предсердной аритмии также очень важно для уменьше- ния дискомфорта пациента и преодоления страха перед процедурой. Идея метода НВКД основана на том, что эти аритмии связаны со спиральными волнами электрическо- го возбуждения в миокарде, поэтому для остановки аритмии необходимо избавиться от спиральных волн. Один из способов сделать это – вызвать дрейф спиральной волны с помощью серии импульсов, подава- емых с одного или нескольких имплантированных электродов. При достаточном приближении кончика спиральной волны к электрической границе миокарда (например, зоне фиброзного кольца, разделяющей и изолирующей предсердия и желудочки сердца) спиральная волна исчезает. Энергетические параметры импульсов при НВКД близки к характеристикам нормальных импульсов плоских волн синусного узла, поэтому НВКД не оказывает повреждающего действия на миокард. Метод НВКД основывается на из- вестном факте из теории автоволн: если в среде сосуществуют два или больше источников возбуждения (автоколебательные источники, спиральные волны и т.д.) с разной частотой, то источник с наибольшей частотой вытесняет все остальные источники. Значит, частота НВКД должна быть несколько выше часто- ты возникшей спиральной волны. Начала теории НВКД для случая предельно плотной спиральной волны изложены в статье [4]. Клинические результаты применения НВКД в случае пароксизмальной тахикардии, когда возникает одна спиральная волна, изложены в статьях [24, 20]. Имеются теоретические результаты в области НВКД спиральной волны, вращающейся вокруг невозбу- димой неоднородности в 2-мерной среде (например, вокруг рубца после инфаркта миокарда) и на 1-мерной модели [19]. В частности, в указанной книге показывается, что стимуляция, генерирующая последова- тельность плоских волн, не может ликвидировать одиночную спиральную волну, вращающуюся вокруг препятствия размером существенно больше размера ядра спирали; получены теоретические оценки макси- 269 мального периода стимуляции, необходимого для устранения спиральной волны, в зависимости от радиуса круглого препятствия. Предложен и другой способ НВКД, основанный на воздействии внешнего электрического поля на весь миокард. В этом случае спиральные волны, вращающиеся вокруг небольших невозбудимых или замедленно возбудимых областей, отделяются от них, начинают дрейфовать и аннигилируют друг с другом или на границе миокарда. Этот метод рассматривается в серии работ [13, 8, 28, 27, 3]. Случай расположения длинного плоского электрода вблизи ядра спиральной волны проанализирован теоретически, включая вычислительный эксперимент, в работе [6]. Стимуляция с электрода, расположен- ного внутри или вблизи ядра спиральной волны, была рассмотрена в работе [29]. Управление автоволнами изучалось не только на примере миокарда, но и в химических средах [14, 30, 9]. Интересен пример создания химических логических элементов «И» и «ИЛИ» с помощью перегородок в ванночке, где происходит реакция [23]. Однако перечисленные работы рассматривали процессы управления с общими моделями возбудимых сред и не учитывали специфических свойств миокарда. В частности, не была проанализирована роль ани- зотропии, нет исследований на реалистичных моделях геометрии желудочков и биофизических моделях клеток миокарда, включающих описание ионных токов и межклеточное взаимодействие (диффузионную компоненту системы реакции-диффузии). Кроме того, никем не была учтена механоэлектрическая обрат- ная связь, оказывающая большое влияние на спонтанный дрейф спиральных волн. Эти факторы являются определяющими для распространения волн в сердце, и неучет их может объяснить несоответствие теории и результатов лабораторных экспериментов, изложенных, в частности, в работах [25, 21, 12]. Клиниче- ские испытания показывают не очень высокую эффективность НВКД, реализованной в существующей аппаратуре. Цель настоящей статьи – изучение дрейфа двумерных спиральных волн, индуцированного внешней сти- муляцией, с помощью математических моделей однородного изотропного неподвижного миокарда. Задачи: • используя модели кардиомиоцита Алиева–Панфилова (AP) [1] и тен Тюшер–Панфилова (TP06) [22], на изотропном квадрате смоделировать динамику спиральных волн без самостоятельного дрейфа спи- ралей при различных значениях натяжения, • сравнить эффект внешнего воздействия в зависимости от количества, расположения, размера элек- тродов: один точечный, два точечных, один линейный электрод, • выяснить, при каких значениях параметров стимуляции происходит «прорыв» внешней стимуляции, • изучить, когда и как происходит взаимодействие плоских волн и основной спиральной волны. 2 Материалы и методы 2.1 Модели, алгоритмы, численные методы Мы воспользовались тремя моделями клетки сердечной мышцы: феноменологической безразмерной мо- делью AP [1] в двух вариантах, APsimple и APmu, и биофизической моделью клетки рабочего миокарда желудочков сердца человека TP06 [22]. Для расчета распространения волн возбуждения в ткани использо- вали монодоменные уравнения реакции-диффузии, описывающие однородную изотропную среду, имеющие следующий общий вид: ∂u ∂~v = D∆u + f (u, ~v ) + Istim (~r, t), = ~g (u, ~v ), ∂t ∂t где u = u(~r, t) – трансмембранный потенциал клетки в точке ~r = (x, y) в момент времени t, D – коэффи- циент диффузии, ∆u = uxx + uyy – Лапласиан на плоскости (x, y), ~v = ~v (~r, t) – вектор остальных фазовых переменных модели, f (u, ~v ), ~g (u, ~v ) – функции, зависящие от модели кардиомиоцита, Istim (~r, t) – внешний стимулирующий ток. В случае модели APsimple система имеет вид: ∂u = D∆u − ku(u − a)(u − 1) − uv + Istim (x, y, t), ∂t ∂v = (u)(ku − v), ∂t 270 ( 1, если u < a; где (u) = η, иначе. При использовании модели APmu система принимает вид: ∂u = D∆u − ku(u − a)(u − 1) − uv + Istim (x, y, t), ∂t   ∂v µ1 v =− + · (v + ku(u − a − 1)), ∂t u + µ2 где k = 8, a = 0.1,  = 0.01, µ1 = 0.2. µ2 = 0.3. Модель TP06 имеет следующий общий вид: ∂u Iion − Istim (x, y, t) = D∆u − , ∂t Cm Iion = IKr + IKs + IK1 + Ito + IN a + IbN a + ICaL + IbCa + IN aK + IN aCa + IpCa + IpK . Здесь внутриклеточные процессы описываются суммой ионных токов Iion = Iion (~r, t); Cm – емкость кле- точной мембраны. Значения параметров модели были взяты из [22]; они соответствуют физиологическим характеристикам кардиомиоцитов. Стимулирующий ток имел величину Ist , подавался на область Ωstim с периодом Tstim импульсами дли- тельности tstim начиная с момента τ0 : ( n o Ist , если (x, y) ∈ Ωstim , t > τ0 , Tt−τ 0 6 Ttstim ; Istim (x, y, t) = stim stim 0, иначе. Величину импульсов заданной длительности мы находили с помощью модели кардиомиоцита. Численно определяли минимальную величину Imin импульса, при которой возникает потенциал действия, и полагали Ist := 4Imin . Момент начала стимуляции выбирали так, чтобы при t = τ0 спиральная волна «контроли- ровала» всю расчетную область. При расчетах на ионной модели клетки длительность стимулирующих импульсов 1.5 мс была выбрана как одна из типичных длительностей импульсов кардиостимуляторов мар- ки Talos (производитель Biotronik). Задачу исследования динамики спиральных волн обычно упрощают, изучая траекторию точки, вокруг которой вращается спираль. Для ее нахождения задают некоторую линию уровня трансмембранного по- тенциала, тогда центр (также используют термины «начало», «кончик», «точка дислокации» [4], англ. tip) спиральной волны соответствует точке перехода переднего фронта волны в задний фронт. Действительно, если следовать по линии уровня, то мы перейдем от точек переднего фронта к точкам заднего фронта, а значит, есть точка перехода между ними. Как было показано в работе [5], эта точка ~rtip может быть аппроксимирована следующей системой уравнений: u(~rtip , t) = u∗ , u(~rtip , t + ∆t) = u∗ , где значения u∗ и ∆t выбираются индивидуально для каждой модели. В нашем случае их значения были следующими: u∗ = 0 мВ для модели ТР06, u∗ = 0.5 для модели АР; ∆t = 10 мс в модели ТР06, 2 модельные единицы (мод. ед.) времени в модели АРsimple, 1 мод. ед. времени в модели АРmu (ниже мы укажем коэффициенты перевода мод. ед. времени T и длины L в мс и мм соответственно). По траектории движения этой точки определяют среднюю скорость дрейфа спиральной волны и тип её динамики. Важной характеристикой спиральной волны является ее (временной) период Tsw . Известно, что период недрейфующей волны совпадает с периодом колебаний фазовых переменных модели вне ядра спирали. Мы находили период спиральных волн как одну десятую времени между максимумами трансмембранного потенциала, разделенными десятью периодами, в точке вне области стимуляции и вне ядра волны. Мы задавали период внешней стимуляции относительно Tsw . Натяжение характеризует склонность спиральной волны к распаду. Мы находили натяжение по мето- дике из работ [15, 16, 17]: решали систему ∂u D ∂u = D∆u + f (u, ~v ) + · , ∂t R ∂x 271 ∂~v = ~g (u, ~v ), ∂t на квадрате x ∈ [0, N ∆r], y ∈ [0, N ∆r], t ∈ [0, T ], N = 600, ∆r = 0.4, при D = 1, радиусе нити вихревого кольца R = 10 (подробнее см. в [17]) с начальными условиями, при которых возникает спиральная волна, и достаточно большим T (400 мод. ед. для APsimple, 500 мод. ед. для APmu, 2 с для ТР06). Находили центр спиральной волны ~rtip (t) и разлагали его в сумму периодической функции (вращения) и линейной функции (дрейф): ~rtip (t) = ~r per tip (t) + V ~ drif t t, ~ drif t = (Vxdrif t , Vydrif t ). Натяжение b2 (обозначение из [2]) находили по формуле где V b2 = −Vxdrif t · R. Мы вычислили натяжение филамента, длительность потенциала действия ДПД-90, скорость одномер- ных волн и временной период спиральной волны во всех использованных в расчетах моделях кардиомио- цита (табл. 1). Мод. ед. времени равнялась T мс, а мод. ед. длины – L мм. Величину T находили, исходя из равенства ДПД-90 в модели AP и модели TP06, а величину L – по равенству скоростей одномерной волны в этих моделях. Натяжение для APsimple, a = 0.10 нам не удалось найти в силу слишком большого радиуса ядра волны. Таблица 1: Параметры и характеристики моделей клетки миокарда Модель, Натяжение, ДПД-90, T, Скорость одно- L, Период параметр a мод. ед. мод. ед. мс мерной волны, мм спиральной времени мод. ед. волны, мод. ед. APsimple: a = 0.03 0.3 6.42 48 1.73 18.9 14.4 a = 0.04 0.25 6.27 49 1.66 20.0 14.6 a = 0.05 0.19 6.14 50 1.60 21.3 15.4 a = 0.06 0.1 6.00 51 1.53 22.7 17.2 a = 0.07 −0.085 5.86 52 1.45 24.4 21.4 a = 0.08 −0.5 5.74 53 1.37 26.3 30.8 a = 0.09 −1.5 5.61 55 1.28 29.3 56 a = 0.10 5.49 56 1.16 32.7 184 APmu a = 0.10 −0.51 25.07 12.2 1.57 5.3 26 TP06 0.6 мм2 /мс 306 мс 1 0.68 мм/мс 1 249 мс Длительность и сила повторных стимулов для использованных нами моделей кардиомиоцитов указаны в таблице 2. Таблица 2: Длительность и сила повторных стимулов Модель клетки Сила стимулов Длительность стимулов APsimple 20 0.1 APmu 3 0.18 TP06 17 1.5 Краевым условием было равенство нулю потока потенциала через границу квадрата. Для создания спиральной волны использовали протокол S1S2. Численный метод – метод сеток на основе явного метода Эйлера. Значение фазовых переменных на следующем шаге по времени находили по формулам " ! # k+1 k uki+1,j − 2uki,j + uki−1,j uki,j+1 − 2uki,j + uki,j−1 k k ui,j = ui,j +dt D + + f (ui,j , ~v i,j ) + Istim ((i dr, j dr), k dt) , dr2 dr2 ~v k+1 v ki,j + dt ~g (uki,j , ~v ki,j ), i,j = ~ где dr – шаг по пространству, dt – шаг по времени, неотрицательные целые числа i = 1, 2, . . . , N – индекс по x, j = 1, 2, . . . , N – по y, k = 0, 2, . . . , M – по t, значения фазовых переменных в фиктивных узлах сетки с i = 0, i = N + 1, j = 0, j = N + 1 находили из краевых условий (uk0,j = uk2,j и т.д.). 272 2.2 Программная реализация Расчеты на модели AP были проведены на разработанной нами в ИММ УрО РАН программе на кластере УрО РАН «УРАН». Для распараллеливания расчетов на несколько ядер одного вычислительного узла мы применили технологию OpenMP. Компилятор – Intel C Compiler «icc». Моделирование низковольтной дефибрилляции включает разработку инструментария, соответствую- щего поставленным в работе задачам. Большое число доступных для варьирования параметров вычисли- тельного эксперимента предполагает выполнение большой серии расчетов с использованием современной высокопроизводительной компьютерной техники. Однако, чаще всего такой подход не предполагает обра- ботку и визуализацию модельных величин в процессе расчета, что оказывается особенно полезным при начальном подборе параметров моделирования. Для визуализации моделируемого процесса во время проведения вычислений мы разработали приложе- ние с графическим интерфейсом (рис. 1). В этой программе задаются параметры используемой модели, численной схемы, начальной и низковольтной стимуляции (через отдельные панели, рис. 2). Одной из ос- новных особенностей приложения является возможность управления параметрами модели и расчета (кроме изменения количества узлов сетки) в процессе вычислений. Вычисления и обработка команд пользователя выполняются в отдельных потоках. В качестве средств разработки мы использовали язык C++11 (версия языка C++) и кроссплатфор- менный фреймворк Qt [18]. Последний позволяет разрабатывать графический интерфейс и предполагает написание объектно-ориентированной архитектуры приложения. Визуализация расчетной сетки осуществ- лялась через интерфейс прикладного программирования OpenGL, для которого в Qt предусмотрена под- держка с возможностью встраивания графической сцены в разрабатываемые приложения. Рис. 1: Графический интерфейс приложения для моделирования эффекта низковольтной дефибрилляции. На панели управления присутствуют поля для ввода параметров, кнопки запуска и остановки вычислений, а также кнопки вызова окна начальной стимуляции и низковольтной стимуляции 273 Рис. 2: Окно параметров протокола начальной стимуляции (слева) и окно параметров низковольтной сти- муляции (справа) 3 Результаты 3.1 Результаты на модели APsimple Следующие параметры модели были постоянными и общими для всех экспериментов: k = 8, η = 0.1, D = 1. Параметры сетки: шаг по времени dt = 0.008, шаг по пространству dr = 0.4, размер сетки 300 × 300. Параметры стимуляции S1S2: стимул S1 подавался на половину квадрата x < 60, а стимул S2 – на половину квадрата y < 60 в момент 10. Внешняя стимуляция подавалась на небольшую область x < 4, y < 4. Начало стимуляции τ0 = 200. Относительный период стимуляции, общее время расчета, моменты начала и окончания дрейфа на границу квадрата и тип ответа спиральной волны на стимуляцию приведены в таблицах 3, 4, 5, 6. Таблица 3: Общее время расчета, тыс. мод. ед. времени, в модели APsimple a Относительный период внешней стимуляции 0.93 0.95 0.97 0.99 1.00 1.01 0.03 10 10 10 20 10 10 0.04 10 10 10 20 10 10 0.05 10 10 20 10 10 10 0.06 10 10 20 10 10 10 0.07 10 10 20 10 10 10 0.08 10 10 10 10 10 10 0.09 10 10 10 10 10 10 0.10 10 10 10 10 10 10 Типичная траектория центра спиральной волны приведена на рис. 3. Воздействие электрода оказывается эффективным при небольших периодах стимуляции: с течением временени его плоские волны занимают все бо́льшую и бо́льшую часть квадрата. В тот момент, когда плоская волна приближается к началу спиральной достаточно близко (на величину около ширины волны от переднего до заднего фронта), начинается индуцированный дрейф спиральной волны. Когда спираль- ная волна почти достигает границы квадрата, она может либо исчезнуть (см. табл. 6, ответ типа А), либо продолжить дрейф вдоль границы квадрата (ответ типа В), либо остановиться и вращаться вокруг непо- движного ядра (ответ типа С). Кроме того, при больших периодах стимуляции (при a < 0.1) и во всех рассмотренных случаях (при a = 0.1), плоские волны электрода не подавляют спиральную волну, то есть внешняя стимуляция оказывается неэффективной (ответ типа D). 274 Таблица 4: Момент начала дрейфа спиральной волны в модели APsimple (мод. ед. времени) a Относительный период внешней стимуляции 0.93 0.95 0.97 0.99 1.00 1.01 0.03 1500 2000 2900 7700 — — 0.04 1500 2000 2800 6700 — — 0.05 1700 2300 3500 — — — 0.06 1600 2300 4000 — — — 0.07 1900 2200 3800 — — — 0.08 1600 2200 3700 — — — 0.09 1700 2300 4100 — — — 0.10 — — — — — — Таблица 5: Момент окончания дрейфа спиральной волны на границу квадрата в модели APsimple (мод. ед. времени) a Относительный период внешней стимуляции 0.93 0.95 0.97 0.99 1.00 1.01 0.03 3000 4100 6500 17000 — — 0.04 3100 4500 7000 16000 — — 0.05 4000 6100 11000 — — — 0.06 5400 8000 15000 — — — 0.07 4700 6500 11500 — — — 0.08 2800 4000 6800 — — — 0.09 2300 3300 6000 — — — 0.10 — — — — — — При снижении периода стимуляции воздействие на спиральную волну становится более эффективным, а стимуляция с периодом 0.99 периода спиральной волны и более неэффективна. Отличие модели APsimple от биофизических моделей кардиомиоцита заключается в том, что эффект Венкебаха [26] (в данном слу- чае – возникновение автоволны при каждой второй стимуляции в силу попадания половины импульсов на период рефрактерности) возникает лишь при существенном снижении периода внешней стимуляции. Для более точного определения наименьшего эффективного периода стимуляции с электрода мы воспользова- лись другой моделью клетки миокарда, а именно – моделью APmu. Рис. 3: Типичная траектория центра спиральной волны при использовании модели АРsimple (a = 0.08, относительный период стимуляции 0.95). Исходное положение начала спиральной волны – около точки (65, 60). Положение стимулирующего электрода показано буквой «э». Индуцирован дрейф вначале влево вниз, а затем вдоль границы квадрата против часовой стрелки. 275 Таблица 6: Тип ответа спиральной волны на внешнюю стимуляцию в модели APsimple a Относительный период внешней стимуляции 0.93 0.95 0.97 0.99 1.00 1.01 0.03 A B B C D D 0.04 B B B C D D 0.05 B B C D D D 0.06 B B C D D D 0.07 B B B D D D 0.08 B B B D D D 0.09 B B B D D D 0.10 D D D D D D Условные обозначения: A – дрейф от электрода и исчезновение спиральной волны, B – дрейф к границе квадрата, затем вдоль границы, C – дрейф к границе квадрата и остановка, D – нет эффекта. 3.2 Результаты на модели APmu Следующие параметры модели были постоянными и общими для всех экспериментов: k = 8, a = 0.1, µ1 = 0.2, µ2 = 0.3,  = 0.01, D = 1. Параметры сетки: шаг по времени dt = 0.008, шаг по пространству dr = 0.4, размер сетки 400 × 400 при использовании точечных электродов и 300 × 300 при использовании линейного электрода. Мы использовали шаг по времени в 3 раза меньше рекомендованного в работе [7] (0.029 мод. ед.), чтобы избежать неустойчивости численного решения задачи. Параметры стимуляции S1S2: стимул S1 подавался на половину квадрата x < 80, а стимул S2 – на половину квадрата y < 80 в момент 44. Период спиральной волны составил 26, общее время расчета – 20 тыс. мод. ед. времени. 3.2.1 Стимуляция с одного точечного электрода Внешняя стимуляция подавалась на небольшую область x < 4, y < 4. Относительный период стимуляции, общее время расчета, моменты начала и окончания дрейфа на границу квадрата и тип ответа спиральной волны на стимуляцию приведены в таблице 7. Таблица 7: Результаты расчетов на модели APmu при внешней стимуляции с одного точечного электрода Номер Отн. период Начало Окончание дрейфа Тип ответа эксперимента стимуляции дрейфа на границу квадрата волны 1 0.7 — — D 2 0.75 — — E 3 0.8 1200 1350 A 4 0.85 1450 1650 A 5 0.9 2000 2200 A 6 0.95 4200 4500 C 7 0.97 7400 7800 C 8 0.99 — — D 9 1.0 — — D Условные обозначения: A – дрейф от электрода и исчезновение спиральной волны, B – дрейф к границе квадрата, затем вдоль границы, C – дрейф к границе квадрата и остановка, D – нет эффекта, E – распад волны. Типичные траектории центра спиральной волны приведены на рис. 4. При слишком низком периоде внешней стимуляции (относительный период 0.7) эффект от нее отсут- ствует; при несколько большем периоде (относительный период 0.75) из плоских волн возникают многочис- ленные спиральные волны. Если период внешней стимуляции еще выше (относительный период 0.8–0.97), происходит требуемое вытеснение спиральной волны плоскими волнами. В экспериментах 3–5 спиральная волна натолкнулась на границу квадрата и исчезла (ответ типа А), а в экспериментах 6 и 7 она, прибли- зившись к границе квадрата, остановилась (ответ типа С). При еще большем периоде внешней стимуляции она вновь теряет эффективность (ответ типа D). 276 Рис. 4: Примеры траекторий центра спиральной волны при использовании модели APmu (относительные периоды стимуляции: слева – 0.85, эксперимент 4; справа – 0.95, эксперимент 6). Исходное положение начала спиральной волны – около точки (105, 85). Положение стимулирующего электрода показано буквой «э». Результат внешней стимуляции: слева – спиральная волна исчезла, дойдя до правой границы квадрата; справа – индуцирован дрейф вправо вниз до границы квадрата, затем вниз до угла, где дрейф спиральной волны остановился. Результаты показывают, что нижняя граница периода эффективной внешней стимуляции составляет около 0.8 периода спиральной волны, а верхняя – около 0.97. При стимуляции с более высокой частотой может возникнуть распад спиральной волны, который в случае волн в миокарде соответствует фибрилля- ции, что являлось бы крайне нежелательным осложнением электротерапии. Стимуляция с более низкой частотой менее эффективна: чем ниже частота, тем больше времени должна занимать внешняя стимуля- ция. Влияние частоты стимуляции на ее эффективность мы изучим более подробно на биофизической модели миокарда TP06 в следующем разделе статьи. 3.2.2 Стимуляция с двух электродов Имплантированные кардиостимуляторы могут иметь не один стимулирующий электрод, а два, поэтому важно изучить особенности индукции дрейфа спиральных волн при стимуляции с двух электродов. Мы рассмотрели два варианта расположения электродов: в смежных вершинах и в серединах противополож- ных сторон квадрата. Результаты расчетов приведены в таблицах 8, 9. Таблица 8: Результаты расчетов на модели APmu, электроды в смежных вершинах квадрата Номер Отн. период Начало Окончание дрейфа Тип ответа эксперимента стимуляции дрейфа на границу квадрата волны 10 0.7 — — D 11 0.75 — — D 12 0.8 920 1150 A 13 0.85 1180 1410 A 14 0.9 1700 2000 C 15 0.95 2800 3100 C 16 0.97 4750 5000 C 17 0.99 18000 18400 B 18 1.0 — — D Условные обозначения: A – дрейф от электрода и исчезновение спиральной волны, B – дрейф к границе квадрата, затем вдоль границы, C – дрейф к границе квадрата и остановка, D – нет эффекта, E – распад волны. Время начала и время окончания дрейфа при стимуляции с двух электродов оказались меньше, чем при использовании одного. Это говорит о том, что два электрода оказываются эффективнее, чем один. Сравнение этих же показателей для разного расположения электродов дает основание сделать вывод о том, что электроды в серединах противоположных сторон квадрата оказываются более эффективны, чем 277 Таблица 9: Результаты расчетов на модели APmu, электроды в серединах противоположных сторон квадрата Номер Отн. период Начало Окончание дрейфа Тип ответа эксперимента стимуляции дрейфа на границу квадрата волны 19 0.7 — — D 20 0.75 — — E 21 0.8 750 1100 A 22 0.85 950 1250 A 23 0.9 1350 1600 A 24 0.95 2550 2850 A 25 0.97 4400 4600 C 26 0.99 14200 14800 B 27 1.0 — — D Условные обозначения: A – дрейф от электрода и исчезновение спиральной волны, B – дрейф к границе квадрата, затем вдоль границы, C – дрейф к границе квадрата и остановка, D – нет эффекта, E – распад волны. в смежных вершинах. В то же время, когда электроды были расположены в смежных вершинах, случаев распада спиральной волны не было, а другое положение электродов вызвало такой распад в одном случае. Эти результаты говорят о том, что, возможно, использование электродов в серединах противоположных сторон более опасно, чем в смежных вершинах, хотя мы не имеем доказательств того, что спиральная волна при расположении электродов в смежных вершинах не распадается. 3.2.3 Стимуляция с длинного линейного электрода Мы рассмотрели также стимуляцию, подаваемую с электрода, занимающего целую сторону квадрата. Этот случай интересен с теоретической точки зрения: в работах [11, 10] приведена формула, позволяющая по периодам спиральной волны и стимуляции предсказать скорость дрейфа спирали как результат взаи- модействия с плоскими волнами. Электрод был расположен у левой стороны квадрата x < 4. Результаты расчетов приведены в таблице 10. Таблица 10: Результаты расчетов на модели APmu, электрод расположен у левой стороны квадрата Номер Отн. период Начало Окончание дрейфа Тип ответа эксперимента стимуляции дрейфа на границу квадрата волны 28 0.775 650 750 A 29 0.8 775 860 A 30 0.825 890 1000 A 31 0.85 1000 1200 A 32 0.875 1070 1250 A 33 0.9 1400 1600 A 34 0.925 1850 2150 C 35 0.95 2900 3150 C 36 0.975 4900 5250 C 37 0.98 6500 6800 C 38 0.9825 7300 7900 B 39 0.985 6600 7000 B 40 0.99 11000 12300 E Условные обозначения: A – дрейф от электрода и исчезновение спиральной волны, B – дрейф к границе квадрата, затем вдоль границы, C – дрейф к границе квадрата и остановка, D – нет эффекта, E – распад волны. Результаты показывают, что наиболее эффективна стимуляция при ее относительном периоде 0.775−0.9, причем, как и в рассмотренных ранее вариантах расположения электродов, время вытеснения спирали растет с увеличением периода стимуляции. При большем относительном периоде стимуляции, до 0.985, 278 спиральная волна не исчезает, а остается прижатой к границе квадрата. Распад волны наблюдается при периоде 0.99. Мы измерили относительную скорость дрейфа спиральной волны без учета вращения вокруг ядра (гра- фик ее компонент V x, V y в зависимости от относительного периода стимуляции на рис. 5). Компонента V x положительна, то есть волна дрейфует в направлении от электрода, и убывает почти до 0 с ростом периода стимуляции (до 0.98), что подтверждает наше наблюдение о снижении эффективности стимуляции с ро- стом периода. Если период стимуляции находится между 0.98 и 0.99, эта продольная компонента дрейфа резко возрастает до 0.1. Поперечная (по отношению к линейным волнам) компонента V y при периодах стимуляции до 0.98 слабо от них зависит, а при повышении периода стимуляции свыше 0.98 резко меняет знак и становится примерно равной 0.1. При этом возрастает неустойчивость всей системы, появляются дополнительные разрывы волн. Рис. 5: Скорость дрейфа на модели APmu с линейным электродом (область стимуляции вблизи границы квадрата x = 0). 3.3 Результаты на модели ТP06 Приведем параметры, использованные нами при расчетах на биофизической модели миокарда. Во всех экспериментах коэффициент диффузии был равен D = 0.154 мм2 /мс (физиологическое значение). Пара- метры сетки: шаг по времени dt = 0.02 мс, шаг по пространству dr = 0.4 мм, размер сетки 400 × 400 (сторона квадрата тогда равна 160 мм – примерный размер окружности экватора желудочков в норме у взрослых). Параметры стимуляции S1S2: стимул S1 подавался на левую часть квадрата x < 32, а стимул S2 – на половину квадрата y < 80 в момент 300. Период спиральной волны составил 249 мс, общее вре- мя расчета – 60 сек. Внешняя стимуляция подавалась на область x < 8, примыкающую к левой стороне квадрата. Относительный период стимуляции и моменты начала и окончания дрейфа на границу квадрата при- ведены в таблице 11. Примеры траекторий центра спиральной волны приведены на рис. 6. Расчеты показали, что вытеснение спиральной волны в модели TP06 происходило по совсем иному механизму, нежели в моделях АР. Внешняя стимуляция приводит к возникновению дополнительных спи- ральных волн, дрейфующих от электрода. В тот момент, когда одна из них приближается к «центральной» 279 Таблица 11: Результаты расчетов на модели TP06 Номер Отн. период Начало Момент исчезновения эксперимента стимуляции дрейфа, с спиральной волны, с 1 0.80 8 10 2 0.84 5 10 3 0.88 4.5 6 4 0.92 7.5 10 5 0.96 5.5 9 6 1.00 7.5 12 7 1.02 8.5 9 8 1.04 10.5 15 спиральной волне, последняя смещается и начинает двигаться в ту или иную сторону (см. графу «Начало дрейфа» в табл. 11). При одном из последующих контактов с другими спиральными волнами «централь- ная» волна аннигилирует (см. графу «Момент исчезновения спиральной волны» в той же табл.). Все остальные спиральные волны дрейфуют и являются довольно короткоживущими, поэтому прекращение внешней стимуляции приводит к исчезновению всех спиральных волн. Интересно отметить, что в моде- ли ТР06, в отличие от обеих вариантов модели АР, внешняя стимуляция была эффективна даже при небольшом превышении ее периода над периодом спиральной волны (эксперименты 7 и 8). Рис. 6: Траектории центра спиральной волны при использовании модели TP06 (относительные периоды стимуляции: A, B – 0.8, эксперимент 1; C, D – 1.02, эксперимент 7). Слева (A, C) – зависимость координат x, y от времени t, x – красные точки, y – синие точки; справа (B, D) – траектории в координатах x, y. Начальное положение спиральной волны – около точки (50, 70). 280 4 Обсуждение результатов, заключение Результаты расчетов на модели APsimple показывают, что внешняя стимуляция при одном и том же относительном периоде стимуляции может быть как эффективной, так и неэффективной в зависимости от параметров 0-мерной модели. При небольшом периоде стимуляции легче всего спиральная волна вы- теснялась тем легче, чем больше по модулю было натяжение. При периоде стимуляции, близком к макси- мальному (т.е. периоду спиральной волны), навязать ритм при натяжении менее 0.2 (мод. ед.) оказалось невозможно. Наблюдался краевой эффект остановки дрейфа спиральной волны у границы квадрата. Дальнейшие расчеты на анатомических моделях сердца покажут, сохраняется ли он в реалистичном трёхмерном случае. Наиболее эффективной в наших экспериментах оказалась стимуляция с наименьшим использованным нами относительным периодом стимуляции для модели APsimple – 0.93 – и периодом 0.8 для модели APmu. Снижение периода стимуляции ниже этого значения приводит к индуцированному распаду спирали, т.е. к ятрогенной фибрилляции, что является тяжелым осложнением лечения, поэтому в практических приложениях метода, по-видимому, нужно избегать минимального периода стимуляции. Большие значения периода, например 0.85 для модели APmu, чуть менее эффективны, но намного более безопасны. Стимуляция с двух электродов, как мы уже отмечали выше, оказалась эффективнее стимуляции с одного электрода. В работах [11, 10] была предложена формула, связывающая относительный период стимуляции плос- кими волнами и относительную скорость вынужденного дрейфа поперек фронтов этих волн линейной зависимостью. Наши результаты на модели APmu хорошо согласуются с ней вплоть до относительного периода стимуляции 0.98. Расчеты на реалистичной биофизической модели TP06 показали, что существует качественно иной путь ликвидации спиральной волны в среде, нежели наблюдаемое нашими предшественниками и нами на других моделях выталкивание спирали от места стимуляции. Наши дальнейшие исследования будут посвящены компьютерному моделированию навязывания ритма с учетом механических деформаций, механоэлектри- ческой обратной связи и анизотропии миокарда. Благодарности Работы по разработке необходимых программ (раздел 2.2) поддержаны субсидией на выполнение го- сударственного задания ИММ УрО РАН (тема № 0387-2014-0043). Остальная часть работы поддержана грантом РНФ 14-11-00702 (ИММ УрО РАН). При проведении расчетов был использован кластер «Уран» ИММ УрО РАН. Список литературы [1] R.R. Aliev and A.V. Panfilov. A simple two-variable model of cardiac excitation. Chaos, Solitons and Fractals, 7(3):293–301, 1996. [2] V.N. Biktashev, A.V. Holden, and H. Zhang. Tension of organizing filaments of scroll waves. Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, 347(1685):611–630, 1994. [3] Bryan J. Caldwell, Mark L. Trew, and Arkady M. Pertsov. Cardiac response to low-energy field pacing challenges the standard theory of defibrillation. Circulation: Arrhythmia and Electrophysiology, 8(3):685– 693, 2015. [4] E.A. Ermakova, V.I. Krinsky, A.V. Panfilov, and A.M. Pertsov. Interaction between spiral and flat periodic autowaves in an active medium. Biofizika, 31(2):318–323, 1986. In Russian. [5] Flavio Fenton and Alain Karma. Vortex dynamics in three-dimensional continuous myocardium with fiber rotation: Filament instability and fibrillation. Chaos: An Interdisciplinary Journal of Nonlinear Science, 8(1):20–47, 1998. [6] Georg Gottwald, Alain Pumir, and Valentin Krinsky. Spiral wave drift induced by stimulating wave trains. Chaos: An Interdisciplinary Journal of Nonlinear Science, 11(3):487–494, 2001. 281 [7] Monica Hanslien, Robert Artebrant, Aslak Tveito, Glenn Terje Lines, and Xing Cai. Stability of two time-integrators for the Aliev–Panfilov system. International journal of numerical analysis and modeling, 8(3):427–442, 2011. [8] D. Hornung, V. N. Biktashev, N. F. Otani, T. K. Shajahan, T. Baig, S. Berg, S. Han, V. I. Krinsky, and S. Luther. Mechanisms of vortices termination in the cardiac muscle. Royal Society Open Science, 4(3), 2017. [9] On-Uma Kheowan, Chi-Keung Chan, Vladimir S. Zykov, Orapin Rangsiman, and Stefan C. Müller. Spiral wave dynamics under feedback derived from a confined circular domain. Phys. Rev. E, 64:035201, Aug 2001. [10] V.I. Krinsky and K.I. Agladze. Interaction of rotating waves in an active chemical medium. Physica D: Nonlinear Phenomena, 8(1):50 – 56, 1983. [11] V.I. Krinsky, A.B. Medvinsky, and A.V. Panfilov. Evolution of autowave vortices. Znanie, Moscow, 1986. (in Russian) = В.И. Кринский, А.Б. Медвинский, А.В. Панфилов. Эволюция автоволновых вихрей. Знание, Москва, 1986. [12] Wenwen Li, Crystal M. Ripplinger, Qing Lou, and Igor R. Efimov. Multiple monophasic shocks improve electrotherapy of ventricular tachycardia in a rabbit model of chronic infarction. Heart Rhythm, 6(7):1020 – 1027, 2009. [13] Stefan Luther, Flavio H. Fenton, and Bruce G. Kornreich et al. Low-energy control of electrical turbulence in the heart. Nature, 475:235–239, 2011. [14] Alexander S. Mikhailov and Kenneth Showalter. Control of waves, patterns and turbulence in chemical systems. Physics Reports, 425(2–3):79 – 194, 2006. [15] A.V. Panfilov, R.R. Aliev, and A.V. Mushinsky. An integral invariant for scroll rings in a reaction-diffusion system. Physica D Nonlinear Phenomena, 36:181–188, June 1989. [16] A.V. Panfilov and A.N. Rudenko. Two regimes of the scroll ring drift in the three-dimensional active media. Physica D Nonlinear Phenomena, 28:215–218, September 1987. [17] A.V. Panfilov, A.N. Rudenko, and V.I. Krinsky. Scroll rings in the three-dimensional active medium with two component diffusion. Biofizika, 31:850–854, 1986. [18] Qt framework. https://www.qt.io/. [19] S. Sinha and S. Sridhar. Patterns in Excitable Media: Genesis, Dynamics, and Control. Taylor & Francis, 2014. [20] Michael O. Sweeney. Antitachycardia pacing for ventricular tachycardia using implantable cardioverter defibrillators:. Pacing and Clinical Electrophysiology, 27(9):1292–1305, 2004. [21] Harikrishna Tandri, Seth H. Weinberg, Kelly C. Chang, Renjun Zhu, Natalia A. Trayanova, Leslie Tung, and Ronald D. Berger. Reversible cardiac conduction block and defibrillation with high-frequency electric field. Science Translational Medicine, 3(102):102ra96–102ra96, 2011. [22] K.H. ten Tusscher and A.V. Panfilov. Alternans and spiral breakup in a human ventricular tissue model. Am. J. Physiol. Heart Circ. Physiol., 291:H1088–1100, 2006. [23] Agota Toth and Kenneth Showalter. Logic gates in excitable media. Journal of Chemical Physics, 103(6):2058–2066, 1995. [24] Mark S. Wathen, Paul J. DeGroot, Michael O. Sweeney, Alice J. Stark, Mary F. Otterness, Wayne O. Adkisson, Robert C. Canby, Koroush Khalighi, Christian Machado, Donald S. Rubenstein, and Kent J. Volosin. Prospective randomized multicenter trial of empirical antitachycardia pacing versus shocks for spontaneous rapid ventricular tachycardia in patients with implantable cardioverter-defibrillators. Circulation, 110(17):2591–2596, 2004. 282 [25] Seth H. Weinberg, Kelly C. Chang, Renjun Zhu, Harikrishna Tandri, Ronald D. Berger, Natalia A. Trayanova, and Leslie Tung. Defibrillation success with high frequency electric fields is related to degree and location of conduction block. Heart Rhythm, 10(5):740 – 748, 2013. [26] K.F. Wenckebach. De Analyse van den onregelmatigen Pols. III. Over eenige Vormen van Allorythmie en Bradykardie, volume 2 of Nederlandsch Tijdschrift voor Geneeskunde. Amsterdam, 1898. (In Dutch). [27] Dan Wilson and Jeff Moehlis. An energy-optimal methodology for synchronization of excitable media. SIAM Journal on Applied Dynamical Systems, 13(2):944–957, 2014. [28] Dan Wilson and Jeff Moehlis. Toward a more efficient implementation of antifibrillation pacing. PLOS ONE, 11(7):1–28, 07 2016. [29] Hong Zhang, Bambi Hu, and Gang Hu. Suppression of spiral waves and spatiotemporal chaos by generating target waves in excitable media. Phys. Rev. E, 68:026134, Aug 2003. [30] V.S. Zykov, A.S. Mikhailov, and S.C. Müller. Controlling spiral waves in confined geometries by global feedback. Phys. Rev. Lett., 78:3398–3401, Apr 1997. 283 Inducing drift of spiral waves in 2D isotropic model of myocardium by means of an external stimulation Sergei F. Pravdin1,2 , Timur V. Nezlobinsky1,2 , Alexander V. Panfilov3 1 – Krasovskii Institute of Mathematics and Mechanics (Yekaterinburg, Russia) 2 – Ural Federal University (Yekaterinburg, Russia) 3 – Ghent University (Ghent, Belgium) Keywords: spiral wave, biophysics, cardiac arrhythmia, pacemaker, low-voltage defibrillation, cardioversion. Spiral waves of electrical excitation appear in many physical, chemical and biological active media, but it is especially important to study spiral waves in the myocardium because they arise only in dangerous arrhythmias and therefore must be eliminated. It is known that spiral waves disappear when moving toward the boundary of the medium. Induction of the drift of the spiral waves towards the boundary is one of the promising strategies for the treatment of cardiac arrhythmias. In this paper, we simulate the dynamics of spiral waves on an isotropic square using phenomenological and biophysical models of cardiac muscle cells. External stimulation is implemented by application of stimulating current from point and linear electrodes. Trajectories of spiral wave drift are found. Two mechanisms for eliminating self-sustaining waves were observed: drift toward the boundary of the medium and annihilation with a new spiral wave. 284