Численное моделирование функционирования ансамбля нейронов коры головного мозга

Аннотация

Д.Ю. Душенин

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

05.11.17 - Приборы, системы и изделия медицинского назначения

05.13.01 - Системный анализ, управление и обработка информации (по отраслям)

Эпилепсия представляет собой хроническое заболевание головного мозга, характеризующееся повторными непровоцируемыми приступами нарушений двигательных, чувствительных, вегетативных, мыслительных или психических функций, возникающих вследствие чрезмерных нейронных разрядов. Перед операцией врач сначала должен определить, где находится зона эпилепсии. Чтобы это сделать, он может провести неинвазивную (ЭЭГ) и инвазивную (ЭКоГ) записи электрической активности кортекса.
Для того чтобы понять в дальнейшем электрическую активность мозга человека, записанную во время приступа, исследователи разработали многочисленные модели. Чтобы смоделировать данные, полученные с одного электрода, необходимо воспроизвести поведение приблизительно 105 отдельных нейронов. К сожалению, используя физиологически точные математические модели эту задачу нельзя выполнить.На основании обзора существующих моделей была выбрана для дальнейшего исследования математическая модель мезоскопической электрической активности мозга человека [1], в основе которой лежат усредненные свойства соседних нейронов. Одна из наиболее полных математических моделей мезоскопической электрической активности нейронов коры – модель Штайн-Росса (Steyn-Ross, 2003) [2].
Модель состоит из системы 8 дифференциальных уравнений, которые описывают распространение ЭЭГ волны. Для того чтобы вызвать приступоподобное состояние в модели были изменены два параметра, связанные с процессом возбуждения [2]. Эти параметры выбраны по двум причинам: во-первых, известно, что повышенное возбуждение провоцирует приступ. Во-вторых, увеличение уровня мембранного потенциала нейронной популяции считается важным фактором, влияющим на появление приступа.
Основная переменная в модели – переменная h_e – усредненный по поверхности соматический мембранный потенциал возбудительных клеток кортекса. В уравнениях используется безразмерный параметр h_e= h~_e·(-70 мВ), где мембранный потенциал нормируется на потенциал покоя. Экспериментальные данные ЭЭГ пропорциональны модельной переменной h_e.
Модельное состояние «приступ» было определено как решение системы обыкновенных дифференциальных уравнений, для которых h_e подвергается устойчивым осцилляциям большой амплитуды. Были вычислены множественные решения системы в диапазоне параметров, используя метод Рунге-Кутты четвертого порядка [3] с временным шагом 0,4 мс в программе MathCAD. Для каждого решения был вычислен спектр мощности для h_e после окончания переходных процессов и установления частоты максимальной мощности. Спектр мощности имеет при этом единственный четко определенный пик.
Для решения дифференциальных уравнений Mathcad имеет ряд встроенных функций, в частности, функцию rkfixed, реализующую метод Рунге–Кутты четвертого порядка с фиксированным шагом. Фактически эта функция предназначена для решения систем дифференциальных уравнений первого порядка. Функция rkfixed(y, x₁, x₂, npoints, D) возвращает матрицу. Первый столбец этой матрицы содержит точки, в которых получено решение, а остальные столбцы – решения и его первые n-1 производные.
Особенностью данной функции является то, что решение возвращается в виде массива с запрошенным при ее вызове количестве строк (рассчитанных точек). Каждая строка содержит значение аргумента и значения рассчитанных в этой точке искомых функций.
На рис 1 показан рассчитанный график временной зависимости потенциала возбудительных клеток кортексаh(t). На рис. 2 приведен алгоритм вычисления спектра, содержащий некоторые особенности: вычитание среднего значения, конкатенация векторов для придания функции «гладкости», выбор числа окон, обнуление спектра.

Рис. 1. Функция решения системы ОДЕ и результат

Рис. 2. Особенности вычисления спектра электрической активности кортекса

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

ЛИТЕРАТУРА
1. Mark A. Kramer, Heidi E. Kirsch and Andrew J. Szeri. Pathological pattern formation and cortical propagation of epileptic seizures // J. R. Soc. Interface. – Springer, 2005.
2. Ильина В.А., Силаев П.К. Численные методы для физиков-теоретиков. т.2. - Москва-Ижевск: Институт компьютерных исследований, 2004. – с. 16-30.