Origin of heterogeneous spiking patterns in spinal dorsal horn neurons (Balachandar & Prescott 2018)

 Download zip file 
Help downloading and running models
Accession:256628
"Neurons are often classified by spiking pattern. Yet, some neurons exhibit distinct patterns under subtly different test conditions, which suggests that they operate near an abrupt transition, or bifurcation. A set of such neurons may exhibit heterogeneous spiking patterns not because of qualitative differences in which ion channels they express, but rather because quantitative differences in expression levels cause neurons to operate on opposite sides of a bifurcation. Neurons in the spinal dorsal horn, for example, respond to somatic current injection with patterns that include tonic, single, gap, delayed and reluctant spiking. It is unclear whether these patterns reflect five cell populations (defined by distinct ion channel expression patterns), heterogeneity within a single population, or some combination thereof. We reproduced all five spiking patterns in a computational model by varying the densities of a low-threshold (KV1-type) potassium conductance and an inactivating (A-type) potassium conductance and found that single, gap, delayed and reluctant spiking arise when the joint probability distribution of those channel densities spans two intersecting bifurcations that divide the parameter space into quadrants, each associated with a different spiking pattern. ... "
Reference:
1 . Balachandar A, Prescott SA (2018) Origin of heterogeneous spiking patterns from continuously distributed ion channel densities: a computational study in spinal dorsal horn neurons. J Physiol 596:1681-1697 [PubMed]
Model Information (Click on a link to find other models with that property)
Model Type:
Brain Region(s)/Organism:
Cell Type(s): Abstract Morris-Lecar neuron; Dorsal Root Ganglion (DRG) cell;
Channel(s): I Na,t; I K; I A; I Potassium;
Gap Junctions:
Receptor(s):
Gene(s):
Transmitter(s):
Simulation Environment: MATLAB;
Model Concept(s): Bifurcation; Activity Patterns; Parameter Fitting; Action Potential Initiation; Action Potentials; Parameter sensitivity;
Implementer(s):
Search NeuronDB for information about:  I Na,t; I A; I K; I Potassium;
% modified_morris_lecar.m
% Arjun Balachandar 2016
% Adapted from Husain Shakil and PLOS comp bio (Prescott & Sejnowski 2008)
%   This program simulates the Morris Lecar model, modified to include
%   Low-Threshold K+ and A-Type K+ currents.
%   The A-Type K+ current is modelled by a modified Connor-Stevens model

%This file is called by Simulate.m, which is in turn called by AutoSim.m
%for each combination of input parameters specified by the user in AutoSim.

