Регуляризация решения обратной задачи ЭЭГ и МЭГ на основе физиологически обусловленных моделей динамики нейрональной активности тема диссертации и автореферата по ВАК РФ 00.00.00, кандидат наук Кузнецова Александра Алексеевна

• Кузнецова Александра Алексеевна
• кандидат наук
• 2022

Оглавление диссертации кандидат наук Кузнецова Александра Алексеевна

Глава 1. Введение

1.1 Обратная задача ЭЭГ и МЭГ

1.2 Анализ межприступных разрядов пациентов с эпилепсией

1.3 Алгоритмы поиска межприступных разрядов в данных

1.4 Кортикальные бегущие волны

1.5 Анализ вызванных потенциалов

1.6 Цели работы

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

1.8 Научная новизна

1.9 Теоретическая и практическая значимость

1.10 Методология и методы исследования

1.11 Основные положения, выносимые на защиту

1.12 Вклад автора в проведенное исследование

1.13 Публикации и апробация работы

Глава 2. Содержание работы

2.1 Решение обратной задачи ЭЭГ/МЭГ в условиях наличия коррелирующих источников с помощью модифицированного бимформера

2.1.1 Основная идея предложенного метода

2.1.2 Ковариационная матрица в пространстве сенсоров

2.1.3 Метод Яее1Р811С08

2.1.4 Метод wReciPSIIC0S

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

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

2.2.1 Модель данных

2.2.2 «Базисные» волны

2.2.3 Оптимальная комбинация бегущих волн

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

2.3 Метод автоматического поиска ирритативных зон в МЭГ данных

пациентов с эпилепсией

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

2.3.2 Валидация найденных событий

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

Глава 3. Заключение

Список литературы

Приложение А. Статья «Modified covariance beamformer for solving

MEG inverse problem in the environment with correlated sources»

Приложение Б. Статья «Анализ локальной динамики

распространения межприступных разрядов с

помощью модели бегущих волн»

Приложение В. Статья «Data-driven approach for the delineation of

the irritative zone in epilepsy in MEG»

Приложение Г. Статья «Correlation of cue-locked FRN and

feedback-locked FRN in the auditory monetary

incentive delay task»

Приложение Д. Статья «Cortical plasticity elicited by acoustically cued

monetary losses: an ERP study»

Глава 1. Введение 1.1 Обратная задача ЭЭГ и МЭГ

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

Более точная диагностика и детальный анализ когнитивных процессов требуют применения методов решения обратной задачи ЭЭГ/МЭГ, которые позволяют по неинвазивным записям электромагнитной активности мозга оценить активность нейронных популяций. Пространственное разрешение МЭГ достаточно высоко, имеет порядок нескольких миллиметров, особенно в участках мозга с высокой кривизной [3]. Пространственное разрешение ЭЭГ ниже, чем у МЭГ, и составляет примерно сантиметр [4]. В силу фундаментальных физических ограничений обратная задача является некорректно поставленной (или недоопределенной) [5] и нуждается во ведении регуляризации для поиска единственного решения [6]. Даже после введения ограничений решение может быть неустойчивым: небольшие погрешности записанных экспериментальных данных могут приводить к значительному изменению решения. Таким образом, пространственное разрешение ЭЭГ/МЭГ в итоге в значительной мере зависит от выбора метода решения обратной задачи.

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

вы, но при этом череп обладает низкой проводимостью. Магнитное поле менее чувствительно к разнице в проводимости тканей [7]. Для решения прямой задачи ЭЭГ, таким образом, наиболее подходящим методом является метод граничного элемента (boundary element method, BEM) [8], который реалистично моделирует разные ткани. Для МЭГ также подходящим и менее вычислительно затратным может быть метод пересекающихся сфер (overlapping spheres) [9].

Пусть записанные ЭЭГ/МЭГ данные для каждого момента времени t представляют из себя вектор x(t)[Mх1], где здесь и далее в квадратных скобках указана размерность вектора (матрицы), и М — количество сенсоров. Пусть решение прямой задачи найдено одним из подходящих методов и хранится в операторе G[мxn], где N — количество источников в кортикальной модели.

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

Среди методов глобальной оптимизации — методов, стремящихся найти такое решение, в котором заданный функционал достигает своего глобального минимума, — мы выделяем метод минимальной нормы (minimum norm estimate, MNE) [10], метод взвешенной минимальной нормы (weighted minimum norm estimate, wMNE) [11], метод минимальных токов (minimum current estimate, MCE) [12], а также методы LORETA [13] и sLORETA [14], FOCUSS [15]. Во всех перечисленных методах целью является нахождение распределенного решения: имеется изначальная информация о расположении источников в кортикальной модели и для каждого из них по записанным ЭЭГ/МЭГ данным необходимо найти амплитуду активации. Задача является недоопределенной, так как количество источников намного больше, чем количество сенсоров: N ^ М.

Метод минимальной нормы, MNE, заключается в поиске такой оценки активации источников s(t)[Nх1], для которой, с одной стороны, разница между проекцией этого вектора на сенсоры и наблюдаемыми данными стремится к минимуму и, с другой стороны, квадратичная норма (L2 норма) вектора тоже стремится к минимуму:

||Gs(i) - x(i)||2 + a||s(i)||2 ^ min ,

s(i)

где || • ||2 — Евклидова норма L2, a — коэффициент регуляризации. Такая регуляризация позволяет получить более устойчивое решение, решая проблему корре-лированности столбцов матрицы прямого оператора G. Кроме того, мы находим самый короткий вектор, который позволяет получить достаточное качество решения с точки зрения квадратичной ошибки. Похожую оптимизационную задачу решает гребневая мультирегрессия [16].

Метод минимальных токов, MCE, по формулировке оптимизационной задачи похож на MNE, но вместо L2 регуляризации использует L\ регуляризацию:

||Gs(i) - x(i)||2 + a||s(i)|| ^ min ,

s (t)

где || • || — L\ норма. Благодаря введению негладкой L\ регуляризации, итоговое решение, полученное таким методом, будет разреженным, так как в результате оптимизации максимальное количество активаций источников будет приравнено к нулю. Одновременно с этим, негладкий регуляризатор делает невозможным нахождение аналитического решения, доступного в методе минимальной нормы, и поэтому решение методом MCE находится численно. Похожую оптимизационную задачу решает LASSO регрессия [17].

Однако чувствительность ЭЭГ/МЭГ сенсоров к источникам, расположенным в разных областях мозга, не одинакова. Электрическое и магнитное поля затухают по мере увеличения расстояния от источника до сенсора, и сенсоры максимально чувствительны к самым близким источникам на поверхности коры. При этом ЭЭГ более чувствительна к глубоким источникам, чем МЭГ [11]. Из-за этого ограничения решения методами MNE и MCE склонны быть смещенными в сторону источников на поверхности и недооценивать глубокие источники. Для того чтобы добавить важность глубинным источникам, используется метод взвешенной минимальной нормы, wMNE, в котором источники перевзвешиваются с соответствующими коэффициентами.

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

