Журнал вычислительной математики и математической физики, 2019, T. 59, № 12, стр. 2102-2110

1. ВВЕДЕНИЕ

В настоящей работе представлен метод численного моделирования ЭКГ с использованием персонифицированной модели туловища пациента. Прямое моделирование ЭКГ в настоящее время используется в различных областях: в фармакологических исследованиях [14], диагностике и планировании лечения сердечно-сосудистых заболеваний [9].

Моделирование электрофизиологии может быть формализовано как полномасштабное изучение электрической активности сердца, от внутриклеточного масштаба до масштаба тканей сердца [7]. Для получения электрокардиограммы необходимо рассчитать потенциалы на поверхности тела, которые в свою очередь определяются в основном распределением электрических потенциалов на поверхности миокарда и анатомическими особенностями пациента, а именно структурой, положением и проводимостями внутренних органов. Таким образом, для моделирования ЭКГ необходимо использовать персонифицированные анатомические модели органов. В данной работе мы предлагаем некоторую технологическую цепочку для моделирования ЭКГ (см. фиг. 1). При персонализированном моделировании ЭКГ рассматриваются как ткани сердца, так и окружающие органы грудной клетки и брюшной полости. Кроме того, предполагается наличие корректной анатомической модели сердца, являющейся результатом ручной сегментации, а также распределения потенциалов на поверхности сердца, полученного в результате решения модономенной или бидоменной задач [12]. Далее в работе эти этапы детально рассматриваться не будут.

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

В разд. 2 представлен метод, являющийся универсальным для автоматической сегментации одновременно нескольких органов. В разд. 3 описан подход к численному моделированию ЭКГ и продемонстрированы результаты численных экспериментов. В разд. 4 приведена и исследована слабая постановка прямой задачи ЭКГ. В разд. 5 описаны результаты моделирования и анализа чувствительности прямой модели ЭКГ к наличию в модели определенных анатомических структур. В разд. 6 предложен ускоренный метод расчета отведений ЭКГ. В разд. 7 сформулированы выводы.

2. ПОСТРОЕНИЕ ПЕРСОНАЛИЗИРОВАННОЙ МОДЕЛИ ТУЛОВИЩА

В работе [1] нами был представлен метод сегментации паренхиматозных органов брюшной полости. В работе [2] предложенная идея была расширена на более широкий класс органов, которые были условно обозначены как органы с равномерной текстурой. Отличительными особенностями этих органов являются небольшой диапазон интенсивностей внутри органов и большая вариативность интенсивностей на границах. Наш подход основан на идее, предложенной в работе [4] для классификации двумерных изображений методом текстурного анализа. В этой работе авторы используют матрицу ${{P}_{{\mathbf{d}}}}$ взаимной встречаемости пикселей ${{p}_{{\mathbf{d}}}}(i,j)$ для характеристики текстурных особенностей изображения. В нашей работе этот подход был расширен для случая трехмерных изображений, а дискретная энтропия использована для определения областей с равномерной текстурой:

Эта дискретная энтропия достигает максимального значения в случае, если число различных (по интенсивности соседей) пар смежных вокселей совпадает с числом пар соседей в рассматриваемой окрестности. Это свойство дискретной энтропии описывает две важные анатомические особенности, используемые в предложенном алгоритме сегментации: воксели внутри органов с равномерной текстурой имеют низкие значения дискретной энтропии, а воксели на границах – высокие значения. Более подробно с предложенным методом можно ознакомиться в работах [1], [2].

Для построения адаптивных неструктурированных расчетных тетраэдральных сеток был модернизирован сеточный генератор на основе открытой библиотеки CGAL [10]. В сеточный генератор был добавлен новый блок для сглаживания границы между разными органами и тканями и учета одномерных особенностей. Новый подход к построению сеток позволяет задавать область с помощью неявных функций (функций уровня), основанных на иерархическом делении пространства на отдельные органы и ткани.

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

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

3. ПРЯМОЕ МОДЕЛИРОВАНИЕ ЭКГ

Обозначим через $\Omega $ расчетную область, соответствующую туловищу или какой-то его части. Подобласть ${{\Omega }_{1}} \subset \Omega $ соответствует миокарду, а подобласть ${{\Omega }_{2}} = \Omega {\backslash }{{\bar {\Omega }}_{1}}$ – окружающим миокард тканям и органам, которые будут приняты во внимание при моделировании ЭКГ. Границу миокарда обозначим через $\Gamma = \partial {{\Omega }_{1}}$. В области $\Omega $ построена тетраэдральная сетка, каждому тетраэдру приписана метка, соответствующая номеру анатомической структуры в воксельной модели, которой он принадлежит. Проводимость некоторых анатомических структур будем считать анизотропной.

