РЕШЕНИЕ ОБРАТНОЙ ЗАДАЧИ МЭГ В МНОГОДИПОЛЬНОЙ МОДЕЛИ Текст научной статьи по специальности «Математика»

ВЕСТН. МОСК. УН-ТА. СЕР. 15. ВЫЧИСЛ. МАТЕМ. И КИБЕРН. 2023. № 3. С. 27-32 Lomonosov Computational Mathematics and Cybernetics Journal

УДК 519.2

Т. В. Захарова1 , А. И. Сабиров2

РЕШЕНИЕ ОБРАТНОЙ ЗАДАЧИ МЭГ В МНОГОДИПОЛЬНОЙ МОДЕЛИ

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

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

ВОТ: 10.55959/М8и/0137-0782-15-2023-47-3-27-32

1. Введение. Магнитоэнцефалография производит большой объем данных, и обработка этих данных с целью восстановления источников сигнала с заданной точностью является чрезвычайно некорректной задачей [1-3].

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

2. Прямая задача. Будем считать, что форма головы достаточна близка к шару, и что во всем объеме головы сохраняется одинаковая проводимость. В таких условиях модель не перестанет отражать фундаментальные свойства, и при этом будет достаточно простой для изучения [4].

Одним из получаемых магнитоэнцефалографией сигналов является нормальная компонента вектора магнитной индукции. Магнитная индукция в точке г, вызванная токами с плотностью ] в объеме V, вычисляется по закону Био-Савара-Лапласа:

в(г) = !£[т х (г 7 ¿У. (1)

4п } |г - г'|3 V

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

Так как первичный ток локализован внутри нейронов и имеет малую длину (много меньше всех других характерных расстояний в задаче), то его можно представить в виде токового диполя с дипольным моментом Q, для которого плотность тока задается с помощью дельта-функции:

Зр(Г) = QS(r - г').

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

Но Q х (г - г')

| r _r'

вд=g -,.з • (2)

1 Факультет ВМК МГУ, доц., к.ф.-м.н., e-mail: [email protected]

2 Факультет ВМК МГУ, студ. магистратуры, e-mail: [email protected]

Из принципа суперпозиции или из той же интегральной формулы закона Био-Савара-Лапласа выводится формула для магнитной индукции, вызванной несколькими диполями:

M д0 M Qj х (r - rj) j=i j=1 |r 'j1

Объемные токи сильно усложняют задачу в общем случае. Составляющая магнитного поля Bv, связанная с объемными токами, может быть определена аналитически только для некоторых случаев. Для сферической модели нормальная компонента Bv точно равна нулю [5], а для эллипсоидной модели Bv была рассчитана в [6].

В сферическом случае получается сильно упростить вычисление нормальной компоненты магнитной индукции, так как в формулу (1) дает вклад только первичный ток.

Наблюдаемая величина к(г) = Bn(r) (нормальная компонента вектора магнитной индукции относительно поверхности головы) приближенно вычисляется с помощью формулы:

M

r

j= |r - rj1

Наблюдаемый сигнал в прямой задаче находится следующим образом:

M

r

K(r,t) = -—з . (Qj(t) X r'3(t)). (5)

j=1 |r - rj(t)|

3. Обратная задача. Обратная задача магнитоэнцефалографии заключается в нахождении числа M источников активности, формирующих магнитное поле на поверхности П головы, их дипольных моментов {Qj(í)}j=i и расположения {rj(í)}j=i по N наблюдениям {к(',t)}N=1.

Количество точек наблюдений ограничено (от нескольких десятков в ЭЭГ до нескольких сотен в МЭГ), при этом часть сигналов теряется при фильтрации сигналов. Поэтому в задаче важно добиться оптимального качества при наименьшем числе наблюдений.

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

Авторами разработаны два метода локализации. Первый метод использует экстремумы магнитного поля [7]. Это точный метод решения, опирающийся на точки минимума и максимума нормальной компоненты магнитной индукции [8]:

в _ № QtqT sin в cos ф ^

4п (r2 + rQ — 2rrQ cos в)3/2

Исходя из формулы при ф = ±п/2 нормаль магнитной индукции равна нулю, так определяется плоскость симметрии. Точки экстремумов можно выразить через следующие уравнения:

cos ф = ±1, (7)

r4 + 14r2rQ + rQ - (r2 + rQ)

cos 0 = —-. (8)

2rrQ v 7

Тогда диполь будет лежать в плоскости симметрии, причем на оси, в которой лежит центр сферы. Модуль радиуса-вектора rQ задается через 9m - средний угол точек экстремумов:

3 - cos2 вт - V9 — 10 cos2 вт + cos4 вт

го = г---. 9

4 2 COS в т.

Известно, что геометрию типичного человеческого мозга можно аппроксимировать как эллипсоид с полуосями 6 х 6, 5 х 9 см [9,10]. Так как две из полуосей примерно равны друг другу, то можно рассмотреть приближение эллипсоида вращения, а именно вытянутого эллипсоида. Таким образом, эллипсоид аппроксимирует человеческий мозг намного точнее, чем сфера. Для эллипсоидальной модели прямая задача, как и в случае сферической геометрии, имеет точное аналитическое решение [6,11], причем качественно изображения магнитного поля в обоих случаях очень похожи: имеется по 2 максимума нормальной компоненты |ВП|. Это послужило основанием для поиска аналогичного решения обратной задачи.

Для диполей, лежащих не слишком глубоко внутри головного мозга, rQ > К/2, первичное магнитное поле доминирует над объемным ВП/ВП < 0,1 [6], поэтому можно приближенно использовать результат сферической модели. Для этого надо провести сечение эллипсоида, проходящее через точки максимумов по нормали к его поверхности, найти радиус кривизны полученного сечения вблизи точек максимумов и воспользоваться формулой (9).

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

Рассмотрим метод локализации через симметрию.

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

г' = (0,0, г)т, э = Сь, 0,зг)т, г = (Гх,Гу,г2)т, (10)

«г) = —~7|3 ' О' х г') = (11)

\г — г \ (1 + г2-2ггг)2

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

2

Используем его для формирования квадратного уравнения:

1 + г2 - 2гГ2г

7=1 + ,2-2^' (13)

,2 - + 1=0. (14)

7 - 1

Введем новый коэффициент, упрощающий вид финального решения:

Шг - Г2г

/3 = - ' , (15)

7 - 1

г = -/3 + 81§п(/3) V/?2 - 1- (16)

Нужно лишь найти эту плоскость и ось. Основываясь на симметрии картин наблюдений, мы

можем получить эти формулы для поиска оси и нормали к плоскости.

Формула для направляющего вектора оси, на которой лежат диполь и центр сферы:

ах = ^2 |я(Гг)|Гг- (17)

г

Нормаль к плоскости, на которой лежат диполь с дипольным моментом и центр сферы:

п = ^2 к(гг)гг. (18)

На практике оказывается, что следующие формулы оценивают параметры лучше:

*—л 1

(19)

аж

п = ^sign(/í(r¿))|/í(r¿)|9r¿. (20)

г

Так как из каждой пары наблюдений получается по одной оценке параметров, всего мы можем построить N (Ж — 1) оценок. И для каждой из них можно построить собственные реализации картин наблюдения и сравнить их с реальными. Мы можем брать такую оценку, которая минимизировала бы разницу с реальными наблюдениями или какую-нибудь взвешенную.

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

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

Предположим, что различные функциональные области головного мозга независимы, и перенесем это свойство на токи. В этом случае мы можем использовать метод независимых компонент. Наблюдения к(*) = {к(тгбудут результатом некоторого преобразования М независимых источников = {т'у(¿)}М=1. Дополнительно предположим, что источники центрированы:

Св(*) = ^(¿М^ = I, Ев(*) = 0 V*, (21)

к(*) = А(в(*)). (22)

Если дополнительно предположить, что т^-(¿) = т^- V* (координаты диполей стационарны), то преобразование А будет линейным:

к(*) = Авф, (23)

= (ААТ )-1Ат к(*) = Вк(*), (24)

С = Е [Вк(*)к(*)т Вт] = ВСКВт, (25)

ВСКВт = I. (26)

Используя псевдообратное преобразование и следствие независимости, мы получаем первый шаг решения — сделать наблюдения сферическими ¿(¿) = £к(*), С = I .И после искать матрицу весов Ш, которая преобразует ¿(¿) к сигналу с независимыми компонентами:

= В = (27)

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

Поэтому, можно использовать "негауссовость" как меру независимости. А в качестве "негаус-совости" — абсолютное значение коэффициента эксцесса:

кигф,) = Ев4(*) — 3Е2 в?(*), (28)

вг(г) = и)тг(Ь), и)тw = |И|2 = 1. (29)

Из сферичности получим

= 4sign(kurt(адтz))№(адт,г:(í))3] - Зад£[(адт,г(£))2]), (30)

Е[(-шт г(г))2] = Е[г(г)г(г)т = |И|2 = 1, (31)

= 481ёп(киЛ(и,тг))(Е[г(и,тг(г))3} - Зад) = 0, (32)

Е[х(г)(ь]тх(г))3] - = 0. (33)

В этом выражении первое слагаемое отвечает за направление вектора а второе — за его длину. И мы можем обобщить его, используя функцию д замены кубической функции над ,штг(Ь) и коэффициент в (нормировочный коэффициент (-3)), при помощи нижеследующих рассуждений:

Е[г(г)д(ь)т г(Щ + = 0. (34)

С помощью метода Ньютона получим выражение

д Р

— =Е1г(Ы1)тд>(и;тгт+(3- (35)

Из сферичности г следует, что Е[х(1)х(1)Тд'(и!Тг(£))] ~ Е[д'(мтг(£))] и

Е|а'(штг(4))]+/)

или

w ^ Е[г(г)д(ь)тг(Щ - Е[д'(ытг(г))№, (37)

ад тт^-. (38)

1М|

Итак, для максимизации эксцесса киЛ(в) с вектором наблюдений г выбирается некоторый начальный вектор w, рассчитывается направление, по которому абсолютное значение эксцесса случайной величины в = wтг растет наиболее быстро, и далее вектор w сдвигается в этом направлении. Решение состоит в построении последовательности векторов w увеличивающих эксцесс киЛ(в), т.е. меру "негауссовости".

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

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

4. Заключение. Настоящая работа посвящена изучению функциональных зон головного мозга человека. Функциональное картирование коры головного мозга является чрезвычайно сложной задачей, возникновение которой обусловлено современным уровнем развития методов неин-вазивного исследования головного мозга. И одним из таких методов является магнитоэнцефа-лография. Методы локализации важны в реальной клинической практике. Во время нейрохирургических вмешательств могут быть повреждены различные участки головного мозга, в том числе и невосполнимые. Поскольку расположение функциональных зон в головном мозге человека индивидуально, врач должен уметь локализовать эти зоны в предоперационном периоде

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

СПИСОК ЛИТЕРАТУРЫ

1. B a i 11 e t S., Mosher J. C., Leahy R.M. Electromagnetic brain mapping / / IEEE Signal Processing Magazine. IEEE, 2001. P. 14-30.

2. Sarvas J. Basic mathematical and electromagnetic concepts of the biomagnetic inverse problem // Physics in Medicine and Biology. 1987. 32. P. 11-22.

3. Fr i s t o n K., Harrison L., D a u n i z e a u J., K i e b e l S., Phillips C., et al. Multiple sparse priors for the MEG inverse problem // Neurolmage. 2008. 39. P. 1104-1120.

4. 11 m o n i e m i R. J., H a m a l a i n e n M.S., K n u u t i l a J. The forward and inverse problems in the spherical model // Biomagnetism: Applications and Theory. Eds by: Weinberg H., Stroink G., Katila T. New York: Pergamon, 1985. P. 278-282.

5. Hamalainen M.S., Hari R., Ilmoniemi R.J., Knuutila J., Lounasmaa O.V. Magnetoencephalography — theory, instrumentation, and applications to noninvasive studies of the working human brain // Reviews of Modern Physics. 1993. 65. P. 413-497.

6. C u f f i n B., Cohen D. Magnetic fields of a dipole in special volume conductor shapes // IEEE Trans. Biomed. Eng. 1977. BME-24. P. 372-381.

7. K a r p o v P. I., Z a k h a r o v a T. V. Magnetoencephalography inverse problem in the spheroid geometry // Journal of Inverse and Ill-Posed Problems. 2019. 27. N 2. P. 159—169.

8. Zakharova T.V., Karpov P. I., Bugaevskii V. M. Localization of the activity source in the inverse problem of magnetoencephalography // Computational Mathematics and Modeling. 2017. 28. N 2. P. 148-157.

9. D a s s i o s G. The magnetic potential for the ellipsoidal MEG problem // J. Comput. Math. 2007. 25. P. 145-156.

10. U i t e r t R., We instein D., Johnson C. Can a spherical model substitute for a realistic head model in forward and inverse MEG simulations? // Proc. 13th Int. Conf. on Biomagnetism. Jena, Germany, 2002. P. 798-800.

11. De Munck J. C. The potential distribution in a layered anisotropic spheroidal volume conductor //J. Appl. Phys. 1988. 64. P. 464-470.

Поступила в редакцию 15.05.23 Одобрена после рецензирования 24.05.23 Принята к публикации 24.05.23