латься на метод адаптивного LCMV бимформера [18], [19]. Задача заключается в том, чтобы для каждого источника r = [xi,yi,Zi] найти такие коэффициенты пространственного фильтра bi[[Mх1], что дисперсия сигнала на выходе b^X минимальна, и при этом удовлетворяется ограничение на единичное усиление сигнала (unit gain) относительно интересующего источника. Оптимизационная задача записывается следующим образом:

min bf XXT Ьг = bf Cxb, bj

rji

subject to b gi = 1,

где Cx — ковариационная матрица в пространстве источников, g^ — столбец матрицы прямого оператора, соответствующий источнику i, (•)т — операция транспонирования. Постановку задачи можно переформулировать так: найти такой пространственный фильтр, который восстановит активность только от целевого источника, так, что её проекция обратно на сенсоры вернет исходные значения, и при этом найденный фильтр должен заглушить активность остальных источников. Из-за фундаментальных ограничений электромагнитной обратной задачи невозможно одновременно подавить вклад от всех нецелевых источников. Поэтому в рамках подхода бимформинга проблема нахождения весов пространственного фильтра формулируется локально [20].

Можно показать, что веса бимформера рассчитываются по следующей формуле:

bf = (gf Cx-1g,)-1gf Cx-1.

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

1.2 Анализ межприступных разрядов пациентов с эпилепсией

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

Эпилепсия — одно из самых распространенных неврологических заболеваний в мире, сопровождающееся для пациента не только наличием приступов, но также и риском возникновения сопутствующих болезней, когнитивных дефицитов, психологических расстройств и неблагоприятных социальных последствий. По данным Всемирной организации здравоохранения, в мире насчитывается более 50 миллионов человек с диагностированной эпилепсией1. Несмотря на то что для большей части пациентов приступы могут быть купированы с помощью правильной комбинации противоэпилептических препаратов, около 30% пациентов имеют фармакорезистентную форму эпилепсии, при которой медикаментозное лечение не позволяет контролировать приступы [21]. В таком случае пациенту может быть назначено нейрохирургическое вмешательство, и примерно в половине случаев такая операция позволяет полностью избавиться от приступов на срок не менее 10 лет, а в 85% случаев хирургическое вмешательство приводит к существенному сокращению частоты приступов и повышению качества жизни пациента [22], [23].

В распространенном случае мультифокальной эпилепсии патологическая активность зарождается в одной компактной части мозга, называемой эпилепто-генной зоной, из которой впоследствии распространяется на другие части мозга, часто затрагивая глубинные структуры, и проецируется на широкие области коры, вызывая у пациента приступ [24]. Локализация эпилептогеной зоны — ключевой практический вопрос при лечении фармакорезистентных форм эпилепсии, и эффективность хирургического вмешательства напрямую зависит от её качества. Нейрохирургическое вмешательство представляет собой удаление ткани в эпилептогенной зоне или рассечение нейронных связей для предотвращения распространения патологической активности. Как правило, для локализации эпилептогенной зоны используются неинвазивные ЭЭГ записи, инвазивная электрокортикограмма (ЭКоГ) или записи глубинных электродов, но все чаще,

1 https://www.who.int/news-room/fact-sheets/detail/epilepsy

и в случае доступности предпочтительнее, используется неинвазивная методика МЭГ [25].

Для определения эпилептогенной зоны различные области головного мозга исследуют на наличие межприступных разрядов — коротких, высокоамплитудных событий длиной 100-200 мс, значительно выделяющихся по сравнению с фоновой ЭЭГ/МЭГ активностью, и, как правило, порождаемых одним или несколькими фокальными источниками. Область мозга, генерирующая межпри-ступные разряды, называется ирритативной зоной, и таких зон у пациента может быть множество. Как правило, одна из ирритативных зон совпадает с эпилепто-генной зоной. Для повышения точности предоперационной диагностики эпилепсии необходима разработка новых методов и подходов для детального анализа межприступных разрядов, причем наиболее ценны неинвазивные методы диагностики, так как они сопряжены с меньшим риском для пациента.

В работе [26] авторы анализируют эпилептиформную активность, вызванную черепно-мозговой травмой у человека и крысы. Авторы отмечают, что неин-вазивные ЭЭГ записи оказались нечувствительны к патологической активности, тогда как инвазивные записи показали ее наличие у 86% пациентов. В то же время в работе [27] было показано, что межприступные МЭГ записи могут содержать значимую информацию, достаточную для правильной локализации эпилеп-тогенной зоны и последующего хирургического вмешательства. Часто проблема поиска эпилептогенной зоны решается с точки зрения анализа распространения активности между областями в масштабах всего головного мозга, однако также есть основания полагать, что эпилептиформная активность имеет распространение и на локальном уровне, но такие наблюдения основываются на инвазивных записях [28], [29] или на модельных данных [30].

1.3 Алгоритмы поиска межприступных разрядов в данных

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

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

Первым этапом анализа обычно является поиск межприступных разрядов в многоканальных МЭГ или ЭЭГ записях. В большинстве эпилептологических медицинских центров анализ межприступных данных производится экспертами вручную. ЭЭГ или МЭГ запись тщательно просматривается на наличие межпри-ступных разрядов, затем найденные события размечаются и, как правило, далее к сегменту данных, содержащих межприступный разряд, подгоняется диполь-ная модель. Наш опыт показывает [31], что для надёжного обнаружения ирри-тативных зон необходимо проанализировать большое количество разрядов. Так как ручная обработка большого объема данных приводит к усталости эксперта и росту вероятности совершения ошибки, как правило, анализ прекращается на минимальном количестве обработанных событий, субъективно считающимся достаточным. Зачастую в результате такого анализа эксперт рассматривает только один итоговый кластер диполей с морфологически устойчивыми временными рядами активаций. Кроме того, часто динамика межприступного разряда формируется одновременной активностью нескольких дипольных генераторов [27], и стандартная процедура подгонки одного диполя к разряду дает смещенную, а подчас и совершенно некорректную картину локализации.

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

Методы, основанные на морфологическом анализе [32], [33], [34] и использовании согласованного фильтра для поиска типичной формы межприступного разряда [35], [36], [37] сталкиваются с проблемой высокой изменчивости морфологии разряда даже для одного пациента, и, тем более, между разными пациентами [38]. Методы поиска по шаблону, активно развивающиеся на данный момент с развитием методов машинного обучения [39], [36], [37], предполагают

предварительную разметку значительного количества записи экспертом, поэтому такие методы нельзя назвать автоматическими и они не решают поставленную задачу.

Методы адаптивной фильтрации оказываются неэффективными и характеризуются низкой специфичностью [40], [41]. Ближе всего к разрабатываемому нами подходу находятся методы, основанные на анализе независимых компонент [42], [43], [44], но даже и в этом случае анализ не является автоматическим, так как выбор нужной компоненты и интерпретация результатов должны быть произведены экспертом.

1.4 Кортикальные бегущие волны

Другой важной концепцией для понимания одного из предложенных в данной работе методов является феномен кортикальных бегущих волн. Первые описания кортикальных бегущих волн — такой активации нейронов, пик которой меняет своё положение в пространстве, — появились ещё в 1930-х годах [45], [46]. Однако большинство современных когнитивных исследований и методов анализа опираются на парадигму Дондерса: представление о пространственно-временной разделенности мозговой активности [47]. С одной стороны, предполагая, что наблюдаемая активность мозга может быть представлена как сумма активаций статических источников, часто моделируемых в виде токовых диполей, можно использовать такие методы анализа как усреднение эпох для выделения интересующей нейронной активности или вычитание эпох в разных условиях. С другой стороны, такое представление не подходит для описания нейронной активности, распространяющейся в пространстве и обладающей высокой изменчивостью между эпохами, и применение перечисленных приемов приведет к уничтожению наблюдаемых явлений в данных [48].

Развитие технологий регистрации электромагнитной активности мозга и возможность многоканальной записи позволили показать, что волновой характер распространения нейронной активности присущ многим сенсорным, моторным и когнитивным системам мозга (для обзора см. [49]). Кортикальные бегущие волны могут появляться спонтанно или могут быть вызваны внешним стиму-

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

Несмотря на то что механизмы генерации кортикальных бегущих волн не до конца ясны [49], существующие работы демонстрируют, что их функции разнообразны. Ниже перечислены несколько примеров.

1. Кортикальные бегущие волны участвуют в координации работы удаленных друг от друга участков мозга. Согласно теории коммуникации через когерентность [51], на макроуровне функциональная интеграция областей мозга осуществляется за счет установления когерентных осцилляций в вовлеченных в процесс областях коры. Однако было показано, что альфа-осцилляции являются бегущими волнами [52], [53], и при этом фаза бегущей альфа-волны определяет мощность активации в гамма-диапазоне [54]. Кроме того, тета-активность также распространяется в виде бегущей волны в неокортексе человека [53], и в гиппокампе [55]. Дальнейшие исследования пространственно-временных характеристик распространения осцилляций позволят точнее понимать механизмы функционального взаимодействия в нейронных сетях головного мозга.

2. Кортикальные бегущие волны вовлечены в механизмы рабочей памяти. Было показано, что испытуемые лучше справлялись с заданием на рабочую память, когда кортикальные шаблоны распространения активности были согласованы [53].

3. Бегущие волны также участвуют в механизмах моторного контроля: скорость реакции в моторной задаче Go / NoGo статистически значимо связана с длительностью распространения волн альфа-активности в престимульном интервале [56].

4. Волны вовлечены в механизмы консолидации долгосрочной памяти во сне. Сонные веретена представляют из себя сферические волны, расходящиеся из одной точки, или радиальные волны с выраженным направлением [57]. Тета-волны во время медленного сна [58] и К-комплексы [59] также имеют сложные волновые шаблоны распространения.

5. Ряд работ показал, что есть основания полагать, что активность источников, порождающих межприступные разряды, также представляет из себя бегущие кортикальные волны [60], [61], [62] и распространяется по коре с определенной скоростью и в определенном направлении.

Остановимся подробнее на последнем пункте. Бегущие волны были описаны в бикукулинной модели эпилепсии [63]. В работе [60] приводятся результаты анализа инвазивных электроэнцефалографических записей приступной активности на разных масштабах: запись производилась одновременно с помощью сетки электродов размером 10x10 см и с помощью микро-сетки 4x4 мм. Авторам удалось показать, что в то время как приступная активность на макроуровне не демонстрировала никакого шаблона распространения, на микроуровне прослеживался отчетливый волновой шаблон. В качестве возможного объяснения авторы предположили, что во время приступа нейронные популяции на макроуровне координируются маленькими нейронными группами с помощью волнового механизма, тем самым еще раз подчеркивая важность анализа локального распространения эпилептической активности.

Другой группой исследователей [61] также в результате анализа инвазив-ных электроэнцефалографических данных было показано, что некоторые меж-приступные разряды демонстрируют устойчивые шаблоны распространения на сенсорах. Распространение активности происходит на участке порядка сантиметра, и такой масштаб может быть разрешим с помощью МЭГ [3]. Наиболее важным заключением из работы является наблюдение, что пациенты с успешным итогом операции демонстрировали статистически значимо более устойчивые траектории распространения межприступных разрядов. В последующей работе авторы подтвердили этот результат, также используя инвазивные данные и проводя анализ в пространстве сенсоров [62].

1.5 Анализ вызванных потенциалов

Одним из распространенных в литературе, благодаря своей простоте и эффективности, подходов к анализу неинвазивных ЭЭГ сигналов является анализ вызванных потенциалов, ВП (evoked potentials, EP). Вызванный потенциал —

усредненный ответ мозга на многократное предъявление стимула, в данном исследовании мы фокусируемся на предъявлении сенсорных слуховых стимулов. Как правило, вызванный потенциал невозможно увидеть при единичном предъявлении стимула, так как целевая активность обладает небольшой амплитудой и скрывается за спонтанной активностью мозга. Первые устойчивые исследования когнитивных ВП относят к 60-м годам, когда технология записи ЭЭГ, а также вычислительные возможности стали достаточными для того, чтобы провести запись и усреднение большого количества испытаний [64]. Аналоги вызванных потенциалов для МЭГ называются вызванными полями (evokedfields, EF).

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

Первым известным когнитивным компонентом ВП стал компонент условного отрицательного отклонения (contingent negative variation, CNV) [65]. В ходе эксперимента испытуемым сначала предъявлялся предупреждающий звуковой сигнал, за которым следовало предъявление визуального стимула. Затем эксперимент повторялся, но испытуемых просили нажимать на кнопку в ответ на предъявление визуального стимула. В этом случае исследователям удалось выделить негативный компонент с латентностью примерно в 500 мс после предъявления предупреждающего сигнала, который был связан с подготовкой к совершению целевого действия.

В данной работе мы используем для локализации компоненты слуховых вызванных потенциалов, поэтому остановимся на них подробнее. Один из первых компонентов N1 можно разделить на несколько частей: через 75 мс появляется первичный фронтально-центральный ответ, ассоциированный с первичной слуховой корой и далее максимума потенциал достигает примерно к 100 мс [66]. Другим важным компонентом является компонент P300, названный по латентно-сти (300 мс после стимула) и полярности (положительный) [67]. В исследовании было показано, что если испытуемый не может предположить, каким будет следующий предъявленный стимул: слуховым или визуальным, то стимул вызывает

высокоамплитудный положительный пик через 300 мс после предъявления. Если испытуемый ожидает получить предъявляемый стимул, то этот компонент имеет значительно меньшую амплитуду [64].

Более сложной экспериментальной парадигмы требует выделение компонента негативности рассогласования (mismatch negativity, MMN). Мы используем слуховую oddball парадигму: испытуемому предъявляется последовательность из повторяющихся звуков с редкими включениями девиантного звука. Предъявление девиантного стимула вызывает негативный ответ мозга, который на ЭЭГ сенсорах присутствует преимущественно в центральных отведениях. Ла-тентность компонента негативности рассогласования может варьироваться, как правило, от 160 до 220 мс после предъявления стимула. Данный компонент наблюдается даже в тех задачах, в которых девиантный стимул не связан для испытуемого с выполнением какого-либо задания. Таким образом, MMN компонент отражает некоторый автоматический процесс соотношения предъявляемого стимула и следа в памяти от предыдущего стимула [68].

Введение диссертации (часть автореферата) на тему «Регуляризация решения обратной задачи ЭЭГ и МЭГ на основе физиологически обусловленных моделей динамики нейрональной активности»

1.6 Цели работы

Как видно из обзора приведенного выше, как правило, выбор априорных предположений о характере активности источников головного мозга для регуляризации обратной задачи определяется не их физиологичностью, а, скорее, техническим удобством, позволяющим, например, получить аналитическое решение для обратного оператора (MNE, wMNE, LORETA), или наличием известного численного алгоритма оптимизации для итерационного поиска решения (MCE, FOCUSS). Целью данной работы является разработка подходов к решению обратной задачи ЭЭГ и МЭГ на основе именно физиологически обусловленных ограничений. В таком случае искомое решение получается более разумным, соответствует физиологической природе процесса и позволяет делать дальнейшие выводы относительно исследуемого явления.

Среди множества методов решения обратной задачи ЭЭГ и МЭГ адаптивные LCMV бимформеры [19], [18], [20] выделяются благодаря своему высокому пространственному разрешению, которое может быть достигнуто, когда наблю-

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

Несмотря на то что было доказано, что для целого ряда процессов головного мозга активность распространяется в форме движущихся по коре волн [49], большинство когнитивных исследований, и, соответственно, методов решения обратной задачи, опираются на предположение о том, что активность мозга может быть представлена как сумма активаций статических источников [48]. Существуют достаточные экспериментальные основания полагать, что локальное распространение межприступных разрядов у пациентов с эпилепсией может быть описано распространением кортикальной волны [60], [61], [62]. Во второй части исследования нашей целью была разработка алгоритма решения обратной задачи ЭЭГ/МЭГ, который использует физиологически обусловленное предположение о пространственно-временной связности искомой активности. Разрабатываемый алгоритм также должен восстанавливать характерные параметры волновой модели, а именно, направление и скорость движения волны.

Анализ МЭГ записей межприступной активности мозга пациентов с эпилепсией позволяет неинвазивно, с минимальным дискомфортом для пациента, с высокой точностью локализовать ирритативные зоны и исследовать пространственно-временную динамику эпилептогенной сети. Идентификация ирритатив-ных зон предполагает поиск и локализацию межприступных разрядов, содержащих ценную диагностическую информацию. Как правило, поиск эпилептиформ-ных событий производится с помощью визуального анализа записей экспертами. Учитывая то, что мы говорим о многоканальных данных с несколькими сотнями каналов, и очень фокальных событиях, такой визуальный анализ — это очень трудозатратная процедура, подверженная смещению со стороны эксперта. Так как ручная обработка большого объема данных приводит к усталости экспер-

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

Подводя итог, цели данной работы:

1. Разработка алгоритмов и методологии решения обратной задачи ЭЭГ/МЭГ для восстановления синхронной активности источников с высоким пространственным разрешением.

2. Исследование свойств предложенных методов в численных экспериментах, а также применение к реальным данным в слуховой парадигме. Анализ вызванных потенциалов в слуховых задачах.

3. Разработка алгоритма решения обратной задачи ЭЭГ/МЭГ, основанного на предположении о волновом распространении активности, исследование его свойств и применение разработанного метода к межприступным записям пациентов с эпилепсией.

4. Разработка алгоритмов автоматической обработки ЭЭГ/МЭГ в межпри-ступном периоде в целях обнаружения зон ирритации в коре головного мозга пациентов с эпилепсией.

5. Разработка методологии обнаружения эпилептогенной зоны на основе математического анализа динамики распространения эпилептической активности в межприступном периоде (алгоритм из пункта 2) и на основе анализа большого количества автоматически обнаруженных и кластеризованных разрядов (алгоритм из пункта 3).

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

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

В первой части исследования мы разработали два новых метода ReciPSIICOS и отбеленный wReciPSIICOS, которые являются модификациями классического адаптивного LCMV бимформера и, с одной стороны, сохраняют его высокое пространственное разрешение, но одновременно с этим, позволяют восстанавливать синхронную активность источников головного мозга. Свойства двух предложенных алгоритмов были исследованы сначала с помощью реалистично смоделированных данных, а затем и на реальных МЭГ данных, записанных в слуховых парадигмах. Мы сравнивали решения, полученные предложенными методами, с классическими подходами MNE и LCMV. Анализ свойств алгоритмов на реальных данных показал схожие результаты с анализом модельных данных. Наши эксперименты показали, что LCMV бимформер значительно чувствителен к присутствию коррелирующих источников в данных, решение подвержено проблеме подавления сигнала (signal cancellation problem), и часто содержит только один из синхронных источников, причем с размазанной картой активации. MNE часто не находит все активные источники или находит слишком смещенные и распределенные активации. При этом ReciPSIICOS и отбеленный wReciPSIICOS показывают высокое качество решения и позволяют находить фокальные билатеральные источники со значительно большим динамическим интервалом активаций.

В слуховых ЭЭГ экспериментах нам удалось показать, что в результате двух последовательных сессий задачи отложенного вознаграждения (monetary incentive delay, MID, в которой с помощью звуковых стимулов кодировались денежные потери, слуховой компонент MMN значительно вырос для тех сигналов, которые предсказывали большие денежные потери. Также мы показали, что FRN компонент модулировался как величиной, так и вероятностью результатов во время слуховой MID задачи, тогда как для dN200 такого эффекта обнаружено не было. Кроме того, dN200 компонент, который связан с обновлением информации о величине предполагаемого выигрыша, коррелирует со стандартным FRN, который связан с отрицательным RPE. Разработанные методы для восстановления активности коррелирующих источников с высоким разрешением позволят провести дальнейший анализ наблюдаемых феноменов.

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

ных волн. Мы представили событие межприступного разряда как суперпозицию предварительно сгенерированных с учетом индивидуальной анатомии бегущих волн. Мы использовали метод LASSO с положительными коэффициентами для того, чтобы оценить оптимальную скорость и направления распространения волн. Мы протестировали работу алгоритма как на модельных данных, так и на реальных МЭГ сигналах и продемонстрировали, что динамика распространения разрядов, записанных в МЭГ, может быть измерена в пространственно-временном масштабе миллиметр/миллисекунда. Несмотря на наличие источников ошибок, часть межприступных разрядов были успешно описаны с помощью волновой модели. Мы заметили, что у всех трех пациентов, данные которых были проанализированы, волновое поведение характерно не для всех межприступных разрядов, причем «волновые» разряды объединяются в хорошо пространственно очерченные кластеры. Более того, для пациентов, у которых были доступны данные об эпилептогенном очаге, эти кластеры совпадают с очагом. На основании этих результатов, которые хорошо согласуются с инвазивными данными [28], [29], мы предполагаем, что анализ межприступных разрядов, записанных в МЭГ, может помочь в локализации эпилептогенного очага.

В третьей части мы предложили метод автоматического поиска межпри-ступных разрядов и их кластеризации для определения ирритативной зоны, состоящий из двух этапов. На первом этапе мы используем разложение данные на независимые компоненты с помощью метода ICA (independent component analysis) и автоматически выбираем «спайковые» компоненты с помощью эвристик. Далее по порогу во временных рядах отобранных компонент мы определяем временные отсчеты, которые являются кандидатами на то, чтобы быть пиками межприступных разрядов. На втором этапе мы валидируем найденные события по результатам пространственно-временной кластеризации записей сенсоров вокруг найденных пиков методом сверточного разреженного кодирования (convolutional sparse coding). Мы валидировали результаты нашего анализа с помощью сравнения с ирритативными зонами, идентифицированными в результате визуальной инспекции, а также сравнения с зоной резекции вместе с известным исходом операции. Автоматически восстановленные ирритативные зоны пересекались с восстановленными по результатам визуального анализа у шести из семи пациентов. У оставшегося пациента и автоматически, и визуально восстановленные зоны, не пересекаются с зоной резекции и для этого пациента исход опера-

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

1.8 Научная новизна

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

Кроме того, мы впервые предложили алгоритм для решения обратной задачи ЭЭГ/МЭГ с использованием предположения о волновом распространении активности (волнового приора). В доступном сейчас программном обеспечении FieldTrip [69], MNE Python [70], Brainstorm [71] реализован широкий спектр методик для решения обратной задачи. В литературе также предлагаются подходы, отказывающиеся от моделирования нейрональной активности с помощью набора токовых диполей, например, в [72] в качестве решения рассматриваются сферические гармоники. Однако несмотря на то что феномен бегущих кортикальных волн за последнее десятилетие приобрел большую популярность и всё больше исследований демонстрируют разнообразие их функции в норме и патологии, ни в одном из предложенных сейчас решений не используется информация о пространственно-временной связности изучаемой активности.

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

1.9 Теоретическая и практическая значимость

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

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

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

1.10 Методология и методы исследования

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

1.11 Основные положения, выносимые на защиту

1. Разработаны два метода-модифицикации ЬСМУ бимформера (Яее1Р811С08, wReciPSIICOS), которые позволяют восстанавливать активность синхронных источников по неинвазивным ЭЭГ и МЭГ данным. Показано превосходство предложенных алгоритмов над классическими МКБ и ЬСМУ бимформером как на модельных, так и на реальных данных.

2. Разработан алгоритм решения обратной задачи ЭЭГ/МЭГ, основанный на предположении о волновом распространении активности. Его свойства исследованы на модельных данных. Полученные результаты для трёх пациентов демонстрируют наличие связи между качеством описания межприступного разряда с помощью волновой модели и принадлежностью разряда к эпилептогенной зоне.

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

1.12 Вклад автора в проведенное исследование

1. В первой части исследования автор: работала над развитием идеи алгоритма, которая изначально принадлежит научному руководителю А. Е. Осадчему, и комплементарна идее, представленной в диссертации Д. И. Алтухова; писала код для реализации всех алгоритмов на языке программирования Matlab2; полностью рассчитала результаты всех вычислительных экспериментов для исследования свойств алгоритма; провела анализ реальных данных для Набора данных 1; создавала визуализации результатов; внесла значимый вклад в написание текста публикации.

2. Во второй части исследования автор: работала над развитием идеи метода, которая принадлежит научному руководителю А. Е. Осадчему и в первой версии была реализована коллегой автора Е. А. Каленковичем, в процессе метод был значительно изменен, по сравнению с первоначальной версией; написала код для реализации метода на языке Matlab и частично на языке Python3; рассчитала результаты всех вычислительных экспериментов; рассчитала результаты на реальных данных пациентов; внесла преимущественный вклад в подготовку публикации.

3. В третьей части исследования автор: продолжая исследование научного руководителя А. Е. Осадчего, внесла значимый вклад в разработку Этапа 1 предложенного алгоритма; написала код на языке Matlab для Этапа 14; протестировала алгоритм на трёх пациентах из второй части исследования.

1.13 Публикации и апробация работы

Выступления на конференциях

2Код можно найти в репозитории https://github.com/kuznesashka/ReciPSIICOS

3 https://github.com/kuznesashka/wave_prior_inverse

4https://github.com/kuznesashka/ASPIRE

1. Международная конференция «The 8-th Mismatch Negativity Conference», тема «A novel beamformer immune to correlated sources and forward model inaccuracies», постерный доклад, Хельсинки, Финляндия, 2018.

2. Международная конференция «The 21-st International Conference on Biomagnetism BIOMAG», тема «MEG-based functional microscopy using traveling wave priors», постерный доклад, Филадельфия, США, 2018.

3. Международная конференция «Tubingen Systems Neuroscience Symposium», тема «MEG-based epilepsy diagnostics using traveling wave priors», постерный доклад, Тюбинген, Германия, 2018.

4. Международная конференция «MEG UK», тема «Traveling wave model for SOZ localization in MEG data», устный доклад, Кардифф, Великобритания, 2019.

5. Международная конференция «Organization for Human Brain Mapping», тема «MEG based functional microscopy using traveling wave priors: a new technology for exploring epilepsy», постерный доклад, Рим, Италия, 2019.

6. Международная конференция «BCI-Samara», тема «MEG based functional microscopy using traveling wave priors», постерный доклад, приз за лучший постер, Самара, Россия, 2019.

7. Международная конференция «Baltic Forum: Neuroscience, Artificial Intelligence and Complex Systems», тема «Local propagation of MEG interictal spikes: source reconstruction with traveling waves priors», устный доклад, Калининград, Россия, 2021.

Публикации повышенного уровня

1. Kuznetsova A., Nurislamova Y., Ossadtchi A. Modified covariance beamformer for solving MEG inverse problem in the environment with correlated sources. NeuroImage. 2021. 228. 117677. DOI: 10.1016/j.neuroimage.2020.117677.

2. Krugliakova E., Klucharev V., Fedele T., Gorin A., Kuznetsova A., Shestakova A. Correlation of cue-locked FRN and feedback-locked FRN in the auditory monetary incentive delay task. Experimental Brain Research. 2018. 23. 141-151.

3. Gorin A., Krugliakova E., Nikulin V., Kuznetsova A., Moiseeva V., Klucharev V., Shestakova A. Cortical plasticity elicited by acoustically cued monetary losses: an ERP study. Scientific Reports. 2020. 10. 21161.

4. Fedele T., Chirkov V., Kryuchkova A., Koptelova A., Stroganova T., Kuznetsova A., Kleeva D., Ossadtchi A. Data-driven approach for the delineation of the irritative zone in epilepsy in MEG. PLOS One. Новая версия рукописи с исправленными замечаниями рецензентов подана в редакцию.

5. Kuznetsova A., Fedosov N., Kleeva D., Ossadtchi A. Interictal traveling cortical waves and their implications for non-invasive MEG based detection of the seizure onset zone. NeuroImage. Рукопись подана в редакцию.

Публикации стандартного уровня

1. Кузнецова А. А., Осадчий А. Е. Анализ локальной динамики распространения межприступных рязрядов с помощью модели бегущих волн. Журнал высшей нервной деятельности. 2022. т. 72, No 3. с. 370—386. DOI: 10.31857/S0044467722030078

Глава 2. Содержание работы

2.1 Решение обратной задачи ЭЭГ/МЭГ в условиях наличия коррелирующих источников с помощью модифицированного бимформера

Пространственное разрешение ЭЭГ/МЭГ и результат восстановления активации источников по записанным сигналам сенсоров критически зависят от подхода, используемого для решения некорректно поставленной обратной задачи. В последние годы всё популярнее становятся решения с помощью адаптивных бимформеров [19], [18], [20]. В случаях активации небольшого количества некоррелированных друг с другом источников решение бимформера оптимально и обеспечивает высокое пространственное разрешение. Однако также известно, что бимформеры склонны ошибаться в тех случаях, когда активации источников коррелируют друг с другом: в результате итоговые временные ряды обладают низким отношением сигнал-шум (signal-to-noise ratio, SNR), а кортикальные карты распределений полученных активаций часто бессмысленны.

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

2.1.1 Основная идея предложенного метода

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

це ковариаций пространства сенсоров. Эта процедура проецирования является дополняющей к разработанной коллегами ранее в методе PSIICOS [74]. Изначально метод PSIICOS был разработан для анализа коннективности по МЭГ данным, в частности, для неинвазивного обнаружения взаимодействий между источниками с околонулевой фазовой задержкой. Проблема, которую решает PSIICOS, заключается в наличии в МЭГ данных артефактов объемной проводимости. Активации источников, которые на самом деле являются независимыми, на сенсорах могут проявить себя как коррелированная активность. Таким образом, цель метода заключается в том, чтобы отстроиться от эффектов объемной проводимости и оценить истинные корреляции источников. PSIICOS использует операцию проецирования, которая применяется к кросс-спектральной матрице пространства сенсоров, представленной как элемент М2-мерного векторного пространства. Было показано [74], что PSIICOS может достаточно хорошо разделить подпространство протечки сигнала и подпространство, содержащее вклад истинно связанных источников.

Однако наличие коррелирующих источников и является причиной проблемы, решаемой в данной работе: взаимодействия источников приводят к эффекту подавления сигнала (signal cancellation problem) при решении обратной задачи с помощью адаптивного бимформера. Решение устроено так, что в случае коррелированных источников можно подобрать для них взаимообратные коэффициенты пространственного фильтра, чтобы таким образом искусственно уменьшить целевой функционал. Применение аналогичного описанному выше подхода, основанного на проекциях, позволяет нам эффективно решать задачу оценивания временных рядов коррелированных источников с помощью адаптивного бим-формера. Мы используем комплементарную версию проекции PSIICOS для ковариационной матрицы сенсоров, чтобы вместо подавления вклада мощностей источников, подчеркнуть их, и, наоборот, уменьшить вклад коррелирующих источников. Такая проекция не удаляет активность коррелированных источников, а скорее выборочно обрабатывает их вклады в ковариационную матрицу и создает достаточно точную аппроксимацию идеальной ковариации данных, которую можно было бы гипотетически наблюдать, если бы эти источники были независимыми. Мы назвали новые методы ReciPSIICOS и отбеленный wReciPSIICOS, потому что предложенные алгоритмы решают проблему, комплементарную той, которую решает PSIICOS.

Список литературы диссертационного исследования кандидат наук Кузнецова Александра Алексеевна, 2022 год

Список литературы

1. Berger H. Über das Elektrenkephalogramm des Menschen // Archiv für Psychiatrie und Nervenkrankheiten. — 1929. — дек. — т. 87, № 1. — с. 527— 570. -DOI: 10.1007/bf01797193. — ÜRL: https://doi.org/10.1007/bfD1797193.

2. Cohen D. Magnetoencephalography: Detection of the Brain's Electrical Activity with a Superconducting Magnetometer // Science. — 1972. — февр. — т. 175, № 4022. — с. 664—666. — DOI: 10. 1126/science. 175.4022.664. — ÜRL: https://doi.org/10.1126/science.175.4022.664.

3. High-resolution retinotopic maps estimated with magnetoencephalography / K. Nasiotis [и др.] // Neurolmage. — 2017. — янв. — т. 145. — с. 107—117. — DOI: 10.1016/j.neuroimage.2016.10.017. — ÜRL: https://doi.org/10.1016/j. neuroimage.2016.10.017.

4. Michel C. M., Brunet D. EEG Source Imaging: A Practical Review of the Analysis Steps // Frontiers in Neurology. — 2019. — апр. — т. 10. — DOI: 10.3389/fneur.2019.00325. — ÜRL: https://doi.org/10.3389/fneur.2019.00325.

5. Helmholtz H. Üeber einige Gesetze der Vertheilung elektrischer Strome in körperlichen Leitern mit Anwendung auf die thierisch-elektrischen Versuche // Annalen der Physik und Chemie. — 1853. — т. 165, № 6. — с. 211—233. — DOI: 10.1002/andp. 18531650603. — ÜRL: https://doi.org/10.1002/andp. 18531650603.

6. Magnetoencephalography—theory, instrumentation, and applications to noninvasive studies of the working human brain / M. Hamalainen [и др.] // Reviews of Modern Physics. — 1993. — апр. — т. 65, № 2. — с. 413—497. — DOI: 10.1103/revmodphys.65.413. — ÜRL: https://doi.org/10.1103/revmodphys. 65.413.

7. A study of dipole localization accuracy for MEG and EEG using a human skull phantom / R. Leahy [и др.] // Electroencephalography and Clinical Neurophysiology. — 1998. — апр. — т. 107, № 2. — с. 159—173. — DOI: 10.1016/s0013-4694(98)00057-1. — ÜRL: https://doi.org/10.1016/s0013-4694(98)00057-1.

8. Hamalainen M., Sarvas /.Realistic conductivity geometry model of the human head for interpretation of neuromagnetic data // IEEE Transactions on Biomedical Engineering. — 1989. — февр. — т. 36, № 2. — с. 165—171. — DOI: 10.1109/10.16463. — URL: https://doi.org/10.1109/10.16463.

9. Huang M. X., Mosher J. C., Leahy R. M. A sensor-weighted overlapping-sphere head model and exhaustive head model comparison for MEG // Physics in Medicine and Biology. — 1999. — янв. — т. 44, № 2. — с. 423—440. — DOI: 10.1088/0031-9155/44/2/010. — URL: https://doi.org/10.1088/0031-9155/44/2/010.

10. Hamalainen M. S., Ilmoniemi R. /.Interpreting magnetic fields of the brain: minimum norm estimates // Medical &amp Biological Engineering &amp Computing. — 1994. — янв. — т. 32, № 1. — с. 35—42. — DOI: 10.1007/ bf02512476. — URL: https://doi.org/10.1007/bf02512476.

11. Assessing and improving the spatial accuracy in MEG source localization by depth-weighted minimum-norm estimates / F.-H. Lin [и др.] // NeuroImage. — 2006. — май. — т. 31, № 1. — с. 160—171. — DOI: 10.1016/j.neuroimage.2005. 11.054. — URL: https://doi.org/10.1016/j.neuroimage.2005.11.054.

12. Uutela K., Hamalainen M., Somersalo E. Visualization of Magnetoencephalographic Data Using Minimum Current Estimates // NeuroImage. — 1999. — авг. — т. 10, № 2. — с. 173—180. — DOI: 10.1006/nimg.1999.0454. — URL: https://doi.org/10.1006/nimg.1999.0454.

13. Pascual-Marqui R., Michel C., Lehmann D. Low resolution electromagnetic tomography: a new method for localizing electrical activity in the brain // International Journal of Psychophysiology. — 1994. — окт. — т. 18, № 1. — с. 49—65. — DOI: 10.1016/0167-8760(84)90014-x. — URL: https://doi.org/10. 1016/0167-8760(84)90014-x.

14. Pascual-Marqui R. Standardized low resolution brain electromagnetic tomography (sLORETA): technical details // Methods and findings in experimental and clinical pharmacology. — 2002.

15. Gorodnitsky I. F., George J. S., Rao B. D. Neuromagnetic source imaging with FOCUSS: a recursive weighted minimum norm algorithm // Electroencephalography and Clinical Neurophysiology. — 1995. — окт. — т.

95, № 4. - с. 231-251. - DOI: 10.1016/0013- 4694(95) 00107 - a. - URL: https://doi.org/10.1016/0013-4694(95)00107-a.

16. Hoerl A. E., Kennard R. W. Ridge Regression: Biased Estimation for Nonorthogonal Problems // Technometrics. — 1970. — февр. — т. 12, № 1. — с. 55—67. — DOI: 10.1080/00401706.1970.10488634. — URL: https://doi.org/ 10.1080/00401706.1970.10488634.

17. Tibshirani R. Regression Shrinkage and Selection Via the Lasso // Journal of the Royal Statistical Society: Series B (Methodological). — 1996. — янв. — т. 58, № 1. — с. 267—288. — DOI: 10.1111/j.2517-6161. 1996.tb02080.x. — URL: https://doi.org/10.1111/j.2517-6161.1996.tb02080.x.

18. Sekihara K., Nagarajan S. Adaptive Spatial Filters for Electromagnetic Brain Imaging. — Springer Berlin Heidelberg, 2008. — DOI: 10.1007/978-3-54079370-0. — URL: https://doi.org/10.1007/978-3-540-79370-0.

19. Localization of brain electrical activity via linearly constrained minimum variance spatial filtering / B. V. Veen [и др.] // IEEE Transactions on Biomedical Engineering. — 1997. — т. 44, № 9. — с. 867—880. — DOI: 10.1109/10.623056. — URL: https://doi.org/10.1109/10.623056.

20. Greenblatt R., Ossadtchi A., Pflieger M. Local linear estimators for the bioelectromagnetic inverse problem // IEEE Transactions on Signal Processing. — 2005. — сент. — т. 53, № 9. — с. 3403—3412. — DOI: 10.1109/ tsp.2005.853201. — URL: https://doi.org/10.1109/tsp.2005.853201.

21. Outcomes of surgical treatment of patients with pharmacoresistant epilepsy / V. V. Krylov [и др.] // Zhurnal nevrologii i psikhiatrii im. S.S. Korsakova. — 2016. — т. 116, 9. Vyp. 2. — с. 13. — DOI: 10.17116/jnevro20161169213-18. — URL: https://doi.org/10.17116/jnevro20161169213-18.

22. The long-term outcome of adult epilepsy surgery, patterns of seizure remission, and relapse: a cohort study / J. de Tisi [и др.] // The Lancet. — 2011. — окт. — т. 378, № 9800. — с. 1388—1395. — DOI: 10.1016/s0140-6736(11)60890-8. — URL: https://doi.org/10.1016/s0140-6736(11)60890-8.

23. Schuele S. U., Lüders H. O. Intractable epilepsy: management and therapeutic alternatives // The Lancet Neurology. — 2008. — июнь. — т. 7, № 6. — с. 514—

524. -DOI: 10.1016/s1474-4422(08)70108-x. — URL: https://doi.org/10.1016/ s1474-4422(08)70108-x.

24. Epilepsy / O. Devinsky [h gp.] // Nature Reviews Disease Primers. — 2018. — Man. — t. 4, № 1. — DOI: 10.1038/nrdp.2018.24. — URL: https://doi.org/10. 1038/nrdp.2018.24.

25. Ictal onset localization of epileptic seizures by magnetoencephalography / C. Tilz [h gp.] // Acta Neurologica Scandinavica. — 2002. — ceHT. — t. 106, № 4. — c. 190—195. — DOI: 10.1034/j.1600-0404.2002.02047.x. — URL: https: //doi.org/10.1034/j.1600-0404.2002.02047.x.

26. A Translational Study on Acute Traumatic Brain Injury: High Incidence of Epileptiform Activity on Human and Rat Electrocorticograms and Histological Correlates in Rats / I. G. Komoltsev [h gp.] // Brain Sciences. — 2020. — aBr. — t. 10, № 9. — c. 570. — DOI: 10. 3390/brainsci10090570. — URL: https://doi.org/10.3390/brainsci10090570.

27. Ictal and interictal MEG in pediatric patients with tuberous sclerosis and drug resistant epilepsy / A. Koptelova [h gp.] // Epilepsy Research. — 2018. — $eBp. — t. 140. — c. 162—165. — DOI: 10.1016/j.eplepsyres.2017.12.014. — URL: https://doi.org/10.1016/j~.eplepsyres.2017.12.014.

28. Microseizures and the spatiotemporal scales of human partial epilepsy / M. Stead [h gp.] // Brain. — 2010. — aBr. — t. 133, № 9. — c. 2789—2797. — DOI: 10.1093/brain/awq190. — URL: https://doi.org/10.1093/brain/awq190.

29. Millimeter-scale epileptiform spike patterns and their relationship to seizures / A. C. Chamberlain [h gp.] // 2011 Annual International Conference of the IEEE Engineering in Medicine and Biology Society. — IEEE, 08.2011. — DOI: 10. 1109/iembs.2011.6090174. — URL: https://doi.org/10.1109/iembs.2011. 6090174.

30. Minimal model of interictal and ictal discharges "Epileptor-2" / A. V. Chizhov [h gp.] // PLOS Computational Biology / nog peg. M. Bazhenov. — 2018. — Man. — t. 14, № 5. — e1006186. — DOI: 10.1371/journal.pcbi.1006186. — URL: https://doi.org/10.1371/journal.pcbi.1006186.

31. Automated interictal spike detection and source localization in MEG using ICA and spatial-temporal clustering / A. Ossadtchi [и др.] // Proceedings IEEE International Symposium on Biomedical Imaging. — IEEE. — DOI: 10.1109/ isbi.2002.1029375. — URL: https://doi.org/10.1109/isbi.2002.1029375.

32. Diambra L., Malta C. P. Nonlinear models for detecting epileptic spikes // Physical Review E. — 1999. — янв. — т. 59, № 1. — с. 929—937. — DOI: 10.1103/physreve.59.929. — URL: https://doi.org/10.1103/physreve.59.929.

33. Gotman J., Gloor P. Automatic recognition and quantification of interictal epileptic activity in the human scalp EEG // Electroencephalography and Clinical Neurophysiology. — 1976. — нояб. — т. 41, № 5. — с. 513—529. — DOI: 10.1016/0013-4694(76)90063-8. — URL: https://doi.org/10.1016/0013-4694(76)90063-8.

34. Gotman J.Automatic Detection of Seizures and Spikes // Journal of Clinical Neurophysiology. — 1999. — март. — т. 16, № 2. — с. 130—140. — DOI: 10. 1097/00004691-199903000-00005. — URL: https://doi.org/10.1097/00004691-199903000-00005.

35. Wilson S. B., Emerson R. Spike detection: a review and comparison of algorithms // Clinical Neurophysiology. — 2002. — дек. — т. 113, № 12. — с. 1873—1881. — DOI: 10.1016/s1388-2457(02)00297-3. — URL: https : //doi.org/10.1016/s1388-2457(02)00297-3.

36. Ko C.-W., Chung H.-W. Automatic spike detection via an artificial neural network using raw EEG data: effects of data preparation and implications in the limitations of online recognition // Clinical Neurophysiology. — 2000. — март. — т. 111, № 3. — с. 477—481. — DOI: 10.1016/s1388-2457(99)00284-9. — URL: https://doi.org/10.1016/s1388-2457(99)00284-9.

37. Automatic Detection of Epileptiform Events in EEG by a Three-Stage Procedure Based on Artificial Neural Networks / N. Acir [и др.] // IEEE Transactions on Biomedical Engineering. — 2005. — янв. — т. 52, № 1. — с. 30—40. — DOI: 10.1109/tbme.2004.839630. — URL: https://doi.org/10.1109/tbme.2004.839630.

38. Flanagan D., Agarwal R., Gotman J. Computer-aided Spatial Classification of Epileptic Spikes // Journal of Clinical Neurophysiology. — 2002. — март. — т. 19, № 2. — с. 125—135. — DOI: 10.1097/00004691-200203000-00003. — URL: https://doi.org/10.1097/00004691-200203000-00003.

39. Acir N., Guzeli§ C. Automatic spike detection in EEG by a two-stage procedure based on support vector machines // Computers in Biology and Medicine. — 2004. — окт. — т. 34, № 7. — с. 561—575. — DOI: 10.1016/j.compbiomed.2003. 08.003. — URL: https://doi.org/10.1016/j.compbiomed.2003.08.003.

40. Faure C. Attributed strings for recognition of epileptic transients in EEG // International Journal of Bio-Medical Computing. — 1985. — май. — т. 16, № 3/4. — с. 217—229. — DOI: 10.1016/0020-7101(85)90056-x. — URL: https: //doi.org/10.1016/0020-7101(85)90056-x.

41. Freeman W. /., Barrie J. M. Analysis of Spatial Patterns of Phase in Neocortical Gamma EEGs in Rabbit // Journal of Neurophysiology. — 2000. — сент. — т. 84, № 3. — с. 1266—1278. — DOI: 10.1152/jn.2000.84.3.1266. — URL: https://doi.org/10.1152/jn.2000.84.3.1266.

42. Isolation of epileptiform discharges from unaveraged EEG by independent component analysis / K. Kobayashi [и др.] // Clinical Neurophysiology. — 1999. — окт. — т. 110, № 10. — с. 1755—1763. — DOI: 10.1016/s1388-2457(99)00134-0. — URL: https://doi.org/10.1016/s1388-2457(99)00134-0.

43. Systematic source estimation of spikes by a combination of independent component analysis and RAP-MUSIC / K. Kobayashi [и др.] // Clinical Neurophysiology. — 2002. — май. — т. 113, № 5. — с. 713—724. — DOI: 10.1016/s1388-2457(02)00046-9. — URL: https://doi.org/10.1016/s1388-2457(02)00046-9.

44. Systematic source estimation of spikes by a combination of independent component analysis and RAP-MUSIC / K. Kobayashi [и др.] // Clinical Neurophysiology. — 2002. — май. — т. 113, № 5. — с. 725—734. — DOI: 10.1016/s1388-2457(02)00047-0. — URL: https://doi.org/10.1016/s1388-2457(02)00047-0.

45. Adrian E. D., Matthews B. H. C. The interpretation of potential waves in the cortex // The Journal of Physiology. — 1934. — июль. — т. 81, № 4. — с. 440— 471. — DOI: 10.1113/jphysiol.1934.sp003147. — URL: https://doi.org/10.1113/ jphysiol.1934.sp003147.

46. ADRIAN E. D, YAMAGIWA K. THE ORIGIN OF THE BERGER RHYTHM // Brain. — 1935. — т. 58, № 3. — с. 323—351. — DOI: 10.1093/brain/58.3.323. — URL: https://doi.org/10.1093/brain/58.3.323.

47. Donders F. On the speed of mental processes // Acta Psychologica. — 1969. — т. 30. — с. 412—431. — DOI: 10.1016/0001-6918(69)90065-1. — URL: https: //doi.org/10.1016/0001-6918(69)90065-1.

48. Alexander D. M., Trengove C., Leeuwen C. van. Donders is dead: cortical traveling waves and the limits of mental chronometry in cognitive neuroscience // Cognitive Processing. — 2015. — июль. — т. 16, № 4. — с. 365— 375. — DOI: 10.1007/s10339-015-0662-4. — URL: https://doi.org/10.1007/ s10339-015-0662-4.

49. Cortical travelling waves: mechanisms and computational principles / L. Muller [и др.] // Nature Reviews Neuroscience. — 2018. — март. — т. 19, № 5. — с. 255—268. — DOI: 10.1038/nrn.2018.20. — URL: https://doi.org/10.1038/nrn. 2018.20.

50. Spiral Wave Dynamics in Neocortex / X. Huang [и др.] // Neuron. — 2010. — дек. — т. 68, № 5. — с. 978—990. — DOI: 10.1016/j.neuron.2010.11.007. — URL: https://doi.org/10.1016/j.neuron.2010.11.007.

51. Fries P. Rhythms for Cognition: Communication through Coherence // Neuron. — 2015. — окт. — т. 88, № 1. — с. 220—235. — DOI: 10.1016/j. neuron.2015.09.034. — URL: https://doi.org/10.1016/j.neuron.2015.09.034.

52. Hindriks R., Putten M. J. van, Deco G. Intra-cortical propagation of EEG alpha oscillations // NeuroImage. — 2014. — дек. — т. 103. — с. 444—453. — DOI: 10.1016/j.neuroimage.2014.08.027. — URL: https://doi.org/10. 1016/j. neuroimage.2014.08.027.

53. Theta and Alpha Oscillations Are Traveling Waves in the Human Neocortex / H. Zhang [и др.] // Neuron. — 2018. — июнь. — т. 98, № 6. — 1269—1281.e4. — DOI: 10.1016/j.neuron.2018.05.019. — URL: https://doi.org/10.1016/j.neuron. 2018.05.019.

54. Propagating Neocortical Gamma Bursts Are Coordinated by Traveling Alpha Waves / A. Bahramisharif [и др.] // Journal of Neuroscience. — 2013. — нояб. — т. 33, № 48. — с. 18849—18854. — DOI: 10.1523/jneurosci.2455-13.2013. — URL: https://doi.org/10.1523/jneurosci.2455-13.2013.

55. Lubenov E. V., Siapas A. G. Hippocampal theta oscillations are travelling waves // Nature. - 2009. - май. - т. 459, № 7246. - с. 534-539. - DOI: 10.1038/nature08010. — URL: https://doi.org/10.1038/nature08010.

56. Human Cortical Traveling Waves: Dynamical Properties and Correlations with Responses / T. M. Patten [и др.] // PLoS ONE / под ред. P. A. Valdes-Sosa. -2012. - июнь. - т. 7, № 6. - e38392. - DOI: 10.1371/journal.pone.0038392. -URL: https://doi.org/10.1371/journal.pone.0038392.

57. Rotating waves during human sleep spindles organize global patterns of activity that repeat precisely through the night / L. Muller [и др.] // eLife. - 2016. -нояб. - т. 5. - DOI: 10.7554/elife.17267. - URL: https://doi.org/10.7554/elife. 17267.

58. Massimini M. The Sleep Slow Oscillation as a Traveling Wave // Journal of Neuroscience. - 2004. - авг. - т. 24, № 31. - с. 6862-6870. - DOI: 10. 1523/jneurosci.1318-04.2004. - URL: https://doi.org/10.1523/jneurosci.1318-04.2004.

59. Distribution, Amplitude, Incidence, Co-Occurrence, and Propagation of Human K-Complexes in Focal Transcortical Recordings / R. A. Mak-McCully [и др.] // eneuro. - 2015. - июль. - т. 2, № 4. - ENEURO.0028-15.2015. - DOI: 10.1523/eneuro.0028-15.2015. - URL: https://doi.org/10.1523/eneuro.0028-15.2015.

60. Human seizures couple across spatial scales through travelling wave dynamics / L.-E. Martinet [и др.] // Nature Communications. - 2017. - апр. - т. 8, № 1. -DOI: 10.1038/ncomms14896. - URL: https://doi.org/10.1038/ncomms14896.

61. Spatiotemporal Mapping of Interictal Spike Propagation: A Novel Methodology Applied to Pediatric Intracranial EEG Recordings / S. B. Tomlinson [и др.] // Frontiers in Neurology. - 2016. - дек. - т. 7. - DOI: 10.3389/fneur.2016. 00229. - URL: https://doi.org/10.3389/fneur.2016.00229.

62. Reproducibility of interictal spike propagation in children with refractory epilepsy / S. B. Tomlinson [и др.] // Epilepsia. - 2019. - апр. - т. 60, № 5. - с. 898-910. - DOI: 10.1111/epi. 14720. - URL: https://doi.org/10.1111/ epi.14720.

63. Huang XSpiral Waves in Disinhibited Mammalian Neocortex // Journal of Neuroscience. — 2004. — ho^6. — t. 24, № 44. — c. 9897—9902. — DOI: 10. 1523/jneurosci.2705-04.2004. — URL: https://doi.org/10.1523/jneurosci.2705-04.2004.

64. Luck S. J. An Introduction to the Event-Related Potential Technique. — 2005.

65. Walter W. Cooper R. A. V. Contingent Negative Variation: An Electric Sign of Sensori-Motor Association and Expectancy in the Human Brain // Nature. — 1964.

66. Naatanen R., Picton T. The N1 wave of the human electric and magnetic response to sound: a review and an analysis of the component structure // Psychophysiology. — 1987. — uronB. — t. 24, № 4. — c. 375—425. — DOI: 10.1111/j.1469-8986.1987.tb00311.x.

67. Evoked-potential correlates of stimulus uncertainty / S. Sutton [h gp.] // Science (New York, N.Y.) — 1965. — 26 ho*6. — t. 150, № 3700. — c. 1187—1188. — DOI: 10.1126/science.150.3700.1187.

68. The mismatch negativity (MMN): towards the optimal paradigm / R. Naatanen [h gp.] // Clinical Neurophysiology: Official Journal of the International Federation of Clinical Neurophysiology. — 2004. — shb. — t. 115, № 1. — c. 140—144. — DOI: 10.1016/j.clinph.2003.04.001.

69. FieldTrip: Open Source Software for Advanced Analysis of MEG, EEG, and Invasive Electrophysiological Data / R. Oostenveld [h gp.] // Computational Intelligence and Neuroscience. — 2011. — t. 2011. — c. 1—9. — DOI: 10.1155/ 2011/156869. — URL: https://doi.org/10.1155/2011/156869.

70. Gramfort A. MEG and EEG data analysis with MNE-Python // Frontiers in Neuroscience. — 2013. — t. 7. — DOI: 10.3389/fnins.2013.00267. — URL: https://doi.org/10.3389/fnins.2013.00267.

71. Brainstorm: A User-Friendly Application for MEG/EEG Analysis / F. Tadel [h gp.] // Computational Intelligence and Neuroscience. — 2011. — t. 2011. — c. 1— 13. — DOI: 10.1155/2011/879716. — URL: https://doi.org/10.1155/2011/879716.

72. Petrov Y. Harmony: EEG/MEG Linear Inverse Source Reconstruction in the Anatomical Basis of Spherical Harmonics // PLoS ONE / под ред. J.-C. Baron. — 2012. — окт. — т. 7, № 10. — e44439. — DOI: 10.1371/journal. pone.0044439. — URL: https://doi.org/10.1371/journal.pone.0044439.

73. Kuznetsova A., Nurislamova Y., Ossadtchi A. Modified covariance beamformer for solving MEG inverse problem in the environment with correlated sources // NeuroImage. — 2021. — март. — т. 228. — с. 117677. — DOI: 10.1016/j. neuroimage.2020.117677. — URL: https://doi.org/10.1016/j.neuroimage.2020. 117677.

74. Ossadtchi A., Altukhov D., Jerbi K. Phase shift invariant imaging of coherent sources (PSIICOS) from MEG data // NeuroImage. — 2018. — дек. — т. 183. — с. 950—971. — DOI: 10.1016/j. neuroimage. 2018. 08 . 031. — URL: https: //doi.org/10.1016/j.neuroimage.2018.08.031.

75. Cortical plasticity elicited by acoustically cued monetary losses: an ERP study / A. Gorin [et al.] // Scientific reports. - 2020. - Vol. 10, no. 21161.

76. Correlation of cue-locked FRN and feedback-locked FRN in the auditory mon" etary incentive delay task / E. Krugliakova [et al.] // Experimental Brain Re" search. - 2018. - Vol. 23. - P. 141-151.

77. Naatanen R., Ilmoniemi R. J., Alho K. Magnetoencephalography in studies of human cognitive brain function // Trends in Neurosciences. — 1994. — янв. — т. 17, № 9. — с. 389—395. — DOI: 10.1016/0166-2236(94)90048-5. — URL: https://doi.org/10.1016/0166-2236(94)90048-5.

78. Кузнецова А. А., Осадчий А. Е. Анализ локальной динамики распространения межприступных рязрядов с помощью модели бегущих волн // Журнал высшей нервной деятельности. — 2022. — т. 72, № 3. — с. 370—386. — DOI: 10.31857/S0044467722030078.

79. Mosher J., Leahy R. Source localization using recursively applied and projected (RAP) MUSIC // IEEE Transactions on Signal Processing. — 1999. — т. 47, № 2. — с. 332—340. — DOI: 10.1109/78.740118. — URL: https://doi.org/10. 1109/78.740118.

80. Data-driven approach for the delineation of the irritative zone in epilepsy in MEG / T. Fedele [et al.] // PLOS One. -.

81. Multivariate Convolutional Sparse Coding for Electromagnetic Brain Signals / T. D. L. Tour [и др.] // NeurIPS Proceedings. — 2018. — май. — URL: http: //arxiv.org/abs/1805.09654.

Приложение А

Статья «Modified covariance beamformer for solving MEG inverse problem in

the environment with correlated sources»

Authors: Aleksandra Kuznetsova, Yulia Nurislamova, Alexei Ossadtchi Abstract: Magnetoencephalography (MEG) is a neuroimaging method ideally suited for non-invasive studies of brain dynamics. MEG's spatial resolution critically depends on the approach used to solve the ill-posed inverse problem in order to transform sensor signals into cortical activation maps. Over recent years non-globally optimized solutions based on the use of adaptive beamformers (BF) gained popularity. When operating in the environment with a small number of uncorrelated sources the BFs perform optimally and yield high spatial resolution. However, the BFs are known to fail when dealing with correlated sources acting like poorly tuned spatial filters with low signal-to-noise ratio (SNR) of the output timeseries and often meaningless cortical maps of power distribution. This fact poses a serious limitation on the broader use of this promising technique especially since fundamental mechanisms of brain functioning, its inherent symmetry and task-based experimental paradigms result into a great deal of correlation in the activity of cortical sources. To cope with this problem, we developed a novel data covariance modification approach that allows for building beamformers that maintain high spatial resolution when operating in the environments with correlated sources. At the core of our method is a projection operation applied to the vectorized sensor-space covariance matrix. This projection does not remove the activity of the correlated sources from the sensor-space covariance matrix but rather selectively handles their contributions to the covariance matrix and creates a sufficiently accurate approximation of an ideal data covariance that could hypothetically be observed should these sources be uncorrelated. Since the projection operation is reciprocal to the PSIICOS method developed by us earlier (Ossadtchi et al., 2018) we refer to the family of algorithms presented here as ReciPSIICOS. We assess the performance of the novel approach using realistically simulated MEG data and show its superior performance in comparison to the classical BF approaches and well established MNE as a method immune to source synchrony by design. We have also applied our approach to the MEG datasets from the two experiments involving two different auditory tasks. The analysis of experimental MEG

datasets showed that beamformers from ReciPSIICOS family, but not the classical BF, discovered the expected bilateral focal sources in the primary auditory cortex and detected motor cortex activity associated with the audio-motor task. In most cases MNE managed well but as expected produced more spatially diffuse source distributions. Notably, ReciPSIICOS beamformers yielded cortical activity estimates with SNR several times higher than that obtained with the classical BF, which may indirectly indicate the severeness of the signal cancellation problem when applying classical beamformers to MEG signals generated by synchronous sources.

elsevier

Contents lists available at ScienceDirect

Neurolmage

journal homepage: www.elsevier.com/locate/neuroimage

Modified covariance beamformer for solving MEG inverse problem in the r environment with correlated sources c«r

Aleksandra Kuznetsova, Yulia Nurislamova, Alexei Ossadtchi*

Center for Bioelectric Interfaces, Higher School of Economics, Moscow 101000, Russia

ARTICLE INFO

Keywords: MEG

Inverse problem Beamforming Synchronous sources Correlated sources Data covariance Spatial resolution

ABSTRACT

Magnetoencephalography (MEG) is a neuroimaging method ideally suited for non-invasive studies of brain dynamics. MEG's spatial resolution critically depends on the approach used to solve the ill-posed inverse problem in order to transform sensor signals into cortical activation maps. Over recent years non-globally optimized solutions based on the use of adaptive beamformers (BF) gained popularity.

When operating in the environment with a small number of uncorrelated sources the BFs perform optimally and yield high spatial resolution. However, the BFs are known to fail when dealing with correlated sources acting like poorly tuned spatial filters with low signal-to-noise ratio (SNR) of the output timeseries and often meaningless cortical maps of power distribution.

This fact poses a serious limitation on the broader use of this promising technique especially since fundamental mechanisms of brain functioning, its inherent symmetry and task-based experimental paradigms result into a great deal of correlation in the activity of cortical sources. To cope with this problem, we developed a novel data covariance modification approach that allows for building beamformers that maintain high spatial resolution when operating in the environments with correlated sources.

At the core of our method is a projection operation applied to the vectorized sensor-space covariance matrix. This projection does not remove the activity of the correlated sources from the sensor-space covariance matrix but rather selectively handles their contributions to the covariance matrix and creates a sufficiently accurate approximation of an ideal data covariance that could hypothetically be observed should these sources be uncorrelated. Since the projection operation is reciprocal to the PSIICOS method developed by us earlier (Ossadtchi et al., 2018) we refer to the family of algorithms presented here as ReciPSIICOS.

We assess the performance of the novel approach using realistically simulated MEG data and show its superior performance in comparison to the classical BF approaches and well established MNE as a method immune to source synchrony by design. We have also applied our approach to the MEG datasets from the two experiments involving two different auditory tasks.

The analysis of experimental MEG datasets showed that beamformers from ReciPSIICOS family, but not the classical BF, discovered the expected bilateral focal sources in the primary auditory cortex and detected motor cortex activity associated with the audio-motor task. In most cases MNE managed well but as expected produced more spatially diffuse source distributions. Notably, ReciPSIICOS beamformers yielded cortical activity estimates with SNR several times higher than that obtained with the classical BF, which may indirectly indicate the severeness of the signal cancellation problem when applying classical beamformers to MEG signals generated by synchronous

1. Introduction

MagnetoencephalogTaphy (MEG) is a noninvasive neuroimaging method hallmarked by the millisecond scale temporal resolution and subcentimeter spatial resolution. As such, this method is very well suited

for studying the fine features of spatio-temporal dynamics exhibited by neural circuits. The high temporal resolution of MEG is concomitant to the nature of the underlying electrophysiological processes in the brain tissue. As to the spatial resolution, it crucially depends on how the ill-posed inverse problem is approached to recover the distribution of neu-

* Corresponding author.

E-mail addresses: [e-mail] (A. Kuznetsova), [e-mail] (Y. Nurislamova), [e-mail] (A. Ossadtchi). https://doi.Org/10.1016/j.neuroimage.2020.117677

Received 17 June 2020; Received in revised form 10 November 2020; Accepted 17 December 2020 Available online 29 December 2020

1053-8119/© 2020 The Author(s). Published by Elsevier Inc. This is an open access article under the CC BY-NC-ND license

(http://creativecommons.Org/licenses/by-nc-nd/4.0/)

ral sources from the magnetic field measurements around the head. Because the MEG inverse problem is inherently ill-posed no universal approach for solving it exists and each method makes specific assumptions about the properties of the neuronal activity.

One assumption common in the MEG literature is that the neuronal substrate that produces the observed data can be approximated by a small number of active focal sources typically represented by a set of equivalent current dipoles (Hamalainen and Hari, 1993). This multiple equivalent current dipoles (ECD) model originally implied the use of the least squares dipole fit to identify the location and orientation parameters of the dipoles (An et al., 2008). In practice the parametric ECD model can be used to represent spatially extended sources which was demonstrated in De Munck et al. (19B8) where dipole model error was studied for extended sources of the simulated visual evoked potentials (VEPs). The difficulties associated with this approach are the unknown number of sources and the need to solve a non-convex optimization task can be gracefully resolved with Multiple Signals Classification (MUSIC) family of methods, e.g. RAP-MUSIC (Mosher and Leahy, 1999). This method allows constructing a constellation of ECDs explaining the observed data produced by several simultaneously active (but not synchronous) sources. Recent extension called tRAP-MUSIC (Makela et al., 201B) furnishes an elegant approach that yields an Improved estimate on the number of true active sources present in the data. An additional attractive feature of MUSIC family of approaches Is that they do not invert data covariance matrix which avoids the need to employ regularization strategies that may bias the results. However, MUSIC methods can only find source location and are not designed for the source timeseries estimation task. Although the absence of data covariance inversion step make MUSIC approaches fairly robust against source timeseries correlations, applying MUSIC family methods to the data generated by synchronous sources (with source timeseries correlation close to 1) requires an explicit scan for synchronized cliques which may be time-consuming and impractical, especially when such cliques comprise more than two sources.

More than two decades ago local linear estimators were introduced as an alternative approach to estimating source distributions and time-series of specific neuronal sources from MEG and EEG data (Van Veen et al., 1997). This technique is based on the adaptive beamforming principle (Borgiotti and Kaplan, 1979) originally developed to detect sources in radar signals. In neuroimaging we employ adaptive linearly constrained minimum variance (LCMV) beamformers (Greenblatt et al., 2005; Sekihara et al., 2001; Vail Veen et al., 1997) both in the scanning mode and for estimation of neuronal source timeseries. To recover activity of a given source the LCMV beamformer finds a spatial filter by solving the optimization task to minimize the output power under the constraint of unit gain in the "direction" of the given source. Spatial distributions delivered by beamformers when operated in optimal conditions tend to be focal, which represents the "most interesting feature" of this approach as described in Borgiotti and Kaplan (1979). However, as a result of the "greedy" optimization linear dependencies between the neuronal source timeseries are utilised by the algorithm to minimize the output power (Sekihara and Nagarajan, 200B). This undesired phenomenon is called signal cancellation. Therefore applying the original LCMV beamforming approach we must assume that the measurements are generated by sources with uncorrelated timeseries.

In reality, in the neuroimaging applications this assumption is often violated which leads to suboptimal performance of the adaptive beam-formers. Thus, cortical sources exhibit transient synchrony (Varela et al., 2001) that manifests ongoing integrative processes (Fries, 2015). Another reason for such synchronization is the time-locking of brain activity (called event-related potentials, ERPs) to task events, such as movement or stimulus onset, for a wide range of cognitive, motor and sensory paradigms (Gascoyne et al., 2016). Synchronous ERPs often occur in both hemispheres at the functionally homologous areas. Bidirectional interaction of neuronal populations also leads to synchronization of their activity.

In this paper we propose a model-based extension of the LCMV beamforming that makes it robust against source synchrony and allows for reliable estimation of both locations and timeseries of synchronous sources. To explain the basic concept of this paper, we first illustrate the two opposite but related problems of (1) Estimating functional connectivity based on the indirect measurements and (2) Estimating source timeseries or source power distribution using adaptive beamforming. The problem of the former is the presence of artifacts of volume conduction: truly independent sources result into correlated sensor activity, The goal is to tune away from the volume conduction and estimate true source correlations. However, the presence of a non-trivial solution to the first problem, i.e. the presence of sources with true dependencies, is the problem of the latter: correlated sources lead to well known and unwanted signal cancellations in adaptive beamfoiming. As already stated, the original adaptive beamformer is optimal if all source activities are uncorrelated, i.e. a situation where any observed functional connectivity (measured as correlation) is purely an artifact of volume conduction. Correlated sources cause a significant reduction in the SNR when adaptive beamformers are used to process the data (Sekihara and Nagarajan, 200B).

Here we focus on the latter problem and develop an improvement for adaptive beamforming to make it robust against source correlations. The key component of an adaptive beamformer is the data covariance matrix that contains information about both the source-power distribution and about the linear source dependencies. The manifestation of linear dependencies present in the data covariance matrix is implicitly used by the LCMV beamformer optimal weights expression to further minimize beamformer output variance. Here we propose a procedure that uses data covariance generative model to selectively handle covariance matrix components to create a close approximation of an ideal data covariance matrix that would hypothetically be observed in the absence of correlation between source timeseries. Importantly, our co-variance modification procedure is data-independent. It relies only on the forward model and can be considered as a deterministic extra step in building an adaptive LCMV-based inverse operator.

This paper Is organized as follows. First we introduce the data model and describe the original adaptive beamforming approach (Van Veen et al., 1997). We then briefly illustrate both analytically and graphically the problem of correlated sources and emphasize the situation when more than two coupled sources are present. Next, having introduced the basic concepts and outlined the problem of correlated source timeseries we review the existing solutions and stress that most of them focus on the idea of nulling the sources identified as coupled to the target one. This approach may become impractical when the number of coupled sources is large and therefore next we introduce two covariance modification methods that allow us to suppress the manifestation of source correlations in the data covariance matrix without the need to identify the synchronized cliques. Using extensive simulations and real data from two paradigms we compare the performance of the adaptive beam-formers built using the original and the modified data covariance matrices. We also match the performance of the beamforming approaches against the minimum-norm estimation technique (Hamalainen and II-moniemi, 1994) that is (1) by design not affected by source timeseries correlation and (2) represents a solid and trusted reference that was also shown to manage well with localizing focal source in real data (Komssi et al., 2004). In the discussion we summarize our findings and oudine the shortcomings and potential directions to further advance the proposed approach.

2. Methods

2.1. Data model

We assume that event-related potentials (ERPs) measured with electrophysiological methods such as electroencephalography (EEG) and magnetoencephalography (MEG) can be represented as a superposition

of the contributions from a finite number of sources. Only the contributions that are sufficiently phase-locked to the task onset moment make it through the averaging procedure and form the evoked response (Luck and Kappenman, 2011). The observation equation linking vector s(r) of the averaged source-space activity with ERP's timeslice x(r) measured by the array of M sensors at time instance I can be written as

R

x(t) = X g,s,(0 + n(t) = Gs(f) + n(<). (1)

where x(i) is [M x 1] vector and s(() = [■.l(i:,.... iA(r)J7 is an [J? x 1 ] vector, as we assume, that this evoked activity is generated by a relatively small number R of focal cortical sources. Here g, Is the topography of the i-th equivalent current dlpole and G=[g1,....gJiJisa matrix of the oriented source topographies. Noise term n(r) represents the sum of the remainders of the induced and task unrelated activity that is supposed to be sufficiently suppressed by the event-related averaging procedure. For compactness in the presentation that follows we will stick to the fixed orientation source model. However, in Section 2.5 we show how to handle the unknown source orientations and apply the proposed methodology to the vector beamformer (Sekihara et al., 2001).

One of the key assumptions implicit in this ERP data model is that the stimulus-locked cortical activation profiles are expected to be focal. The number of active stimulus-related sources R is usually unknown, but it is assumed to be several orders of magnitude less than the total number of cortical sources. Note that in this paper we suggest an improvement of the standard beamformer without making additional assumptions about the spatial distribution of the sources as compared to beamformers themselves. Thus (1) represents basically the same source model as the one in (Van Veen et al., 1997).

It is customary to consider the MEG data in the space of virtual sensors obtained by calculating the projection coefficients onto a subset of principal directions of the forward model matrix. The formulation presented here accommodates this point of view without any conceptual changes. The value of M corresponding to the number of physical sensors will then reflect the number of virtual sensors corresponding to the number of principal directions capturing a specific percentage of variance in the original forward model. Correspondingly, the forward matrix G will be replaced by the product with the matrix containing the coordinates of principal components In the original space.

The LCMV beamforming approach extended here does not formally require the spatial whiteness of the additive noise term (Van Veen et al., 1997) and therefore we do not make assumptions regarding its spatial covariance structure. Whole data covariance is assumed to have full rank.

The active source locations are unknown, and finding them is the goal of the EEG and MEG inverse problem solving. We approach this problem with the knowledge of the forward model that matches every i-th location of the dipolar sources with a topography vector gt that contains the weights for the contribution of the i-th unit dipole to the sensor measurements. Our goal is then to identify the grid nodes containing active dipoles. Sufficient accuracy of computing g(- is a strong requirement to ensure adequate performance of the inverse solvers, including adaptive beamformers.

2.2. Adoptive LCMV beamformer

Adaptive linearly constrained minimum variance (LCMV) beam-former (Sekihara et a)., 2001; Van Veen et al., 1997) is a local linear estimator with a unit gain weight constraint (Greenblatt et a)., 2005). Over the recent years, this approach gained popularity as an efficient inverse solver of the MEG inverse problem (Darvas et al., 2004).

Due to the fundamental limitations of the electromagnetic inverse problem, it is impossible to globally suppress contribution from all nontarget sources. Therefore, within the beamforming approach, the problem of finding the spatial filter weights is formulated locally as mini-

mizing the spatial filter output variance under the unit gain constraint with respect to the source of Interest.

Beamformers can be used to estimate the activity of a specific ROI or in a scanning mode to assess the distribution of activity over the entire cortex. Various forms of beamformers exist that can be classified based on the source-space and sensor-space norms (Greenblatt et al., 2005). Unlike global estimators (MNE, wMNE, MCE, etc.), beamformers tuned to different cortical locations do not depend on each other and their summed output projected back to the source space is generally not supposed to be equal to the measured data X.

Good spatial resolving power is one of the most attractive features of the adaptive beamforming technique (Borglotti and Kaplan, 1979). It is achieved using the data covariance matrix that conveys the information about the subset of active sources to the beamformer. This information is then used by the beamformers to efficiently distribute the available degrees of freedom to suppress only these active sources.

2.2.1. Adaptive beamforming principle

Consider an elementary cortical dipolar source with free orientation at r, = [x„y,,z,]. To reconstruct activity via vector LCMV beamformer from the EEG or MEG data, one uses a spatial filter b, tuned to this dipolar source. These weights are calculated as the solution to the following optimization problem (Sekihara et al., 2001)

minimize b'Crb.

; a)

subject to b' g, = 1,

where (■)' is the transpose operator, C, = £{x(;)x(Orl is the sensorspace covariance matrix, matrix g, = g(r,) is a topography of the i-th source. Using Lagrange multipliers method, it can be shown that the optimal solution is

bf = [grc^i-'gic;1 <3)

The calculated spatial filters b, could then be used to reconstruct the source timeseries vector estimates as

3,(0 = »?*(/)■ (4)

It is also possible to employ the beamformer in the scanning mode and compute power distribution profile af = Var(r;) for the entire set of N cortical locations r„i= 1,...,N which can be done without the explicit computation of (Sekihara and Nagarajan, 2008):

of = b?'CIb, = [gfC;1g(r1. (5)

As follows from Eqs. (3) and (5), given the forward model, the co-variance matrix fully determines the beamformer weights and the output power of source estimates when applied to the data x(f).

The described approach does not introduce any assumptions on the number of active sources or their spatial distribution. However, the LCMV adaptive beamformer assumes that the measured neural activity is produced by a small number of focal cortical sources (Borgiotti and Kaplan, 1979).

2.2.2. Adaptive beamforming in the environment with correlated sources Despite the high localization efficiency and high spatial resolution,

adaptive LCMV beamformers suffer from the signal cancellation problem when operating in the suboptimal environment. The optimal spatial filters calculated In (3) reflect activity of intrinsic sources correctly only when the measurements are generated by a small number of focal sources with uncorrelated timeseries.

These assumptions, however, are rarely met in the real experimental conditions because the brain integrative mechanisms result In a high level of synchrony across different areas (Fries, 2015). Thus, synchronous, stimulus-locked responses occur in many cortical areas, such as, for instance, ERPs occurring in bilateral, functionally homologous locations (Gascoyne et al., 2016).

Such correlation of sources causes significant reduction in the SNR of sources timeseries estimated with adaptive beamformers. For example,

-1 -0.5 0 0.51

Fig. 1. LCMV beamformer reconstruction in the case of three active sources with pairwise correlations p^r p1^ and = 0.9. The estimated variance ah i = 1,2,3 for each of the three sources as a function of pairwise correlation p"2, /r^ is color-coded into the heat-maps. In cach plot, the axis reflect coupling of the fiist source, r to the other two sources r and r. quantified by the correlation coefficients and p"y The impossible combinations of correlation coefficient values given unit variance sources are shown in white.

it can be shown (Sekihara and Nagarajan, 2008), that in the case of two interacting sources the estimated power of each source decreases as a quadratic function of the correlation coefficient of their timeseries and in the case of unit source variance, the output power can be expressed

8f = l-(/i")2, ( = 1,2

(6)

As follows from (6) in the case of complete synchrony, the adaptive beamformer will simply output zero. Intuitively this can be understood as follows. In order to meet the constraint, the beamformer provides unit gain for the target source. The functional optimized by the adaptive beamformer requires minimization of the output power. In the presence of another source with correlated activity, the beamformer adjusts the weight vector in such a way that, on the one hand, the gain of unity constraint is met, and on the other hand, the activity of the correlated sources is subtracted from the target activity to minimize the output power. Therefore, in the case of perfect correlation, the beamformer produces a zero SNR with respect to the activity of the target source.

The situation is aggravated in the environment with larger number of correlated source. Even moderate correlation of source timeseries causes significant reduction in the SNR of adaptive beamformer based estimates. Consider an environment with three active unit variance sources r,, r2, r3 with activation timeseries characterized by pairwise correlation coefficients

It is relatively straightforward to show, that in this case the adaptive LCMV estimated source power S,, ( = 1,2,3 depends on the pairwise correlation coefficients as

F F F

(pS) -1

where

where F = (p^)

W -1

(0ÏÏ)

■ 1

(7)

Panels A and B of Figure 1 show color-coded estimated power of the three sources r3, r2 and r3 for various degrees of coupling between the first source and the other two sources when r2 and r3 are strongly coupled, to = 0.9. As expected, the output power of the second and the third sources is reduced primarily due to their strong mutual coupling, Fig. l.B. In this case low SNR in the estimates of r2 and r3 complicates their detection and makes it problematic to implement the strategies suggested in Dalai et al. (2006); Popescu et al. (200B) to ameliorate the problems caused by source correlation. Additionally, we also observe a rapid reduction in the estimate of rj power with the growth of coupling between this first source and either of the two remaining sources. Fig. l.A.

In real life setting this situation may occur, for example, in an auditory-motor experimental paradigm, where the subject is required to

perform a motor action (e.g. press the button) In response to a deviant auditory stimulus. According to our indexing scheme source rj models sensory-motor response and sources r, and r3 are located bilaterally in the primary auditory cortex and respond synchronously to the auditory stimulus.

Another likely scenario is when the source distribution is represented in the form of two groups of sources with significant synchrony within each such group but no correlation between the groups. In this case the within group synchrony will affect beamformer's performance and will result in timeseries with underestimated amplitudes and erroneous source distributions.

2.2.3. Existing solutions

Several approaches have been developed to improve beamforming in the presence of correlated sources. (Dalai et al., 2006) suggested that the entire region that may potentially include a source correlated to the activity in the region of Interest (ROI) should be suppressed. This idea can be implemented using an SVD derived constraint based on the topography of the cortical patch containing the interfering source. This approach requires an apriori knowledge of the locations of correlated sources. Nulling the activity of multiple regions (or spatially-extended ones) with this method reduces the number of degrees of freedom, that otherwise could be used for suppressing the interfering sources.

Brookes et al. (2007) suggested building a beamformer based on the constraints that are calculated from the topographies of correlated sources using their linear combination. An "amplitude optimization" routine was suggested to compute the optimal mixing coefficients for this procedure. Application of this method to real data requires explicit scanning over all possible pairs using a coarse grid, which is not time efficient and prone to errors if the seed source is set opriori.

Two beamformers that allow to overcome correlated sources issue were evaluated in Popescu et al. (2008): (1) a linearly constrained minimum variance beamformer with partial sensor coverage (LCMV-PSC), and (2) a multiple constrained minimum-variance beamformer with coherent source region suppression (MCMV-CSRS). It was demonstrated that the latter exhibits precise localization and minimal amplitude and phase distortion for a broad range of relative positions of the interfering source within the suppression region. With this method, again, degrees of freedom are consumed because of the assumption regarding the location of the interfering source that maintains the regional zero constraint.

Quraan and Cheyne (2010) compared various solutions available to the date of that publication to cope with correlated sources in beamforming. When prior information about the location of correlated sources is available, the method of beamforming with coherent source region suppression described in Dalai et al. (2006) appeared to be the most effective, including the case of closely located (3 cm apart) correlated sources. The authors concluded that this solution, when carefully exer-

cised, can significantly improve localization accuracy, but does not fully solve the amplitude bias problem.

Diwakar et al. (2011) introduced a dual-core beamforming idea, which is an extension of Sekihara's vectorized LCMV approach (Sekihara et al., 2001). Dual-core beamformer is built using the constraint created from the topographies of two spatially disjoint regions. This approach allows for finding the pairs of highly correlated sources and eliminates the computationally expensive search for topographies mixing parameter and optimal dipole orientation required in the approach of Brookes et al. (2007).

The approaches based on finding pairs of correlated sources would, for example, fall to detect a hub coupled simultaneously to more than one additional source. Further, as our simulations show, synchrony between more than two sources has a complex effect on the suppression of the beamformer output power. Therefore, in a number of practical situations, the beamformers limited to considering only a pair of sources could fail. Moiseev et al. (2011) presented a detailed treatment of the multiple constrained minimum-variance beamformers with coherent source region suppression (MCMV-CSRS) and offered a set of practical solutions and scanning statistics to be used for identifying cortical regions with correlated activity.

The approaches for coping with source synchrony problem that we have described so far are conceptually similar and spin around the idea of suppressing the activity of the sources correlated to the current ROI. Additional insights into the problem of adaptive beamfoming in the environment with correlated sources could be gained from the data covariance matrix, which contains the information about source synchrony. Formulation of the minimum variance adaptive beamforming optimization problem implies the absence of correlation of the underlying sources which leads to a specific structure of the data covariance matrix. When this structure is violated the adaptive beamforming algorithm results in significant suppression of the SNR of the estimated source timeseries, and in the case of perfect synchrony the adaptive beamformer would completely cancel the signal.

The approach we propose here is based on the analysis of the structure of the data covariance matrix. Kimura et ai. (2007) previously developed a method that is conceptually close to ours. They used a forward model and the least squares approximation approach to find an estimate of the source-space covariance matrix that corresponded to the interaction of a small number of sources. They then nulled the off-diagonal elements of this matrix and projected it back to the sensor space. This new decorrelated matrix was used for building the conventional beam-formers. The performance of this approach depends on the estimated number of active sources chosen to comprise the source-space covariance matrix. The selection of the number of sources to be represented in the source-space covariance matrix is a non-trivial problem that affects the performance of this technique.

Here we propose an alternative way to modify data covariance matrix that makes adaptive LCMV beamforming robust against correlation of neural sources activity but does not require the full blown source modelling. In contrast to the approach described in Kimura et al. (2007), our covariance modification procedure does depend on the data. The proposed covariance modification procedure creates a close approximation of an ideal data covariance that could hypothetically be observed in the absence of correlation of source timeseries. It relies only on the whole head forward model and can be considered as a deterministic extra step in building an adaptive LCMV-based inverse operator.

In our approach, we consider the sensor-space data covariance matrix Cx as an element of Ai2-dimensional vector space and derive its expression as a function of active source topographies scaled by intrinsic source timeseries variance and covariance values. We then employ a data independent projection procedure to mitigate the contribution explained by coupling of intrinsic sources. This way the modified covariance matrix Cx approximates the one that could have been obtained should the same but non-Interacting constellation of sources is mea-

sured. We then apply the standard beamforming approach but use this modified data covariance matrix to compute the beamformer weights.

2.3. ReciPSIICOS beamforming

In this section we introduce the new beamformer modification, which is immune to the contributions from correlated sources in the data. The key component of the proposed algorithm is the projection procedure complementary to the PSIICOS (Ossadtchi et al., 201B) introduced by us earlier.

We originally designed PSIICOS projection technique for connectivity analysis, particularly for non-invasive detection of true zero-phase interactions between sources which addressed problem (1) Estimating functional connectivity outlined in the introduction. The procedure employs a projection operation applied to the sensor space cross-spectral matrix treated as an element of MI-dimensional vector space. We showed that PSIICOS could sufficiently well disentangle the subspace containing spatial leakage power from the subspace containing the contributions from the true zero-phase coupling of sources having different locations.

Here we address problem (2) Estimating source timeseries or source power distribution using adaptive beamforming and use a similar projection based approach for solving this complementary problem. Instead of suppressing the variance in the subspace filled with source power, we use the complementary version of PSIICOS projection to emphasize it and instead suppress the contributions to the sensor data covariance matrix the spatial components modulated by the source-space timeseries covariance. We refer to the new method as ReciPSIICOS because the proposed pipeline solves the problem reciprocal to that solved by PSIICOS.

2.3.1. General idea and vectorhed sensor-space covariance matrix

MxM sensor-space covariance matrix C, = £{x(i)xT(0} plays a pivotal role in adaptive beamforming. Vectorized covariance uec(Cx) of the evoked response can be expressed in terms of the elements c**, i,j = 1,.... R of the source-space covariance matrix of R active sources and their topographies g(- as

R

vec(Cx) = vec(E(xV)xT(_!)}) = £ i>ec(g,B?>"

i-I

R R

+ 2 S ""te^ + gjgf)tf + vec(C„) (8)

j-i+i

where C„ is the noise covariance matrix.

Eq. (8) demonstrates that matrix Cx can be decomposed into two types of additive terms: auto-terms corresponding to the source power Sj=i "ectgyg^e" and also the palrwise cross-products of source topographies weighted with source covariance coefficients J,*, ZjL,+i(g,gJ' + gjgf )cjj. These cross-terms are present in the covariance matrix due to the non-zero off-diagonal elements of the source-space covariance matrix.

The very existence of the cross-terms in the source-space covariance matrix leads to the undesired performance of the adaptive beamformer. This happens because the weight vectors are formed in such a way that the correlated sources are utilized to minimize the beamformer output power.

To dwindle the contribution of source covariance terms to the sensor-space covariance matrix we propose to use a projector operating in the M2-dimensional space. This processing pipeline is shown in Fig. 2. We apply this projector P to the vectorized data covariance matrix to weaken the contribution of the cross-terms. The resulting postprojection matrix when reshaped back approximates the sensor-space covariance that could have been obtained if no coupling was present. The obtained covariance matrix is then used in designing the LCMV beamformer according to the standard approach.

Next we propose two different approaches to building the M1 x M2 projection matrix P dubbed ReciPSIICOS and whitened ReciPSIICOS.

1. Vectorize sensor-space covariance matrix

c -

C11 C12

CA/l CM 2

C1M

cmm

Cll Cl2

cim

cmm

[a/2 x 1]

X "ecf&gj' + BjsF)<ij

2. Project and keep within PDM manifold

C-r = vec

Cll

C12

ClAi

cmm

(

Cll C12

ClAi

cmm

Projection

3. Apply usual beamformer with the projected matrix

minimize b,

bTc"rh*b,

subject to b; g; = 1

Fig. 2. The main steps of the proposed approach. Wc manipulate with a the sensor-space covariance matrix Cx 1. We consider Cx as an element of a linear vector space and first vectorize it to enable standard operations defined in the linear vector spaces. The vectorized sensor-space covariance matrix if <:■;€.) is shown with the red dot in the coordinate axes corresponding the trivial basis of the linear space of matrices. As a vector and based on the generative model (8) it can be viewed as a sum of two non-orthogonal vectors. The first vector is a linear contribution of sensor-space covariance auto-terms modulated by source power The second vector is a sum of the pairwise cross-products g gj °f source topographies weighted with source space covariance coefficients . -f gj&J)cij-

Our goal is to suppress the contribution of this second vector to C, 2. We do this using the special prccomputed projection matrix P that we apply to r^riCj; the result is then reshaped back (uec-1) into MxM matrix and care is taken to make sure is a positive definite matrix. 3. Finally, adaptive beamformer spatial filters are calculated as usual but using the modified sensor-space data covariance.

2.3.2. ReciPSIICOS projection

The first approach to suppressing the manifestation of source activity correlation in the sensor-space covariance is based on projecting the vectorized data covariance uec(Cx) onto the source power subspace defined here as the X-dimensional principal subspace of collection of the vectorized auto-products of source topographies uecfg.g^),

i = [1.....IV], where N is the total number of sources in the forward

model and K is the user defined parameter. See Fig. 2 and Projection matrix in Fig. 3. A rigorous criterion used to determine K is described in Section 2.4.

Importantly, unlike in Kimura et al. (2007) our approach for modification of the sensor-space covariance aimed at eliminating the manifestation of source activity correlation does not require estimating the the sources present in the data, neither their count nor locations. Rather, we use a projector matrix P informed by the sensor-space covariance generative model (8) and built using the forward model alone. To construct P we estimate the principal source power subspace Sffw/ of some fixed dimension K capturing the largest proportion of variance that could have been potentially conveyed to the vectorized covariance matrix via pec(g.-g^), i = [1.....JV] under the assumption of uniform distribution of activity across all JV sources represented in the whole brain forward model matrix. This is exactly reciprocal to what is done in (Ossadtchi et al., 2018) where we solve the complementary problem of mitigating volume conduction effects and project the vectorized data

covariance away from the principal source power subspace Spulr. The projection operator P obtained this way is directly complementary to that of PSIICOS.

More specifically, to build the projector we use the following sequence of steps:

1. Construct matrix Gpwr, where columns correspond to the vectorized auto-products of topographies for all available JV sources in the forward model. For compactness we are considering here the fixed orientation case and treat arbitrary source dipole orientations later in Section 2.5. We will denote the vectorized auto-product of topogra-

. can be created

phies as q,, = vec(g,gj) = g, ® gf. Then, matrix Gpl by stacking qi; as columns, i.e.

Gpwr = [in.....1/,.....q.VA'l

(9)

Compute the singular value decomposition of matrix Gfwr obtained in the previous step, Gpwr = VpwrSpwr{\pwr)T. and create the projector onto the principal source power subspace S*

P = IT

(v* Y

<r\ pwr f '

(10)

where U*, is formed from the first K left singular vectors. K is the projection rank and the only parameter which is needed to be chosen manually according to the guidelines given in Section 2.4. Note that this projection is complementary to that used

Preprocessing

Correlation matrix On

Bandpass filter, artifact rejection

Epoching

ERPs

uYy'

Projection matrix

G = [gl,. •• ,gj, - ,gjv] 1

Gpwr = [qib----q«, ■ • •, q/\w], Qu = t'cci g.g/ )

I SVD decomposition

p = u* fuft' V

Gpwr - Upior Spiuî- Vj^yj.

Projection rank

Output

Сж = uec 1 (P • vec(Cxj)

Tabs Jx

Tabs :n

Ё|Л|Е5

Use Cat>s instead of Ст in

x л

LCMV beamformer

Fig. 3. A schematic presentation of the ReciPSHCOS algorithm. Preprocessing: We first epoch the data and compute evoked response fields. Usually, at this step to account for MEG data not being numerically full rank, we project the data into a space of virtual sensors derived from the forward model matrix and capturing user-specified percentage of its variance. Projection matrix: This is the main ingredient of our method. We are creating projection matrix P aimed at emphasizing the source power component (and at the same time suppressing source correlation component] of the sensor-space covariance, see Fig. 2 and expression [81 (green

, IV into a single matrix Gr„., see

box]. We operate in the product space of outer product of source topographies and collect all auto outer products q;f = gfg; Eq. (13). Next, we perform SVD of this matrix to identify the subspace capturing the largest proportion of variance for any given dimension K. Then, using the first

K left singular vectors u^kt1.....upwrK of G,,, r we form the projection matrix P to be applied to the vectorized regular sensor-space covariance matrix computed

during the Preprocessing step. The projection matrix needs to be formed only once for a given forward model. Output: Lasdy, we modify sensor-space covariance matrix C, by projecting its vectorized version into the A-dimensional subspace capturing most of the variation in the columns of Gpwr. We then reshape the projected vector back into the matrix form and apply the spectral flip approach (Duin and Pekalska, 2010] to ensure non-negative definiteness of the resulting matrix Cf1. We then use C'r'r instead of C . in (3) and (5) to compute the adaptive beamformer weights and output power correspondingly.

in Ossadtchi et al. (2018) where our goal was to suppress the manifestation of volume conduction.

3. Apply the obtained projection matrix P to the vectorized sensorspace covariance matrix vec(Cx) in order to project it onto the source power subspace and reduce the contribution of the components modulated by the off-diagonal elements of the source-space covariance matrix. The projected matrix is then

Cx = чес1(Р-вес(С^)). (11)

4. The projection procedure followed by reshaping returns a symmetric matrix, see below. However, the projection procedure does not guarantee the positive definiteness of the resulting matrix C,, and we ensure this property by replacing its negative eigenvalues with their absolute values. This operation is inline with spectral flip approach introduced in Duin and Pekalska (2010) which has to be used carefully. See Section 3.2.3 for further analysis. Thus, we calculate the final data covariance matrix with suppressed source covariance component contribution from the source-covariance of the correlated sources as

Cf = E|A|E

(12)

where E and A are the matrices containing eigenvectors and eigenvalues of the matrix after the projection C, and | ■ | denotes element-wise absolute value. In Section 3.2.3 using the examples of real data covariance matrices we explore In greater detail the issue with nonpositive eigenvalues in the post-projected matrix.

5. Use projected covariance matrix C°bl in (3) Instead of the original C, to compute spatial filters b;, in (4) to reconstruct the source ac-

tivation timeseries and in (5) to estimate the source power a, distribution.

We can prove the symmetry of C , if we realize that M2-dimensionai vector uec(Cx) projected onto the if-dimensional principal subspace of Gpwr = (^^(EffB7)}' n = ■ is a linear combination of K left singular vectors u, of C,pic,. Since each singular vector is simply a linear combination of the columns of Gpwr we can write that P ■ rcciC^ } = Vakuk = 2f=1 at Pl; ■ Reshaping this expression back into M x

M matrix form results into a summation of scaled symmetric matrices that are symmetric Itself. Similar arguments can be made to prove that the Whitened ReciPSHCOS projection described next also leads to a symmetric matrix. Therefore, the eigenvalues of this matrix are real numbers.

In vector diagrams in Fig. 2 we show the decomposition of data co-variance Cx as a summation of components from S^wr and S^or subspaces of the M2-dimensional space. These principal subspaces are not orthogonal. After projection performed at Step 2 the contribution from Sfor to Cx is supposed to be reduced and the components projecting into Sjfwr are to be left Intact.

The ReciPSHCOS procedure of source reconstruction Is schematically shown in Fig. 3 for the case of ERP analysis. We use MEG data projected into the forward model principal subspace defined as the subspace capturing 90% of the energy in the forward operator matrix. This allows reducing the number of dimensions to about 60 - 80 virtual sensors. Following the preprocessing, we split the data into the epochs according to the timestamps of the stimulus onsets. Next, we calculate average

evoked response (ERP). We then use these average ERPs to compute sensor-space covariance matrix. At this step, we have to make sure that the ERP curve is of sufficient length and that our data covariance matrix is full rank.

We follow the steps listed in the beginning of this section to prepare projection operator based on the forward model matrix corresponding to the virtual sensors. This has to be done only once for a given forward model matrix. If the forward model matrix needs to be recomputed for some reason, the ReciPSIICOS projector has to be recalculated.

We then apply this projection operator to the vectorized data co-variance matrix and reshape the result of the projection to a square matrix and heurlstically fix negative eigenvalue problem as described above. Our experience shows that contribution of the eigendirections modulated by the negative eigenvalues is quite slim and does not exceed 20% of the total energy in the eigenvalue spectrum of the projected matrix, (see section 3.2.3) In our view the proposed heuristics provides an adequate balance of efficiency and simplicity as compared to the full blown consideration of Riemannian geometry of the correlation matrices, which nevertheless constitutes an interesting direction for further improvement of our method. In this section we have described a simple projection into the power subspace. This method delivers a reasonable performance and provides adaptive beamforming scheme with immunity to the source-space correlations inevitably present in the MEG and EEG data. See Section 3 for description of simulations and real-data analysis results where we refer to this method as ReciPSIICOS.

Because and non-trivially overlap, the described projection procedure depletes power from both subspaces, but does so at a faster pace for S*r subspace. The ultimate balance between the power in the two subspaces achieved with specific value of projection rank K may serve as a quality metric for the procedure. We use these considerations in Section 2.4 where we suggest a strategy for choosing optimal projection rank.

In what follows we will describe another, slightly more complicated procedure, that allows achieving a better balance between the depletion ratio of power in the S^wr and subspaces.

2.3.3. Whitened ReciPSIICOS projection

Projector P built according to the second approach directly projects the vectorized sensor-space covariance away from the subspace that captures dominant variance modulated by the source-space correlations. This is done by first estimating the K dimensional principal correlation subspace S* that captures most of the power in the collection of "ec^g.gj + gjgf)' M = P. ■ ■ ■. vectors and then projecting the covariance matrix away from However, in order to maximally spare the source power terms from this projection we operate this projection in the space whitened with respect to In what follows we refer to this method as Whitened ReciPSIICOS. Step-by-step algorithm of building and applying the whitened projector is presented below as well as the schematic representation in Fig. 4.

1. Construct matrix Gcor whose columns span correlation subspace Scttr. The columns of this matrix contain vectorized sums of symmetric source topographies outer-products for all ordered pairs of

sources. Using qtj = vec ^g gj ) = g, ® gj shortcut form matrix

G»r= [-,qii +q,h-] (13)

and calculate matrix C„„ = GcorG^,r.

2. Construct matrix Gpur whose columns span power subspace Spwr Columns of this matrix are the vectorized sums of source topography vectors as in the previous projector

(14)

.Ijtf»]

Then calculate the matrix Cpulr = GpwrG^r. 3. Using eigendecomposition of Cpwr. compute the whitening operator W

where Epwr is the matrix of eigenvectors of Cpwr and diagonal matrix \pwr contains the corresponding eigenvalues.

4. Apply whitening transformation to Cc0, to obtain C^ = w r W7^

pWt tot pu:r'

5. Extract the principal subspace of by means of eigenvalue decomposition as

(16)

6. Form matrix projecting away from the source correlation subspace S, and operating in the space whitened with respect to Spwr

p = W"1 il-E"

pwr \ ft

W

(17)

where 1 is identity matrix, E"'^ is the matrix of the first K eigenvectors of matrix Wf(4, is the whitening matrix computed earlier and K is the projection rank and the only parameter in the introduced algorithm that needs to be chosen manually.

7. Apply the obtained matrix P to the vectorized sensor-space covariance matrix i.'Cf: (C, j in order to project it away from the principal source correlations subspace and thus reduce the contribution of the components modulated by the off-diagonal elements of the source-space covariance matrix. The projected matrix is then

C^imr'tP-weiC,)) (18)

8. Again, as with constructing the first projection it can be proved that the projected matrix is symmetric. Also, since this whitened projection procedure does not guarantee the positive definiteness of the resulting matrix Cx, we ensure this by replacing its negative eigenvalues with their absolute values, In other words, we calculate the final data covariance matrix that has the contribution of the correlated sources suppressed as

qf =E|A|E'

(19)

pwr pwr pwr'

(15)

where E and A are the matrices containing eigenvectors and eigenvalues of C^.

9. Use the projected covariance matrix C°bs in (3) instead of the original Cx to compute the spatial filters b, to reconstruct the source activation timeseries (4) and to estimate the source power trf distribution

in (5).

Note that this projection is still conceptually complementary to the original PSHCOS (Ossadtchi et al., 201B), as it projects the vectorized covariance away from the AC-dimensional principal subspace S* computed using the columns of Gcor.

Both projection operations will inevitably affect both subspaces. However, as it has been demonstrated earlier in Ossadtchi et al. (2018), the use of spectral value decomposition procedure (SVD) allows identifying a low dimensional subspace of or capturing most of

the power of the auto-product terms uec(g.g^), i = [1.....N] and the

cross-product terms uec^gj + g^gf J, i,j = [1.....JV].

As in the previous case, this projection procedure is designed without taking into account the Riemannian geometry of the manifold of the correlation matrices. Interestingly, this approach is based on projecting the original vectorized covariance matrix away from Stlir yields even smaller fraction of negative eigenvalues. In our experience, this whitened projection operator consistently results Into less than 20% of the total energy in the eigenvalue spectrum of the projected matrix.

2A. Optimal projection rank

In the previous sections we introduced two projection procedures aimed at suppressing the contribution of the components modulated by the off-diagonal elements of the source-space covariance matrix. The projection procedures require only one parameter, the projection rank, that has to be preset by the operator. Since the projection depends only on the forward model, a possible solution is to select the projection rank K based on the simulations like the ones described in Section 3.1.

Whitening operator

G pwr — [qi 11 * • •) Qti) • • -,Qiviv],q« = vec(gigf)

Eigendccomposition

Gpim-Gpim- — EpujrApmrEp^,. ~Wpwr — ^pwr^-pu/r ^pw

Correlation matrix C

Projection matrix

Gcor = [..., qij + qji, • - ■] , <hj = vec(gigf )

Whitening transform

C'" = W G G W

^cor — VT pwr^cor^'cor v* pwr

Eigendecomposition

K

Projection rank

Fig. 4. Schematic presentation of the Whitened ReciPSIICOS algorithm.

Output

Cx = vec-1 (P ■ vec (Cx))

Caxbs = E|A[ET

Use S instead of C LCMV beamformer 'a; in

The other approach is based on the following considerations. We operate in the Mz-dimensional space and consider Eq. (8) as the observation equation of the vec(Cx). This equation represents vec(Cx) as a superposition of activity of sources with the topographies qs = vec(g, gf), i = l.....iiandh,, = uec(gigJ+gj.gf),/ = 1.....R, j = (i + 1).....it. Vectors q(i and h; . mainly project into Spult and S*r subspaces correspondingly. To calculate maximum possible energy in each of the subspaces we can assume that the "sources" with "topographies" qd and h,y are activated with unit independent activations and then the total energy in both subspaces can be measured as the trace of the M1 x M1 correlation matrices Cpu.r and Ccor introduced earlier in Section 2.3.3. After the projection procedure of rank k is performed by means of matrix Pt, the subspace correlation matrices of the observable vec(Cx) will read as VkCpwtP^ and PkCc„Pl The fraction of power left in both subspaces can be then computed as

lrace( C )

and

trace(PtCcorPT)

(roce(Cc„)

Based on these equations we can plot the curves (PpuJr

(20) (21)

№),

parameterized by projection rank parameter k. An example of such a plot is shown in Fig. 5 that depicts the fraction of past projection power in both subspaces parameterized by the projection rank k.

With the reduction of projection rank in the ReciPSIICOS procedure that performs projection onto the power subspace, more power is depleted from both the subspaces Spilr and S^. Conversely, the Whitened ReciPSIICOS implements the projection away from the correlation subspace, so the power Is suppressed while the projection rank grows.

Our method works because there exists a range of projection rank values, where the power Is depleted faster from the correlation subspace than from the power subspace, i.e. Sp = - > 0. However, when projection rank k reaches a certain value K*, the rate of depletion of power from the two subspaces become comparable. Consequently, for the values of k > K*, the difference Sp changes its sign to negative, i.e. the power subspace starts to lose power faster than the correlation subspace. This projection rank value K* can be considered

as optimal. Our simulations show that the projection operation remains fairly stable for a broad range of projection rank values k. Yet, these values are affected by the number of virtual sensors used and would significantly differ for different MEG probes.

Fig. 5 demonstrates the power suppression curves for the two proposed methods and for the two MEG arrays: 204-gradiometer array of Neuromag Vectorview-306 (panel A) and 275 gradiometer array of CTF (panel C). In both cases, the number of virtual sensors was selected so that their leadfiled matrix captured 99% of variance present in the original leadfield matrix which resulted into MNms = 50 and MCTF = 42 number of sensors. Panels B and D show the logarithm of marginal gain, defined as log(dPe0r/dk) - log(dPpwJdk), as a function of projection rank. According to the imposed assumptions, the optimal projection rank is the value which sets the logarithm of marginal gain in power depletion to zero, which also means that the angle of the tangent to power suppression curve equals to 45° at this point.

This approach demonstrates the clear superiority of Whitened ReciPSIICOS method over plain ReciPSIICOS, as for each value of preserved power in Spwr Whitened ReciPSIICOS suppresses more power in

2.5. HandJing unknown source orientations

Anatomically, the dipole orientations coincide with the direction of apical dendrites of the pyramidal neurons and therefore are predominantly orthogonal to the cortical mantle.

Modern tools of MRI data analysis allow for a very accurate extraction and precise parametrization of the cortical surface with the number of nodes on the order of several hundreds of thousands, which in turn results in a reasonable accuracy of the orientation specification. However, in practice significantly sparser grids of tens of thousands of nodes are used which leads to the increased uncertainty in the orientation parameter. To compensate for this it is customary to place a triplet of dlpoles in each node to accommodate any orientation to be learnt from the data. The local forward matrix Is then represented by a triplet of topographies corresponding to the three dipoles. As shown In Ahlfors et al. (2010), even with realistically shaped forward models the third dimension in the node forward model conveys to the data median 6% of energy of the local source with the most powerful orientation. Based on this, we follow the common practice in MEG inverse modelling and reduce local

Expanding these auto-products of topographies of the oriented

dipolesin terms of the topographies g*, gf, I = [1.....N] oriented along

x and y axis of the local tangential plane we obtain

««(g f'(gf'f) = g-'®gf' =(G,e,.)®(G,.e,.)

= ® s' + efflf (ef ® ef ® B*) + ® sf

(22)

In the above for compactness we have replaced the operation of vectorized outer product of two vectors with an equivalent operation of Kronecker product of these two vectors since vec(&, b) = a ® b.

Therefore, in order to accommodate the arbitrary orientation constraint, elements g, ® g, in Eqs. (9) and (14) have to be replaced by

rg/ ® gj'. gf®gf + gf ®gf, gf®gf].

Arbitrary orientations lead to a slightly more cumbersome manipulation for the cross-products of topographies. Similarly to the auto-products step sums of vectorized outer-products of topographies of different arbitrarily oriented sources present in Eq. (13) can be expanded using Kronecker product notation as

wc(i?'(ïj'>1') + uec(g"J(g"')T) = g,' ® g/ +g/ = (6,8,) ® (G,8,)

+(G;8,) ® (G,e,) = e*sj(;gf ® gj + gj ® gp + vf-t/J/tf ® gj

+gj ® gf) + eJfffCgf ® g* + g'j ® gf) + 8>e>(gf ® gj + gj ® gf). (23)

Therefore, elements nee ^f'(B,J )r) + ^(g^Cgf')7" j in Eq. (13) are to be replaced with

In order to compare the proposed algorithms with other relevant source reconstruction methods we ran several Monte Carlo simulations. To simulate the MEG signals, we used cortical surface model with 20000 vertices reconstructed from anatomical MR1 data using FreeSurfer software (Flschl, 2012). The forward model matrix G for freely oriented dipolar sources was computed with Brainstorm software (Tadel et al., 2011) using the overlapping spheres procedure. Each location on the cortical grid was served by two topography vectors confined to the locally tangential plane as determined by the first two right singular vectors of the local [M x 3] forward matrix.

In different experiments, we studied the interaction of two or three sources. For each Monte Carlo trial, a random set of pairs or triplets of dipolar sources was picked as the target stimulus-related sources. One hundred epochs were generated and then averaged to obtain the event related field (ERF). The activations of target sources s(t) were modeled with 10-Hz sinusoids.

The simulated data mimicked two experimental conditions: the first with highly correlated source activation timeseries and the second with uncorrelated ones. To model these two extreme cases, the phase difference between activation timeseries was set to zero for the correlated condition and to | for the uncorrelated condition. Each trial onset was jittered with respect to the task onset by adding a random shift generated from zero-mean noimal distribution with the standard deviation corresponding to ^ phase difference.

We modeled task irrelevant activity with 1000 task-unrelated cerebral sources whose locations and timeseries varied with each epoch. Source locations were mapped on the nodes of the high resolution cortical grid (20000 vertices). The activation timeseries were narrow-band signals obtained via zero-phase filtering of the realizations of Gaussian (pseudo)random process by the fifth order band-pass IIR filters in the bands corresponding to theta (4-7 Hz), alpha (8-12 Hz), beta (15-30 Hz)

and gamma (30-50 Hz, 50-70 Hz) activity. Their relative contributions were scaled in accordance with the well-known 1 // characteristic of the MEG spectrum. We scaled the brain noise components to match typical signal-to-noise ratio of real-life recordings. To project these sources into the sensor space, the corresponding columns of the forward matrix was computed for the high resolution source grid. We simulated 100 epochs of ERF data and for each epoch a new randomly picked set of noisy sources was chosen and new noisy timeseries with approximately 1 // spectrum were generated. We defined the SNR in our simulated data in the sensor space as the ratio of Frobenius norms of data matrices for the induced and brain noise components filtered in the band of interest (0.5-7 Hz), corresponding to the ERF response.

We did not perform whitening of the data (and the model) with respect to noise covariance C„. In the ERP setting this covariance could be estimated from a small number of samples over the pre-stimulus interval. Whitening with respect to such a noisy covariance estimate would require regularization for which there is no universally working recipe (Engemann and Gramfort, 2014). Also, the whitening changes the forward operator which would then require recomputing ReclPSHCOS and Whitened ReciPSIICOS projections. The former is fast and the latter is way slower and takes about 5-10 minutes on a regular PC and normally needs to be done once per forward model. However, In the case of 500 Monte-Carlo simulations per condition with new noise samples on every MC trial this would result into a very significant computational burden that we decided to avoid. Given that the proper whitening may Improve the source localization accuracy, the results reported here for the low SNR cases can be conservative.

The high resolution grid of 20,000 vertices was used only for data simulation. For the source reconstruction process we employed a 4 times sparser cortex grid with 5000 vertices. Both projectors (10) and (17) were computed using the sparser grid. We ran 500 simulations and compared two versions of the proposed ReciPSIICOS method against the classical adaptive vector LCMV beamformer (Sekihara et al., 2001) and MNE techniques.

We have also projected the data into the principal subspace capturing 95% of variance in the forward model matrix. This helps us to improve the condition number of data covariance matrices. In addition, we follow a safe strategy and augment matrix inversions with a simple Tikhonov regularization, however, the regularization parameters are usually very small and fall into 10~3 - 10~2 range.

2.7. Performance metrics

In our Monte Carlo simulation analysis, goodness of localization was estimated using three metrics: source localization bias, point spreading radius and ratio of successful detection. While the data were simulated using dense cortical model with 20,000 sources, source reconstruction procedure was based on the sparse cortical model with 5000 sources to make the reconstruction procedure more realistic. We considered separately the scenarios with two and three target sources. In order to run the comparative analysis for each Monte Carlo trial we repeated the following steps: