📂動力学

ダイナミカルシステムとしてのメムリスターヒンドマーシュ・ローズニューロンモデル

概要

ヒンドマッシュ-ローズモデル^{Hindmarsh-Rose model}はニューロン^neuronにおけるイオンチャネル^{ion channel}に関連した電圧の変化を説明するモデルである。メムリスター^memristerとはメモリー^memoryとレジスター^resisterを合わせた概念の素子として、ヒンドマッシュ-ローズモデルにメムリスターを追加したシステムについて説明する。

モデル ¹

メムリスター基盤のヒンドマッシュ-ローズ 回路と等価な 微分方程式の無次元^{dimensionless}システムは以下のような 非自律的 ノンスムースシステムで表すことができる。 $$ \begin{align*} \dot{x} =& y - ax^{3} + bx^{2} + kxz + I \\ \dot{y} =& c - dx^{2} - y \\ \dot{z} =& \alpha g (z) + \beta x \end{align*} $$

ここで $I = f \cos \left( \omega t \right)$ は外部電流^{external current}であり、$g(z)$ は以下のように定義された 関数である。 $$ \begin{align*} g(z) =& \sign \left( z + 1 \right) + \sign \left( z - 1 \right) - z \\ =& \begin{cases} -2 -z & , \text{if } z < -1 \\ -z & , \text{if } -1 \le z \le 1 \\ 2 - z & , \text{if } z > 1 \end{cases} \end{align*} $$

変数 ²

• $x(t)$: $t$ 時点で膜電位^{membrane potential}に関連した変数である。
• $y(t)$: $t$ 時点でイオンチャネルのソディウム^sodiumとポタッシウム^potassiumの移動に関連したスパイキング^spiking変数である。
• $z(t)$: $t$ 時点で適応電流^{adptation current}に関連した変数である。

パラメータ

• $a=1$, $b=3$, $c=1$, $d=5$: 抵抗に関連する値である。
• $\omega = 1$: 外部電流の周期である。
• $f = 0.1$: 外部電流の大きさである。
• $k = 0.9$: カップリングの強度^{coupling strength}を示す。
• $\alpha = 0.1$, $\beta = 0.8$: メムリスターシナプス^{memristor synapse}パラメータである。

説明

メムリスターヒンドマッシュ-ローズニューロンモデル、以下HRモデルは上記のようにドメインを二つの 平面 $z = \pm 1$ で分けるノンスムースシステムで、三つのサブシステム^subsystemを持つのが特徴である。

バイフィケーション

$f$ を バイフィケーションパラメータとしてバイフィケーションダイアグラム](../1006)を描くと、上記のように $f^{\ast} \approx 0.21$ 以下の値で カオティックな行動を示す。

コード

全体コード

以下は ジュリアでHRモデルのバイフィケーションダイアグラムを再現するコードである。十分な事前反復^preiteration後にサブシステムが変わる時点の $x$ 値を集めて、点図を描くと良い。

using CSV, DataFrames, ProgressMeter

function RK4(f::Function, v::AbstractVector, h=1e-2, nonsmooth=0.0)
    V1 = f(v, nonsmooth)
    V2 = f(v + (h/2)*V1, nonsmooth)
    V3 = f(v + (h/2)*V2, nonsmooth)
    V4 = f(v + h*V3, nonsmooth)
    return v + (h/6)*(V1 + 2V2 + 2V3 + V4), V1
end

const _a = 1.0
const _b = 3.0
const _c = 1.0
const _d = 5.0
const _k = 0.9
# const _f = 0.1
const _ω = 1.0
const _α = 0.1
const _β = 0.8 
function factory_hrnm(idx::Int64, _f::Number; ic = [0.0, 0.0, 0.0, 0.1], tspan = [0, 100], dt = 1e-3)
    function HR(txyz::AbstractVector, nonsmooth::Real)
        t,x,y,z=txyz
    
        ṫ = 1
        ẋ = y - _a*x^3 + _b*x^2 + _k*x*z + _f*cos(_ω*t)
        ẏ = _c - _d*x^2 - y
        ż = _α*nonsmooth + _β*x
        return [ṫ, ẋ, ẏ, ż]
    end

    
    t_ = first(tspan):dt:last(tspan)

    ndatapoints = count(first(tspan) .< t_ .≤ last(tspan))

    len_t_ = length(t_)
    x = ic; DIM = length(x)
    traj = zeros(2DIM, len_t_+1)
    traj[1:DIM, 1] = x

    
    for tk in 1:len_t_
        nonsmooth = sign(x[4]+1) + sign(x[4]-1) - x[4]
        # if x[4] < -1
        #     nonsmooth = -2 -x[4]
        # elseif x[4] > 1
        #     nonsmooth = 2 - x[4]
        # end
        x, dx = RK4(HR, x, dt, nonsmooth)
        if tk+1 ≥ (len_t_ - ndatapoints)
            traj[        1:DIM , tk+1] =  x
            traj[DIM .+ (1:DIM), tk  ] = dx
        end
    end
    traj = traj[:, (end-ndatapoints):(end-1)]'

    return traj
end
factory_hrnm(T::Type, args...; ic = [0.0, 0.0, 0.0, 0.1], tspan = [0, 100], dt = 1e-3) = 
DataFrame(factory_hrnm(args...;  ic, tspan, dt), ["t", "x", "y", "z", "dt", "dx", "dy", "dz"])

hrzn = []
vrtc = []
@showprogress @threads for dr = eachrow(schedules)
        filename = "bifurcation/hrnm/$(lpad(dr.idx, 5, '0')).csv"
        data = factory_hrnm(DataFrame, dr.idx, dr.f, tspan = [0, 1500]); data = data[1000(nrow(data) ÷ 1500):end , :]
        # CSV.write(filename, data)

        idx_sampled = abs.(diff(data.dz)) .> 0.1
        sampledx = data[Not(1), :x][idx_sampled]
        append!(hrzn, fill(dr.f, length(sampledx)))
        append!(vrtc, sampledx)
end
scatter(hrzn, vrtc, ms = 1, legend = :none, msw = 0, ma = 0.1)