Комплексная диагностика паркинсонизма по данным магнитной энцефалографии

Дергузов А.В. (arkadid@list.ru) (1), Махортых С.А. (1), Семечкин Р.А. (2)

(1) Институт математических проблем биологии РАН (2) Московский государственный университет, факультет вычислительной математики и

кибернетики

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

В последнее время интенсивно развивается подход к распознаванию, использующий спектральное представление сигнала ортогональными функциональными разложениями [3, 4]. При этом элементами признакового пространства являются векторы коэффициентов Фурье.

В работе [6, 7] предложен метод классификации данных магнитной энцефалографии на основе разложений пространственного распределения магнитного поля по поверхности головы от спонтанных и вызванных источников, возникающих в объеме мозга в процессе нормальной и патологической активности (в работе [6, 7] рассматривались случаи заболевания паркинсонизмом). На вход предлагаемого метода подается меняющийся во времени 148-мерный вектор измерений, на выходе получается вектор спектральных компонент по сферическим гармоникам. Коэффициенты разложения используются в качестве признаков. В упомянутой статье получены критерии надежности распознавания патологий. Однако достичь 100 % распознавания моментов времени с патологией не удается в силу отмеченной выше сложности изучаемой системы и наличия внешних и внутренних помех. Для проверки выводов, получаемых в рамках метода [6, 7], предлагается использовать дополнительные алгоритмы классификации на основе 1) решения обратной задачи локализации источников поля; 2) исследования размерностных характеристик сигнала на участках нормальной и патологической активности (будут рассмотрены случаи болезни Паркинсона - Tinnitus, периодическая генерация тремора и дистония).

G = 1 -

Исходные экспериментальные данные получены с помощью 148-канального измерительного стенда Magnes 2500 WH в Медицинской школе Нью-Йоркского университета.

Обратные задачи магнитной энцефалографии локализации источников биомагнитной активности мозга

Основное внимание в настоящем работе будет уделяться задаче пространственной локализации источников патологического сигнала, связанного с разновидностями болезни Паркинсона с использованием спектральных характеристик пространственного распределения поля магнитной индукции (МИ) по поверхности головы. Измеряемый сигнал представляет собой пространственно-временную структуру: 148-мерный вектор измерений в 148 точках на поверхности головы, развернутый во временной ряд с частотой опроса датчиков 500 Гц.

Вся работа по анализу получаемых данных может быть разбита на следующие этапы:

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

2. Выбор моментов времени для решения обратной задачи локализации токовых источников по пространственной картине поля на поверхности головы в эти моменты времени (с использованием результатов [6, 7]).

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

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

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

Если известны приближенные положение и направление токового диполя, в этом случае, используя выражения (9)-(11) можно вычислить величину магнитного поля на всей поверхности головы. Меру точности оценки в этом случае можно определить так [8]:

tw (Б, -B0)2. (4)

Рис. 1. Исходные экспериментальные данные пациента с болезнью Паркинсона.

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

Здесь в0 - теоретическая оценка, Bi - измеренные значения, wi - вес измерительного канала, N - число каналов. Функция невязки записывается в виде

f = TLW (B. -B0)2 min. (5)

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

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

Решение обратной задачи с использованием спектральных разложений

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

Спектральный подход пространственной локализации источников биомагнитной активности записей МЭГ протестирован на клинических записях, содержащих аномальную активности, связанную с болезнью Паркинсона и ее разновидностями (синдромы Tinnitus, акинезия), генерация тремора, Dystonia.

Рис. 2. Сечения ЯМР-томограммы головного мозга с локализованным источником магнитной

активности.

На рис. 3 приведена локализация токовых дипольных источников записей МЭГ как с паркинсонической патологией, так и с нормальным типом активности, число дипольных источников полагалось равным единице. При этом использовалась методика визуального анализа записей.

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

[6, 7].

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

На предварительном этапе на основе визуального анализа выделялись временные участки в записи МЭГ с наличием патологии (по амплитуде сигнала) и сигналом в норме. Для выбранных записей соответствующим различным типам активности решалась обратная задача локализация участков биомагнитной активности согласно методу минимизации функции невязки МИ (4) с использованием программы MRIAN В.В. Сычев, ИМПБ РАН.

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

В результате применения данной программы на рис.2 представлены два взаимно перпендикулярных сечения ЯМР-томограммы, проходящие через источник (токовый диполь). Кружком отмечено местоположение диполя, а цветом - величина момента, меняющаяся от максимального (черный кружок) до нулевого (белый кружок). Отрезок, выходящий из кружка, указывает направление момента диполя.

Для исключения влияния изменений амплитуды источника на картину распределения МИ по поверхности черепной коробки проводилась нормализация записей масштабированием (приведение всех значений МИ в каналах записей к усредненному абсолютному значению по каналам):

B]norm (t.) = 148 B1 (t.)

zlBk (t. )

где - B1 (tt) исходные значения МИ в канале . в момент времени t;-. Процедура распознавания типа активности подробно описана в работах [6, 7].

- 200 mm

-150 -100 -50 0 50

Рис 3. Распределение локализованных токовых дипольных источников записей МЭГ для

сигнала в норме - х и патологии - ж.