Мы продемонстрируем метод численного решения прямой задачи ЭКГ на примере работы [6]. Для этого мы будем использовать сеточную модель туловища, предоставленную исследовательской группой Института экспериментальной сердечно-сосудистой медицины Фрайбургского университета (Германия). В нашей работе мы дополним поставку задачи из [6] необходимыми граничными условиями и другими допущениями. Это позволит нам перейти к слабой постановке задачи и провести анализ ее корректности. Расчетная сетка для туловища $\Omega $ содержит 431 449 вершин и 2 435 642 тетраэдра (см. фиг. 3). Подобласть сердца ${{\Omega }_{1}}$ содержит 196 614 вершин и 950 763 тетраэдра (см. фиг. 4 слева).

Входными данными для рассматриваемой задачи являются:

1. Значение трансмембранного напряжения ${{V}_{m}} \in {{H}^{1}}({{\Omega }_{1}})$.

2. Усредненный по объему тензор внутриклеточной проводимости ${{\sigma }_{i}}$, заданный для всех материалов из подобласти ${{\Omega }_{1}}$.

3. Усредненный по объему тензор внеклеточной проводимости ${{\sigma }_{e}}$, заданный для всех материалов из области $\Omega $.

4. Ориентация векторов анизотропии u, заданная для каждого тетраэдра (фиг. 4 справа).

Тензоры ${{\sigma }_{i}}$ и ${{\sigma }_{e}}$ являются симметричными положительно определенными, их компоненты – ограниченные функции. Для расчета трансмембранного напряжения ${{V}_{m}}$ в клетках миокарда используется “бидоменная модель без бассейна” [12], в которой процессы, происходящие в ткани миокарда, моделируются комбинацией процессов, происходящих в двух наложенных друг на друга областях, одна из которых моделирует внутриклеточное пространство, а другая – внеклеточное, при этом не учитывается влияние окружающих органов (“бассейна”) на электрическую активность сердца. При выполнении условия непротекания тока через границу миокарда $\Gamma $ и предположения о том, что волокна тканей миокарда на его поверхности направлены вдоль поверхности, т.е. вдоль поверхности направлен и вектор анизотропии, выполняется следующее соотношение:

Неизвестной функцией является ${{\Phi }_{e}}$ – потенциал вне клетки, потенциал внутри клетки складывается из внеклеточного потенциала и трансмембранного напряжения: ${{\Phi }_{i}} = {{\Phi }_{e}} + {{V}_{m}}$. В данной работе были использованы значения трансмембранного напряжения ${{V}_{m}}$ для области ${{\Omega }_{1}}$, предоставленные рабочей группой Института экспериментальной сердечно-сосудистой медицины Фрайбургского университета (Германия).

Для нахождения внеклеточного потенциала необходимо решить следующую краевую задачу [6]:

где $\Phi _{e}^{{(k)}}$ и $\sigma _{e}^{{(k)}}$ – сужение функции ${{\Phi }_{e}}$ и тензора ${{\sigma }_{e}}$ на область ${{\Omega }_{k}}$, $k = 1,2$, соответственно, ${{{\mathbf{n}}}_{\Gamma }}$ – единичная нормаль к $\Gamma $ (для определенности будем считать, что нормаль является внешней по отношению к ${{\Omega }_{1}}$), ${\mathbf{n}}$ – единичная нормаль к $\partial \Omega $. Третье уравнение в (3) обеспечивает непрерывность внеклеточного потенциала. Четвертое уравнение получается из следующих соображений: во-первых, нормальная компонента внеклеточного тока через границу миокарда должна быть непрерывной, т.е. $\sigma _{e}^{{(1)}}\nabla \Phi _{e}^{{(1)}} \cdot {{{\mathbf{n}}}_{\Gamma }} - \sigma _{e}^{{(2)}}\nabla \Phi _{e}^{{(2)}} \cdot {{{\mathbf{n}}}_{\Gamma }} = 0$ на $\Gamma $, во-вторых, внутриклеточное пространство изолировано от окружающих миокард тканей, т.е. ${{\sigma }_{i}}\nabla {{\Phi }_{i}} \cdot {{{\mathbf{n}}}_{\Gamma }} = {{\sigma }_{i}}\nabla (\Phi _{e}^{{(1)}} + {{V}_{m}}) \cdot {{{\mathbf{n}}}_{\Gamma }} = 0$ на $\Gamma $. Складывая последние два уравнения и учитывая (2), приходим к четвертому уравнению в (3). Пятое уравнение обеспечивает непротекание тока через границу туловища.

