CC BY
48
УДК 681.3:612.8
doi:10.15217/issn1684-8853.2017.5.12
РЕГУЛЯРИЗАЦИЯ ОБРАТНОЙ ЗАДАЧИ ЭЭГ/МЭГ ЛОКАЛЬНЫМ КОРТИКАЛЬНЫМ ВОЛНОВЫМ ПАТТЕРНОМ
А. А. Горшкова'6, аспирант
А. Е. Осадчий6'в, PhD
А. Л. Фрадков3'6, доктор техн. наук, профессор
аСанкт-Петербургский государственный университет, Санкт-Петербург, РФ
6Институт проблем машиноведения РАН, Санкт-Петербург, РФ
вЦентр нейроэкономики и когнитивных исследований, Национальный исследовательский университет
«Высшая школа экономики», Москва, РФ
Постановка проблемы: пространственное разрешение электроэнцефалографии/магнитоэнцефалографии зависит от метода решения обратной задачи, которая в силу фундаментальных физических причин является некорректно поставленной и имеет бесконечно большое количество решений. В последние несколько лет появились новые свидетельства о том, что нейрональная активность распространяется по коре в соответствии с волновым паттерном, характеризуемым некоторым направлением и скоростью распространения волны. Новые данные о пространственно-временной динамике распространения нейрональной активности по коре головного мозга требуют пересмотра существовавшей долгое время парадигмы, в которой пространственная структура активности рассматривалась независимо от временной динамики. Цель исследования: разработка нового метода локализации электрической активности головного мозга, который позволяет достаточно точно восстановить исходную активность по имеющимся данным на сенсорах в предположении, что эта активность имеет волновую структуру. Результаты: разработан новый математический аппарат регуляризации обратной задачи, ограничивающий решения множеством пространственно-временных динамик, удовлетворяющих двумерному волновому уравнению, определенному на нерегулярной сетке узлов, аппроксимирующих кортикальную поверхность. Это новый метод, реализованный в соответствии с хорошо себя зарекомендовавшей общей методологией регуляризации некорректно поставленных задач на основании минимизации Q-нормы решения. Сравнение на модельных данных с двумя наиболее распространенными методами решения обратной задачи показало, что новый метод, в отличие от них, сохраняет волновую структуру, обеспечивает наибольшую точность оценки моделируемой активности. Параметры регуляризующей волны рассчитываются в соответствии с минимизацией относительной невязки решения в пространстве сенсоров.
Ключевые слова — обратная задача электроэнцефалографии, волновое уравнение, регуляризация, разностный метод, управление.
Введение
Магнитоэнцефалография (МЭГ) и электроэнцефалография (ЭЭГ) — неинвазивные технологии регистрации электрической активности головного мозга, обладающие высоким временным разрешением и позволяющие изучать быстро протекающие нейрональные процессы [1], не нарушая целостность тканей [2, 3]. Пространственное разрешение этих технологий зависит от метода решения обратной задачи, которая в силу фундаментальных физических причин является некорректно поставленной (ill-posed) и имеет бесконечно большое количество решений [4, 5]. Для решения таких задач используют методы регуляризации, по сути соответствующие добавлению дополнительных свойств, которым должно обладать решение [6, 7]. Самым простым методом является тихоновская регуляризация, заключающаяся в поиске решения с минимальной энергией. Именно этот тип регуляризации лежит в основе широко распространенного на сегодняшний день метода MNE (Minimum Norm Estimates) [2] решения обратной задачи в ЭЭГ и МЭГ. Байесовская интерпретация этого подхо-
да показывает, что дополнительные требования к решению соответствуют априорному распределению вероятностей конфигурации нейрональ-ной активности [8]. Например, регуляризация по Тихонову соответствует предположению о независимости и одинаковой мощности активности во всех узлах сетки, аппроксимирующей кору головного мозга. Более реалистичные предположения о пространственной гладкости решения реализованы в методе LORETA (Low Resolution Electromagnetic Tomography) [9], в котором в качестве нормы используется мощность пространственно-высокочастотной компоненты, минимизация которой дает гладкие решения.
В последние несколько лет появились дополнительные свидетельства о наличии кортикальных волн разных масштабов [10-12]. Речь идет о том, что нейрональная активность распространяется по коре в соответствии с волновым паттерном, характеризуемым некоторым направлением и скоростью распространения волны [12, 13]. В соответствии с вышесказанным феноменологические волновые свойства распространения ней-рональной активности могут быть использованы в качестве априорных представлений об искомой
динамической активации. Принципиальным является тот факт, что при исследовании волновой структуры отсутствует возможность разделить временную и пространственную компоненты носителя, на котором определена регистрируемая в пространстве сенсоров и искомая в пространстве источников активность [13]. В настоящий момент не существует метода решения обратной задачи ЭЭГ/МЭГ на основе регуляризации в соответствии с априорно ожидаемой волновой пространственно-временной динамикой искомой кортикальной активности.
Новые данные о пространственно-временной динамике распространения нейрональной активности по коре головного мозга требуют пересмотра существовавшей долгое время парадигмы, в которой пространственная структура активности рассматривалась независимо от временной динамики [13]. Нами был разработан математический аппарат регуляризации обратной задачи, ограничивающий решения множеством пространственно-временных динамик, удовлетворяющих двумерному волновому уравнению, определенному на нерегулярной сетке узлов, аппроксимирующих кортикальную поверхность. Это новый метод, реализованный в соответствии с хорошо себя зарекомендовавшей общей методологией регуляризации некорректно поставленных задач, широко применяемой в различных областях науки, в том числе и при решении обратной задачи ЭЭГ/МЭГ [8]. Принципиально новым является именно использование волновой динамики для формирования априорного ограничения, накладываемого на искомое решение.
Обратная задача ЭЭГ и ее решение
Рассмотрим уравнение наблюдения для Т дискретных моментов времени:
Y(t) = Gx J(t) + e(t), 0 < t < T - 1,
(1)
где Y(t) e мени t; G
lNe — вектор измерений в момент вре-= RNe xNv — матрица прямой модели от
NV источников к Ne сенсорам, i-й столбец которой получается решением прямой задачи для заданной конфигурации сенсоров и диполем, помещенным в i-й источник на сетке коры; J(t) е RNv — неизвестный вектор активаций тока в источниках в момент времени t; e(t) е RNe — вектор измерений ошибок (шума).
Для распределенных источников матрица G необратима, так как Nv >> Ne, и задача является некорректно поставленной. Для ее решения приходится использовать метод регуляризации Тихонова:
min
J
н
GJ-Y +аJ
(2)
где
|GJ-YI]„ =(GJ-Y)i S^GJ-Y) — норма
исходной модели, т. е. насколько мы доверяем полученным данным; || = JTQJ — специально заданная ©-норма матрицы априорных предположений о модели, где Q — симметрическая матрица весов; а> 0 — параметр регуляризации.
Тогда решением задачи (2) является
J = (GtE;:1G + aQ)+ GTSg1Y,
(3)
где + означает псевдообратную матрицу. Формула (3) в дальнейшем в работе будет использована как решение обратной задачи.
Особую роль играет выбор матрицы весов Q априорных предположений о модели. Если Q = 1^, то метод решения по сути является методом наименьших квадратов ММЕ. Если Q задается весами оператора Лапласа, то это метод ЬОКЕТЛ, который порождает решения максимальной гладкости. Заметим, что в ЬОКЕТЛ гладкость учитывается только по пространственной координате. Суть нового метода заключается в рассмотрении такого ограничения, которое может навязывать решению пространственно-временной паттерн.
Параметр регуляризации а можно выбирать разными способами. В данной работе использовался метод Ь-еигув, суть которого заключается в следующем. Строится ^-^-график ©-нормы М© регуляризованного решения относительно соответствующей нормы ошибки уравнения , ,2
GJ-YI . Получающаяся кривая имеет форму
II ИЕ6
буквы Ь (рис. 1), на которой достаточно легко найти точку, минимизирующую обе величины. Таким образом, лучший выбор а соответствует угловой точке получаемой кривой в непрерывном случае и точке, ближайшей к угловой, — в дискретном.
50
40
30
|||J 20
ад
lo
10
0
-10
10
15 20 log||Y - GJ||2
25
30
Рис. 1. Кривая в методе L-curve Fig. 1. The curve in «L-curve» method
5
Предположение о волне
Новый метод заключается во введении пространственно-временного ограничения решения и основывается на предположении, что распространение электрической активности источника представляет собой некоторую волну с параметрами, описываемую волновым уравнением с диссипацией
d2u
dt
2 du
— = с Au - p—, * di
(4)
где c — скорость волны; p — параметр затухания.
Предполагается, что волна имеет неизвестное направление распространения 6 вдоль вытянутой сетки размером N х M и что начальное распределение локализовано рядом с центром сетки. Благодаря наличию временной характеристики в действительности рассматривается расширенное уравнение
Y =GJ + s,
где значения Y(t) и J(t) для t = 0, 1, ..., T - 1 объединены в один вектор:
(G 0 ... О ^ О G ... О
(
Y =
Y (0)
Y (T -1)
G =
0 0 ... G
J =
J(0) ï ( 8(0)
J (T -l)
; e =
E(T - il
Таким образом, задача решается сразу для всего отрезка времени:
min
j
1 /> /> . 2 1 ~ 1
GJ- YI + а J1
(5)
где матрица Q будет построена как специально заданная матрица весов волнового уравнения.
Моделирование волны. Структура матрицы Q
Для моделирования волны в двумерном пространстве использовался разностный метод, который широко распространен в численном моделировании [14]. Параметры к и т соответственно обозначают шаг равномерной сетки в двумерном пространстве и шаг по времени. Дискретный промежуток времени Т выбирался исходя из феноменологических свойств нейрональной активности. В качестве начального условия рассматривалась
4 2
0 100
100
20
0 0
Рис. 2. Начальное условие распространения волны Fig. 2. The initial condition of the wave propagation
гауссова «шапочка», локализованная в центре сетки, а граничные условия брались свободными (рис. 2), поскольку такая начальная активация наиболее соответствует свойствам изучаемой ней-рональной активности, представляющей собой кратковременно локализованное в пространстве колебание потенциала действия в форме пика.
Матрица Q строится «упаковкой» коэффициентов численного моделирования волны. Для этого сначала составляется линейная система уравнений на для каждого временного слоя:
AkUk = 0,
где Uk — вектор значений всех точек на сетке в момент времени Матрица Ak (размера ЫМ х ЫМ) имеет плохую обусловленность и блочный трех-диагональный вид:
( С I О I С I
Ak =
О ^ О
о о
I С I О I с
где матрицы ^ O, I имеют порядок Ы, причем O является нулевой матрицей, в матрице C число Кг есть число соседей ¿-го узла сетки (от 1 до 4):
с=4 h2
Г Ki -1
о о
-1 о
К2 -1
... -1
... о
Г-1 О
0 -1
о ...
о ...
^MN-1 -1
о ... о ...
-1
о
о о
-1
о ï о
о -1
Поскольку Ak не зависит от то переобозначим ее как A. Теперь составим линейную систему уравнений на все , т. е. объединим слои:
AU = 6,
где U — вектор значений всех точек за все время. Тогда матрица А (размера ТШЫ х ТШЫ) имеет следующий вид:
A =
' A" F О О О О Ï
F' А' F О О о
О F' А' F О о
О О F' А' F о
О О О F' А' о
О О О О О о
где F = \INM ; F'= I
NM ;
A' = A -
2 + pi
NM ;
,_ А 1 + рг т А ---2— "
Матрица А представляет собой априорную матрицу для пространственно-временной структуры волны, описываемой уравнением (4). Теперь определим симметрическую матрицу Q следующим образом:
Q:=ATA.
Таким образом, мы определяем ©-норму как норму пространственно-временной структуры волны по заданным параметрам сетки, скорости и затухания для основной задачи (5). Если подставить J в вытянутый вдоль времени вектор всей волны ^ то получим следующее выражение:
Шц = UTQU = ^А^ = 0,
которое доказывает, что в J Ц-норма достигает своего минимума.
Отображение активности на коре головного мозга на сигналы сенсоров. Применение прямой модели
Обратная задача решается на локализованном участке, т. е. заранее известен его центр (х0, у0, г0) — предположительное место активации. Для того чтобы решать обратную задачу на прямоугольной сетке, необходимо выполнить перенос на область вокруг этого центра и вычислить топографии токовых диполей, находящихся в узлах новой сетки (рис. 3).
0,16 0,12 0,08 0,04 0
Рис. 3. Сетка всей коры и прямоугольная сетка Fig. 3. The whole cortex grid and the rectangular grid
Один из способов построения преобразования заключается в рассмотрении двух наборов соответствующих между собой реперных точек в пространстве и задания аппроксимирующего отображения для них. Мы пользовались достаточно известным методом Approximating Thin-Plate Splines [15].
При переносе прямоугольной плоской сетки на кору для удобства был использован сглаженный вариант коры, так как он полностью изоморфен исходной коре, т. е. каждый узел на сглаженной коре можно перенести на соответствующий узел на исходной. Для заданного центра (х0, y0, z0) на сглаженной коре итеративный перенос заключается в следующем.
— Строится сетка рядом с центром ортогонально направлению фиксированного диполя (рис. 4).
— Итеративно используется трехмерный аналог алгоритма, в качестве реперных точек pt выступают узлы сдвигаемой регулярной сетки. Для них строятся проекции на кору и находятся ближайшие к ним узлы сетки коры. В качестве реперных точек q^ берутся середины векторов между p^ и найденными точками на сетке коры. Параметр гладкости X берется в основном в пределах от 0,1 до 10. Итерации происходят до достижения необ-
Рис. 4. Построение ортогональной сетки рядом с центром
Fig. 4. The construction of the orthogonal grid near the center
о> 4Î
Рис. 5. Поверхность коры вокруг центра и регулярная сетка после применения алгоритма и наложения на поверхность почти совпадают
Fig. 5. The cortex surface around the center and the regular grid after applying the algorithm and imposition on the surface almost match
ходимого порога ошибки аппроксимации репер-ных точек. На рис. 5 показан результат применения такого метода.
— Из прямой модели G для всей коры берутся только те источники, которые участвовали в качестве реперных точек на последней итерации, а значения на них интерполируются по значениям на ближайших узлах сетки.
Оценки параметров волны
Параметры исходной волны, такие как скорость с и направление распространения 9, можно оценить исходя из относительной ошибки
Ие* (9, с) = ||У - GJ(Q(0, с))|| / ||У||.
Сканирование параметров заключается в следующем: проверяются заданные заранее направления распространения волны, т. е. используются разные возможные углы 6 (рис. 6), а также для каждого направления проверяется набор скоростей от 0 до максимально возможной, так как численный метод моделирования накладывает условие устойчивости.
0,04 0,03 0,02 0,01 0
-0,03 -0,02
-0,01
0,01
Рис. 6. Направления сеток Fig. 6. Grids directions
Исходная скорость = 0,35355
-1,5 -1 -0,5 0 0,5 1
1,5
0,1 0,2 0,3 0,4 0,5 0,6 0,7 Скорость
0,04 0,03 0,02 0,01 0
Рис. 7. Кодированная цветом таблица отклонения от исходного угла и скорости
Fig. 7. The color coded table of the initial angle and speed deviation
Далее ищутся минимумы относительных ошибок полученных решений и по ним сужаются интервалы сканирования направлений и скорости. Процесс повторяется итеративно до необходимой точности нахождения параметров. В результате получается таблица значений относительной ошибки (рис. 7), по которой можно определить параметры исходной волны.
Результаты моделирования и анализ неизвестных параметров волны
В качестве модели головы как объемного проводника использовалась модель Т^опаЮТЕ (приблизительно 16 х 13 х 12 см) с 15 016 узлами (14 993 после сглаживания). Фиксировался 2010-й узел в качестве центра распространения волны.
Исходная волна J0 моделировалась на небольшой сетке 20 х 6 (приблизительно 5x1 см) в течение 40 отсчетов (при разрешении 598 Гц, что составило 70 мс) с использованием явной разностной схемы, при этом скорость волны равнялась половине от максимально допустимой, а именно ^^ (0,795 мм/мс), а коэффициент затухания р = 0,01. Затем применялся алгоритм переноса сетки на кору, при этом порог ошибки равен 1 х 10-5, параметр гладкости Х= 5. После 10 итераций на коре было получено 94 узла. Мы интерполировали в них значения со 120 узлов на сетке и применили к ним нужные элементы прямой модели G. Обратная задача решалась тремя способами: ММЕ, ЬОИЕТЛ и новым методом в соответствии с выражением (3). Параметр регуляризации а для всех методов мы искали, пользуясь Ь-еигуе методом, и получили следующие значения: ММЕ — а = 2,10 х 1021, ЬОИЕТЛ — а = 2,01 х 10-19, для нового метода а = 6,83 х 10-8.
0
■ Результаты моделирования тремя методами решения обратной задачи ЭЭГ по различным относительным ошибкам
■ The results of modelling by three methods of solution of the inverse EEG problem along different relative errors
Норма MNE LORETA Новый метод
Без шума 10 %-й шум Без шума 10 % -й шум Без шума 10 %-й шум
||Y - GJ||2 / ||Y||2 0,1911 0,3462 2,5199 х 10-8 0,2820 8,7113 х 10-12 0,2963
I|J||Q/||J||2 0,0020 0,0695 0,0001 26,8721 2,1139 х 10-24 2,2466 х 10-13
||Jq - JII2 / II JII2 1,8983 1,8915 0,0833 1,1039 5,5671 х 10-8 0,6017
Исходная активация источников
1,2 1
0,5 0,4 0,3 0,2 0,1 0
-0,01
0,2 0 -0,2 -0,4
Решение методом MNE Решение методом LORETA Решение новым методом
Рис. 8. Визуализация результатов моделирования тремя методами решения обратной задачи ЭЭГ: начальная активация; T = 1; шум 10%
Fig. 8. Visualization of the results of modelling by three methods of solution of the inverse EEG problem: initial activation; T = 1; noise 10%
Результаты моделирования с шумом и без него приведены в таблице, а также на рис. 8 и 9. В таблице показаны относительные погрешности уравнения, относительные ошибки ц-нормы и относительные погрешности решения для нового метода, ММЕ и ЬОИЕТЛ. В случае отсутствия шума имеется многократное улучшение точности оценки волны J0 по всем параметрам, при 10%-м шуме новый метод все равно дает более точные оценки волны J0, чем ММЕ и ЬОИЕТЛ. Однако при визуализации начальной активации источника Т = 1 становится понятно, что только новый метод показывает локализацию центра распространения волны J0, другие методы не справляются с этой задачей (см. рис. 8). Дальнейшая визуализация направления распространения волны также показывает,
что с этой задачей справляется только новый метод (см. рис. 9). Таким образом, при 10%-м шуме использование волнового уравнения в качестве ограничителя исходного решения обратной задачи ЭЭГ дает значительное преимущество по сравнению с известными методами при локализации электрической активности головного мозга и поиску ее направленности.
На рис. 7 и 10 показаны результаты сканирования параметров исходной волны — скорости с и направления распространения 9. Видно, что когда скорость становится больше искомой, относительная ошибка растет пропорционально модулю разности скоростей. Для углов наблюдается та же закономерность в любую сторону от искомого направления. Поэтому оценка параметра скорости происходит перебором возможных
ТЕОРЕТИЧЕСКАЯ И ПРИКЛАДНАЯ МАТЕМАТИКА
T = 1
1,2 1 0,8 0,6 0,4 0,2
T = 16
0,1 0 -0,1 -0,2
T = 24
Исходная активация источников
T = 40
0,05 0
-0,05
0,05
-0,05
T = 1
T = 16
T = 24
T = 40
■
— ■ -R
д- -i
| У
0,8
0,6
0,4
0,2
0
-0,2
-0,4
1 0,15
0,1
0,05
0
I -0,05
L
Решение новым методом
п 0,15
0,1
0,05
0
J -0,05
J ч»
У
0,06 0,04 0,02 0
0,02 0,04 0,06
Рис. 9. Визуализация результатов моделирования новым методом решения обратной задачи ЭЭГ: распространение волны; 1 < T < 40; шум 10 % Fig. 9. Visualization of the results of modelling by new method of solution of the inverse EEG problem: wave propagation; 1 < T < 40; noise 10 %
0
■ Рис. 10. Пространственное представление таблицы
отклонений. Минимум значений соответствует наилучшему выбору параметров скорости и направления волны
■ Fig. 10. The space representation of the deviation
table. Minimum value corresponds to the best choice of the wave and speed parameters
значений от максимальной к минимальной, а оценка параметра направления происходит по всем значениям.
Заключение
В работе описан новый математический аппарат регуляризации обратной задачи, ограничивающий решения множеством пространственно-временных динамик, удовлетворяющих дву-
мерному волновому уравнению. На модельных данных показано, что новый метод дает более точные оценки, чем два других наиболее распространенных метода. Также был разработан способ оценки латентных параметров мозговой активности — скорости с и направления распространения волны 6 — с помощью относительной ошибки решения в предположении, что волна распространяется только в конкретном направлении. Кроме того, был сделан анализ зависимости ошибок решения от различного уровня шума в прямой модели и в пространстве сенсоров, а также от различного размера накладываемой сетки.
Таким образом, разработанный новый метод локализации электрической активности головного мозга позволяет достаточно точно восстановить исходную активность по имеющимся данным на сенсорах в предположении, что эта активность имеет волновую структуру.
Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований, грант № 15-29-01344, и Санкт-Петербургского государственного университета, грант № 6.38.230.2015. Сравнительный анализ указанных в работе методов проведен в ИПМаш РАН при исключительной поддержке Российского научного фонда, грант № 14-2900142.
Литература
1. Jousmaki V. Tracking Functions of Cortical Networks on a Millisecond Timescale// Neural Networks.
2000. Vol. 13. P. 883-889.
2. Hanialaiiieii M., Hari R., Ilmoniemi R. J., Knuutila J., Lounasmaa O. V. Magnetoencephalography — Theory, Instrumentation, and Applications to Noninvasive Studies of the Working Human Brain// Rev. Mod. Phys. 1993. Vol. 65 (2). P. 413-497.
3. Nunez P. L. Electric Fields of the Brain. — N. Y.: Oxford, 1981. — 626 p.
4. Gloor P. Neuronal Generators and the Problem of Localization in Electroencephalography: Application of Volume Conductor Theory to Electroencephalogra-phy// J. Clin. Neurophysiol. 1985. Vol. 2. P. 327-354.
5. Sarvas J. Basic Mathematical and Electromagnetic Concepts of the Biomagnetic Inverse Problem// Phys. Med. Biol. 1987. Vol. 32. P. 11-22.
6. Baillet S., Mosher J. C., Leahy R. M. Electromagnetic Brain Mapping// Signal Processing Magazine. IEEE,
2001. Vol. 18(6). P. 14-30.
7. Jeffs B., Leahy R., Singh M. An Evaluation of Methods for Neuromagnetic Image Reconstruction// IEEE Trans. Biomed. Eng. 1987. Vol. 34. P. 713-723.
8. Grech R. Rewiew on Solving the Inverse Problem in EEG Source Analisys// Journal of NeuroEngineering and Rehabilitation. 2008. Vol. 5. N 25. doi:10.1186/ 1743-0003-5-25
9. Pascual-Marqui R. D. Low Resolution Electromagnetic Tomography: A New Method for Localizing
Electrical Activity in the Brain// International Journal of Psychophysiology. 1994. Vol. 18 (1). P. 49-65.
10. Patten T. M., Rennie C. J., Robinson P. A., Gong P. Human Cortical Traveling Waves: Dynamical Properties and Correlations with Responses// PLOS ONE. 2012. Vol. 7 (6). e38392.
11. Bahramisharif A., van Gerven M. A. J., Aarnout-se E. J., Mercier M. R., Schwartz T. H., Foxe J. J., Ramsey N. F., Jensen O. Propagating Neocortical Gamma Bursts are Coordinated by Traveling Alpha Waves// The Journal of Neuroscience. 2013. Vol. 33 (48). P. 18849-18854.
12. Alexander D. M., Jurica P., Trengove C., Nikolaev A. R., Gepshtein S., Zvyagintsev M., Mathiak K., Schulze-Bonhage A., Ruescher J., Ball T., van Leeuwen C. Traveling Waves and Trial Averaging: The Nature of Single-Trial and Averaged Brain Responses in Large-Scale Cortical Signals// Neuroimage. 2013. Vol. 73. P. 95-112.
13. Alexander D. M., Trengove C., van Leeuwen C. Don-ders is Dead: Cortical Traveling Waves and the Limits of Mental Chronometry in Cognitive Neuroscience // Cogn Process. Nov. 2015. Vol. 16 (4). P. 365-375.
14. Калиткин Н. Н. Численные методы. — СПб.: БХВ-Петербург, 2011. — 592 с.
15. Rohr K., Stiehl H. S. Landmark-Based Elastic Registration using Approximating Thin-Plate Splines// IEEE Transactions on Medical Imaging. June 2001. Vol. 20. N 6. P. 526-534.
UDC 681.3:612.8
doi:10.15217/issn1684-8853.2017.5.12
Regularization of EEG/MEG Inverse Problem with a Local Cortical Wave Pattern
Gorshkov A. A.a>b, Post-Graduate Student, andrey.a.gorshkov@gmail.com
Ossadtchi A. E.b>c, PhD, aossadtchi@hse.ru
Fradkov A. L.a,b, Dr. Sc., Tech., Professor, fradkov@mail.ru
aSaint-Petersburg State University, 7/9, Universitetskaya Emb., 199034, Saint-Petersburg, Russian Federation bInstitute of Problems of Mechanical Engineering of RAS, 61, Bol'shoi Pr. V. O., 199178, Saint-Petersburg, Russian Federation
cCentre for Cognition and Decision Making, National Research University Higher School of Economics, 46b, Volgogradsky Pr., 109316, Moscow, Russian Federation
Introduction: The space resolution of EEG/MEG depends on the inverse problem solution method. This problem is ill-posed because of fundamental physical issues and hence has an infinite number of solutions. Recent research has brought new evidence that cortical waves propagate through the cortex according to a wave pattern characterized by a certain direction and speed. New data about the space-time dynamics of neuronal activity propagation through the cortex require the revision of the long-standing paradigm in which a spatial activity structure is considered independently of the temporal dynamics. Purpose: We develop a new method for localizing electrical activity in the brain, which would allow us to fairly accurately restore the initial activity using the data from sensors under the assumption that this activity has a wave structure. Results: A new mathematical technique has been developed for the inverse problem regularization, which restricts the solutions by a set of space-time dynamics satisfying the two-dimensional wave equation defined on an irregular node grid approximating the cortical surface. This new method is implemented in correspondence with the tried and true methodology of general regularization of ill-posed problems on the base of minimizing the Q-norm of a solution. Model data comparison with two other common methods of inverse problem solution has shown that they do not retain the wave structure, whereas the new method does. It also provides the best estimation accuracy of the modeled activity. The regularization wave parameters are calculated in correspondence with the minimization of a relative error in the space of the sensors. Keywords — Inverse EEG Problem, Wave Equation, Regularization, Difference Method, Control.
References
1. Jousmaki V. Tracking Functions of Cortical Networks on a Millisecond Timescale. Neural Networks, 2000, vol. 13, pp. 883-889.
2. Hamalainen M., Hari R., Ilmoniemi R. J., Knuutila J., Lounasmaa O. V. Magnetoencephalography — Theory, Instrumentation, and Applications to Noninvasive Studies of the Working Human Brain. Rev. Mod. Phys, 1993, vol. 65(2), pp. 413-497.
3. Nunez P. L. Electric Fields of the Brain. New York, Oxford, 1981. 626 p.
4. Gloor P. Neuronal Generators and the Problem of Localization in Electroencephalography: Application of Volume Conductor Theory to Electroencephalography. J. Clin. Neuro-physiol, 1985, vol. 2, pp. 327-354.
5. Sarvas J. Basic Mathematical and Electromagnetic Concepts of the Biomagnetic Inverse Problem. Phys. Med. Biol., 1987, vol. 32, pp. 11-22.
6. Baillet S., Mosher J. C., Leahy R. M. Electromagnetic Brain Mapping. Signal Processing Magazine, IEEE, 2001, vol. 18(6), pp. 14-30.
7. Jeffs B., Leahy R., Singh M. An Evaluation of Methods for Neuromagnetic Image Reconstruction. IEEE Trans. Bi-omed. Eng., 1987, vol. 34, pp. 713-723.
8. Grech R. Rewiew on Solving the Inverse Problem in EEG Source Analisys. Journal of NeuroEngineering and Rehabilitation, 2008, vol. 5, no. 25, doi:10.1186/1743-00 03-5-25.
9. Pascual-Marqui R. D. Low Resolution Electromagnetic Tomography: A New Method for Localizing Electrical Activity
in the Brain. International Journal of Psy-chophysiology, 1994, vol. 18, pp. 49-65.
10. Patten T. M., Rennie C. J., Robinson P. A., Gong P. Human Cortical Traveling Waves: Dynamical Properties and Correlations with Responses. PLOS ONE, 2012, vol. 7(6), e38392.
11. Bahramisharif A., van Gerven M. A. J., Aarnoutse E. J., Mercier M. R., Schwartz T. H., Foxe J. J., Ramsey N. F., Jensen O. Propagating Neocortical Gamma Bursts are Coordinated by Traveling Alpha Waves. The Journal of Neuroscience, 2013, vol. 33(48), pp. 18849-18854.
12. Alexander D. M., Jurica P., Trengove C., Nikolaev A. R., Gepshtein S., Zvyagintsev M., Mathiak K., Schulze-Bonhage A., Ruescher J., Ball T., van Leeuwen C. Traveling Waves and Trial Averaging: The Nature of Single-Trial and Averaged Brain Responses in Large-Scale Cortical Signals. Neuroimage, 2013, vol. 73, pp. 95-112.
13. Alexander D. M., Trengove C., van Leeuwen C. Donders is Dead: Cortical Traveling Waves and the Limits of Mental Chronometry in Cognitive Neuroscience. Cogn Process, November 2015, vol. 16(4), pp. 365-375.
14. Kalitkin N. N. Chislennye metody [Numerical Methods]. Saint-Petersburg, BKhV-Peterburg Publ., 2011. 592 p. (In Russian).
15. Rohr K., Stiehl H. S. Landmark-Based Elastic Registration Using Approximating Thin-Plate Splines. IEEE Transactions on Medical Imaging, June 2001, vol. 20, no. 6, pp. 526534.
УВАЖАЕМЫЕ АВТОРЫ!
Научная электронная библиотека (НЭБ) продолжает работу по реализации проекта SCIENCE INDEX. После того как Вы зарегистрируетесь на сайте НЭБ (http://elibrary.ru/ defaultx.asp), будет создана Ваша личная страничка, содержание которой составят не только Ваши персональные данные, но и перечень всех Ваших печатных трудов, имеющихся в базе данных НЭБ, включая диссертации, патенты и тезисы к конференциям, а также сравнительные индексы цитирования: РИНЦ (Российский индекс научного цитирования), h (индекс Хирша) от Web of Science и h от Scopus. После создания базового варианта Вашей персональной страницы Вы получите код доступа, который позволит Вам редактировать информацию, помогая создавать максимально объективную картину Вашей научной активности и цитирования Ваших трудов.
