📂Dynamics

Memristive Hindmarsh-Rose Neuron Model as a Dynamical System

Overview

The Hindmarsh-Rose model is a model that explains the changes in voltage related to ion channels in neurons. The term memristor is a concept combining memory and resistor. This document describes a system that adds a memristor to the Hindmarsh-Rose model.

Model ¹

The dimensionless system of the memristor-based Hindmarsh-Rose circuit and the equivalent differential equations can be represented as the following non-autonomous non-smooth system. $$ \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*} $$

Here, $I = f \cos \left( \omega t \right)$ denotes the external current, and $g(z)$ is a function defined as follows: $$ \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*} $$

Variables²

• $x(t)$: This is the variable associated with membrane potential at the point $t$.
• $y(t)$: This is the spiking variable related to the movement of sodium and potassium ion channels at the point $t$.
• $z(t)$: This is the variable related to the adaptation current at the point $t$.

Parameters

• $a=1$, $b=3$, $c=1$, $d=5$: These values are related to resistance.
• $\omega = 1$: This is the period of the external current.
• $f = 0.1$: This is the magnitude of the external current.
• $k = 0.9$: This denotes the coupling strength.
• $\alpha = 0.1$, $\beta = 0.8$: These are memristor synapse parameters.

Explanation

The memristor Hindmarsh-Rose neuron model, hereinafter referred to as the HR model, is characterized by having three subsystems, dividing the domain into two planes $z = \pm 1$ as a non-smooth system.

Bifurcation

If one sets $f$ as the bifurcation parameter and plots a bifurcation diagram, as shown above, it demonstrates characteristic behaviors at values below $f^{\ast} \approx 0.21$.

Code

Full Code

The following is the code to reproduce the bifurcation diagram of the HR model in Julia. Collect the values of $x$ at the points where the subsystems change after sufficient preiteration, and plot them on a scatter plot.

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)