4. СЛАБАЯ ПОСТАНОВКА ЗАДАЧИ

Для численного решения задачи (3) будем использовать метод конечных элементов. Заметим, что если ${{\Phi }_{e}}$ – решение задачи (3), то и ${{\Phi }_{e}} + c$, где $c$ – произвольная постоянная функция, тоже является решением этой задачи. Чтобы обеспечить единственность решения, необходимо наложить дополнительное условие на ${{\Phi }_{e}}$, а именно, положим, что ${{\Phi }_{e}} \in {{\tilde {H}}^{1}}(\Omega )$, где

т.е. в качестве решения будем искать функцию, ортогональную постоянной функции в смысле скалярного произведения в ${{L}^{2}}(\Omega )$. Перейдем к слабой постановке исходной задачи: умножим правую и левую части первых двух уравнений (3) на произвольную функцию $\psi \in {{H}^{1}}(\Omega )$ и проинтегрируем по ${{\Omega }_{1}}$ и ${{\Omega }_{2}}$ соответственно. Применяя формулу Грина и принимая во внимание однородное условие Неймана, получаем Во втором уравнении учтено, что ${{{\mathbf{n}}}_{\Gamma }}$ является внутренней нормалью по отношению к ${{\Omega }_{2}}$. Складывая эти уравнения и используя третье и четвертое уравнения из (3), приходим к следующей задаче: найти функцию ${{\Phi }_{e}} \in {{\tilde {H}}^{1}}(\Omega )$ такую, что Здесь введены обозначения

Функция ${{\Phi }_{e}}$, являющаяся решением задачи (4), называется слабым решением задачи (3). Если функция является решением задачи (3) в классическом смысле, то она, очевидно, является и слабым решением. Если слабое решение обладает достаточной гладкостью, то оно является решением задачи (3) в классическом смысле. Подставляя $\psi \equiv 1$ в (4) и используя (5), получаем необходимое условие существования слабого решения:

Нетрудно видеть, что благодаря (2) функция $f$ из (5) удовлетворяют данному условию: В ходе работы была доказана теорема о существовании и единственности решения в слабой постановке.

Отметим, что полученные на практике значения ${{V}_{m}}$ являются лишь приближенным решением бидоменной задачи, и в отдельных случаях условие (6) может нарушаться на величину, близкую к машинной точности или точности решения бидоменной задачи. Однако при решении прямой задачи ЭКГ мы также ищем лишь приближенное решение задачи (3), поэтому нарушение условия (6) в пределах допустимой погрешности на практике не приводит к проблемам.

5. РЕЗУЛЬТАТЫ

Фиг. 5 демонстрирует исходные значения трансмембранных напряжений на желудочках сердца и результаты расчета распределения электрических потенциалов на поверхности туловища ($t = 320$ мс). Для расчета электрокардиографических отведений необходимы значения электрических потенциалов, соответствующих положениям электродов на поверхности тела (см. фиг. 3 справа).

1. Эйнтховен I: ${{\Phi }_{e}}(LA) - {{\Phi }_{e}}(RA)$.

2. Эйнтховен II: ${{\Phi }_{e}}(LL) - {{\Phi }_{e}}(RA)$.

3. Эйнтховен III: ${{\Phi }_{e}}(LL) - {{\Phi }_{e}}(LA)$.

В данной работе было исследовано влияние различных органов, включенных в воксельную модель туловища, на корректность моделирования ЭКГ. В численных экспериментах проводимости отдельных органов последовательно заменялись на среднюю проводимость туловища. Кроме того, проводилось сравнение с моделью, в которой проводимости всех анатомических структур были заменены на усредненные. Фигура 6 демонстрирует результаты: печень и легкие в наибольшей степени влияют на корректность моделирования ЭКГ.

6. БЫСТРЫЙ СПОСОБ РАСЧЕТА ЭЛЕКТРИЧЕСКИХ ПОТЕНЦИАЛОВ В НЕКОТОРЫХ ТОЧКАХ ПОВЕРХНОСТИ ТЕЛА