Для рассматриваемых данных размер обучающей выборки эталонных распределений МИ с патологической активностью составлял 50 записей Bltk , k < 50 - номер записи, 0 < j < 147 -

номер канала. При исследовании реальных данных пациента с диагнозом болезни Паркинсона на выборке 30000 моментов времени, длительностью сеанса 1 минуты было выявлено 67

контрольных записей миограмм (кривых тремора, снимаемых с конечностей испытуемого пациента), 3) с использованием результатов спектральной классификации [6, 7].

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

180 mm

Рис. 4 Локализация дипольных источников для найденных записей МЭГ с патологией на основе спектральных разложений, х - положения диполей эталонных записей, ж - положения диполей найденных патологических записей МЭГ, при этом также локализация патологии дополнительно найденных записей осуществляется в точке (-10, 9, -10).

На рис. 4 приведена локализация дипольных источников найденных записей на основании предварительной спектральной классификации [6, 7]. Для обучающей выборки из 50 записей было найдено 180 точек с хорошо локализующимися в одной области источниками

(рис. 4).

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

моментов времени, удовлетворяющих условию (7) и принадлежащих классу сигналов с патологической активностью.

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

Рис 5. Локализация токовых диполей для сигнала в норме - х, патологии - ж, при предположении наличия одного источника активности. (По результатам визуального анализа).

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

На рис. 5 приведен результат для предполагаемого одного источника активности. Как видно в данном предположении достигается лучшая локализация.

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

Рис. 6 Локализация двух токовых диполей (а - первого, б - второго источников) на основании результатов спектрального анализа, х - положения диполей эталонных записей, ж -положения диполей найденных патологических записей МЭГ, при этом локализация патологии осуществляется в точках а - (19, -17, 26), б - (23, -25, 50) .

Рис 7. Локализация одного токового диполя на основании результатов спектрального анализа для сигнала в норме - х, патологии - ж, при этом локализация патологии осуществляется в

точке (21, -33, 59).

4. Анализ стохастической динамики сигнала

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

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

Для динамических систем принятым представлением развития процесса во времени является построение портрета в фазовом пространстве (в пространстве, координатами которого являются переменные, описывающих поведение системы). Нелинейная стохастическая динамическая система характеризуется странным аттрактором -

i, 1=1

С (s) = lim - YjQ(s-P {хг, xj

[1,<х > 0

о) = < xi - точки в фазовом пространстве; р - расстояние.

[0,ос < 0

Таким образом, размерность DC определяется значением корреляционного интеграла C(s ), характеризующим относительное число пар точек Xi,Xj, удаленных на расстояния ri j = р (xi, xj) < s .

Для вычисления таких статистических средних, как размерность, энтропия, спектр показателей Ляпунова и других характеристик аттрактора, необходимо иметь множество точек, определенных в фазовом пространстве размерности n и принадлежащих аттрактору. Число точек M в расчетах конечно, но обязано быть достаточно большим.

M >Mmin = 102+04D

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

притягивающим множеством в фазовом пространстве, в котором расположены траектории. При анализе таких объектов - сложных, обладающих тонкой структурой - плодотворным является подход, разработанный в теории фракталов.

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

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

Корреляционную размерность можно представить в виде [9]:

D, = lim( > 1 (8)

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

Пусть имеется временной ряд экспериментальных данных, представляющий собой отсчеты некоторой физической величины: {sk }M-1. Если известен шаг по времени At, то время

t=k At. Предполагается, что физическая величина s является одной из переменных динамической системы. Система находится в стационарном режиме, т. е. фазовая траектория проходит внутри аттрактора. Для восстановления аттрактора Такенсом предложен метод временной задержки координат. В n-мерном фазовом пространстве строится последовательность точек вида:

Xk = (sk , Sk+T Sk+(n-1)x ) (9)

k = 0, m -1 = M - (n - 1)t .

Здесь т - временная задержка, n - размерность вложения. Основной результат Такенса состоит в следующем. ЕслиM оо, то множество точек xk е Rn задает вложение исходного аттрактора

почти при любом выборе наблюдаемой переменной, если n не меньше удвоенной размерности исходного аттрактора. Для оценки характеристик реального исследуемого аттрактора можно вычислять характеристики восстановленного аттрактора. С целью уменьшения ошибки, обусловленной конечностью набора экспериментальных точек {sk }M-1, необходимо проводить

расчеты при нескольких различных значениях M и n и добиваться независимости получаемых оценок характеристик от M и n в пределах заданной точности.

Анализ данных МЭГ. Исходные данные МЭГ представляют собой матрицу A е R79x20346 , где последний, 79-ый столбец а79;., представляет собой булеву переменную, принимающую 1/0

в зависимости от наличия аномальной активности мозга в данный момент. Метод решения поставленной задачи следующий: выбирается наиболее информативный столбец матрицы A , т.е. находится min (D[a1 .],D[a2 . ],...,D[a78 . ]). Затем полученный вектор-столбец разбивается

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

Все вышеописанные процедуры были проведены при помощи СУБД MS SQL Server 2000 и программы FRACTAN, разработанной Сычевым В.В. в Институте математических проблем биологии РАН (г. Пущино).