Нейронное напряжение
Как дифференциальные уравнения помогают понять работу клеток головного мозга
Чтобы понять, как работает человеческий мозг, недостаточно усилий одних только биологов. Нейроны работают с помощью, в том числе, электрических сигналов, а значит, можно построить электрическую схему отдельно взятого нейрона и описать ее математически. И здесь нам снова понадобятся дифференциальные уравнения — системы уравнений Фитцхью-Нагумо. О том, как именно математики исследуют нервную систему и мозг, читайте в нашем четвертом материале о самых интересных дифференциальных уравнениях (предыдущие можно прочитать здесь, здесь и здесь).
Термин «вычислительная нейробиология» (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 и проследить за поведением системы.
Юлия Ровник
Математические футболисты оборонялись почти как живые, хотя на штрафные не собирались
Аргентинские физики предложили подход на основе изменяющихся во времени двудольных сетей близости для описания оборонительной тактики футбольных команд. Авторы не только смогли найти приемлемый способ характеризации персональной и зональной опеки, но и на основе данных, собранных с реальных матчей, построили механическую модель футбольной команды, которая сохранила основные статистические особенности обороны своего прототипа, хотя и не смогла воспроизвести штрафные. Исследование опубликовано в Physical Review E. В последние годы спортивные состязания, в особенности футбол, все чаще становятся объектом сбора большого объема данных и их анализа новыми статистическими методами. Исследователям удается извлекать из них все больше полезной информации, что вызывает интерес к этой области со стороны представителей спортивной индустрии. Мы уже рассказывали, как математика и статистика помогает предсказывать успешность пасов, оценивать талант футболистов и отличать опытных игроков от новичков во время тренировки, а недавно руководство одного из клубов решило собирать данные о болельщиках. Одним из самых мощных вычислительных методов в спортивной науке остается сетевой анализ, зарекомендовавший себя и в исследовании других аспектов человеческой деятельности, например, коррупции, работы городских транспортных сетей и даже эволюции мемов. Применительно к футболу этот метод, как правило, основан на построении сети игроков, взаимодействующих друг с другом. Чаще всего ученые рассматривают взаимодействие в рамках одной команды, чтобы количественно охарактеризовать командную работу, однако такой подход не дает оценить иные аспекты футбольного матча, например, эффективность защиты. На решение этой проблемы направили свои усилия аргентинские физики во главе с Андресом Чакомой (Andrés Chacoma) из Национального университета Кордовы. На основе данных, собранных с реальных футбольных матчей, они с помощью сетевого подхода охарактеризовали то, как команды используют тактику опеки. Смоделировав футболистов с помощью системы динамических уравнений, ученые смогли воссоздать основные статистические паттерны, характерные для реальных игроков. Опекой в футболе называют оборонную тактику, связанную с преследованием защитниками и полузащитниками нападающих чужой команды. Преследование выражается в постоянной близости к опекаемому, что можно описать через дистанцию между игроками. Для этого, в свою очередь, необходимо знать координаты всех игроков на поле в любой момент матча. Такие данные были выложены в публичный доступ компанией Metrica Sports. Они были собраны с трех матчей с помощью обработки видеозаписей и обезличены, поэтому неизвестно, о каких командах и играх шла речь. Данные представляли собой информацию о координатах всех 22 игроков с пространственным разрешением 10 сантиметров и частотой 25 кадров в секунду. Для сглаживания шумов, физики усредняли кадры до массива с шагом в одну секунду. Авторы характеризовали динамику матча с помощью временных двудольных сетей близости, то есть таких графов, в которых каждый игрок одной команды связывается лишь с игроками другой команды. При этом связь возникает только тогда, когда расстояния между футболистами меньше некоторого порога. Ученые тщательно исследовали связность и распределения кластеров в таких сетях в зависимости порогового значения, а также то, как эта зависимость меняется со временем. Особый интерес для них представляли события, названные лавинами, когда сеть становится максимально связной, то есть игроки группируются в некоторой точке поля. Такая ситуация возникает как в случае активной опеки, так и во время пробития угловых или во внеигровые моменты. Распределения интенсивности и длительности лавин свидетельствовали об их самоподобном характере, а их параметры можно использовать для количественного описания оборонной тактики. На следующем этапе работы физики описали движения футболистов с помощью второго закона Ньютона. Ускорение каждого игрока определялось аналогом силы вязкости, не дающей ему разгоняться бесконечно, аналогом силы упругости, привязывающей его к точке на поле, определяемой выбранной тренером тактической схемой, а также некоторой силой взаимодействия, зависящей от расстояния до других игроков с соответствующими коэффициентами. Для определения всех свободных параметров модели авторы брали конфигурацию матча в некоторый момент игры и минимизировали разницу между реальными и вычисленными скоростями спустя какое-то время. Чтобы понять, насколько построенная таким способом динамическая модель хорошо воспроизводит статистику близости, физики добавляли в уравнения стохастический шум и запускали симуляцию. Оказалось, что оборона и лавины в виртуальном матче организована примерно по тем же распределениям, что и в реальном матче. Модель не смогла воспроизвести лишь совсем экстремальные лавины, не связанные напрямую с тактической обороной: пробитие штрафных и угловых, а также группировка футболистов вне игры. Подробнее о том, как большие данные и машинное обучение помогают большому спорту, читайте в материале «Победа по расчету».