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

 Download zip file 
Help downloading and running models
"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. ... "
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]
Citations  Citation Browser
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:
Simulation Environment: MATLAB;
Model Concept(s): Bifurcation; Activity Patterns; Parameter Fitting; Action Potential Initiation; Action Potentials; Parameter sensitivity;
Search NeuronDB for information about:  I Na,t; I A; I K; I Potassium;
% Simulate.m
% Arjun Balachandar 2016
% Simulate neuron: inputs and monitor stimulation parameters

%AutoSim.m calls this file, inputting parameters i (i_stim), i_off
%(pre-stim stimulation to alter resting Vm), time of stimulation,
%(ATP-related input parameters not used here) and if this specific
%neuron-model is one of many, or is the only one being simulated by AutoSim
%(i.e. multipleNeurons: if 0, then display single-neuron-specific data,
%e.g. V-t traces)
function [firing_info] = Simulate(i,i_off,gsub,g_A, tim, multipleNeurons)
    check = @modified_morris_lecar; %uses modified_morris_lecar model from 

    %initialize variables
    type = 1;
    type_str = 'R'; %stores firing-pattern outputted by model neuron after simulation
    i_stim = i;
    gA = g_A;
    g_sub = gsub;
    time = tim;
    i_offset = i_off;
    %run modified morris lecar model using input parameters, and store
    %output variables (V-t, currents-t, conductances-t data etc..)
    [V,currents,conductances,spike,numAPs,t] = check(i_stim, i_offset, g_sub, gA, time);

    %I-t data after simulation (i.e. current through each channel)
    INa = currents(:,1);
    IK = currents(:,2);
    IgA = currents(:,3);
    Igsub = currents(:,4);
    Inet = currents(:,5);
    %instantaneous (i.e. current) value of g (conductance) at each
    gA_current = conductances(:,1);
    gsub_current = conductances(:,2);

    %%Determine Firing Pattern
    spike_times = zeros(numAPs,1);
    curAP = 1; %index of AP currently being analyzed
    for i=1:length(spike) %spike stores spike-time of output response after simulation
        if spike(i) > 0
            spike_times(curAP) = t(i);
            curAP = curAP + 1;

    rate = numAPs/time; %calculate firing-rate
    %Classify neuron firing-pattern type
    %type corresponds to numrical classification system, while type_str is the
    %corresponding string-name
    if numAPs == 1
        type = 1;
        type_str = 'SS'; %if one spike, then single-spike neuron (SS)
    elseif numAPs < 1
        type_str = 'R'; %if no spikes, than reluctanct/non-spiking neuron (R)
        type = 0;
    elseif numAPs >= 3 %if more than 2 spikes, than can be of delayed-onset, gap- (Gap) or repetitive- (R - i.e. tonic in paper) spiking
        ISI2 = spike_times(3) - spike_times(2); %inter-spike interval between spikes 3 and 2
        ISI1 = spike_times(2) - spike_times(1); %ISI between t=0 and spike 1
        if 0.5*spike_times(2) > spike_times(1) %if the first spike does not occur late enough, then not considered a delayed-onset neuron
            %^specifically, if 0.5x the time for the second spike occurs after
            %the first spike occurs, then first spike did not occur late
            %enough to be 'delayed-onset'
            if ISI1 > 1.5*ISI2
                type_str = 'Gap';
                %If not delayed-onset, and if time between first and second
                %spike is more than 1.5x the time between 2nd and 3rd spikes,
                %then the initial gap between 1st and 2nd spikes is
                %considered long enough to classify neuron as 'gap-spiking'
                type = 3;
                type_str = 'RF'; %if not a gap neuron, then considered RF
                type = 4;
            type_str = 'DO'; %delayed-onset'
            type = 2;
    %store all pertinent information in firing_info, for eventual use when
    %analyzing simulation data (Sim_Analysis.m)
    firing_info = [];
    firing_info(1) = type;
    firing_info(2) = rate;
    if multipleNeurons == 0
        %Plots --> if this neuron-model is the only one being simulated by
        %AutoSim.m, then display all pertinent output data to user
        %Conductance plots:
        figure('name','gA-t Plot');
        figure('name','gA-t (normalized) Plot');
        axis([0 tim*1000 0 1]);
        figure('name','gsub-t Plot');
        figure('name','gsub-t (normalized) Plot');
        axis([0 tim*1000 0 1]);
        figure('name','Conductances Plot');
        %Plot showing activation of each channel (normalized to the maximum
        %value of g)
        figure('name','Conductances (normalized) Plot');
        axis([0 tim*1000 0 1]);
        figure('name','gA-V plot');
        %Plot showing (semi-real-time, to give impression of dynamic plot)
        %how the instantaneous conductance gA corresponds to V, over time:
%         figure('name','gA-V plot (real-time)');
%         hold on;
%         for k = 1:size(V)-1
%             plot(V(k:k+1),gA_current(k:k+1)*(1.0/gA));
%             pause(0.0001);
%         end
%         hold off;
        %V-t plots:
        figure('name','V-t Plot');
        axis([0 tim*1000 -80 50]);
        file_type = 'V-t_';