Вступление
Сразу сообщу, что данная заметка не имеет отношения к перцептронам, сетям Хопфилда или любым другим искусственным нейронным сетям. Мы будем моделировать работу «настоящей», «живой», биологической нейронной сети, в которой происходят процессы генерации и распространения нервных импульсов. В англоязычной литературе такие сети ввиду их отличия от искусственных нейронных сетей называются spiking neural networks, в русскоязычной же литературе – нет устоявшегося названия. Кто-то называет их просто нейронными сетям, кто-то – импульсными нейронными сетями, а кто-то – спайковыми.
Вероятно, большинство читателей слышали о проектах Blue Brain и Human Brain, спонсируемых Европейским Союзом, под последний проект правительство ЕС выдало около миллиарда евро, что говорит о наличии большого интереса к этой области. Оба проекта тесно связаны и пересекаются друг с другом, даже руководитель у них общий, Генри Маркрам, что может создать некоторую путаницу в том, чем же они отличаются друг от друга. Если кратко, то конечной целью обоих проектов является разработка модели работы целого мозга, всех ~86 миллиардов нейронов. Blue Brain Project – это вычислительная часть, а Human Brain – это больше фундаментальная часть, где работают над сбором научных данных о принципах работы мозга и созданием единой модели. Чтобы прикоснуться к этой науке и попробовать самим сделать нечто подобное, хотя и в значительно меньших масштабах, была написана эта заметка.
На хабре уже было несколько интересных и информативных статей по нейробиологии, что очень радует.
1. Нейробиология и искусственный интеллект: часть первая — ликбез.
2. Нейробиология и искусственный интеллект: часть вторая – интеллект и представление информации в мозгу.
3. Нейробиология и искусственный интеллект: часть третья – представление данных и память
Но в них не рассматривались вопросы вычислительной нейробиологии, или по-другому вычислительной нейронауки, включающей в себя компьютерное моделирование электрической активности нейронов, поэтому я решил восполнить этот пробел.
Немного биологии
Рис. 1 — Схематическое изображение строения нейрона.
Прежде чем приступим к моделированию, нам нужно ознакомиться с некоторыми азами нейробиологии. Типичный нейрон состоит из 3-х частей: тела (сомы), дендритов и аксона. Дендриты принимают сигнал от других нейронов (это input нейрона), а аксон передает сигналы от тела нейрона к другим нейронам (output). Место контакта аксона одного нейрона и дендрита другого нейрона называется синапсом. Сигнал, принимаемый с дендритов, суммируется в теле и если он превышает определённые порог, то генерируется нервный импульс или по-другому спайк. Тело клетки окружено липидной оболочкой, являющейся хорошим изолятором. Ионные составы цитоплазмы нейрона и межклеточной жидкости различаются. В цитоплазме концентрация ионов калия выше, а концентрация натрия и хлора ниже, в межклеточной же жидкости все наоборот. Это связано с работой ионных насосов, которые постоянно перекачивают определенные типы ионов против градиента концентрации, потребляя при этом энергию, запасенную в молекулах АденозиноТриФосфата (АТФ). Самым известным и изученным из таких насосов является натрий-калиевый насос. Он выводит 3 иона натрия в наружу, а внутрь нейрона забирает 2 иона калия. На рисунке 2 изображен ионный состав нейрона и отмечены ионные насосы. Благодаря работе этих насосов в нейроне образуется равновесная разность потенциалов между внутренней стороной мембраны, заряженной отрицательно, и внешней, заряженной положительно.
Рис. 2 — Ионный состав нейрона и окружающей среды
Кроме насосов на поверхности нейрона есть ещё ионные каналы, которые при изменении потенциала или при химическом воздействии могут открываться или закрываться, тем самым увеличивая или уменьшая токи определённого типа ионов. Если мембранный потенциал превышает некоторый порог, открываются натриевые каналы, а так как снаружи больше натрия, то возникает электрический ток направленный внутрь нейрона, что ещё больше увеличивает мембранный потенциал и ещё сильнее открывает натриевые каналы, происходит резкое увеличение мембранного потенциала. Физики назовут это положительной обратной связью. Но, начиная с какого-то значения потенциала, более высокого чем пороговый потенциал открытия натриевых каналов, открываются и калиевые каналы, благодаря чему ионы калия начинают течь в наружу, уменьшая мембранный потенциал и тем самым возвращая его к равновесному значению. Если же первоначальное возбуждение меньше порога открытия натриевых каналов, то нейрон вернётся к своему равновесному состоянию. Что интересно, амплитуда генерируемого импульса слабо зависит от амплитуды возбуждающего тока: либо импульс есть, либо его нет, закон «всё или ничего».
Кстати, именно принцип «всё или ничего» и вдохновил Мак-Каллока и Питтса на создание моделей искусственных нейронных сетей. Но область искусственных нейросетей развивается по своему, и главной её целью является наиболее оптимальное решение практических задач, безотносительно к тому, насколько это соотносится с процессами обработки информации в живом мозге. В то время как спайковые нейронные сети – это модель работы настоящего мозга. Можно собрать спайковую сеть для распознования визуальных образов, но для практического применения лучше подойдут классические нейронные сети, они проще, считаются на компьютере быстрее и для них придуманно множество алгоритмов для обучения под конкретные практические задачи.
Принцип «всё или ничего» наглядно изображён на рисунке 3. Внизу изображён входной ток, направленный к внутренней стороне мембраны нейрона, а вверху – разность потенциалов между внутренней и внешней стороной мембраны. Поэтому согласно доминирующей ныне концепции в живых нейронных сетях информация кодируется во временах возникновения импульсов или, как сказали бы физики, – путем фазовой модуляции.
Рис. 3 — Генерация нервного импульса. Внизу изображен подаваемый внутрь клетки ток в пкА, а вверху мембранный потенциал в мВ
Возбудить нейрон можно, например, воткнув в него микроэлектрод и подав ток внутрь нейрона, но в живом мозге возбуждение обычно происходит путем синаптического воздействия. Как уже было сказано, нейроны соединяются друг с другом с помощью синапсов, образующихся в местах контакта аксона одного нейрона с дендритами другого. Нейрон, от которого идет сигнал, называется пресинаптическим, а тот к которому идет сигнал – постсинаптическим. При возникновении импульса на пресинаптическом нейроне, он выделят в синаптическую щель нейротрансмиттеры, которые открывают натриевые каналы на постсинаптическом нейроне, а дальше происходит цепь описанных выше событий, приводящих к возбуждению. Кроме возбуждения нейроны могут и тормозить друг друга. В случае если пресинаптический нейрон тормозный, то он выделят в синаптическую щель тормозный нейротрансмиттер открывающий хлорные каналы, а так как снаружи хлора больше, то хлор течет внутрь нейрона из-за чего отрицательный заряд на внутренней стороне мембраны увеличивается (не забываем, что ионы хлора в отличии от натрия и калия заряжены отрицательно), вгоняя нейрон в ещё более неактивное состояние. В таком состоянии нейрон труднее возбудить.
Математическая модель нейрона
На основе описанных выше динамических механизмов работы нейрона может быть составлена его математическая модель. На данный момент созданы различные как относительно простые модели, вроде «Inregrate and Fire», в которой нейрон представляется в виде конденсатора и резистора, так и более сложные, биологически правдоподобные, модели, вроде модели Ходжкина-Хаксли, которая гораздо сложнее как в вычислительном плане так и в плане анализа её динамики, но она гораздо точнее описывает динамику мембранного потенциала нейрона. В данной же статье мы будем использовать модель Ижикевича [1], она представляет из себя компромисс между вычислительной сложностью и биофизической правдоподобностью. Несмотря на свою вычислительную простоту, в этой модели можно воспроизвести большое количество явлений, происходящих в настоящих нейронах. Модель Ижикевича задается в виде системы дифференциальных уравнений (Рисунок 4).
Рис. 4 — Модель Ижикевича
Где a, b, c, d, k, Cm различные параметры нейрона. Vm — это разность потенциалов на внутренней и внешней стороне мембраны, а Um — вспомогательная переменная. I — это внешний постоянный приложенный ток. В данной модели наблюдаются такие характерные для нейронов свойства как: генерация спайка в ответ на одиночный импульса внешнего тока и генерация последовательности спайков с определённой частотой при подаче на нейрон постоянного внешнего тока. Isyn — сумма синаптических токов от всех нейронов, с которыми связан этот нейрон.
В случае если на пресинаптическом нейроне генерируется спайк, на постсинаптическом происходит скачок синапического тока, который экспоненциально затухает с характерным временем.
Переходим к кодингу
Итак, мы приступаем к самому интересному. Пора закодить на компьютере виртуальный кусок нервной ткани. Для этого будем численно решать систему дифференциальных уравнений, задающих динамику мембранного потенциала нейрона. Для интегрирования будем использовать метод Эйлера. Кодить будем на С++, рисовать с помощью скриптов написанных на Python с использованием библиотеки Matplolib, но у кого нет Питона могут рисовать с помощью Exel.
Нам понадобятся двумерные массивы Vms, Ums размерности Tsim*Nneur для хранения мембранных потенциалов и вспомогательных переменных каждого нейрона, в каждый момент времени, Tsim это время симуляции в отсчетах, а Nneur количество нейронов в сети.
Связи будем хранить в виде двух массивов pre_con и post_con размерности Ncon, где индексами является номера связей, а значениями являются индексы пресинаптических и постсинаптических нейронов. Ncon — число связей.
Так же нам понадобится массив для представления переменной, модулирующей экспоненциально затухающий постсинаптический ток каждого синапса, для этого создаем массив y размерности Ncon*Tsim.
const float h = .5f; // временной шаг интегрирования в мс
const int Tsim = 1000/.5f; // время симуляции в дискретных отсчетах
const int Nexc = 100; // Количество возбуждающих (excitatory) нейронов
const int Ninh = 25; // Количество тормозных (inhibitory) нейронов
const int Nneur = Nexc + Ninh;
const int Ncon = Nneur*Nneur*0.1f; // Количество сязей, 0.1 это вероятность связи между 2-мя случайными нейронами
float Vms[Nneur][Tsim]; // мембранные потенциалы
float Ums[Nneur][Tsim]; // вспомогательные переменные модели Ижикевича
float Iex[Nneur]; // внешний постоянный ток приложенный к нейрону
float Isyn[Nneur]; // синаптический ток на каждый нейрон
int pre_conns[Ncon]; // индексы пресинаптических нейронов
int post_conns[Ncon]; // индексы постсинаптических нейронов
float weights[Ncon]; // веса связей
float y[Ncon][Tsim]; // переменная модулирующая синаптический ток в зависимости от спайков на пресинапсе
float psc_excxpire_time = 4.0f; // характерное вермя спадания постсинаптического тока, мс
float minWeight = 50.0f; // веса, размерность пкА
float maxWeight = 100.0f;
// Параметры нейрона
float Iex_max = 40.0f; // максимальный приложенный к нейрону ток 50 пкА
float a = 0.02f;
float b = 0.5f;
float c = -40.0f; // значение мембранного потенциала до которого он сбрасываеться после спайка
float d = 100.0f;
float k = 0.5f;
float Vr = -60.0f;
float Vt = -45.0f;
float Vpeak = 35.0f; // максимальное значение мембранного потенциала, при котором происходит сброс до значения с
float V0 = -60.0f; // начальное значение для мембранного потенциала
float U0 = 0.0f; // начальное значение для вспомогательной переменной
float Cm = 50.0f; // электрическая ёмкость нейрона, размерность пкФ
Как уже было сказано, информация кодируется во временах возникновения импульсов, поэтому создаем массивы для сохранения времен их возникновения и индексов нейронов где они возникли. Далее их можно будет записать в файл, с целью визуализации.
float spike_times[Nneur*Tsim]; // времена возникновения спайков
int spike_neurons[Nneur*Tsim]; // индексы нейронов на которых происходят спайки
int spike_num = 0; // номер спайка
Разбрасываем случайно связи и задаем веса.
void init_connections(){
for (int con_idx = 0; con_idx < Ncon; ){
// случайно выбираем постсипантические и пресинаптические нейроны
pre_conns[con_idx] = rand() % Nneur;
post_conns[con_idx] = rand() % Nneur;
weights[con_idx] = (rand() % ((int)(maxWeight - minWeight)*10))/10.0f + minWeight;
if (pre_conns[con_idx] >= Nexc){
// если пресинаптический нейрон тормозный то вес связи идет со знаком минус
weights[con_idx] = -weights[con_idx];
}
con_idx++;
}
}
Установка начальных условий для нейронов и случайное задание внешнего приложенного тока. Те нейроны для которых внешний ток превысит порог генерации спайков, будут генерировать спайки с постоянной частотой.
void init_neurons(){
for (int neur_idx = 0; neur_idx < Nneur; neur_idx++){
// случайно разбрасываем приложенные токи
Iex[neur_idx] = (rand() % (int) (Iex_max*10))/10.0f;
Isyn[neur_idx] = 0.0f;
Vms[neur_idx][0] = V0;
Ums[neur_idx][0] = U0;
}
}
Основная часть программы с интегрированием модели Ижикевича.
float izhik_Vm(int neuron, int time){
return (k*(Vms[neuron][time] - Vr)*(Vms[neuron][time] - Vt) - Ums[neuron][time] + Iex[neuron] + Isyn[neuron])/Cm;
}
float izhik_Um(int neuron, int time){
return a*(b*(Vms[neuron][time] - Vr) - Ums[neuron][time]);
}
int main(){
init_connections();
init_neurons();
float expire_coeff = exp(-h/psc_excxpire_time); // для экспоненциально затухающего тока
for (int t = 1; t < Tsim; t++){
// проходим по всем нейронам
for (int neur = 0; neur < Nneur; neur++){
Vms[neur][t] = Vms[neur][t-1] + h*izhik_Vm(neur, t-1);
Ums[neur][t] = Ums[neur][t-1] + h*izhik_Um(neur, t-1);
Isyn[neur] = 0.0f;
if (Vms[neur][t-1] > Vpeak){
Vms[neur][t] = c;
Ums[neur][t] = Ums[neur][t-1] + d;
spike_times[spike_num] = t*h;
spike_neurons[spike_num] = neur;
spike_num++;
}
}
// проходим по всем связям
for (int con = 0; con < Ncon; con++){
y[con][t] = y[con][t-1]*expire_coeff;
if (Vms[pre_conns[con]][t-1] > Vpeak){
y[con][t] = 1.0f;
}
Isyn[post_conns[con]] += y[con][t]*weights[con];
}
}
save2file();
return 0;
}
Полный текст кода можно скачать тут. Код свободно компилируется и запускается хоть под виндой с Visual Studio 2010 Express Edition или MinGW, хоть под GNU/Linux c GCC. После завершения работы программа сохраняет растр активности, времена и индексы возникновения нервных импульсов, и осциллограммы среднего мембранного потенциала в файлах rastr.csv и oscill.csv, соответственно. Можно их прямо в Exelе и нарисовать. Либо с помощью прикрепленных мною питоновских скриптов, но для их работы нужна библиотека Matplotlib. Для тех у кого GNU/Linux установить из репозиториев пакет python-matplotlib не составит труда, пользователям же Windows придется вручную скачать отсюда последовательно пакеты numpy, scypy, pyparsing, python-dateutil и matplotlib.
Рисование растра активности — plot_rastr.py
Рисование среднего мембранного потенциала — plot_oscill.py
Результаты
Рис. 5 — Активность нейронной сети. Вверху по оси y отложены номера нейронов, а моменты времени когда на нейроне генерируется спайк отмечены точкой, внизу – среднее количество спайков в 1 мс времени.
Рис. 6 — Средний мембранные потенциал, «электроэнцефалограмма».
Вот что получается при моделировании 1 секунды активности 125 нейронов. Мы наблюдаем периодическую активность на частоте ~3 Гц, похожую на дельта-ритм.
[1] E.M. Izhikevich, Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting, USA, MA, Cambridge: The MIT Press., 2007
