CC BY
16
УДК 517.958
И.А. Павельчак1
ЧИСЛЕННЫЙ МЕТОД ОПРЕДЕЛЕНИЯ ЛОКАЛИЗОВАННОГО НАЧАЛЬНОГО УСЛОВИЯ В МОДЕЛЯХ ФИТЦ-ХЬЮ-НАГУМО И АЛИЕВА-ПАНФИЛОВА*
Рассматривается обратная задача для моделей Фитц-Хью-Нагумо и Алиева-Панфилова, описывающих распространение волн в возбудимых средах. Задача состоит в определении локализованного начального условия по измерениям, проводимым на внешней границе плоской области. Предлагается численный метод решения обратной задачи, приводятся результаты численного решения обратной задачи для областей, близких к различным сечениям сердца.
Ключевые слова: обратная задача, модель Фитц-Хью-Нагумо, модель Алиева-Панфилова, локализованное начальное условие.
1. Введение. В настоящее время одну из наиболее важных сфер применения методов математического моделирования образуют проблемы медицины и биологии. Среди разнообразных явлений, изучаемых математическим методом, значительный интерес представляют процессы распространения волн в возбудимых средах, в частности распространение нервных импульсов в миокарде. Для математического описания процессов распространения импульсов широкое распространение получили модели Фитц-Хью 11агумо [1-3] и Алиева-Панфилова [4], представляющие собой начально-краевые задачи для эволюционных уравнений в частных производных. В этой работе рассматривается численный метод решения обратной задачи для модели Фитц-Хью 11агумо и Алиева-Панфилова, состоящей в определении локализованного начального условия по измерениям, проводимым на внешней границе области.
2. Постановка прямой и обратной задач. Рассмотрим модель Фитц-Хью-Нагумо [1-3] в случае двух пространственных переменных. Требуется найти функции и(х,у,1), го(ж, у, являющиеся решением начально-краевой задачи
щ = DAu — и(и — а) (и — 1)
Wt = /Зи — jw,
^(x,y,t) = О,
и(х,у, 0) = (р(х,у), , у, 0) = 0,
w,
где G
(x,y)eG, ie[0,T], (1)
(x,y)eG, ie[0,T], (2)
(x,y)€ г, t € [0,T], (3)
(x,y) G G, (4)
(x,y) G G, (5)
заданные положительные постоянные,
ограниченная область с границей Г, D, а, /3, j заданная функция. Функция и(х, у, t) представляет собой трансмембранный потенциал; функция w(x, у, t) — медленную восстанавливающую переменную, связанную с ионными токами [1-3]. Задача (1)-(5) достаточно хорошо описывает процесс распространения возбуждения в миокарде [5]. Приведем пример численного решения задачи (1)-(5), иллюстрирующий эффект распространения исходного локального возбуждения. Пусть граница области Г представляет собой эллипс с большой полуосью 100 и малой полуосью 80, D = 1, а = 0.15, ¡3 = 0.005, 7 = 0.025, а функция (р(х,у) = +у описывает локализованное начальное возбуждение. На рис. 1 изображены
графики функции и(х, у, t) для четырех моментов времени. На них видно, как формируется и распространяется возбуждение в форме волны, а также, как оно проходит через границу области.
Недостатком модели Фитц-Хью-Нагумо является то, что в ряде случаев она дает отрицательные значения потенциала u(x,y,t) и неточно воспроизводит его форму [4, 5]. Рассмотрим модель
1 Факультет ВМК МГУ, асп., e-mail: pavelchakiQgmail.com
* Работа выполнена при финансовой поддержке РФФИ, проект № 11-01-00259.
Рис. 1. Прямая задача для модели Фитц-Хью-Нагумо: а — £ = 0; б — Ь = 50; в —
Ь = 150; г — Ь — 200
Алиева-Панфилова в случае двух пространственных переменных. Требуется найти функции и(х, у, £), являющиеся решением начально-краевой задачи
щ = О А и
Щ = - ( ^о +
ки(и — а) (и — 1)
и + Ц2
ит,
(гу + — а — 1)),
^ / ч
и(х,у, 0) = <р(ж,у), = 0,
(ж, у) € С, (ж, у) € С,
(ж, у) 6 Г,
I 6 [0,Т],
< е [0,Т],
< е [0,т],
(ж, у) € С, (ж, у) € С,
(6)
(7)
(8)
(9) (10)
где О — ограниченная область с границей Г, £), /г, а, £о, Мъ М2 — заданные положительные постоянные, (р(х,у) — заданная функция. Как и в модели Фитц-Хью-Нагумо, функция и(х,у^) представляет собой трансмембранный потенциал, функция £) — медленную восстанавливающую перемен-
ную, связанную с ионными токами.
Модели Фитц-Хью-Нагумо и Алиева-Панфилова используются для анализа различных процессов, происходящих в миокарде. На их основе было проведено моделирование электрической активности сердца, в том числе исследованы случаи возникновения и распространения спиральных волн, например в [4, 6, 7]. В ряде работ для моделей Фитц-Хью-Нагумо и Алиева-Панфилова рассматривались обратные задачи, состоящие в определении входящих в уравнения параметров и функций по измерениям на границе области, например [5, 8]. Мы рассмотрим задачу определения неизвестного начального возбуждения (р(х,у) по измерениям функции и(х^у^) на границе области Г или ее части Г].:
и{х,у,г) = ф(х ,у,£),
(х,у) егх С Г.
(И)
Эта обратная задача имеет плохую обусловленность, в том числе из-за некоторых особенностей рассматриваемых моделей. Параметр а в уравнении (1) или (6) представляет собой порог возбуждения для мышечных клеток. Как показывают анализ моделей и вычислительные эксперименты, при начальном потенциале ср(х,у) < а У(х,у) Е С? потенциал и(х,у^) быстро стремится к нулю с ростом £ на всей области (7, т.е. возбуждения не происходит, а следовательно, граничное условие (11) не
содержит существенной информации о начальном возбуждении. Если же для некоторых (ж, у) € О <р(х, у) > а, то происходит возбуждение среды и распространение решения в виде волны, но максимум волны при некотором > 0 не зависит от максимума начального распределения потенциала. Так, например, начальные распределения потенциала
ф, у) = 0.5е-((®-зо)2+у2)/1б? у) = е-((х-ж?+у*)Ц^ ф у) = 2е-((ж-зо)Ч?,2)/1б
при Б = 1, а = 0.15, ¡3 = 0.005, 7 = 0.025 дают решения задачи (1)-(5), аналогичные показанным на рис. 1,б-г и отличающиеся малым сдвигом волны, но не ее максимумом. Учитывая некорректность поставленной обратной задачи и возможную сложность геометрии области С?, в которой решается задача (1)-(5) или (6)—(10), рассмотрим обратную задачу в предположении, что в начальный момент времени значение потенциала в области О определяется следующим образом:
<р{х,у,Цх,у,о) = Ае-((ж"4)2+(^)2)/'72,
где А — заданная постоянная. Таким образом, обратная задача состоит в нахождении параметров функции начального распределения по измерениям потенциала на части границы области. При решении обратной задачи величину А будем брать равной максимуму потенциала М = тах(0(ж, у, ¿)) на наблюдаемой границе IV Итак, обратная задача для модели Фитц-Хью Нагумо состоит в нахождении параметров ж, у, а по заданной функции ф(х,у,1) (11) и при известных параметрах I), а, /3, 7. Аналогично ставится обратная задача для модели Алиева-Панфилова.
3. Численный метод решения обратных задач. Рассмотрим численный метод решения сформулированных обратных задач. Рассмотрим модель Фитц-Хью Нагумо. Пусть для точной функции ф(х, у, ¿), (ж, Ь) е Гь г е [0, Т], существует решение обратной задачи <р(х, у, Ц ж, у, ст), определяемое параметрами ж, у, д. Однако функция ф(х,у,1) неизвестна, а задана функция ф$(х,у,1), такая, что
т
ф(ж,у,¿) — ф$(ж,у,Щ2(И(И ^ б2.
0 Г1
Обозначим через и(ж, у, ж, у, а), го(ж, у, Ц ж, у, ст) решение задачи (1)-(5) для заданных значений ж, у, а, определяющих функцию ср(х,у;х,у,а). В качестве приближенного решения обратной задачи будем рассматривать значения ж, у, а, для которых
г
[и(ж, у, ж, у, ст) — фв(х, у, ¿)]2ШсЙ ^ 52.
0 Г1
Таким образом, решение обратной задачи сводится к минимизации функции
г
Ф(ж, у, с) = J ! (и(х, у, ж, у, ст) — ф(х, у, 1))2ё1сИ
0 Г!
с естественными ограничениями на переменные ж, у, а: (ж, у) € О, 0 < а ^ стах, где стах определяется диаметром рассматриваемой области О. Для минимизации Ф(ж, у, о) будем использовать метод градиентного спуска. В связи с этим рассмотрим вопрос о вычислении градиента функции Ф(ж, у, а). Для удобства записи переобозначим ж = р\, у = Р2, о = Рз- Частные производные функции Ф(р1,р2;рз) находятся следующим образом:
г
д /* /* ди
ф-Ф(РьРг,Рз) = I I 2—(х,у,Р,р1,р2,рг)(и(х,у,Р,р1,р2,рг) - Ф(х, у, ¿)) (¿Ш,
0 Г1
где Рг — одна из переменных ж, у, а, а
ди
— (х,у,1;р1,р2,рз) = Ц'г(х,у,р,р1,р2,рг)
ОРг
является решением системы
т т
дШг
т
т
дп
= П\1- I 'ДЗ»1' - 2(1 + а)и + а] - И7* = Риг-ч\Уг,
(х,у,1) = О,
Щх,У, 0) = -^(ж,у), от
1Щж,у,0) = 0,
(ж,у)еС, [О,Т],
(ж,у)еС, ге[0,т],
(ж,у)еГ, ¿6 [О, Г],
(ж, у) € С, (ж, у) € С.
С помощью вычисленного таким образом градиента производится переход от (хп, уп, оп) к (жп+1, Уп+1, сп+1). Итерационный процесс останавливается как только выполняется неравенство Ф(жп,Уп,0п) <
Для модели Алиева-Панфилова обратная задача решается аналогично с использованием следующей системы для |р(ж,у,¿;р1,р2,р3) = ^¿(ж,у,ПРИ вычислении частных производных функции Ф(рьр2,Рз):
= />Д/- /,•/ 'ДЗ»1' - 2(1 + п) // + а] - ;/И', -
а*
аЖ;
дЬ \ и + [м 2
+ /Х2) - /¿1«^
= - Ио
//1 го
(и + /х2):
(И^ + Щ(2«-а - 1))-
(го + — а — 1)),
т
дп
ж, у, 0) = О,
^¿(ж,у,0) = ^(ж ,у),
дрг
[Щж,у,0) = 0,
(ж,у)еС, ¿е[0,т],
(ж,у)еС, ге[0,Т], (ж, у) € г, ге[о,т],
(ж, у) € С, (ж, у) € С.
Рассмотрим вопрос о выборе начального приближения для неизвестных значений ж, у. Из рис. 1 видно, что задача (1)-(5) с локализованным источником описывает распространение кольцевой волны. Опишем алгоритм выбора начального приближения в предположении, что область С?, в которой решается задача, ограничена снаружи выпуклой кривой Г1 и, возможно, имеет внутреннюю границу, состоящую из не более чем двух замкнутых выпуклых кривых. В качестве первого приближения возьмем значения ж0, уо, полученные следующим образом. Обозначим через М наибольшее значение функции у, I) при (ж, у) € Г1, £ е [О, Т1]. Введем вспомогательную функцию г (ж, у), (ж, у) € Г1, принимающую значения наименьшего времени для которого в данной точке границы у, I) ^ 0.9М, и равную Т, если фй(х,у,£) < 0.9М в точке (ж,у) Ш € [0,Т]. Зададим на внешней границе Г1 точки (ж1,у1), (ж2,у2), (жз, уз), такие, что они делят ее на равные по длине дуги. Найдем точку (ж4,у4), такую, что для нее значение функции т(ж, у) будет минимальным. Выберем из (жг,уг), г = 1,... ,4, три точки (аг,Ьг), г = 1,2,3, следующим образом. Если расстояние между (ж4,у4) и одной из трех точек (ж¿,Уг); г = 1,2,3, меньше где Ь — длина внешней границы, возьмем (а*,^) = (ж¿,Уг), г = 1,2,3. В противном случае возьмем в качестве (щ, Ь^), г = 1,2, 3, точку (Ж4, у 4) и две ближайшие к ней точки из (ж¿,Уг), г = 1,2,3. Положим в качестве первого приближения к неизвестным ж, у
хо_1( г(а2,Ыа1 +т(а1,Ь1)а2 г(а3,Ь3)а1 + т(а!,Ь1)а3 т(а3,Ь3)а2 + т(а2,Ь2)а3\
3 \ т(а1, 61) + т(а2, Ь2) т(аь Ъг) + т(а3, Ь3) г(а2, Ь2) + г(а3, Ь3) /'
0 _ 1 /т(а2,Ь2)Ь1 + т(а1, Ь\)Ь2 т(а3,Ь3)Ь1 + т(а1,Ь1)Ь3 т(а3,Ь3)Ь2 + т(а2,Ь2)Ь3\
У 3 \ т(а1, 61) + т(а2, Ь2) т(аь 61) + т(а3, Ь3) г(а2, Ь2) + г(а3, Ь3) /'
Далее это приближение улучшается минимизацией функции
3 _ 2
F(x, у, X) = Y^ bi) - у/(х - di)2 + (у - 6г)2) ,
г = 1
где (di,bi) — выбранные ранее три точки, а переменная А представляет собой приближение скорости распространяющейся волны. Минимизация проводится методом градиентного спуска, в качестве первого приближения А0 берется его оценка для моделируемой среды. Итерации проводятся до заранее заданной точности либо до заданного их максимального числа. Полученная в результате точка (^start? 2/start) берется в качестве начального приближения для градиентного метода решения обратной задачи xq = xstart? У о = 2/start- В качестве начального приближения для а выбирается некоторое значение сто из интервала (0,сгтах).
4. Вычислительные эксперименты. Описанный численный метод решения обратных задач был применен для определения локализованного начального условия в моделях Фитц-Хью-Нагумо и Алиева-Панфилова. Во всех вычислительных экспериментах параметры модели Фитц-Хью-Нагумо были равны D = су = 0.15, ß = 0.005, 7 = 0.025. Параметры модели Алиева-Панфилова были равны D = 1, к = 8, а = 0.15, бо = 0.002, /¿х = 0.2, ¿¿2 = 0.3. Задачи (1)-(5) и (6)—(10) решались в областях G, приближенных к различным сечениям сердца (см. рис. 2-4), с помощью метода конечных элементов; для программной реализации использовалась библиотека deal.II*. Число конечных элементов при расчетах бралось порядка 150000. Метод хорошо проявил себя для выпуклых областей с вырезами в них. На рис. 2-4 показаны точки на границе, которые участвовали в отыскании первого приближения, и точные начальные координаты ж, у.
Рис. 2. Область без вырезов Рис. 3. Область с одним вырезом Рис. 4. Область с двумя вырезами
Рассмотрим схему первого вычислительного эксперимента по решению обратной задачи. Прямая задача для модели Фитц-Хью-Нагумо (1)-(5) решалась в области без внутренних границ, изображенной на рис. 2, с заданной начальной функцией (ж, у) = )/36. В результате вычислялась i/j(x,y,t) на границе (ж, у) Е Гх, t Е [0,Т]. В нее вносилась погрешность и получалась £). Затем с этой функцией решалась обратная задача с использованием описанного численного метода. Описанная схема применялась и для решения обратной задачи для модели Алиева-Панфилова. В табл. 1 приведены результаты решения обратной задачи для обеих моделей в случае ¿2 = o.oiHV'll2.
Более интересны эксперименты в областях, имеюших внутренние границы. Для таких областей измерения u(x,y,t) = i/j(x,y,t) по-прежнему проводятся только на внешней границе Гх- В табл. 2 приведены результаты решения обратной задачи для обеих моделей в области, изображенной на рис. 3 с заданной функцией (р(х,у) = е~((ж+35) +(у+45) )/36 и погрешностью 52 = 0.01||^||2.
Приведем примеры решения обратной задачи в области с двумя вырезами. В табл. 3 приведены результаты решения обратной задачи для обеих моделей в области, изображенной на рис. 4, с заданной функцией ср(х,у) = е~(х +(^~20) )/36 и погрешностью 52 = 0.01||^||2.
В табл. 4 приведены результаты решения обратной задачи для обеих моделей в области, изображенной на рис. 4, с заданной функцией ср(х,у) = е~(х +(^-20) )/36 и погрешностью 52 = 0.1||^||2-
*А Finite Element Differential Equations Analysis Library (http://www.dealii.org/).
Табл и ца 1
Результаты в области без вырезов
Модель Фитц-Хью-Нагумо Модель Алиева-Панфилова
X У а X У а
Точные параметры 15 5 6 15 5 6
Первое приближение 20.5955 7.72449 10 21.6981 8.26497 10
Результат 15.154 5.06008 6.02941 14.9258 4.98226 6.04139
Таблица 2
Результаты в области с одним вырезом
Модель Фитц-Хью-Нагумо Модель Алиева-Панфилова
X У а X У а
Точные параметры -35 -45 6 -35 -45 6
Первое приближение -33.9044 -51.4028 10 -37.5176 -51.8835 10
Результат -35.7705 -46.5982 5.83098 -35.0571 -45.5898 6.28534
Табл и ца 3
Результаты в области с двумя вырезами
Модель Фитц-Хью-Нагумо Модель Алиева-Панфилова
X У а X У а
Точные параметры 0 20 6 0 20 6
Первое приближение 1.45972 24.9571 10 11.7676 26.0682 10
Результат 0.371937 19.9782 5.9758 0.144125 19.9833 6.0797
Приведенные примеры показывают, что увеличение погрешности в исходной информации не дает существенного увеличения погрешности в решении обратной задачи. Это объясняется тем, что наиболее информативной частью исходных данных для решения обратной задачи является время прихода волны в данную точку границы.
Таблица 4
Результаты в области с двумя вырезами
Модель Фитц-Хью-Нагумо Модель Алиева-Панфилова
X У а X У а
Точные параметры 0 20 6 0 20 6
Первое приближение 2.0381 24.2992 10 10.9539 25.8233 10
Результат 0.44535 19.9716 6.00481 0.299884 19.9964 6.00032
Приведенные результаты вычислительных экспериментов показывают, что предложенный численный метод решения обратной задачи дает достаточно хорошие результаты как в областях с простой геометрией границы, так и в областях достаточно сложной формы.
СПИСОК ЛИТЕРАТУРЫ
1. FitzHugh R. Mathematical models of threshold phenomena in the nerve membrane // Bull. Math. Biophysics. 1955. N 17. P. 257-278.
2. FitzHugh R. Impulses and physiological states in theoretical models of nerve membrane // Biophysical J. 1961. N 1. P. 445-466.
3. Nagumo J., ArimotoS., Yoshizawa S. An active pulse transmission line simulating nerve axon//Proc. IRE. 1962. N 50. P. 2061-2070.
4. Aliev R.R., Panfilov A.V. A simple two-variable model of cardiac excitation // Chaos Solutions and Fractals. 1996. 7. N 3. P. 293-301.
ВЕСТН. МОСК. УН-ТА. СЕР. 15. ВЫЧИСЛ. МАТЕМ. И КИБЕРН. 2011. № 3
13
5. Sundnes J., Lines G.T., Cai X. et al. Computing the Electrical Activity in the Heart. Springer, 2006.
6. Елькин Ю.Е., Москаленко А.В., СтармерЧ.Ф. Спонтанная остановка дрейфа спиральной волны в однородной возбудимой среде // Математическая биология и биоинформатика. 2007. 2. № 1. С. 73-81.
7. МедвединскийА.Б., Русаков А. В., Москаленко А. В. и др. Исследование автоволновых механизмов вариабельности электрокардиограмм во время высокочастотных аритмий: результаты математического моделирования // Биофизика. 2003. 48. № 2. С. 314-323.
8. Не Y., К eyes D. Е. Reconstructing parameters of the FitzHugh-Nagumo system from boundary potential measurements //J. Comput. Neurosci. 2007. 23. N 2. P. 251-264.
Поступила в редакцию 14.12.10
ВЕСТН. МОСК. УН-ТА. СЕР. 15. ВЫЧИСЛ. МАТЕМ. И КИБЕРН. 2011. № 3
47
NUMERICAL METHOD FOR DETERMINING THE LOCALIZED INITIAL CONDITIONS IN FITZ-HUGH-NAGUMO AND ALIEV-PANFILOV MODELS
Pavelchak I. A.
We consider the inverse problem for FitzHugh-Nagumo and Aliev-Panfilov models. These models describe wave propagation in excitable medium. The stated inverse problem is to determine localized initial conditions through measurements on the boundary of 2D domain. We present a numerical method for the stated problem and show results of numerical simulations for domains that are similar to some profiles of heart.
Keywords: inverse problem, FitzHugh-Nagumo model, Aliev-Panfilov model, localized initial conditions.
