CC BY
25
УДК 519.612:632.4
Е. В. Захаров, Р. Е. Зимоздра2
О ГРАНИЦАХ ПРИМЕНИМОСТИ СФЕРИЧЕСКОЙ МОДЕЛИ ДЛЯ РЕШЕНИЯ ЗАДАЧ ЭЛЕКТРОЭНЦЕФАЛОГРАФИИ
Работа посвящена развитию методов локализации нейронных источников мозга, чей потенциал регистрируется в виде электроэнцефалограммы (ЭЭГ). Рассматривается краевая задача для уравнения Пуассона. Выведено граничное интегральное уравнение Фредгольма II рода. Предложена методика численного решения интегрального уравнения, приводятся результаты некоторых вычислительных экспериментов.
Ключевые слова: ЭЭГ, прямая задача, интегральные уравнения.
1. Введение. Головной мозг представляет собой проводящую среду, в объеме которой распределено множество электрически активных элементов (нейронов и глиальных клеток), являющихся первичными источниками электроэнцефалограммы (ЭЭГ), регистрируемой на поверхности головы. Регистрация ЭЭГ — один из неинвазивных способов получения физиологической и клинико-диагностической информации о деятельности головного мозга [1]. Метод ЭЭГ позволяет судить о наличии очаговых поражений, общемозговых расстройств, а также их характере. Одним из важных приложений ЭЭГ является диагностика эпилепсии. Известно, что за генерацию эпилептического припадка отвечают локальные группы нервных клеток, так называемые очаги. На основе ЭЭГ-данных неврологи могут локализовать эпилептогенный очаг [2].
Целью данной работы является определение границ применимости основных математических моделей, применяемых в ЭЭГ, а также ответ на вопрос: насколько влияют неоднородности внутреннего строения головы на точность решения прямых задач.
2. Постановка задачи. Прямая задача ЭЭГ состоит в определении потенциала электрического поля на поверхности головы. Делая различные предположения о распределении и количестве источников, можно добиться совпадения экспериментального потенциала с модельным. Прямая задача должна учитывать структуры тканей мозга с различными проводимостями и другие геометрические параметры для правильного вычисления индуцированного потенциала.
Обратная задача ЭЭГ состоит в определении положений источников по измеренным потенциалам на поверхности головы. При этом для решения каждой обратной задачи требуется решение большого числа прямых задач. Известно, что существует бесконечное число конфигураций внутренних источников тока, которые могут создавать любое заданное на поверхности головы распределение потенциала. Это означает, что обратная задача ЭЭГ в общей постановке не имеет единственного решения. Чтобы сделать решение задачи единственным, приходится вводить некоторую модель источников тока.
Одной из них является модель дипольных источников. В первых работах по ЭЭГ диполи использовались из соображений простоты вычислений и сокращения машинного времени. Однако, как выяснилось, дипольные источники моделируют большинство наблюдаемых электрофизиологических явлений.
Впрочем, большинство авторов чаще всего используют представление не о диполе как таковом, а об эквивалентном диполе. Источник мозговой активности охватывает большую группу нейронов, однако, если считать, что эта область мала по сравнению с расстоянием до отводящих электродов, ее можно представить в виде "эквивалентного диполя" [3]. Неизвестных в задаче в этом случае шесть: три пространственные координаты расположения диполя и три компоненты вектора тока, характеризующего направление и величину момента диполя.
Перейдем к постановке задачи: рассмотрим трехмерную область П, ограниченную поверхностью дО,. В квазистатическом приближении электрическое поле головного мозга описывается потенциалом
1 Факультет ВМК МГУ, проф., д.ф.-м.н., e-mail: zspecQcs.msu.ru
2 Факультет ВМК МГУ, асп., e-mail: bigzimQyandex.ru
и(М), который в области П есть решение краевой задачи:
А = МеП, = 0 Ред п
ст о оп
где / — плотность источника, сто — проводимость, п — внешняя нормаль к поверхности.
Исходя из электропроводности тканей мозга, можно показать, что время распространения сигнала очень мало. Таким образом, задача решается на одном срезе по времени. При выполнении принципа непрерывности сигнала по времени можно предположить, что найденные решения будут близки в близкие моменты времени [2].
3. Сферическая модель: метод интегральных уравнений. Общепринятой моделью в ЭЭГ считается так называемая сферическая модель головы, где в качестве основного объекта берется обладающая сферической симметрией область, заполненная однородной (или неоднородной) средой с помещенными в нее источниками дипольного типа. Выясним, насколько данная модель отражает анатомические и геометрические особенности строения головы.
3.1. Вывод интегрального уравнения. Представим потенциал и(М) в виде суммы и(М) = = щ(М) + у(М), где ио — потенциал, создаваемый источниками, а V — неизвестный потенциал индуцированного поля, который создается из-за наличия границ и неоднородной электропроводности. Если известны положения диполей и их моменты, тогда потенциал щ, создаваемый диполями, известен и требуется найти индуцированный потенциал V.
В областях однородности проводимости для функции у(М) имеем задачу Неймана для уравнения Лапласа
Ау(М) = О, М е П,
с неоднородным граничным условием на скальпе (граница с непроводящей средой):
ду(Р) дщ{Р)
=--= *Ь(Р), Р € дО.
Представим V в виде потенциала простого слоя:
„(Р)=/М^,,
J У Нмр
дп
продифференцируем по п и, воспользовавшись теоремой о скачке, получим для плотности потенциала ср(М) интегральное уравнение Фредгольма второго рода:
2жср(Р) + // ^ ё(™ = (1)
дп
где М и Р — произвольные точки области П (Р — точка наблюдения, М — точка интегрирования), ду(Р)
*Ь(Р) = —г--первичное поле источника (диполя).
оп
3.2. Сведение к СЛАУ. Приближенное решение ф(М) интегрального уравнения (1) представим некоторой функцией ф(М, а,..., сп), зависящей от параметров с^, к = 1,..., п, определяемых в процессе решения. При определении свободных параметров с^ воспользуемся невязкой уравнения (1):
г{Р) = 2тгф{Р) + Л ^ {н~) Ф{-М) ё(™ ~ дп
При этом точное решение <р(М) обращает невязку г в нуль. Потребуем, чтобы при подстановке в невязку г приближенного решения ф(М, с\,..., сп) она обращалась в нуль в заданной системе точек ] = 1,... ,п, поверхности дО,. Наиболее удобно приближенное решение ф(М) представить в виде линейной комбинации
п
ф{М) = ^скфк{М),
к=1
где фк(М) — характеристические функции. Чтобы ввести их, разобьем поверхность dil на п конечных элементов д£1к и положим
fl, если МедПк, v ; [0, если М i дПк.
Для определения совокупности коэффициентов {с*;} получаем следующую систему линейных алгебраических уравнений:
Y^ ckAjk = bj, Ajk = 2тгSjk I jj ( д^) d(JM, k=1 an
Sjk — дельта Кронекера, bj = FQ(Pj).
3.3. Выбор конечных элементов. Большое значение имеет число конечных элементов, которое должно быть не меньше числа искомых параметров: 6 для однодипольной модели, 12 для двухди-польной и т. д. Однако, учитывая неустойчивость решения, связанную с шумами, желательно иметь большее число элементов — оптимально 20 для одно- и двухдипольной моделей.
Аппроксимируем сферическую поверхность dil описанным многогранником с достаточно большим числом равновеликих граней д£1к. Среди правильных многогранников наибольшим является икосаэдр, у которого 20 треугольных граней. В качестве точек Pj выберем их центры. К сожалению, 20-гранная аппроксимация сферы выглядит довольно грубой, но есть простой способ улучшить ее точность. Представим, что икосаэдр вписан в сферу и увеличим число граней в четыре раза путем подразбиения каждого ребра на две равные части. Новые вершины окажутся внутри сферы, поэтому сдвинем их к поверхности. Процесс измельчения может быть повторен неоднократно.
4. Сферическая модель: метод сферических функций. Человеческий мозг окружен четырьмя основными слоями, которые значительно отличаются по электропроводности и соответственно влияют на измерение потенциалов на поверхности головы. Речь идет о спинномозговой жидкости, твердой мозговой оболочке, кости черепа и коже скальпа. Значения электропроводности G указанных слоев, по данным Гнездицкого [3], составляют: для спинномозговой жидкости G = 1 (Ом-м)"1; для черепа 0.04 (Ом-м)-1; для скальпа 0.28 (Ом-м)-1 (почти такая же проводимость, как у мозговой ткани). Средние размеры толщины слоев колеблются в районе 2 мм (для твердой мозговой оболочки), 8 мм (для черепа) и 4 мм (для скальпа) при радиусе кривизны головы 8-9 см.
Наличие вышеуказанных слоев вводит электрические неоднородности, которые могут привести к изменению потенциального поля. Присутствие тканей с высокой проводимостью приводит к смазыванию и усреднению картины, а слои с низкой проводимостью вызывают ослабление поля при регистрации на поверхности.
4.1. Однородная модель. Пусть диполь, имеющий радиальную и тангенциальную компоненты тг и mt, помещен внутрь однородной сферы радиуса R и проводимости а. Без ограничения общности предположим, что диполь находится на оси z на расстоянии г от центра, а дипольный момент располагается в верхней полуплоскости xz (любая другая позиция или ориентация диполя может быть получена вращением).
Тогда поверхностная плотность потенциала (р(а,/3) в сферической системе координат задается формулой, выведенной в [4]:
(р(а,Р) = -— ^-Ьп-1 (птгРп( cosa) + mtPn (cos a) cos /3) , (2)
n=1
где b = yi — эксцентриситет диполя в однородной модели, Рп(cos а) и P^(cos a) — полиномы Лежандра и присоединенные полиномы Лежандра соответственно.
4.2. Трехслойная модель. Теперь рассмотрим слоистую модель, состоящую из однородной сферы нервной ткани радиуса ri, окруженной концентрическим слоем внешнего радиуса г2, представляющим череп, и другим концентрическим слоем внешнего радиуса R, представляющим скальп. Пусть а
обозначает проводимость нервной ткани и скальпа (их проводимости принимаются равными, исходя из медицинских данных), проводимость черепа обозначим as.
Тогда поверхностная плотность потенциала ф в неоднородной модели задается следующей формулой:
1 + £(2n + I)2 ф(а,/3) = -—У,-Ьп~ ——-— (nfhrPn (cos а) + fhtPn (cos а) cos /3) , (3)
471(7 ТЬ (^ J-)
п= 1 4 '
dn = ((п + 1)С + п) + (1 - О СГ+1 - /Г+1) +1) - п(1 - с)2 (I)2n+1,
i = - f 1 = - f2 = -
S ! Л j-, 1 J J- j-, '
a К R
Всюду ниже используются следующие значения параметров: £ = 0.0125, fi = 0.87, /2 = 0.92.
Оценим, какой вклад вносит неоднородность. Для радиального (mt = 0) диполя с известным эксцентриситетом b в трехслойной модели вычислим потенциал по формуле (3) и, многократно применяя формулу (2), найдем соответствующий ему диполь в однородной модели, который генерирует поле с наименьшим отклонением от заданного.
Минимизация табличной функции одной переменной была осуществлена с помощью численных методов. Полученные результаты приведены в следующей таблице:
ь 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8
ь 0.0615 0.12325 0.1853 0.2476 0.3097 0.3705 0.4276 0.4781
ъ/ъ 0.615 0.61625 0.6177 0.619 0.6194 0.6175 0.611 0.597625
Похожую задачу ставили авторы в [4]; оказалось, что в обоих случаях зависимость между эксцентриситетами аппроксимируется линейной функцией. Следовательно, любое положение источника в неоднородной модели может быть предсказано исходя из однородной модели.
5. Модификации сферической модели.
5.1. Модель "полусферы". Поскольку мозговая ткань сосредоточена в верхней части головы, логичным шагом будет использование модели, связанной с "полусферой". Для введения конечных элементов достаточно взять икосаэдр и удалить в нем нижнюю вершину вместе с пятью примыкающими к ней гранями. Полученное геометрическое тело (состоящее из 15 треугольных граней и пятиугольного основания, допускающего триангуляцию) аппроксимирует усеченную сферу, причем усечение проводится ниже центра, — поэтому слово "полусфера" в названии модели взято в кавычки.
При локализации источников в верхней части головы указанная модель не дает существенного преимущества над сферической.
5.2. Модель эллипсоида. Выясним, как сказывается на решении несферичность черепа. Для этого выберем декартову систему координат, связанную с реальным объектом (головой). Ось х проведем в направлении от выступающей костной точки затылочного бугра к переносице, ось у — от правого к левому слуховому проходу, ось z — от самой нижней до самой верхней точки. Согласно Гнездицкому [3], средние размеры полуосей, найденные по результатам измерений у двадцати испытуемых, составляют: Ох = 9.5 см, Оу = 8 см, Oz = 12 см. Таким образом, голова лучше аппроксимируется эллипсоидом с указанными длинами полуосей.
При сравнительной локализации источника наихудшие результаты достигались на малых эксцентриситетах (,b0z < 0.3, Ь0 х < 0.2, Ьоу < 0.5) по каждой из осей. Также для диполей, выбранных в модели эллипсоида с большим эксцентриситетом (boz > 0.7) по z, не удалось подобрать соответствующий диполь в сферической модели, так что в данных случаях сферическая модель оказывалась неприменима.
Ьэл 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
Ьсф 0Oz — — 0.225 0.391 0.543 0.696 0.888 — —
Ъсф Ох 0.283 0.219 0.299 0.386 0.474 0.5625 0.647 0.725 0.786
Ъсф 0.387 0.437 0.489 0.545 0.601 0.658 0.724 0.797 0.877
6. Заключение. Суммируя вышеизложенное, приходим к выводу, что однородная сферическая модель довольно проста, не требует громоздких вычислений и применима для локализации большинства дипольных источников. Таким образом, ослабление потенциалов черепом и эффекты "смазывания" поля не являются решающими факторами, влияющими на точность локализации.
В то же время деформация сферической поверхности (уплощение вдоль осей х ж у) сильнее всего влияет на потенциальное поле. Отклонение черепа от сферы оказывает наибольшее влияние при локализации источника в теменной доле головного мозга. При ином расположении источника локализация гораздо меньше зависит от используемой модели.
СПИСОК ЛИТЕРАТУРЫ
1. Захаров Е.В., Коптелов Ю.М. О решении одной задачи математической обработки электроэнцефалографических данных // ДАН СССР. 1987. 292. № 3. С. 576-581.
2. Попов А. М. Решение обратной задачи ЭЭГ с помощью стохастических методов распознавания образов // Вестн. Моск. ун-та. Сер. 15. Вычисл. матем. и киберн. 2006. № 3. С. 91-100.
3. Гнездицкий В. В. Обратная задача ЭЭГ и клиническая электроэнцефалография (картирование и локализация источников электрической активности мозга). М.: МЕДПресс Информ, 2004.
4. Ary J. Р., Klein S.A., Fender D.H. Location of sources of evoked scalp potentials: corrections for skull and scalp thickness // IEEE Trans, on Biomed. Eng. 1981. 28. N 6. P. 447-452.
Поступила в редакцию 19.02.14
ABOUT APPLICABILITY OF THE SPHERICAL MODEL FOR SOLVING THE ELECTROENCEPHALOGRAPHY PROBLEMS
Zakharov E. V., Zimozdra R. E.
This paper concerns the evolution of the cerebral sources localization techniques based on electroencephalogram (EEG) registration. Neumann problem for the Poisson's equation is considered. Eredholm integral equation of the second kind is established. Numerical solution technique is offered, and some numerical results are shown.
Keywords: EEG, forward problem, integral equations.
