logo

Memristive Hindmarsh-Rose Neuron Model as a Dynamical System 📂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 1

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. x˙=yax3+bx2+kxz+Iy˙=cdx2yz˙=αg(z)+βx \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=fcos(ωt)I = f \cos \left( \omega t \right) denotes the external current, and g(z)g(z) is a function defined as follows: g(z)=sign(z+1)+sign(z1)z={2z,if z<1z,if 1z12z,if z>1 \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*}

Variables2

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

Parameters

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

Explanation

alt text

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=±1z = \pm 1 as a non-smooth system.

Bifurcation

alt text

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

Code

Full Code

alt text

The following is the code to reproduce the bifurcation diagram of the HR model in Julia. Collect the values of xx 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)

  1. Fuhong Min, Zhi Rui; Boundary dynamics of a non-smooth memristive Hindmarsh–Rose neuron system. Chaos 1 October 2022; 32 (10): 103117. https://doi.org/10.1063/5.0107067 ↩︎

  2. https://en.wikipedia.org/wiki/Hindmarsh%E2%80%93Rose_model ↩︎