Импульсные нейронные сети

Это введение в работу с импульсными нейронными системами с помощью пакета Julia DifferentialEquations. Мы рассмотрим три разных модели: накопления и возбуждения с утечкой, Ижикевича и Ходжкина-Хаксли. Наконец, мы ознакомимся с двумя механизмами, имитирующими синаптические входы как в реальных нейронах: альфа-синапсом и синапсом Цодыкса-Маркрама. Начнем с модели накопления и возбуждения с утечкой (LIF).

Модель накопления и возбуждения с утечками

Модель LIF (накопления и возбуждения с утечками) является расширением модели накопления и возбуждения (IF). В то время как в модели IF входные сигналы просто накапливаются, пока не произойдет возбуждение, в модели LIF входные сигналы также накапливаются, но их воздействие убывает до равновесного потенциала. Это означает, что если несколько сигналов быстро поступают один за другим, вероятность клеточного импульса гораздо выше, чем при поступлении импульсов с большими промежутками между ними. LIF является более реалистичной нейронной моделью, чем IF, так как из реальных наблюдений известно, что нейронные импульсы сильно зависят от частоты поступления входных импульсов.

В модели LIF есть пять параметров, gL, EL, C, Vth, I, и мы определяем ее в функции lif(u, p, t).

using DifferentialEquations
using Plots
gr()

function lif(u, p, t)
    gL, EL, C, Vth, I = p
    (-gL * (u - EL) + I) / C
end
lif (generic function with 1 method)

Наша система описывается одним дифференциальным уравнением: (-gL*(u-EL)+I)/C, где u — это напряжение, I — входной сигнал, gL — проводимость утечки, EL — равновесный потенциал проводимости утечки, а C — электроемкость мембраны. Как правило, любое изменение напряжения задерживается (фильтруется) электроемкостью мембраны. Вот почему мы делим все уравнение на C. При отсутствии входного сигнала извне напряжение всегда стремится к EL. Если u больше EL, u уменьшается, пока не достигнет EL. Если u меньше EL, u увеличивается, пока не достигнет EL. Единственное, что еще может повлиять на напряжение, — это внешний входной сигнал I.

Наша функция lif требует определенной структуры параметров, так как она должна быть совместима с интерфейсом DifferentialEquations. Сигнатура входных данных имеет вид lif(u, p, t), где u — это напряжение, p — коллекция параметров, описывающих уравнение, а t — время. У вас может возникнуть вопрос: почему в уравнении нет времени, хотя требуется вычислить изменение напряжения относительно времени? Дело в том, что за время отвечает решатель ОДУ. Одно из преимуществ использования решателя ОДУ вместо вычисления изменения u в цикле for заключается в том, что алгоритмы многих решателей ОДУ могут динамически регулировать временной шаг для максимальной эффективности и точности.

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

function thr(u, t, integrator)
    integrator.u > integrator.p[4]
end

function reset!(integrator)
    integrator.u = integrator.p[2]
end

threshold = DiscreteCallback(thr, reset!)
current_step = PresetTimeCallback([2, 15], integrator -> integrator.p[5] += 210.0)
cb = CallbackSet(current_step, threshold)
CallbackSet{Tuple{}, Tuple{DiscreteCallback{DiffEqCallbacks.var"#78#81"{Vector{Int64}}, DiffEqCallbacks.var"#79#82"{Main.var"#1#2"}, DiffEqCallbacks.var"#80#83"{typeof(SciMLBase.INITIALIZE_DEFAULT), Bool, Vector{Int64}, Main.var"#1#2"}, typeof(SciMLBase.FINALIZE_DEFAULT)}, DiscreteCallback{typeof(Main.thr), typeof(Main.reset!), typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(SciMLBase.FINALIZE_DEFAULT)}}}((), (DiscreteCallback{DiffEqCallbacks.var"#78#81"{Vector{Int64}}, DiffEqCallbacks.var"#79#82"{Main.var"#1#2"}, DiffEqCallbacks.var"#80#83"{typeof(SciMLBase.INITIALIZE_DEFAULT), Bool, Vector{Int64}, Main.var"#1#2"}, typeof(SciMLBase.FINALIZE_DEFAULT)}(DiffEqCallbacks.var"#78#81"{Vector{Int64}}([2, 15]), DiffEqCallbacks.var"#79#82"{Main.var"#1#2"}(Main.var"#1#2"()), DiffEqCallbacks.var"#80#83"{typeof(SciMLBase.INITIALIZE_DEFAULT), Bool, Vector{Int64}, Main.var"#1#2"}(SciMLBase.INITIALIZE_DEFAULT, true, [2, 15], Main.var"#1#2"()), SciMLBase.FINALIZE_DEFAULT, Bool[1, 1]), DiscreteCallback{typeof(Main.thr), typeof(Main.reset!), typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(SciMLBase.FINALIZE_DEFAULT)}(Main.thr, Main.reset!, SciMLBase.INITIALIZE_DEFAULT, SciMLBase.FINALIZE_DEFAULT, Bool[1, 1])))