Для построения электрокардиограммы по результатам математического моделирования необходимо решить систему (3) для разных моментов времени сердечного цикла (для каждого момента времени на вход подаются новые значения для трансмембранного напряжения ${{V}_{m}}$). При этом требуется найти значения потенциала ${{\Phi }_{e}}$ лишь в точках поверхности, соответствующих положениям электродов. Приведем быстрый способ поиска решения системы $A{\mathbf{x}} = {\mathbf{f}}$ только в $K$ заданных узлах сетки, где $A$ – симметричная положительно определенная матрица порядка $N$. Будем искать ${{{\mathbf{x}}}_{s}} \in {{\mathbb{R}}^{K}}$ – вектор, содержащий решение в выбранных узлах сетки. Пусть ${{M}_{s}}$ – матрица размера $K \times N$, составленная из $K$ строк матрицы ${{A}^{{ - 1}}}$, соответствующих выбранным узлам сетки, тогда ${{{\mathbf{x}}}_{s}} = {{M}_{s}}{\mathbf{f}}$. Через ${{{\mathbf{m}}}_{i}}$ обозначим $i$-ю строку матрицы ${{A}^{{ - 1}}}$. Для того, чтобы найти ${{{\mathbf{m}}}_{i}}$, нужно решить систему ${{{\mathbf{m}}}_{i}}A = {\mathbf{e}}_{i}^{{\text{т}}}$, где ${{{\mathbf{e}}}_{{\mathbf{i}}}}$ – базисный вектор-столбец. Поскольку $A = {{A}^{{\text{т}}}}$, имеем $A{\mathbf{m}}_{i}^{{\text{т}}} = {{{\mathbf{e}}}_{i}}$, т.е. матрица ${{M}_{s}}$ строится по $K$ решениям исходной системы.

Как видно из (5), если в начале вычислений сгенерировать матрицу ${{A}_{{in}}}$ конечно-элементной дискретизации оператора $\nabla \cdot {{\sigma }_{{\mathbf{i}}}}\nabla $, то на каждом шаге по времени для вычисления правой части ${\mathbf{f}}$ достаточно умножать вектор значений трансмембранного напряжения ${{V}_{m}}$ на эту матрицу. В табл. 1 приводится сравнение обычного и ускоренного алгоритмов. Численные эксперименты проводились для $K = 3$.

На вход алгоритмам подается сетка, тензоры проводимости и массив индексов узлов сетки, соответствующих расположению электродов. В строках 2–4 генерируются матрицы $A$ и ${{A}_{{in}}}$ и строится ILU предобуславливатель для $A$. В строках 5–8 ускоренного алгоритма вычисляется матрица ${{M}_{s}}$, строки которой суть строки матрицы ${{A}^{{ - 1}}}$, соответствующие выбранным узлам. Для этого $K$ раз решается система линейных уравнений размерности $N$. В строках 9–16 реализуется цикл по времени, состоящий из ${{N}_{t}}$ шагов. В строках 10–11 считывается вектор значений трансмембранного напряжения ${{V}_{m}}$ для текущего временного шага и на его основе генерируется правая часть. В строке 12 метода без ускорения решается система уравнений размерности $N$, тогда как в ускоренном методе производится умножение матрицы размерности $K \times N$ на вектор. В строках 13–15 вычисляются требуемые разности потенциалов на текущем шаге по времени. В строке 17 осуществляется вывод рассчитанной ЭКГ.

Таким образом, в методе без ускорения требуется ${{N}_{t}}$ раз решать СЛАУ с матрицей порядка $N$, а в ускоренном методе требуется решать систему с этой же матрицей только $K$ раз и ${{N}_{t}}$ раз умножать матрицу размерности $K \times N$ на вектор. Пусть алгоритмическая сложность решения СЛАУ с матрицей $A$ с точностью $\varepsilon $ есть ${{C}_{\varepsilon }}(N)$. Тогда алгоритмическая сложность метода без ускорения есть ${{N}_{t}}{{C}_{\varepsilon }}(N)$, a ускоренного метода составляет $K({{C}_{\varepsilon }}(N) + N) \approx K{{C}_{\varepsilon }}(N)$. В данной задаче $K < 10$, в то время как число моментов времени сердечного цикла, на которых необходимо найти решение системы, составляет несколько сотен. Приведенный способ расчета значительно ускоряет вычислительный процесс.

Алгоритм 1

Алгоритм решения прямой задачи ЭКГ и вариант его ускорения для $K = 3$ (${{N}_{t}}$ – число шагов по времени, $idE$ – массив индексов электродов, используемых для расчета отведений, ${{E}_{1}}$, ${{E}_{2}}$, ${{E}_{3}}$ – отведения I, II, III)

Таблица 1 демонстрирует среднее время работы алгоритмов технологической цепочки решения прямой задачи ЭКГ.

7. ВЫВОДЫ

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