Как дифференциальные уравнения помогают понять работу клеток головного мозга
Чтобы понять, как работает человеческий мозг, недостаточно усилий одних только биологов. Нейроны работают с помощью, в том числе, электрических сигналов, а значит, можно построить электрическую схему отдельно взятого нейрона и описать ее математически. И здесь нам снова понадобятся дифференциальные уравнения — системы уравнений Фитцхью-Нагумо. О том, как именно математики исследуют нервную систему и мозг, читайте в нашем четвертом материале о самых интересных дифференциальных уравнениях (предыдущие можно прочитать здесь, здесь и здесь).
Термин «вычислительная нейробиология» (computational neuroscience) появился совсем недавно — в 1985 году его предложил профессор Бостонского университета Эрик Шварц (Eric Schwartz), и это как раз та ситуация, когда наука появилась намного раньше — в середине XX века, — чем ее название.
К тому моменту ученым уже давно было известно, что мозг и нервная система построены из электрически возбудимых клеток — нейронов, предназначенных для приема извне, обработки, хранения, передачи и вывода вовне информации с помощью электрических и химических сигналов. Нейроны могут соединяться один с другим, формируя сети; в головном мозге человека насчитывается около 90–95 миллиардов нейронов.
Описать работу нейронов проще всего через математические системы с заданными параметрами. Этим и занимаются ученые, работающие в области вычислительной нейробиологии, с помощью моделей, предложенных математиками почти семьдесят лет назад.
В 1952 году Алан Ллойд Ходжкин (Alan Hodgkin) и Эндрю Хаксли (Andrew Huxley) создали математическую модель для описания электрических механизмов, обуславливающих генерацию и передачу нервного сигнала в гигантском аксоне кальмара. Физиологической основой нервного импульса выступают потенциалы действия (их также называют спайками). Они генерируются благодаря разнице в концентрации ионов внутри и вне клетки.
На мембране клетки находятся ионные каналы — белковые молекулы, образующие в мембране поры, через них ионы могут проходить с внутренней стороны мембраны на наружную и наоборот. Мембрана клетки деполяризуется, реполяризуется, возвращается к исходному уровню поляризации и так далее. Нейрон реагирует на сигналы.
Вот электрическая схема, соответствующая этой модели:
В этой схеме каждый компонент возбуждаемой клетки имеет свой биофизический аналог. Внутреннему липидному слою клеточной мембраны соответствует электроемкость (Cm). Потенциал-зависимые ионные каналы отвечают за нелинейную электрическую проводимость (Gn, где n — отдельный вид ионных каналов).
Эта составляющая системы реализуется благодаря белковым молекулам, образующим потенциал-зависимые ионные каналы, каждый из которых отмечен некоторой вероятностью открытия, чья величина зависит от электрического потенциала, то есть напряжения мембраны клетки. Каналы мембранных пор отвечают за пассивную проводимость (GL, где L означает leak — «течь, утечка»).
Электрохимический градиент побуждает ионы к движению через мембранные каналы, он показан с помощью источников напряжения с соответствующей электродвижущей силой (En/EL), величина которой определяется реверсивным потенциалом для соответствующего вида иона. Ионные транспортеры соответствуют источникам тока (Ip).
Через некоторое время после публикации работы Ходжкина и Хаксли ученый из биофизической лаборатории в Мэриленда Ричард Фитцхью (Richard FitzHugh) занялся анализом математических свойств их модели.
В то время электронные компьютеры не были в свободном доступе, поэтому ему пришлось пользоваться аналоговым. Для решения уравнений использовалась огромное количество оборудования.
В итоге в 1961 году ученый создал упрощенную двухмерную версию модели Ходжкина-Хаксли (оригинальная модель более точно описывала генерацию потенциалов действия, но ее фазовые кривые были в четырехмерном пространстве, что заставляло постоянно использовать проекции).
Решения модели Фитцхью гораздо проще рассматривать, а также с их помощью проще находить геометрические объяснения важных биологических феноменов. В следующем году подобная система была предложена Дзин-ити Нагумо (Jin-Ichi Nagumo) с соавторами.
Рассмотрим запись системы и входные данные:
где v — потенциал мембраны, w — переменная восстановления, Iext — величина стимулирующего тока (эта переменная изображает экспериментальное подключение тока к мембране).
Мы пытаемся смоделировать генерацию спайка в гигантском аксоне кальмара. Из экспериментов мы знаем, что следующие утверждения верны:
1) изначально нейрон находится в невозбужденном состоянии;
2) если мы экспериментально изменим потенциал на небольшую величину, нейрон вернется в невозбужденное состояние;
3) если возбуждение будет сильнее, чем некая граничная величина, потенциал увеличится до очень большого значения, то есть сгенерируется спайк;
4) после спайка нейрон возвратится в состояние покоя.
В нашей модели состояние покоя нейрона равно положениям равновесия для переменной v. Так как система должна возбуждаться и возвращаться в невозбужденное состояние, нам необходима вторая переменная восстановления w, которая движется медленнее (для этого есть параметр τ, отвечающий за «масштаб времени») и возвращает систему в состояние покоя (для этого минус w).
Добавление константы a ко второму уравнению позволило Фитцхью сдвигать положение покоя по кубической кривой. С помощью нее можно влиять на устойчивость состояния покоя. Во втором уравнении также присутствует константа b, из-за которой положение равновесия при I = 0 лежит на правой возрастающей ветке и является устойчивым.
Электрофизиологические исследования также показывают, что действие умеренного тока на мембрану приводит к периодическим спайкам. Если внешний ток слишком велик, спайки блокируются. Чтобы смоделировать подобное поведение, необходимо уравнение третьей степени для потенциала мембраны.
Построим график, показывающий, как траектории v и w меняются в зависимости от значений значениях I, растущих от 0 до 1,5. Синяя линия отвечает за переменную w, красная — за переменную v.
При этом используются следующие параметры: a = −0,7; b = 0,8; τ = 12,5. При таких параметрах и I = 1 траектории наиболее напоминают пульсацию.
Теперь введем еще несколько терминов. Если мы приравняем левую часть дифференциального уравнения y’ = f(x, y) к p, то получим изоклину — кривую на плоскости, вдоль которой поле, задаваемое этим уравнением, будет иметь один и тот же наклон p.
Если же y’ = p = 0, то такую изоклину можно будет назвать нульклиной. Наша система состоит из двух уравнений, поэтому, чтобы найти ее нульклины, нужно приравнять сначала dv/dt а потом dw/dt к 0. Видно, что с возрастанием I нульклина для v поднимается.
Мы получили нульклину для v (синяя кривая, похожая на букву n) и нульклину для w, красную прямую. Видно, что с возрастанием I нульклина для v поднимается; нульклина для w не зависит от I (следующий график).
На пересечении этих нульклин находятся положения равновесия системы, то есть состояния потенциала покоя, в которых клетка находится в невозбужденном состоянии. Типы положений равновесия модели меняются в зависимости от параметров, и эти изменения симулируют биологические процессы.
Сначала построим фазовый портрет системы, позже рассмотрим, как меняются типы положений равновесия.
На картинке выше изображены некоторые типичные траектории. Область «no man’s land» очень нестабильна.
Движение v-нульклины (изменения потенциала мембраны) отражает, как мембрана гиперполяризуется перед приходом к состоянию покоя и деполяризуется перед активным спайком. Вот как выглядит этот процесс в динамике:
Теперь нужно найти положения равновесия. Они находятся на пересечениях нульклин. Можем найти все вещественные корни с помощью numpy.roots (при этом если n — количество корней, а m — степень уравнения, то мы знаем, что n ≤ m.).
Мы можем определить типы всех положений равновесия с помощью вычисления собственных значений матрицы Якоби. Автоматически покрасим положения равновесия на графиках в такие цвета: устойчивый узел — синий, неустойчивый узел — оранжевый, седло — фиолетовый, устойчивый фокус — красный, неустойчивый фокус — зеленый. Желто-оранжевым цветом нарисуем также траекторию, проходящую через такое положение.
Обобщим сведения о поведении системы, которые мы получили, когда строили графики выше. Для этого разметим график таким образом: все до первого перегиба v-нульклины назовем областью 1, все после второго перегиба v-нульклины назовем областью 3, а то, что лежит между ними, — областью 2.
В области выше v-нульклины dv/dt < 0, а в области ниже — dv/dt > 0. В области справа от w-нульклины dw/dt > 0, а в области слева — dw/dt < 0.
Рассмотрим, что происходит, если положение равновесия находится в области 3. Поскольку v движется быстрее, чем w, v быстро увеличивается, пока w остается практически неизменной. На рисунке видно, что это приводит к почти горизонтальной части кривой v-w.
Когда кривая v-w приближается к v-нульклине, скорость изменения v замедляется, и w становится заметнее. Так как все dw/dt все еще положительно, w увеличивается и кривая движется вверх. Затем положение равновесия притягивает кривую, и движение заканчивается в нем (пример — график 8).
Если фиксированная точка находится в области 2, то почти все, что мы обсуждали до сих пор, верно. Разница в том, что как только кривая выходит за пределы правого колена v-нульклины, dv/dt становится отрицательной и v быстро затухает. При движении влево кривая пересекает w-нульклину справа налево.
С этого момента, хотя v и w уменьшаются, эволюция v доминирует, и кривая снова становится горизонтальной. Это продолжается до тех пор, пока кривая не достигнет левой части v-нульклины.
Кривая начинает приближаться к v-нульклине и медленно двигаться вниз. Когда она касается левого колена v-нульклины, она быстро движется к правой части v-нульклины. Это движение никогда не достигает фиксированной точки и поэтому продолжает повторяться (например, графики 3 и 4).
Это оставляет нам один последний случай для обсуждения — когда положение равновесия находится в области 1. Результат выглядит как на графике 1, рассуждения подобны.
Ниже вы можете выбрать значение I и проследить за поведением системы.
Юлия Ровник