автореферат диссертации по информатике, вычислительной технике и управлению, 05.13.18, диссертация на тему:Численные методы решения обратных задач для математических моделей возбуждения сердца

Московский государственный университет им. М.В. Ломоносова. Факультет вычислительной математики и кибернетики

Численные методы решения обратных задач для математических моделей возбуждения

сердца

05.13.18 - Математическое моделирование, численные методы и комплексы

программ

АВТОРЕФЕРАТ диссертации на соискание ученой степени кандидата физико-математических наук

На правах рукописи

1

Павельчак Иван Алексеевич

2 2 НОЯ 2012

005055357

Москва - 2012

005055357

Работа выполнена на кафедре математической физики факультета вычислительной математики и кибернетики Московского государственного университета имени М.В. Ломоносова.

Защита состоится 12 декабря 2012 г. в 15 ч. 30 мин. на заседании диссертационного совета Д 501.001.43 при Московском государственном университете имени М.В. Ломоносова, расположенном по адресу: 119991, ГСП-1, Москва, Ленинские горы, МГУ, 2-й учебный корпус, факультет вычислительной математики и кибернетики, ауд. 685.

С диссертацией можно ознакомиться в библиотеке факультета ВМК МГУ имени М.В. Ломоносова.

Автореферат разослан « Д- » WDJzf^ А 2012 г.

Ученый секретарь

диссертационного совета,

доктор физико-математических наук,

Научный руководитель:

Официальные оппоненты:

Ведущая организация:

доктор физико-математических наук, профессор,

Денисов Александр Михайлович доктор физико-математических паук, профессор,

Братусь Александр Сергеевич; доктор физико-математических наук, профессор,

Ягола Анатолий Григорьевич Институт проблем передачи информации им. A.A. Харкевича РАН

профессор

Захаров Е.В.

Общая характеристика работы

Актуальность работы. В настоящее время методы математического моделирования активно используются в медицине. Одной из важных сфер применения этих методов являются проблемы кардиологии. Математические методы и компьютерные технологии позволяют анализировать различные процессы сердечной активности и совершенствовать диагностику кардиологических заболеваний.

Методы математического моделирования играют большую роль в исследовании электрофизиологических процессов, происходящих в сердце, и выявлении различных нарушений сердечной деятельности. Электрофизиологиче-скис процессы в сердечной мышце характеризуются изменением во времени трансмембранного потенциала. Для описания процесса возбуждения сердца в терминах трансмембранного потенциала предложен ряд математических моделей, см. обзор в [9]. Широкое распространение получили монодоменные модели, представляющие собой начально-краевые задачи для квазилинейных эволюционных систем уравнений в частных производных, рассматриваемых в областях с достаточно сложной геометрией [10|. К числу монодоменных моделей относятся модели Фитц-Хыо-Нагумо |11-13] и Алиева-Панфилова [14], активно используемые для анализа различных процессов возбуждения сердца.

Важным направлением в применении математических методов и компьютерных технологий в кардиологии является разработка численных методов и программного обеспечения для решения различных задач диагностики заболеваний сердца. Многие методы вычислительной диагностики базируются на решении обратных задач для математических моделей возбуждения сердца. Они могут быть использованы для диагностики различных болезней, например таких широко распространенных заболеваний, как аритмия и инфаркт

миокарда. Разработка численных методов решения обратных задач для математических моделей возбуждения сердца., их программная реализация и применение в электрофизиологии сердца безусловно являются актуальными.

Цель диссертационной работы

• разработка численных методов решения обратных задач математических моделей возбуждения сердца

• программная реализация предложенных численных методов решения обратных задач

• проведение вычислительных экспериментов с целью анализа точности решения обратных задач и возможности применения предложенных методов для диагностики кардиологических заболеваний

Научная новизна. Рассмотрены новые постановки обратных задач для математических моделей возбуждения сердца. Разработаны численные методы решения поставленных обратных задач. Создан программный комплекс, реализующий предложенные численные методы. Проведены численные эксперименты, показавшие достаточно хорошую точность решения обратных задач предложенными методами.

Практическая значимость. Предложенные численные методы и созданный комплекс программ могут быть использованы для разработки методов и средств диагностики кардиологических заболеваний.

На защиту выносятся следующие основные результаты и положения:

• численные методы решения задач определения локализованного источника в математических моделях возбуждения сердца;

• численные методы решения задач определения параметров математических моделей возбуждения сердца и определения области, утратившей

4

способность к возбуждению;

• программный комплекс, реализующий разработанные численные методы; результаты вычислительных экспериментов решения обратных задач, показывающие возможность использования предложенных методов для диагностики кардиологических заболеваний.

Апробация работы. Основные результаты диссертации докладывались

• V международной конференции "Математические идеи П.Л. Чебышева и их приложение к современным проблемам естествознания"(Обнинск, 14-18 мая 2011 года)

• научной конференции "Тихоновские чтения"(Москва, МГУ им. М.В. Ломоносова, 14 июня 2011 года)

• научном семинаре лаборатории обработки биоэлектрической информации в Институте проблем передачи информации им. A.A. Харкевича РАН

• научной конференции "Ломоносовские чтения" (Москва, МГУ им. М.В. Ломоносова, 16 25 апреля 2012 года)

• научном семинаре "Обратные задачи математической физики"под руководством профессоров A.B. Бакушинского, A.B. Тихонравова и А.Г. Яголы в Научно-исследовательском вычислительном центре МГУ им. М.В. Ломоносова

• научно-исследовательском семинаре кафедры математической физики факультета ВМиК МГУ им. М.В. Ломоносова