В нашем случае условием является thr(u,t,integrator). Оно срабатывает, когда integrator.u > integrator.p[4], где p[4] — пороговый параметр Vth. В результате выполнения условия вызывается функция reset!(integrator). Она сбрасывает u в значение равновесного потенциала p[2]. Затем мы заключаем как условие, так и его результат в обратный вызов DiscreteCallback, называемый порогом. Есть еще один особенно полезный обратный вызов, который называется PresetTimeCallback. Он увеличивает значение входного параметра p[5] при t=2 и t=15 на 210.0. Оба обратных вызова затем объединяются в CallbackSet. Моделирование системы почти завершено. Осталось лишь задать числовые значения начального напряжения и параметров.

u0 = -75
tspan = (0.0, 40.0)
# p = (gL, EL, C, Vth, I)
p = [10.0, -75.0, 5.0, -55.0, 0]

prob = ODEProblem(lif, u0, tspan, p, callback = cb)
ODEProblem with uType Int64 and tType Float64. In-place: false
timespan: (0.0, 40.0)
u0: -75

Начальное напряжение будет u0 = - 75, что равно равновесному потенциалу, то есть мы начинаем с устойчивой точки. Затем мы зададим временной интервал, на котором будет производиться моделирование. Согласно определению модели LIF масштаб времени измеряется миллисекундами. Затем мы зададим следующие параметры: p = [10.0, -75.0, 5.0, -55.0, 0]. Напомним, что gL, EL, C, Vth, I = p. Наконец, мы заключаем все это в вызов ODEProblem. Не забываем о CallbackSet. На этом наша модель готова. Теперь осталось получить решение, просто вызвав solve.

sol = solve(prob)
retcode: Success
Interpolation: specialized 4th order "free" interpolation, specialized 2nd order "free" stiffness-aware interpolation
t: 153-element Vector{Float64}:
  0.0
  9.999999999999999e-5
  0.0010999999999999998
  0.011099999999999997
  0.11109999999999996
  1.1110999999999995
  2.0
  2.0
  2.6300346520855156
  2.922604932362957
  ⋮
 38.34157807913922
 38.78215044450299
 38.78215044450299
 39.22272279235999
 39.22272279235999
 39.663295181066125
 39.663295181066125
 40.0
 40.0
u: 153-element Vector{Float64}:
 -75.0
 -75.0
 -75.0
 -75.0
 -75.0
 -75.0
 -75.0
 -75.0
 -59.978080290282385
 -57.329991819757744
   ⋮
 -75.0
 -50.40489536606269
 -75.0
 -50.40489597428093
 -75.0
 -50.40489455510494
 -75.0
 -54.419257979386686
 -75.0

Во-первых, в выходных данных функции solve сообщается, получилось ли вообще решить систему. В данном случае мы видим, что получилось, так как код возврата (retcode) имеет значение Success (Успешно). Затем мы получаем числовые значения временных шагов и решений для u. По самим числам понять что-либо не очень легко, поэтому давайте построим график решения.

plot(sol)

Как мы видим, пока нет входного сигнала, модель остается на уровне -75. В момент t=2 входной сигнал увеличивается на 210 и начинается пульсация модели. Пульсация начинается не сразу, так как входной сигнал должен сначала насытить электроемкость мембраны. Обратите внимание: как только пульсация начинается, она быстро становится крайне регулярной. Когда входной сигнал снова усиливается в момент t=15, возбуждение, как и следовало ожидать, возрастает, но по-прежнему остается крайне регулярным. Это одна из особенностей LIF. Частота возбуждения является регулярной при постоянном входном сигнале и представляет собой линейную функцию от его силы. Модели LIF можно сделать менее регулярными. Например, на входе можно добавить шум определенных типов. Кроме того, можно создать множество моделей LIF и соединить их синаптическими связями. Вместо рассмотрения этих вопросов мы перейдем к модели Ижикевича, которая известна тем, что способна генерировать разнообразные кривые нейронных импульсов при постоянном воздействии.

Модель Ижикевича

Модель Ижикевича — это двухмерная модель нейронных импульсов. Она была разработана на основе анализа бифуркации кортикального нейрона. Из-за двухмерного характера эта модель позволяет генерировать гораздо более сложные кривые импульсов, чем модель LIF. Динамика импульсов зависит от четырех параметров и входного сигнала a, b, c, d, I = p. Все принципы действия такие же, как и в предыдущем случае, за исключением небольших изменений в определениях функций, которые отвечают за второе измерение.

#Модель Ижикевича
using DifferentialEquations
using Plots

function izh!(du, u, p, t)
    a, b, c, d, I = p

    du[1] = 0.04 * u[1]^2 + 5 * u[1] + 140 - u[2] + I
    du[2] = a * (b * u[1] - u[2])
end
izh! (generic function with 1 method)

Это наша модель Ижикевича. Здесь можно заметить два важных изменения. Во-первых, обратите внимание на дополнительный входной параметр du. Это последовательность разностей. du[1] соответствует напряжению (первому измерению системы), а du[2] — второму измерению. В оригинальной работе Ижикевича второе измерение обозначается как u, что усложняет понимание математической записи. В данном руководстве мы в целом будем придерживаться соглашений, принятых в Julia и DifferentialEquations, а не соглашений определенных моделей, и поэтому везде будет применяться обозначение du. Мы не будем самостоятельно определять du вне функции, но этот параметр будет использоваться внутри решателя ОДУ. Еще одно изменение — наличие знака ! после имени функции. Он означает, что параметр du будет предварительно размещаться в памяти перед интегрированием, а затем обновляться на месте, что позволяет существенно сэкономить время. Нам осталось лишь добавить обратные вызовы, которые будут обрабатывать импульсы и увеличивать входной сигнал.

function thr(u, t, integrator)
    integrator.u[1] >= 30
end

function reset!(integrator)
    integrator.u[1] = integrator.p[3]
    integrator.u[2] += integrator.p[4]
end

threshold = DiscreteCallback(thr, reset!)
current_step = PresetTimeCallback(50, integrator -> integrator.p[5] += 10)
cb = CallbackSet(current_step, threshold)
CallbackSet{Tuple{}, Tuple{DiscreteCallback{DiffEqCallbacks.var"#78#81"{Int64}, DiffEqCallbacks.var"#79#82"{Main.var"#3#4"}, DiffEqCallbacks.var"#80#83"{typeof(SciMLBase.INITIALIZE_DEFAULT), Bool, Int64, Main.var"#3#4"}, typeof(SciMLBase.FINALIZE_DEFAULT)}, DiscreteCallback{typeof(Main.thr), typeof(Main.reset!), typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(SciMLBase.FINALIZE_DEFAULT)}}}((), (DiscreteCallback{DiffEqCallbacks.var"#78#81"{Int64}, DiffEqCallbacks.var"#79#82"{Main.var"#3#4"}, DiffEqCallbacks.var"#80#83"{typeof(SciMLBase.INITIALIZE_DEFAULT), Bool, Int64, Main.var"#3#4"}, typeof(SciMLBase.FINALIZE_DEFAULT)}(DiffEqCallbacks.var"#78#81"{Int64}(50), DiffEqCallbacks.var"#79#82"{Main.var"#3#4"}(Main.var"#3#4"()), DiffEqCallbacks.var"#80#83"{typeof(SciMLBase.INITIALIZE_DEFAULT), Bool, Int64, Main.var"#3#4"}(SciMLBase.INITIALIZE_DEFAULT, true, 50, Main.var"#3#4"()), SciMLBase.FINALIZE_DEFAULT, Bool[1, 1]), DiscreteCallback{typeof(Main.thr), typeof(Main.reset!), typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(SciMLBase.FINALIZE_DEFAULT)}(Main.thr, Main.reset!, SciMLBase.INITIALIZE_DEFAULT, SciMLBase.FINALIZE_DEFAULT, Bool[1, 1])))

