CC BY
8
УДК 519.632
С. И. Соловьев^, С. Р. Туйкина2
МЕТОД ЧИСЛЕННОГО РЕШЕНИЯ ОБРАТНОЙ ЗАДАЧИ ДЛЯ МОДЕЛИ ЭЛЕКТРИЧЕСКОЙ АКТИВНОСТИ СЕРДЦА*
Для двумерной модифицированной математической модели Фитц-Хью-Нагумо рассматривается обратная задача, состоящая в определении зависящего от пространственных переменных коэффициента системы уравнений в частных производных. Дополнительные динамические измерения потенциала проводятся на всей внутренней границе области. В работе предлагается численный метод решения этой обратной задачи и приводятся результаты вычислительных экспериментов.
Ключевые слова: математическая модель Фитц-Хью-Нагумо, обратная задача, численный метод.
1. Введение. При разработке и совершенствовании диагностики кардиологических заболеваний широкое распространение получили методы математического моделирования электрофизиологии сердца. Математические модели Фитц-Хью 11агумо. Алиева-Панфилова, бидоменная модель описывают процесс возбуждения сердца в терминах трансмембранного потенциала. Они представляют собой начально-краевые задачи для системы квазилинейных уравнений в частных производных [1-4].
При разработке неинвазивных способов диагностики заболеваний сердца возникает необходимость решать обратные задачи определения параметров этих математических моделей [5-13]. Численные методы решения некоторых обратных задач для моделей возбуждения сердца были предложены в работах [1, 5, 8-10, 12, 13].
В данной работе для модифицированной математической модели Фитц-Хью-Нагумо рассматривается задача определения области сердца, пораженной инфарктом миокарда. Эта обратная задача состоит в определении зависящего от пространственных переменных коэффициента системы уравнений в частных производных по дополнительной информации, соответствующей измерениям потенциала, проводимым катетерами на границе желудочков. Задача решается в двумерной области, представляющей собой сечение сердца и его желудочков горизонтальной плоскостью.
В работе предлагается численный метод решения этой обратной задачи и приводятся результаты вычислительных экспериментов, позволяющие оценить эффективность предложенного метода. Численные методы решения обратных задач, близких по постановке к рассматриваемой задаче, предложены в [1, 7, 8, 10, 13]. Основное отличие данной работы состоит в том, что математическая двумерная область соответствует реальной геометрии сердца и его желудочков, динамические измерения потенциала проводятся на всей внутренней, а не внешней границе области сердца.
2. Постановка обратной задачи и численный метод ее решения. Модифицированная модель Фитц-Хью-Нагумо имеет следующий вид:
ди
— = ИАи — д(х, у)и(и — а)(и — 1) — го, (ж, у)б<2) ^ €Е (0, СГ*], ди]
— = ри-'уш, (ж,у)б<2, ге(о,т],
ди
—(ж,у,г) = о, (ж,у)ег, ге(о,т],
и(х,у,0) =р(х,у), Ц](х,у, 0) = 0, (ж,у)€<э,
1 Факультет ВМК МГУ, доц., к.ф.-м.н., e-mail: solQcs.msu.su
2 Факультет ВМК МГУ, доц., к.ф.-м.н., e-mail: tuikQcs.msu.su
* Работа выполнена при финансовой поддержке РФФИ, проект № 14-01-00244.
где функция и(х, у, 1) это трансмембранный потенциал, функция у, £) связана с ионными токами, функция р(х, у) описывает локализованное начальное возбуждение миокарда. Г граница области Функция д(х, у) € С1^) описывает область сердца /, пораженную инфарктом. Эта функция такова, что д(х, у) и 0 в области /€(,), и у) и 1 в Тогда нелинейный коэффи-
циент модели, описывающий способность миокарда к возбуждению, ^ = д(х,у)и(и — а) (и — 1) и О в области /, что соответствует потери в I способности к возбуждению среды, а, /3, 7 реактивные коэффициенты, I) коэффициент электропроводности (I), а, /3, 7 положительные постоянные).
Будем предполагать, что область I задается п параметрами А1,...,АН. Обозначим А = = (А1,..., А„). Функция д зависит от этих параметров: д = д(х, у; А). Рассмотрим следующую обратную задачу. Пусть известны коэффициенты I), а, /3, 7. Требуется определить непрерывную функцию д(х, у; А), если задана дополнительная информация о решении задачи (1), соответствующая р(х, у), а именно, задана функция
(р(х, у, г) = и(х, у, (ж, у) € Г1 и Г2,
где и(х,уЛ) решение задачи (1), соответствующее функции р(х,у), Г1, Г2 внутренние границы области <3 (рис. 1).
Рис. 1
Пусть для функции д = д(х, у; А) и локализованного начального возбуждения миокарда р(х, у) прямая задача (1) имеет решение на границе Г1 и Г2, равное 1р(х,уЛ). Пусть дополнительная информация 1р(х,уЛ) нам известна с погрешностью е, т.е. задана функция у, £), такая, что
т
о Г^иГо
Будем минимизировать невязку
т
I I (<Ре(Х:У:*) ~ <р(х,уЛ))2 (11 (И < £2.
<5(А) = I I («(.т. у, А) — (ре(х, у. (11 (1,1,
О ПиГо
используя градиентный итерационный метод с критерием <5(А) < е2 для окончания процесса минимизации. Для этого найдем выражение для градиента функции <5(А).
Пусть функции д(х, у; А) соответствует решение задачи (1) {и(х, у. 1; А). и;(х, у, Ц А)}, а д(х, у; А + <5А) решение {и(х, у, 1; А + <5А), и;(х, у, А + ¿А)}. Введем обозначения:
6и(х, у. Ц А. 8Х) = и(х, у, 1; А + 8А) — и(х, у. 1; А). 8гп(х, у, 1; А, 8А) = у, А + ¿А) — у, А), /(«) = и(и — а) (и — 1), <7^ = дя/д^'г
Так как имеет место равенство
/(и(х, у, А + 6\))д(х, у; X + 5\) -/(и(х, у, Ц \))д(х, у; А) =
1Ы^2д\л(х,у;Х)6Хз + /Ци)6и д(х,у;Х) + Я,
3 = 1
где К = О (5и2 + 5А2), то функции би, бы являются решениями задачи дби
т
дб'ш
~дГ
дби дп
= Бби - би; - ¡(и)^2дх^х,у;Х)6Ху - /¡1(и)6ид(х,у;Х) - К, 3 = 1
= /Зби — 7бт, (ж, у) € ¿6 (О, Т], (ж,у,г) = о, (ж,у)еГ, ¿е(0,Т],
(2)
би(х, у, 0) = 0, бт(х, у, 0) = 0, (х,у)€.(^. Введем а(х,у,1), Ь(х,у,1) как решение сопряженной начально-краевой задачи да
т
дЬ
= —БАа — РЬ + а/и(и)д(х, у; Л), (х,у) €. ¿6 [0,Т),
= а + 'уЬ, (®,у)€(Э, ¿€[0,Т),
т
да
0—(х,у,г) = 2(и-<р), (х,у)&т1ит2, ¿е [О,т], до
1>—(ж,у,г) = 0, (ж, у) б Г\(Г! иГ2), ¿6 [О,Г], дп
(3)
а(ж, у, Т) = 0, Ь(х, у, Т) = 0, (ж, у) € С?. Учитывая, что функции и являются решениями (2), получим
г
3 =
о Я
дби
~Ж
{ д5и)
Б А би + б'ш + ¡'и(и)бид( ж, у; Л) ) + Ь I -
{ Ы
/Зби + ^уб'ш
да \ (дЬ
би ( — + Б А а Н~ /5Ь — /^(и)ад( ж, у; Л) ) + би; ( — — а — 'уЬ
йх <1у (Й =
г
д
— (аби + Ьби;) — (БаАби — БбиАа)
Ы
о Я
Используя формулу Грина, начальные и граничные условия для би, бы, а и Ь, найдем
г г
йх <1у (Й.
t=Q
йх <1у
дби да \
Ба—--Оби— ) Ш сЙ =
дп дп)
О Г 0 Г^иГа
,1 = JJ (аби + Ьб'ш) Я
С другой стороны, это выражение равно
г
3 = — J II (х,у;Х)бХз + К)ёхёу<]Л
о Я
2би(и — (р) Ш (Й.
3 = 1
Так как приращение 6Б = Б(Х + <5А) — <5(Л) функции <5(Л) есть
т
5S =
(и(х, у, t; А + <5А) - <р€)2 - (и(х, у, t; А) - <р€)2
di dt =
О TiUro
т
[2(и — <ре)6и + <S?í2] di dt.
то с учетом полученного выражения для J будем иметь
т
О Г^Го
т
5S =
а ^ /(»,) ^^ д\. (х, у; \)5\¡ + Rj dx d/y dt + j j 5u2 di dt.
о д ■>-1 о г
Пренебрегая величинами второго порядка малости, для градиента получим выражение
т
а1(11)яА. (х~ у; А) ¿х ¿у сИ, 1 ^ ^ п.
0S
DX¡
о Q
С помощью вычисленного таким образом градиента производится переход от Ага к Аш+1. Итерационный процесс останавливается при выполнении неравенства ¿>(Аш) ^ е2.
3. Результаты вычислительных экспериментов. Описанный численный метод решения обратной задачи был применен для определения областей I специального вида. Для известной функции г(х, у; А): г(х, у; А) < 0. (х, у) € I, и г(х, у; А) > 0. (х, у) € Q\I, положим
1 1
д(х, у; А) = - + - arctg(02r(x, у, А)).
¿j 7Т
В проведенных вычислительных экспериментах использовалась функция
((х - Ai) eos А5 - (у - A2)sinA5^2 ((х - Ai)sinA5 + (у - Х2)сояХ5\2 г{х,у,л) — - + - — 1. (4J
А,
А4
соответствующая области I для случая п = 5. Для моделирования начального возбуждения сердца использовалась функция
р(х, у) = ехр
а
где (ж*.у*) координаты точки возбуждения, а а некоторый задаваемый параметр.
Прямые задачи для модифицированной модели Фитц-Хью Нагумо (1) решались в области ) (см. рис. 2), представляющей собой сечение сердца и его желудочков горизонтальной плоскостью, с помощью метода конечных элементов. На рис. 2 показаны: зоны возбуждения сердца (зона 1); границы, на которых производятся измерения (зона 2)\ область инфаркта (зона 3).
Рис. 2
В результате решения прямой задачи вычислялась Тр(х, у, £) на внутренней границе (х, у) € € Г1 и Г2, I € [0,Т], вносилась погрешность эксперимента е и находилась функция (ре(х,уЛ), такая, что
т
j j у, t; А) — (ре(х, у. t))2 dl dt ^ е2.
О ПиГо
В ходе вычислительных экспериментов решались обратные задачи по восстановлению функции г(х, у. А).
Начальное приближение А° = (А?, А®, А®. А®. Ад) находилось следующим образом. Фиксировались значения А® = А® = 1, Ад = 0. Область С} при помощи триангуляции Делоне разбивалась на конечное число треугольников, из которых произвольным образом выбиралось т треугольников. Координаты центра тяжести Ат-го треугольника обозначались А 5 к и к- Для Хк = (А5 к, Х?2 к, А®. А^.Ад) вычислялись значения невязки
т
Sk
S( Хк)= j j (ii(x,y,t; Х^.) — <pe(x,y,t))2 dl dt, 1 ^ к ^ m.
0 riur2
Затем находилось А:*, соответствующее минимальному значению Sk (<Smm = тт^), и полагалось
к
Х° = Хк . Далее с использованием описанного выше метода решалась обратная задача с функцией (pf и этим значением А0.
Расчеты проводились для модели с параметрами: D = 1, а = 0.15. /3 = 0.005. 7 = 0.025 и а = 4. Результаты вычислительного эксперимента приведены на рис. 3, 4.
Рис. 3
Рис. 4
Выбор начального приближения для А изображен на рис. 3. Отмечены значения Бк(А) для наборов А^ к, к, А® = = 1, Ад = 0, 1 ^ А: ^ т = 9. Результат точка с координатами А^ к = 11.38 и А^ = 123.54, в которой <5тт(А^) — 10.3.
На рис. 4 показан результат восстановления функции (4). Пунктирной чертой показана тестовая область инфаркта для набора А1 = 10. А2 = 125. Аз = 6, А4 = 4, Ад = 7г/3, а сплошной результат для А1 = 9.7. А2 = 123.61. Аз = 5.92, А4 = 3.25. Ад = 0.71. Значение невязки ¿>(А) в данном случае равно 0.004.
Проведенные вычислительные эксперименты показали, что в рамках рассматриваемой постановки обратной задачи, положение и форма области, пораженной инфарктом, в случае дополнительных измерений на всей внутренней границе может быть восстановлена по данным одного эксперимента.
СПИСОК ЛИТЕРАТУРЫ
1. Sundries J., Lines G. Т., Cai X., et al. Computing the Electrical Activity in the Heart. Berlin: Springer, 2006.
2. FitzHugh R. Mathematical models of threshold phenomena in the nerve membrane // Bull. Math. Biophysics. 1955. N 17. P. 257-278.
3. FitzHugh R. Impulses and physiological states in theoretical models of nerve membrane // Biophysical J. 1961. N 1. P. 445-466.
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.
5. He Y., К eyes D.E. Reconstructing parameters of the FitzHugh-Nagumo system from boundary potential measurements //J. Comput. Neuroscience. 2007. 23. N 2. P. 251-264.
6. Denisov A.M., Zakharov E. V., Kalinin A.V., Kalinin V. V. Numerical method for solving an inverse electrocardiography problem for a quasi stationary case //J. Inverse and 111 Posed Problems. 2012. 20. N 4. P. 501-512.
7. Pavel'chak I. A., Tuikina S.R. Numerical solution method for the inverse problem of the modified FitzHugh-Nagumo model // Computational Mathematics and Modeling. 2012. 23. N 2. P. 208-215.
8. Pavel'chak I. A., Tuikina S.R. Numerical solution of an inverse problem for the modified Aliev-Panfilov model // Computational Mathematics and Modeling. 2013. 24. N 1. P. 14-21.
9. Denisov A.M., Zakharov E. V., Kalinin A. V. Numerical solution of the localized inverse problem of electrocardiography // Computational Mathematics and Modeling. 2015. 26. N 2. P. 168-174.
10. Solov'eva S. I., Tuikina S.R. Numerical solution of the inverse problem for the mathematical model of cardiac excitation // Computational Mathematics and Modeling. 2016. 27. N 2. P. 162-171.
Н.Денисов A.M., Калинин В. В. Обратная задача для математических моделей возбуждения сердца // ЖВМ и МФ. 2010. 50. № 3. С. 539-543.
12. Денисов A.M., Захаров Е.В., Калинин А.В. Метод определения проекции точечного очага аритмии на поверхность сердца на основе решения обратной задачи электрокардиографии / / Математическое моделирование. 2012. № 4. С. 22-30.
13. Денисов A.M., Павельчак И. А. Численный метод определения локализованного начального возбуждения для некоторых математических моделей возбуждения сердца // Математическое моделирование. 2012. № 7. С. 59-66.
Поступила в редакцию 14.09.16