%[V,currents,conductances,w,z,a,b,spike,numAPs,t,ATP,ATP_spike]
function [V,currents,conductances,spike,numAPs,t] = modified_morris_lecar(i_stim, i_off, g_sub, gA, time)
    tspan = time*1000; %time is in seconds, tspan is in milliseconds
    dt = 0.1;
    loop = ceil(tspan/dt); % no. of iterations of euler

    %Initialize constants:
    C = 2;
    numAPs = 0;
    %Fast Inward Current:
    gNa = 20;
    ENa = 50;
    %Delayed Rectifier Current (IKdr or recovery variable w)
    gK = 20;
    EK = -100;
    phi_w = 0.15;
    %Slow subthreshold Potassium Current
    phi_z = 0.15;
    %Leak Current (Il):
    gl = 2;
    El = -70;
    %AHP Current:
    g_ahp = 0;
    

    %Allocate memory for variable vectors:
    t = (1:loop)*dt;
    V = zeros(loop,1);
    w = zeros(loop,1);
    z = zeros(loop,1);
    a = zeros(loop,1);
    b = zeros(loop,1);
    q = zeros(loop,1);
    currents = zeros(loop,5);
    conductances = zeros(loop,2);
    INa = zeros(loop,1);
    IK = zeros(loop,1);
    IgA = zeros(loop,1);
    Igsub = zeros(loop,1);
    spike = zeros(loop,1);
    ref = zeros(loop,1);
    
    %gahp2*q2*(v-Vk)

    %Calculate steady-state initial conditions (current = 0) using Euler method:
    for i=1:loop-1 %loop/4
          %dV/dt = (i_stim-gna*minf(V)*(V-Vna)-gk*w*(V-VK)-gl*(V-Vl)-gsub*z*(V-VK)-gA*a^4*b*(V-VK))/c
          dV_dt = (i_off - gNa*m_inf(V(i))*(V(i)-ENa) - gK*w(i)*(V(i)-EK) - gl*(V(i)-El) - gA*((a(i))^4)*b(i)*(V(i)-EK) - g_sub*z(i)*(V(i) - EK) - g_ahp*q(i)*(V(i) - EK))/C;
          V(i+1) = V(i) + dt*dV_dt; %forward Euler equation

          dw_dt = phi_w*(w_inf(V(i))- w(i))/tau_w(V(i));
          w(i+1) = w(i) + dt*dw_dt; %forward Euler equation

          dz_dt = phi_z*(z_inf(V(i)) - z(i))/tau_z(V(i));
          z(i+1) = z(i) + dt*dz_dt; %forward Euler equation
          
          da_dt = (a_inf(V(i)) - a(i))/tau_a(V(i));
          a(i+1) = a(i) + dt*da_dt; %forward Euler equation
          
          db_dt = (b_inf(V(i)) - b(i))/tau_b(V(i));
          b(i+1) = b(i) + dt*db_dt; %forward Euler equation
          
          %AHP: (not used)
          dq_dt = (q_inf(V(i)) - q(i))/tau_q(V(i));
          q(i+1) = q(i) + dt*dq_dt;
          
          INa = gNa*m_inf(V(i))*(V(i)-ENa);
          IK = gK*w(i)*(V(i)-EK);
          IgA = gA*((a(i))^4)*b(i)*(V(i)-EK);
          Igsub = g_sub*z(i)*(V(i) - EK);
          Igahp = g_ahp*q(i)*(V(i) - EK);
          Inet = INa + IK + IgA + Igsub + Igahp + gl*(V(i)-El);
          currents(i,:) = [INa, IK, IgA, Igsub, Inet];
          
          gA_current = gA*((a(i))^4)*b(i);
          gsub_current = g_sub*z(i);
          conductances(i,:) = [gA_current, gsub_current];
          

          spike(i) = (V(i) > 0).*(~ref(i));
          ref(i+1) = (V(i) > 0);
          
    end
    
    % Set initial-conditions
    V(1) = V(loop-2); %V(loop/2 - 1)
    V_rest = V(1);
    display(V_rest);
    w(1) = w(loop-2);
    z(1) = z(loop-2);
    a(1) = a(loop-2);
    b(1) = b(loop-2);
    q(1) = q(loop-2);
    
    conductances(1,:) = conductances(loop-3,:);
    currents(1,:) = currents(loop-3,:);
    
    spike(1) = spike(loop-2); % =0 %loop -2
    ref(1) = ref(loop-2); % =0
    
    
    %Euler method:
    for i=1:loop-1
    %for i=(loop/6+1):(loop-1)
          if i>loop*0.8
              i_stim = 0.0;
          end
          %dV/dt = (i_stim-gna*minf(V)*(V-Vna)-gk*w*(V-VK)-gl*(V-Vl)-gsub*z*(V-VK)-gA*a^4*b*(V-VK))/c
          dV_dt = (i_off + i_stim - gNa*m_inf(V(i))*(V(i)-ENa) - gK*w(i)*(V(i)-EK) - gl*(V(i)-El) - gA*((a(i))^4)*b(i)*(V(i)-EK) - g_sub*z(i)*(V(i) - EK)  - g_ahp*q(i)*(V(i) - EK))/C;
          V(i+1) = V(i) + dt*dV_dt; %forward Euler equation

          dw_dt = phi_w*(w_inf(V(i))- w(i))/tau_w(V(i));
          w(i+1) = w(i) + dt*dw_dt; %forward Euler equation

          dz_dt = phi_z*(z_inf(V(i)) - z(i))/tau_z(V(i));
          z(i+1) = z(i) + dt*dz_dt; %forward Euler equation
          
          da_dt = (a_inf(V(i)) - a(i))/tau_a(V(i));
          a(i+1) = a(i) + dt*da_dt; %forward Euler equation
          
          db_dt = (b_inf(V(i)) - b(i))/tau_b(V(i));
          b(i+1) = b(i) + dt*db_dt; %forward Euler equation
          
          %AHP:
          dq_dt = (q_inf(V(i)) - q(i))/tau_q(V(i));
          q(i+1) = q(i) + dt*dq_dt;
          
          INa = gNa*m_inf(V(i))*(V(i)-ENa);
          IK = gK*w(i)*(V(i)-EK);
          IgA = gA*((a(i))^4)*b(i)*(V(i)-EK);
          Igahp = g_ahp*q(i)*(V(i) - EK);
          Inet = INa + IK + IgA + Igsub + Igahp + gl*(V(i)-El);
          
          gA_current = gA*((a(i))^4)*b(i);
          gsub_current = g_sub*z(i);
          
          conductances(i,:) = [gA_current, gsub_current];
          currents(i,:) = [INa, IK, IgA, Igsub, Inet];
              

          spike(i) = (V(i) > 0).*(~ref(i));
          ref(i+1) = (V(i) > 0);
          
          if spike(i) > 0
              numAPs = numAPs + 1;
          end
    end
    