Одна из ключевых особенностей модели Ижикевича заключается в том, что каждый импульс увеличивает второе измерение u[2] на предварительно заданную величину p[4]. Между импульсами u[2] убывает до значения, которое зависит от p[1], p[2] и равновесного потенциала p[3]. В остальном код не слишком отличается от модели LIF. Нам снова нужно задать параметры, и моделирование можно запускать.

p = [0.02, 0.2, -50, 2, 0]
u0 = [-65, p[2] * -65]
tspan = (0.0, 300)

prob = ODEProblem(izh!, u0, tspan, p, callback = cb)
ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 300.0)
u0: 2-element Vector{Float64}:
 -65.0
 -13.0
sol = solve(prob);
plot(sol, vars = 1)

Такой тип импульсов называется «пачечной активностью». Периоды возбуждения перемежаются с паузами. Обратите внимание, что входной сигнал начинается в момент t=50 и остается постоянным на всем протяжении моделирования. Одним из механизмов, обуславливающих такой вид пульсации, является гиперполяризация, индуцированная импульсами. Ее возникновение объясняется вторым измерением, поэтому давайте рассмотрим эту переменную внимательнее.

plot(sol, vars = 2)

Второе измерение u[2] увеличивается с каждым импульсом. Когда оно становится слишком большим, система не может сгенерировать следующий импульс, пока u[2] не упадет до достаточно низкого значения, чтобы пульсация могла возобновиться. Этот процесс повторяется. В этой модели пульсация уже не является регулярной, как в случае с LIF. Здесь имеются две частоты: частота в течение периода пульсации и частота этих периодов. Модель LIF обладала одной частотой, являющейся функцией от силы входного сигнала. Давайте посмотрим, можно ли сгенерировать пульсацию другого типа, изменив параметры.

p = [0.02, 0.2, -65, 8, 0]
u0 = [-65, p[2] * -65]
tspan = (0.0, 300)

prob = ODEProblem(izh!, u0, tspan, p, callback = cb)
sol = solve(prob);
plot(sol, vars = 1)

Такой тип называется регулярной пульсацией, и мы получили его, просто уменьшив p[3] и увеличив p[4]. Обратите внимание: хотя этот тип пульсации называется регулярным, регулярный характер наступает не сразу. Кратковременное значение частоты в начале выше. Это явление называется адаптацией частоты пульсации и свойственно реальным нейронам. Генерировать можно также множество других типов пульсации. Ознакомьтесь с оригинальной работой Ижикевича и создайте нейрон по своему вкусу!

Модель Ходжкина — Хаксли

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

using DifferentialEquations
using Plots

# Функции скорости передачи по калийному ионному каналу
alpha_n(v) = (0.02 * (v - 25.0)) / (1.0 - exp((-1.0 * (v - 25.0)) / 9.0))
beta_n(v) = (-0.002 * (v - 25.0)) / (1.0 - exp((v - 25.0) / 9.0))

# Функции скорости передачи по натриевому ионному каналу
alpha_m(v) = (0.182 * (v + 35.0)) / (1.0 - exp((-1.0 * (v + 35.0)) / 9.0))
beta_m(v) = (-0.124 * (v + 35.0)) / (1.0 - exp((v + 35.0) / 9.0))