Публикации. Материалы диссертации опубликованы в 8 печатных работах, из них 4 статьи в журналах списка ВАК |1-4|, 2 статьи [5, б| и 2 тезиса докладов [7, 8).

Структура и объем диссертации Диссертация состоит из введения, обзора литературы, 3 глав, заключения и библиографии. Общий объем диссертации 91 страница, из них 83 страницы текста, включая 17 рисунков. Библиография включает 52 наименования на 7 страницах.

Содержание работы

Во введении обоснована актуальность диссертационной работы, сформулирована цель и аргументирована научная новизна исследований, показана практическая значимость полученных результатов, представлены выносимые на защиту научные положения, сделан обзор работ по математическим моделям возбуждения сердца и обратным задачам для этих моделей.

В первой главе диссертации разрабатываются численные методы решения задач определения источника в моделях возбуждения сердца.

В §1.1 приведены рассматриваемые модели возбуждения сердца. Это две монодоменные модели, описывающие изменение трансмембранного потенциала в сердечной ткани. Наибольший интерес с практической точки зрения представляет изучение этих моделей в трехмерных областях, соответствующих геометрии сердца. Однако, учитывая теоретическую и вычислительную сложность соответствующих трехмерных задач, в большом числе работ исследуются двумерные модели, позволяющие учесть многие эффекты, характерные для процессов возбуждения сердца, см. например [15-17]. Поскольку обратные задачи имеют существенно большую численную трудоемкость, чем прямые, в диссертации рассматриваются обратные задачи для двумерных моделей возбуждения сердца.

Модель Фитц-Хыо-Нагумо [11—13] является наиболее часто употребляемой при исследовании процессов возбуждения сердца. Она представляет собой начально-краевую задачу:

щ = ОАи-и(и-а){и- 1)-к>, (х,у)еС, ¿б[0,Т], (1)

ин=?и~г», (х,у)еС,. <6[0,П, (2)

= (х.у)ег, 4 е [0,21, (3)

дп

и(х, у, 0) = #>(*. У)> (х, у) 6 С, (4)

т(т, у., 0) = 0, {х,у)еС. (5)

Здесь в - ограниченная область с границей Г, Д а, ¡3,7 - заданные положительные постоянные. Функция и(х,уЛ) - это трансмембранный потенциал; функция - медленная восстанавливающая переменная, связанная с

ионными токами. Постоянная Б представляет собой коэффициент электропроводности среды; а - коэффициент порога возбуждения среды; параметры 3 и 7 определяют свойства бегущей волны. Модель (1)-(5) позволяет качественно описывать процесс распространения возбуждения в миокарде и дает хорошую точность моделирования таких наблюдаемых характеристик, как продолжительность импульса и скорость его распространения |9, 15, 18-20]. Но некоторые другие характеристики процесса, такие как форма импульса и восстанавливающие свойства среды, модель Фитц-Хью-Нагумо описывает не точно |9. 14]. Модель Алиева-Панфилова [14] более точно описывает форму наблюдаемых импульсов в миокарде. Модель'Алиева-Панфилова записывается гак:

щ = ЯДи - *и(и - *)(« - 0 - ии(х'у) 6 Т]' (6)

щ = _ (£о + (ш + ки(и-а- 1)), {X,у) е 6\ ¿€ [0,Т], (7)

V и+ р.2 )

= 0, (х,у)ет, £ е 10,71, (8)

дп

9.9

/1

Л

N

ы

£Э

(а) Фптц-Хьет-Нагумо

(Ь) Алиева-Панфилова

Рис. 1: Профили волн для моделей

и(х,у,0) = <р(х,у), и>(х,у, 0) = О,

(х,у)€в, (9) (х,у)еа. (Ю)

Здесь (3 - ограниченная область с границей Г, О, к, а, е0, Мь № ~ заданные положительные постоянные, р(х, у) - заданная функция. Как и в модели Фнтц-Хью- Нагумо, функция и(х,у,г) - это трансмембранный потенциал, функция ги(х, у, I) - медленная восстанавливающая переменная, связанная с ионными токами. Постоянные коэффициенты «г0, Рь № уравнения (7) позволяют точно приблизить форму моделируемого импульса к экспериментальным данным [14]. Различия формы фронта распространяющегося импульса показаны на рисунке 1. В работе рассматриваются две модели - часто употребляемая Фитц-Хью-Нагумо, на основе которой уже проведено множество различных исследований электрической активности сердца, и более точная Алиева Панфилова, лучше описывающая форму фронта распространяющегося возбуждения.

Также в §1.1 предложены численные методы решения прямых задач для рассматриваемых моделей, основанные на методе конечных элементов. Прямые задачи решаются в областях двумерного пространства, приближенных к некоторым сечениям сердца. На рисунке 2 приведен пример решения нря-

Рис. 2: Решение прямой задачи для модели Фитц-Хью-Нагумо: (а) 0; (б) 4-50: (в) 1-150; (г) /.-200.

мой задачи для модели Фитц-Хью-Нагумо с помощью описанного численного метода. На нем показано формирование и распространение бегущей волны, соответствующей изменению трансмсмбранного потенциала при локализованном начальном возбуждении <р(х,у) (рис. 2 а.))

В §1.2.1 рассматривается задача определения локализованного начального возбуждения <р(т,у) для моделей Фитц-Хью-Нагумо и Алиева-Панфилова, связанная с диагностикой такого заболевания, как аритмия. Обратная задача состоит и определении функции <р(х,у) по измерениям потенциала и(х,уЛ) на внешней границе Г! области й

и(т, у, 0 = Ф{т, у, <), (-Т, у) € Г! с Г, I е [0, т].

Начальное возбуждение считается локализованным, а именно предполагается, что 1р(х. у) имеет вид

В §1.2.2 предлагается численный метод решения обратных задач определения начального возбуждения, основанный на минимизации функции невязки

Ф(х,у,а)

(и(х, у, t; х, у, а) - %]'(х, у, t.)fdldt.

о г,

Предложенный метод использует априорную информацию об искомых характеристиках и, таким образом, базируется на принципах решения некорректно поставленных задач, предложенных в фундаментальных работах А. Н. Тихонова [21, 22], М. М. Лаврентьева [23] и В. К. Иванова [24]. Для вычисления частных производных функции Ф решаются начально-краевые задачи для уравнений в частных производных. Предложен алгоритм вычисления начального приближения для градиентного метода минимизации, основанный на анализе времени прихода импульса на границу.

В §1.2.3 приводятся результаты вычислительных экспериментов с применением предложенного метода. Результаты решения обратных задач имеют хорошую точность. В табл. 1 показаны результаты решения обратной задачи для обеих моделей в области, имеющей два выреза, с функцией ip(x,y) = e-((x-23)4(y-44)2)/36 и погрешностью равной S2 = O.Olll^ll2.

Таблица 1: Результаты в области с двумя вырезами

Модель Фитц-Хью Нагумо Модель Алиева-Панфилова

X У а X У а

Точные знач. 23 44 6 23 44 6

Результат 22.6153 44.0489 5.96315 22.6689 44.0211 5.92909

В §1.3.1 рассматривается задача определения локализованного источника в уравнении для моделей Фитц-Хыо-Нагумо и Алиева-Панфилова. В этом случае начально-краевая задача для модели Фитц-Хью-Нагумо имеет вид щ = DAw - и{и - а)[и - 1) - w + 1{х, у, í; х, у, i, а. в), (х, у) € G,t е [О, Г], wt = ¡3u- 710, (i,y) € G,í € [0,Т],

—(x,s,,i) =0, е r,í е [о, т],

дп

и(х,у, 0) = 0, (z,2/)€G,

iü(a:,t/;0) = 0,

где / - функция локализованного источника вида

/(*, y, í; 5, у, i, <х, 0) = e-to-V+to-vWe-^.

Требуется определить параметры функции источника I по.измерениям потенциала и(х, у, t) на внешней границе Ti области G

u{x,y,t)=il>{x,y,t), (г,у) е Tj С Г, ÍG[0,T].

В отличии от постановки задачи, сформулированной во втором параграфе, здесь к неизвестным параметрам добавляются í и 0. Аналогично ставится задача для модели Алиева-Панфилова.

В §1.3.2 предлагается численный метод решения обратных задач, основанный на минимизации функции невязки F. Для вычисления частных производных функции F записывается сопряженная задача для уравнений в частных производных. Для модели Фитц-Хью- Нагумо сопряженная задача имеет следующий вид:

да —DAa — ,3b + n(3u2 — 2(1 + а)и + a), (x,y)eG, íe[0,T],

а + 7Ь, (x,y)eG, t € [0, T],

D^{x,y¡t) = 2(u-v), (х,у)€Г1г t е [0,Т],

дп

Sí дЬ

dt.

(х,у,<).= о, (.т!у)ег\г1! t. е [о,т],

an

а{х. у, Т) = О, (.X, у) € G,

6(х,з/,Т)=0, (i,2/) eG,

при этом частиыс производные функции F вычисляются по формуле

т

Ja(x,y,t)—dadt,

ÖF 9Aj

о с

где Аj - один из параметров х,уЛ,а,в. В §1.3.3 приводятся результаты вычислительных экспериментов с применением разработанного метода.

В §1.4.1 ставится задача определения локализованного начального условия для модифицированных моделей, в которых рассматривается изменение двух потенциалов в областях, соответствующих сердцу и торсу. Рассмотрим модель Фитц-Хью-Нагумо:

ut= Dku-u{u-a)(u-l)-w, *€[0,Г},

w, = ßu-itu, (т, у)вН, L e [0,T],

$i(i.y.i) = o, {x,i/)едн, t € [о,т],

an

u{x, y, 0) = <p(x, у), (x, у) e H,

ш(х,у, 0) = 0, (x.y)eH, к которой добавляются уравнения

Av = Au, (х,у)еН, te[0,T},

Av = о, (x,y)eG\H, t e [о,т],

^-(x,y,i) =0, (x,y)edG, i 6 [0,71,

ort

1'дс G - ограниченная область с границей dG, такая, что Н С G. Функция u{x,y,t) описывает изменение потенциала электрического поля, инициированного изменением трансмембранного потенциала u(x,y,t) [25]. Потенциал v(x, у, ¿), определенный вне области Я, может быть измерен на границе обла-

12

сти о, которая интерпретируется как торс человека. Ставится сле,лующая обратная задача. Требуется определить функцию <р(х, у; х. у) = е~1и'х) ]/" если задана

Аналогично ставится задача дгтя модели Алиева-Панфилова.

В §1.4.2 предлагается численный метод решения обратных задач, основанный на минимизации функции невязки Ф. Для вычисления частных производных функции Ф решаются начально-краевые задачи для уравнений в частных производных. В §1.4.3 приводятся результаты вычислительных экспериментов с использованием предложенного метода.

Результаты первой главы опубликованы в работах [1, 4, 6, 7|.

Во второй главе диссертации разрабатываются численные методы определения параметров в моделях возбуждения сердца.

В §2.1.1 рассматривается задача определения коэффициентов моделей возбуждения сердца. Для модели Фитц-Хью-Нагумо это О, а, ¡5,7, а для модели Алиева-Панфилова - О. к, а, е0, А«1, Задача состоит в определении некоторого набора параметров модели по измерениям функции и(х,у.1) на границе области Г:

и{х,у,1) = и{х,у,1), {х,у) е г, г е [о,т].

Некоторые другие постановки обратных задач, состоящих в определении входящих в уравнения моделей Фитц-Хью-Нагумо и Алиева-Панфилова параметров по измерениям на границе области, рассматривались в ряде работ, см. например [9, 26-28]. Основное отличие исследуемых в диссертации обратных задач состоит в том, что неизвестными могут считаться любое количество параметров модели.

В §2.1.2 предлагается численный метод решения обратной задачи, основанный на минимизации функции невязки Т. Для вычисления частных производных функции Т решаются начально-краевые задачи для уравнений в частных производных. В §2.1.3 приводятся результаты вычислительных экспериментов с применением разработанного метода. Коэффициенты моделей восстанавливаются с высокой точностью по одиночке. Неизвестные пары находятся с несколько меньшей точностью. Три и более неизвестных коэффициента находятся с заметно меньшей точностью, при этом требуе тся существенно более близкое первое приближение к искомому набору.

В §2.2.1 решается задача, моделирующая проблему определения области, утратившей способность к возбуждению, например, в результате инфаркта миокарда. Рассматривается начально-красная задача для модифицированной модели Фитц-Хью-Нагумо

щ = ОАи-х{х,у)и(и-а){и-1)-и', (х,у)ев, 1в[0,Т], (11)

ги( = ¿и — 7 и>, (х,у)ес, * 6 [0,71,(12)

= о, (1,у)бГ, г е [0,71, (13)

дп

и{х,у,0) = <р{х,у), {х,у)£С, (14)

го(х, у, 0) = 0, {х,у) е й. (15)

Функция х(х<У) <= С1 (С) принимает значения, близкие к нулю в области Н С С, и значения, близкие к единице в области С\Н. Нелинейный источник и(и — о)(и — 1) в уравнении (1) определяет способность среды к возбуждению, а в модифицированной модели нелинейный источник вида х(х,у)и(и — а)(и — 1) характеризует среду, способную к возбуждению в области й \ Я и пс способную к возбуждению в области Я. Таким образом, математическая модель (11)—(15) может быть использована для описания процессов возбуждения в сердце, часть которого (область Я) потеряла способность к возбуждению. Подобный подход для других моделей возбуждения

сердца рассматривался в [9|. Будем считать, что граница области Н задастся п параметрами Ль ..., А„. Положим функцию х(х,у; Ль ■ ■ ■ > А") равной

Х'(х, у: Ль ..., А„) = ^ + ^ arctan(02£(x, у; Аь ..., А„));

где (](х, у; Ai,..., А„) - известная функция, принимающая значения д{х.у; Ai, • • ■, А„) < 0, (х, у) е Н и д{х,у\ Аь ..., А„) > 0, (x,y)eG\H, а б - заданная постоянная. Сформулируем обратную задачу для модифицированной модели (11)—(15). Требуется найти функцию х(х, у: Ль ..., А„), определяющую границу области Я, если на внешней границе Г^ области G заданы решения задачи (11)—(15)

щ(х,у,1)=1ч(*,уЛ (*.!/) е Г,, Í 6 [0,71, i = l,...,m,

соответствующие различным начальным условиям и;(х, у, 0) = ¿¡(х,у). Коэффициенты моделей и функции <Pi(x,y), (х,у) е G,i = l,...,m, заданы. Аналогично ставится задача для модели Алиева-Панфилова.

В §2.2.2 предлагается численный метод решения обратной задачи, основанный на минимизации функции невязки Y. Для вычисления частных производных функции Y записываются сопряженные задачи для уравнений в частных производных. Для модели Фитц-Хью-Нагумо сопряженные задачи имеют следующий вид:

^ =-DA0i - ¿Ь,-+ aix(3u?-2(1+ a)u¡ +a), (x,y)eG, Í6[0,T], dt.

= + (x,y)eG, Í6[0,T],

dt

D~^{x,y,t) = 2(¡i, — (х,у)еГг, íe[0,T],

on

—í(x, y, t) = 0. (i,») €Г\Г,, í € [0, T],

dri

а,(x,y,T) = o,

б,(х,у,Г)=0, {x,y)eG.

15

Частные производные функции У вычисляются по формуле

щ{х,у,1)щ(ич - а){щ - 1)—1-<1,0(11. /=1 о с"

В §2.2.3 приводятся результаты вычислительных экспериментов с применением разработанного метода. В них искалась область Н. потерявшая способность к возбуждению, по известным решениям одной и двух прямых задач. На рис. 3 показан результаты одного такого вычислительного эксперимента для модели Алиева-Панфилова по двум решениям прямых задач. На рисунке заштрихованы внутренние вырезы области й: крестами обозначены точки локализации начальных распределений для прямых задач: простой линией обозначена точная искомая область, жирной линией - полученный предложенным методом результат. Относительная погрешность исходных данных равна 2%.

-200 -150 -100 -50 0 50 100 150 200

Рис. 3: Результат восстановление области по 2 решениям

Результаты второй главы опубликованы в работах [2, 3. 5. 8].

В третьей главе диссертации описываются программы комплекса, включающего себя 18 программ, реализующих предложенные численные методы решения прямых и обратных задач. Программы представляют собой консольные приложения, написанные па языке программирования С - . Они прини-

16

0Y ЗА,

мают на вход файл с параметрами моделирования и записывают в выходные файлы результаты расчетов. В них используется открытая библиотека deul.ii 1 для работы с разреженными матрицами. Кроме того, в комплексе используются программы gmsh 2 для разбиения областей на конечные элементы и gnuplot 3 для построения графиков. Также в третьей главе приводятся результаты некоторых тестовых расчетов.

В §3.1 приведены программы решения различных прямых задач для моделей возбуждения сердца. Это программы решения прямых задач с функцией локализованного начального распределения для моделей Фитц-Хью-Нагу-мо и Алиева-Панфилова; программы решения прямых задач с локализованной функцией источника в уравнении; программы решения прямых задач с дополнительно введенной внешней областью; программы решения прямых задач для модифицированных моделей, описывающих процессы возбуждения в сердце, часть которого потеряла способность к возбуждению.

В §3.2 описаны программы решения обратных задач, рассмотренных в главе 1. Это программы решения задач определения локализованной функции источника, программы решения задач определения локализованной функции начального распределения, программы решения задач определения локализованной функции начального распределения для моделей с дополнительно введенной внешней областью.

В §3.3 приведены программы решения обратных задач, рассмотренных в главе 2. Это программы решения задач идентификации параметров моделей Фитц-Хыо-Нагумо и Алиева-Панфилова и программы решения задач для модифицированных моделей по определению области, потерявшей способность к возбуждению.

1 A Finite Element Differential Equations Analysis Library (http: 'www.dratti.org )

2 ht.tp://www.geuz.arg/guish/

:t http://www.giiuplot.mfo/

Результаты вычислительных экспериментов на основе разработанного комплекса программ опубликованы в работах |1-8|

В заключении сформулированы основные результаты, полученные в диссертационной работе.

Основные результаты

• Для математических моделей возбуждения сердца созданы и реализованы численные методы решения задач определения локализованного источника возбуждения, определяемого либо функцией начального распределения потенциала, либо функцией источника в уравнении.

• Предложены и разработаны численные методы решения задач идентификации параметров математических моделей и определения области, потерявшей способность к возбуждению, в модифицированных математических моделях возбуждения сердца.

• Создан программный комплекс, реализующий предложенные численные методы. На его основе проведены вычислительные эксперименты, показавшие хорошую точность определения неизвестных характеристик математических моделей возбуждения сердца и возможность использования разработанных методов в медицинской диагностике.

Список публикаций

1. Павельчак И. А. Численный метод определения локализованного начального условия в моделях Фитц-Хыо-Нагумо и Алиева-Панфилова / ! Вестник МГУ. Вычислительная математика и кибернетика. 2011. № 3. С. 7-13.

2. 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. Vol. 23, no. 2. P. 208-215.

3. Павельчак И. А. Численный метод определения параметров в моделях Фитц-Хыо-Нагумо и Алиева-Панфилова // Вычислительные методы и программирование. 2012. Т. 13. С. 172-176.

4. Денисов А. М., Павельчак И. А. Численный метод определения локализованного начального возбуждения для некоторых математических моделей возбуждения сердца // Математическое моделирование. 2012. № 7. С. 59-66.

5. Павельчак И. А., Туйкина С. Р. О численном решении одной обратной задачи для модифицированной математической модели Алиева-Панфилова // Прикладная математика и информатика. 2012. Л' 40. С. 20-28.

6. Павельчак И. А. Метод численного решения задачи определения источника в модели Фитц-Хью-Нагумо / / Прикладная математика и информатика. 2012. № 40. С. 29-37.

7. Павельчак И. А. Восстановление локализованного начального условия в математических моделях процессов возбуждения сердца // V международная конференция "Математические идеи П.Л. Чебышева и их приложение к современным проблемам естествознания", Обнинск. 14-18 мая 2011 г. Тезисы докладов. С. 82.

8. Павельчак И. А., Туйкииа С. Р. О численном решении обратной задачи для модифицированной модели Фитц-Хыо-Нагумо .- ' Научная конференция Тихоновские чтения", Москва, 14 июня 2011 г. Тезисы докладов. С. 62.

Цитированная литература

9. Sundries J.. Lines G. Т.. Cai X. et al. Computing the Electrical Activity in the Heart. Berlin and Heidelberg and New York: Springer, 2006. P. 311. ISBN: 3540334327.

10. Keener J., Sneyd J. Mathematical Physiology. Springer, 1998.

11. FitzHugh R. Mathematical models of threshold phenomena in the nerve membrane // Bull. Math. Biophysics. 1955. no. 17. P. 257-278.

12. FitzHugh R. Impulses and physiological states in theoretical models of nerve membrane // Biophysical J. 1961. no. 1. P. 445-466.

13. Nagumo J., Arimoto S., Yoshizawa S. An active pulse transmission line simulating nerve axon // Proc. IRE. 1962. no. 50. P. 2061-2070.

14. Aliev R. R., Panfilov A. V. A simple two-variable model of cardiac excitation // Chaos Solutions and Fractals. 1996. Vol. 7, no. 3. P. 293-301.

15. Winfree A. T. Varieties of spiral wave behavior: An experimentalist's approach to the theory of excitable media //' Chaos. 1991. Vol. 1, no. 3. P. 303-334.

16. Davidenko J. M., Salomonsz R., Pertsov A. M. et al. Effects of pacing reentrant activity. Theoretical and experimental study // Circ. Res. 1995. no. 77. P. 1166-1179.

17. Елькин Ю. E., Москаленко А. В., Стармер Ч. Ф. Спонтанная остановка дрейфа спиральной волны в однородной возбудимой среде // Математическая биология и бионнформатика. 2007. Т. 2, № 1. С. 73-81.

18. Rinzel J., Keller J. В. Travelling wave solution of a nerve coduction equation // Biofis. J. 1973. no. 13. P. 1313-1337.

19. Pertsov A. M., Ermakova E. A.. Panfilov A. V. Rotating spiral waves in и modified FitzHugh-Nagumo model // Physica D. 198-1. no. 14. P. 117-124.

20. Courtemanche M., Skaggs W., Winfree A. T. Stable three-dimensional action potential circulation in the FitzHugh-Nagumo. model // Physica D. 1990. no. 41. P. 173-183.

21. Тихонов A. H. Об устойчивости обратных задач // Докл. АН СССР. 1943. Т. 39, № 5. С. 195-198.

22. Тихонов А. Н., Арсенин В. Я. Методы решения некорректных задач. Москва: Наука, 1974. С. 288.

23. Лаврентьев М. М., Романов В. Г., Шишатский С. Т. Некорректные задачи математической физики и анализа, Москва: Наука, 1980. С. 286.

24. Иванов В. К., Васин В. В., Танана В. П. Теория линейных некорректных задач и ее приложения. Москва: Наука, 1978. С. 206.

25. Титомир JL, Кнеппо Л. Математическое моделирование биоэлектрического генератора сердца. Москва: Наука, 1999. С. 447.

26. Не Y., Keyes D. Е. Reconstructing parameters of the FitzHugh-Nagumo system from boundary potential measurements //J. Comput. Neurosci. 2007. Vol. 23, no. 2. P. 251-264.

27. Moreau-Villeger V., Delingette H., Sermesant M. ., et al. Building maps of local apparent conductivity of the epicardium with a 2 D electrophysiological model of the heart // IEEE Transactions on Biomedical Engineering. 2006. no. 53. P. 1457-1466.

28. Cox S. J., Wagner A. Lateral overdetermination of the FitzHugh-Nagumo system / - Inverse Problems. 2004. no. 20. P. 1639-1647.

Напечатано с готового оригинал-макета

Подписано в печать 29.10.2012 г. Формат 60x90 1/16. Усл.печ.л. 1,0. Тираж 100 экз. Заказ 417

Издательство ООО "МАКС Пресс" Лицензия ИД N 00510 от 01.12.99 г. 119992, ГСП-2, Москва, Ленинские горы, МГУ им. М.В. Ломоносова, 2-й учебный корпус, 527 к. Тел. 8-495-939-3890. Тел./факс 8-495-939-3891.

Московский государственный университет им. М.В. Ломоносова. Факультет вычислительной математики и кибернетики

На правах рукописи

04201350383

Павельчак Иван Алексеевич

Численные методы решения обратных задач для математических моделей возбуждения

сердца

05.13.18 - Математическое моделирование, численные методы и комплексы

программ

ДИССЕРТАЦИЯ на соискание ученой степени кандидата физико-математических наук

Научный руководитель

д. ф.-м. н.; проф.

Денисов Александр Михайлович

Москва - 2012

Содержание

Введение 3

Обзор литературы 7

Глава 1. Задачи определения источника возбуждения 11

1.1. Математические модели возбуждения сердца. Численные методы их решения............................11

1.2. Обратная задача определения локализованного начального условия ................................................................15

1.3. Обратная задача определения функции источника в уравнении 2-5

1.4. Обратная задача определения начального возбуждения, рассматриваемая в двух областях...................37

Глава 2. Задачи определения параметров моделей 46

2.1. Обратная задача поиска параметров моделей...........46

2.2. Обратная задача поиска области, потерявшей способность к возбуждению ...............................54

Глава 3. Программный комплекс 68

3.1. Программы решения прямых задач для моделей возбуждения сердца.................................68

3.2. Программы решения задач определения локализованного источника возбуждения ........................74

3.3. Программы решения обратных задач определения параметров моделей................................79

Заключение 84

Литература 85

Введение

Актуальность работы. В настоящее время методы математического моделирования активно используются в медицине. Одной из важных сфер применения этих методов являются проблемы кардиологии. Математические методы и компьютерные технологии позволяют анализировать различные процессы сердечной активности и совершенствовать диагностику кардиологических заболеваний.

Методы математического моделирования играют большую роль в исследовании электрофизиологических процессов, происходящих в сердце, и выявлении различных нарушений сердечной деятельности. Электрофизиологические процессы в сердечной мышце характеризуются изменением во времени трансмембранного потенциала. Для описания процесса возбуждения сердца в терминах трансмембранного потенциала предложен ряд математических моделей, см. обзор в [1]. Широкое распространение получили монодоменные модели, представляющие собой начально-краевые задачи для квазилинейных эволюционных систем уравнений в частных производных, рассматриваемых в областях с достаточно сложной геометрией [2]. К числу монодоменных моделей относятся модели Фитц-Хью-Нагумо [3-5] и Алиева-Панфилова [6], активно используемые для анализа различных процессов возбуждения сердца.

Важным направлением в применении математических методов и компьютерных технологий в кардиологии является разработка численных методов и программного обеспечения для решения различных задач диагностики заболеваний сердца. Многие методы вычислительной диагностики базируются на решении обратных задач для математических моделей возбуждения сердца. Они могут быть использованы для диагностики различных болезней, например таких широко распространенных заболеваний, как аритмия и инфаркт

миокарда. Разработка численных методов решения обратных задач для математических моделей возбуждения сердца, их программная реализация и применение в электрофизиологии сердца безусловно являются актуальными.

Цель диссертационной работы

• разработка численных методов решения обратных задач математических моделей возбуждения сердца

• программная реализация предложенных численных методов решения обратных задач

• проведение вычислительных экспериментов с целью анализа точности решения обратных задач и возможности применения предложенных методов для диагностики кардиологических заболеваний

Научная новизна. Рассмотрены новые постановки обратных задач для математических моделей возбуждения сердца. Разработаны численные методы решения поставленных обратных задач. Создан программный комплекс, реализующий предложенные численные методы. Проведены численные эксперименты, показавшие достаточно хорошую точность решения обратных задач предложенными методами.

Практическая значимость. Предложенные численные методы и созданный комплекс программ могут быть использованы для разработки методов и средств диагностики кардиологических заболеваний.

На защиту выносятся следующие основные результаты и положения:

• численные методы решения задач определения локализованного источника в математических моделях возбуждения сердца.

• численные методы решения задач определения параметров математических моделей возбуждения сердца и определения области, утратившей

способность к возбуждению;

• программный комплекс, реализующий разработанные численные методы; результаты вычислительных экспериментов решения обратных задач, показывающие возможность использования предложенных методов для диагностики кардиологических заболеваний.

Апробация работы. Основные результаты диссертации докладывались на:

• V международной конференции "Математические идеи П.Л. Чебышева и их приложение к современным проблемам естествознания"(Обнинск, 14-18 мая 2011 года)

• научной конференции "Тихоновские чтения"(Москва, МГУ им. М.В. Ломоносова, 14 июня 2011 года)

• научном семинаре лаборатории обработки биоэлектрической информации в Институте проблем передачи информации им. A.A. Харкевича РАН

• научной конференции "Ломоносовские чтения"(Москва, МГУ им. М.В. Ломоносова, 16-25 апреля 2012 года)

• научном семинаре "Обратные задачи математической физики "под руководством профессоров А.Б. Бакушинского, A.B. Тихонравова и А.Г. Яголы в Научно-исследовательском вычислительном центре МГУ им. М.В. Ломоносова

• научно-исследовательском семинаре кафедры математической физики факультета ВМиК МГУ им. М.В. Ломоносова

Публикации. Материалы диссертации опубликованы в 8 печатных работах. из них 4 статьи в журналах списка ВАК [7-10], 2 статьи [11, 12] и 2 тезиса докладов [13, 14].

Структура и объем диссертации Диссертация состоит из введения, обзора литературы, 3 глав, заключения и библиографии. Общий объем диссертации 91 страница, из них 83 страницы текста, включая 17 рисунков. Библиография включает 52 наименования на 7 страницах.

Обзор литературы

Методы математического моделирования электрических процессов, происходящих в тканях живых организмов, начали развиваться во второй половине XX века. В 1952 году Аланом Ллойдом Ходжкином и Эндрю Хаксли была предложена модель, описывающая процессы генерации и передачи электрических сигналов в гигантских аксонах кальмара. С тех пор было предложено много моделей, описывающих процессы электрического возбуждения как на клеточном уровне, так и на уровне ткани. Эти модели можно разделить на несколько категорий.

Первую категорию представляют модели ионных токов. К ним относятся модель Ходжкина-Хаксли [15-19] и некоторые её модификации [20]. В них описываются процессы течения ионных токов (Са+, К+) через клеточные мембраны. Они включаются в себя большое количество переменных и параметров. Модели ионных токов содержат большое количество уравнений и имеют очень высокую вычислительную сложность.

Во вторую категорию входят эйкональные модели [21]. Они очень просты и описывают только время, за которое деполяризационная волна распространяется от одной точки пространства до другой. Эти модели не дают достаточно точного описание процессов реакций, происходящих при электрическом возбуждении тканей.

Третью категорию образуют феноменологические модели, описывающие возникновение потенциала действия и его распространение вдоль клеточных мембран. Эти модели делятся на монодоменные, описывающие трансмембранный потенциал, и бидоменные, моделирующие внутри- и внешнеклеточные потенциалы. В моделях ионных токов внутри- и внешнеклеточные потенциалы считаются определенными на непересекающихся областях, разделенных клеточными мембранами. В феноменологических моделях эти потенциалы

считаются определенными на всей рассматриваемой возбудимой области как среднее значение по малому объему. Трансмембранный потенциал представляет собой разницу усредненных внутри- и внешнеклсточных потенциалов. Феноменологические модели возбуждения сердца лучше других подходят для исследования обратных задач, так как они с вычислительной точки зрения не столь сложны как модели ионных токов и, с другой стороны, достаточно точно моделируют процесс изменения потенциалов на уровне ткани.

В диссертационной работе исследуются обратные задачи для двух феноменологических монодоменных моделей возбуждения сердца. Первая рассматриваемая модель - Фитц-Хью-Нагумо [3-5]. На ее основе было проведено множество исследований процессов распространения возбуждения в миокарде, см. например [22-24]. Модель Фитц-Хыо-Нагумо дает хорошую точность моделирования таких наблюдаемых параметров, как продолжительность импульса и скорость его распространения [1, 25-28]. Но некоторые другие характеристики процесса, такие как форма импульса и восстанавливающие свойства среды, модель Фитц-Хью-Нагумо описывает не достаточно точно [1, 6]. Модель Алиева-Панфилова [6] более точно описывает форму наблюдаемых импульсов в миокарде. Она также использовалась в различных исследованиях процессов распространения возбуждения в миокарде [29-33]. Вопрос существования и единственности решения прямых задач для моделей Фитц-Хью-Нагумо и Алиева-Панфилова был рассмотрен в ряде работ, см. например [34-36].

Наибольший интерес с практической точки зрения представляет изучение этих моделей в трехмерных областях, соответствующих геометрии сердца. Однако, учитывая теоретическую и вычислительную сложность соответствующих трехмерных задач, в большом числе работ исследуются двумерные модели, позволяющие учесть многие эффекты, характерные для процессов возбуждения сердца, см. например [22. 28, 37]. Поскольку обратные задачи

имеют существенно большую численную трудоемкость, чем прямые, в диссертации рассматриваются обратные задачи для двумерных моделей возбуждения сердца.

Важным направлением исследований математических моделей возбуждения сердца является разработка методов решения обратных задач, ориентированных на совершенствование методов диагностики. Обратные задачи для математических моделей возбуждения сердца рассматривались в работах [1, 38-45] На настоящее время существует ряд работ, в которых рассматривались обратные задачи для моделей Фитц-Хью-Нагумо и Алиева-Панфилова. В них исследуется возможность восстановления входящих в уравнения функций или параметров по известным значениям потенциала либо на всей границе, либо в некотором количестве точек на ней.

Важным классом обратных задач для моделей возбуждения сердца являются задачи определения источника возбуждения. В работах [38-40] задача определения источника рассматривалась для стационарных моделей, в которых потенциал поля удовлетворяет уравнению Пуассона. В диссертации рассматриваются методы решениях задач определения источника возбуждения для эволюционных моделей возбуждения сердца, таких как Фитц-Хью-Нагу-мо и Алиева-Панфилова.

Ряд работ посвящен поиску неизвестных параметров моделей [41-43]. Параметры моделей могут варьироваться от пациента к пациенту, и, таким образом, задачи определения параметров имеют большое практическое значение для компьютерной диагностики в медицине. В большинстве работ рассматривались задачи определения одного неизвестного параметра. В диссертации рассматриваются постановки задач, в которых определяется как один, так и несколько параметров.

Большое значение для медицинской диагностики также имеют задачи определения участка сердца, пораженного инфарктом. Постановки таких за-

дач рассматривались в работах [1, 44, 45] для бидоменных моделей. В диссертации рассматривается постановка задачи определения области, пораженной инфарктом, для монодоменных моделей Фитц-Хью-Нагумо и Алиева-Панфилова.

Численные методы решения обратных задач, рассматриваемых в диссертации, используют априорную информацию об искомых характеристиках и, таким образом, базируются на принципах решения некорректно поставленных задач, предложенных в фундаментальных работах А. Н. Тихонова [46, 47], М. М. Лаврентьева [48] и В. К. Иванова [49].

Глава 1

Задачи определения источника возбуждения

1.1. Математические модели возбуждения сердца. Численные методы их решения

Одной из наиболее распространенных математических моделей возбуждения сердца является модель Фитц-Хью-Нагумо [3-5], представляющая собой начально-краевую задачу для системы эволюционных уравнений в частных производных. Запишем её в случае двух пространственных переменных:

Здесь С - ограниченная область с границей Г, Б.а./3,7 - заданные положительные постоянные, (р(х,у) -заданная функция. Функция и(х,уЛ) представляет собой трансмембранный потенциал: функция ги(х.уА) - медленную восстанавливающую переменную, связанную с ионными токами [3-5]. Постоянная И представляет собой коэффициент электропроводности среды: а -коэффициент порога возбуждения среды; /3 и 7 определяют свойства бегущей волны. Модель (1.1)—(1.5) позволяет качественно описывать процесс распространения возбуждения в миокарде и дает хорошую точность моделирования таких наблюдаемых характеристик процесса, как продолжительность импульса и скорость его распространения [1, 25-28]. Вместе с тем некоторые другие характеристики процесса, например такие, как форма импульса

щ — ОАи — и{и — а) (и — 1) — 'ш,

Wt = ¡Зи — 7 ги, Вн.

{х./у)ес, (1.1)

(.х,у)ес, *е[о,Т], (1.2) и) е- г. 1:е'().Г: (1.3)

и{х,у, 0) = (р{х,у), 'ш{х., у, 0) = 0,

(х,у)ес.

(1.4)

(1.5)

(а) Фитц-Хью-Нагумо (b) Алиева-Панфилова

Рис. 1.1: Профили волн для моделей

и восстанавливающие свойства среды, модель Фитц-Хью-Нагумо описывает не точно [1, 6]. Модель Алиева-Панфилова [6] более точно описывает форму наблюдаемых импульсов в миокарде. На рисунке 1.1 показано различие профилей волн для моделей Фитц-Хью-Нагумо и Алиева-Панфилова. Модель Алиева-Панфилова представляет собой начально-краевую задачу для системы эволюционных уравнений в частных производных

щ = БДи- ки(и- а){и- 1) -ит, {х,у)еС, * е [О,Т], (1.6)

- ) {w + ки{и- а- 1)), (x,y)€G, íg[0,T], (1.7)

/ ¡llW

Щ = ~ £0 + -;—

V и + №,

ди

—(x,y,t) = о, (х,у)ет, t е [о,т], (1.8)

- ■,:(.!:, у), (.х,у) е G, (1.9)

w(x, у, 0) = 0, (x,y)eG. (1.10)

Здесь G - ограниченная область с границей Г, D, к, а, со, Ц\, fio ~ заданные положительные постоянные, (р{х,у) - заданная функция. Как и в модели Фитц-Хью-Нагумо, функция и(х, у, t) представляет собой трансмембранный потенциал, функция w(x, у, t) - медленную восстанавливающую переменную, связанную с ионными токами.

Рассмотрим численный метод решения начально-краевых задач, опре-

деляющих модели Фитц-Хью-Нагумо и Алиева-Панфилова. Проведем дискретизацию но времени: 1п = ¿п_1 + т, to = 0. Обозначим через ип(х.у) = и(х,у^п), гип(х,у) = и>(ж, т/, ¿п). Рассмотрим следующую задачу:

т?<_1

—-— = D{9Avn + (1 - 6)Аип~1) - fn~\ (1.11)

г

= д»-\ (1.12)

Т

dvn

= 0, (1.13)

¿777, г

= (/?, (1.14)

u>° = 0, (1.15)

где функции fn{x,y) и дп(х.у) для модели Фитц-Хью-Нагумо определяются следующим образом:

р = ип{ип - о/)(ип - 1) + wn, (1.16)

дп — f3un — 7'u/1, (1.17)

а для модели Алиева-Панфилова так:

fn = кип{и'1 - а){ип - 1) + unwv, (1.18)

,п

9п=-(ео + -^j— ) (wn + кип{ип - а - 1)). (1.19)

\ ll,n + ¡lo J

Преобразуем уравнения (1.11), (1-12):

(1 - r6DA)un = (1 + r(l - 0)DA)un~1 - Tfn~\ (1.20)

wn = wn~l+ тдп~1. (1.21)

Функции wn вычисляются по формулам (1.15). (1.21). Для решения краевой задачи (1.13). (1.20) применим метод конечных элементов [50]. Домно-жим (1.20) на базисную функцию ф(х.у) и проинтегрируем по области G:

ипф - 9rDAu1^)às = {ип~1ф + т( 1 - в)ОАип'}ф - rfn~^)ds. (1.22)

G G

Заметим, что Аипф6.з = — _[ V'urгVфds + _[ что с учетом нулевого

с с ос п

граничного условия (1.13) дает

G G

('ип~гф - т( 1 - 6)ОЧуп-{Чф - r/n"V)ds. (1.23)

Пусть область пространства С разбита на четырехузловые конечные элементы с узлами (х3,у3). Вычисляемую функцию ип(х,у), определенную на (.г, у) 6 С, аппроксимируем линейной комбинацией базисных функций: ип(х.у) ~ ^и:гфг{х,у), }'г'{х,у) ~ ^ Р"ф7(х, у). Тогда 1/п является реше-

г г

нием следующей системы линейных алгебраических уравнений:

[М + тОБА)ип = (М - т{ 1 - в)ЭА)ип-1 - гМГ"1, (1.24)

где

Mtj = 0,)

фгф3 ds, (1.25)

AtJ = (V&. \7ф3) =

ШхШх + {Фг)У{Ф3)У)й8. (126)

G

На рисунке 1.2 приведен пример решения прямой задачи для модели Фитц-Хыо-Нагумо с помощью описанного численного метода. На нем показано формирование и распространение бегущей волны, соответствующей изменению трансмембранного потенциала при локализованном начальном возбуждении ip(x,ij) (рис. 1.2 а)).

Для программной реализации использовались внешние свободные компоненты: программа gm,sh 1 для разбиения произвольной области на четырехузловые конечные элементы и библиотека deal и 2 для работы с разреженными матрицами. Более подробно программная реализация описана в главе 3.

1 littp //www gcuz.oig/gmsli/

2 A Finite