end


%Steady-state and Decay functions for gating variables:
%% Fast Sodium Channel:
function [minf] = m_inf(V)
    beta_m = -1.2;
    gamma_m = 18;
    minf = 0.5*(1+tanh((V-beta_m)/gamma_m));
end

%% Delayed Rectifier Potassium Channel:
function [winf] = w_inf(V)
    beta_w = -10;
    gamma_w = 10;
    winf = 0.5*(1+tanh((V-beta_w)/gamma_w));
end

function [tauw] = tau_w(V)
    beta_w = -10;
    gamma_w = 10;
    tauw = 1/cosh((V-beta_w)/(2*gamma_w));
end

%% AHP Potassium Channel
function [qinf] = q_inf(V)
    vhalfq = 0;
    vslopeq = -5;
    qinf = 1/(1+exp((V-vhalfq)/vslopeq));
end

function [tauq] = tau_q(V)
    tauq = 20;
end
%           q2_inf = 1/(1+exp((v-vhalfq2)/vslopeq2))
%             param gahp2=50
%          param vhalfq2=0
%             param vslopeq2=-5
%             param tau_q2=200
%             q2(0)=0

%% Slow Threshold Potassium Current (I_sub):
function [zinf] = z_inf(V)
    beta_z = -21;
    gamma_z = 15;
    zinf = 0.5*(1+tanh((V-beta_z)/gamma_z));
end

function [tauz] = tau_z(V)
    beta_z = -21;
    gamma_z = 15;
    tauz = 1/cosh((V-beta_z)/(2*gamma_z));
end

%% A-Type Potassium Current (I_A) - Modified Connor-Stevens model:
function [ainf] = a_inf(V)
    ainf = 1/(1+exp(-(V+60)/8.5));
end

function [binf] = b_inf(V)
    binf = 1/(1 + exp((V+78)/6));
end

function [taua] = tau_a(V)
    taua = 1/(exp((V + 35.82)/19.69) + exp(-(V + 79.69)/12.7)) + 0.37;
end

function [taub] = tau_b(V)
    r = V+63;
    s = -V-63;
    taub = 19*heav(r)+(1/((exp((V+46.05)/ 5)+exp(-(V+238.4)/37.45))))*heav(s);
end

%%
function [hs] = heav(x)
    if x > 0
        hs = 1;
    elseif x < 0
        hs = 0;
    else
        hs = 0.5;
    end
end
    











Loading data, please wait...