📂동역학

동역학계로써의 멤리스터 힌드마쉬-로즈 뉴런 모델

개요

힌드마쉬-로즈 모델^{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$ 를 바이퍼케이션 파라미터으로 잡고 바이퍼케이션 다이어그램을 그려보면 위와 같이 $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)