alpha_h(v) = 0.25 * exp((-1.0 * (v + 90.0)) / 12.0)
beta_h(v) = (0.25 * exp((v + 62.0) / 6.0)) / exp((v + 90.0) / 12.0)

function HH!(du, u, p, t)
    gK, gNa, gL, EK, ENa, EL, C, I = p
    v, n, m, h = u

    du[1] = (-(gK * (n^4.0) * (v - EK)) - (gNa * (m^3.0) * h * (v - ENa)) -
             (gL * (v - EL)) + I) / C
    du[2] = (alpha_n(v) * (1.0 - n)) - (beta_n(v) * n)
    du[3] = (alpha_m(v) * (1.0 - m)) - (beta_m(v) * m)
    du[4] = (alpha_h(v) * (1.0 - h)) - (beta_h(v) * h)
end
HH! (generic function with 1 method)

У нас имеются три разных вида ионной проводимости: калиевая, натриевая и проводимость утечки. Калиевая и натриевая проводимости потенциалозависимые. Они увеличиваются или уменьшаются в зависимости от напряжения. Открытые ионные каналы могут переходить в закрытое состояние, а закрытые — в открытое. Пожалуй, проще всего будет начать с калиевого тока, описываемого выражением gK * (n^4.0) * (EK - v). Здесь gK — общая проводимость, которой можно было бы достичь при всех открытых калиевых каналах. Если бы все каналы были открыты, n было бы равно 1, чего обычно не случается. Переход из открытого состояния в закрытое моделируется в alpha_n(v), а переход из закрытого состояния в открытое — в beta_n(v). Так как калиевая проводимость потенциалозависимая, эти переходы зависят от v. Числовые значения в alpha_n; beta_n были рассчитаны Ходжкиным и Хаксли на основе большого количества экспериментов с гигантским аксоном кальмара. Они также установили, что для корректного моделирования количества открытых каналов n необходимо возвести в четвертую степень.

С натриевым током ситуация почти не отличается, но у него вместо одной управляющей переменной две: m, h. У проводимости утечки gL управляющих переменных нет, так как она не потенциалозависимая. Перейдем к параметрам. Если вы хотите изучить модель ХХ во всех подробностях, отличное описание можно найти здесь.

current_step = PresetTimeCallback(100, integrator -> integrator.p[8] += 1)

# устойчивое состояние n, m и h
n_inf(v) = alpha_n(v) / (alpha_n(v) + beta_n(v))
m_inf(v) = alpha_m(v) / (alpha_m(v) + beta_m(v))
h_inf(v) = alpha_h(v) / (alpha_h(v) + beta_h(v))

p = [35.0, 40.0, 0.3, -77.0, 55.0, -65.0, 1, 0]
u0 = [-60, n_inf(-60), m_inf(-60), h_inf(-60)]
tspan = (0.0, 1000)

prob = ODEProblem(HH!, u0, tspan, p, callback = current_step)
ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 1000.0)
u0: 4-element Vector{Float64}:
 -60.0
   0.0007906538330645917
   0.08362733690208038
   0.41742979353768533

Для модели ХХ требуется только один обратный вызов — PresetTimeCallback, который активирует входной ток. Сбрасывать напряжение по достижении порога не нужно, так как в модели ХХ имеется собственный механизм реполяризации. Это калиевый ток, который активируется при высоких напряжениях и делает напряжение более отрицательным. Найти подходящие начальные значения управляющих переменных помогают три функции: n_inf; m_inf; h_inf. Эти функции сообщают нам, что устойчивые управляющие значения требуются для начального напряжения. Параметры были выбраны так, чтобы модель приблизительно представляла кортикальную пирамидальную клетку, а не гигантский аксон, с которым изначально работали Ходжкин и Хаксли.

sol = solve(prob);
plot(sol, vars = 1)

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

plot(sol, vars = [2, 3, 4], tspan = (105.0, 130.0))