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]
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;
%Arjun Balachandar 2017
%This program takes the observed firing-pattern of a given neuron at
%various stimulation intensities i_stim, and uses that in conjunction with the
%generated parameter-space slices to determine the possible underlying
%conductance values of gA and gsub characterizing the neuron.
%Speficially, this is done by loading 2 data sets containting 2 ranges of
%values of i_stim (and the FP of the neurons at each slice of i_stim for each combination of gA and gsub) and
%finding areas (conductance combinations) of the parameter-space where the neuron matches both sets of data.
%If no such area found, the algorithm returns no area.
%(See methods).

%***USER MODIFIABLE- refers to variables to be changed by user****

clear all;

%load simulation data

grid_size = 20;

xmax = grid_size;
xmin = 0;
ymax = grid_size;
ymin = 0;

del = 0.1;
dx = del;
dy = del;

num_x = (xmax-xmin)/dx + 1;
num_y = (ymax-ymin)/dy + 1;
FP_domain = zeros(num_istim,num_y,num_x);

best_matches1 = zeros(num_x*num_y,2);
best_matches2 = zeros(num_x*num_y,2);
num_best1 = 0; num_best2 = 0;
perc_match_best1 = 0; perc_match_best2 = 0;

%***USER MODIFIABLE: user inputs the FP of the neuron at each level of i_stim
%(see line below), which the algorithm uses to find areas of best match.
%I_stim order: [50,55,60,65,70,75,80]; from the cumulative data loaded from both data
%FP numbering system: [R, S, D, G, T] = [0, 1, 2, 3, 4] (i.e. [reluctant-,
%single, delayed-, gap- and tonic-spiking] )
target_tot = [0,0,2,3,3,3,3];
%Other example combinations of FPs for a neuron:
%[0,0,0,1,1,1,1]; ; %[0,2,2,2,3,3,3]; 
%[0,0,1,1,3,3,4]; [0,0,1,3,3,3,4]; [0,1,1,1,3,3,4]; [0,0,0,1,1,3,3];
%[0,2,2,3,3,3,3]; [0,2,3,3,4,4,4]; %[0,2,3,3,3,4,4]; [0,2,3,3,3,3,4]; [0,2,3,3,3,3,3];
target = target_tot(1:4); %take first part of target_tot initially (to compare to first data-set), then second part when comapring to second data-set

for t=1:2
    if t==2 %Load data for higher i_stims
        load('AutoSim_distim5_ioff0_dgA0'); %Contains i=70
        target = target_tot(5:7);
    values = zeros(num_istim);
    values = values(:);
    %value above which perc_match must be in order to be considered as a match:
    match_thres = (num_istim-1)/num_istim;
    if match_thres < 0.7
        match_thres = 0.7;
    ind = 1;
    %location of best match (OUTPUT OF THE ALGORITHM):
    perc_match_best = 0.0;
    x_best = 0.0; y_best = 0.0;
    num_xTot = (max_gsub-min_gsub)/dx + 1; %for modified domain size
    num_x = (xmax-xmin)/dx + 1;
    num_y = (ymax-ymin)/dy + 1;
    x_domain = linspace(xmin,xmax,num_x);
    y_domain = linspace(ymin,ymax,num_y);
    [X, Y] = meshgrid(xmin:dx:xmax, ymin:dy:ymax);
    FP_domain_slice = zeros(num_y,num_x);
    match_level = zeros(num_y,num_x);
    slice_size = num_x*num_y;
    param_array_slice = zeros(slice_size, 6 + nSpikeATPs);
    for k=1:num_istim
        %calculate each slice from i_stim = min_istim to i_stim max_istim,
        %and plot them for the user.
        value = min_istim + d_istim*(k-1);
        cur_ind = 1;
        for i=1:(num_istim)*(num_x)*(num_y)
            if param_array(i,ind) == value %current i_stim value
                for j=1:(6+nSpikeATPs)
                    param_array_slice(cur_ind,j) = param_array(i,j);
                cur_ind = cur_ind + 1;
        for i=0:num_x-1
            for j=0:num_y-1
                FP_domain(k+4*(t-1),j+1,i+1) = param_array_slice(j*num_xTot + i + 1,4);
                FP_domain_slice(j+1,i+1) = param_array_slice(j*num_xTot + i + 1,4);
        x_size = num_x;
        y_size = num_y;
        gA_axis = linspace(ymin,ymax,num_y);
        gsub_axis = linspace(xmin,xmax,num_x);
        [aC aC] = contourf(FP_domain_slice.',5);
        title(['istim = ',int2str(value),', ioff = ',int2str(i_off)]);
        set(gca, 'XTick', 1:x_size); % Change x-axis ticks
        set(gca, 'XTickLabel', gsub_axis); % Change x-axis ticks labels.
        set(gca, 'YTick', 1:y_size); % Change x-axis ticks
        set(gca, 'YTickLabel', gA_axis); % Change x-axis ticks labels.
        FP_domain_slice = [];
        FP_domain_slice = zeros(num_y,num_x);
    num_best = 0;
    best_matches = [];
    %finds areas (combinations of cnductances) such that the neuron
    %firing-pattern at each i_stim (for that combination of conductances)
    %matches the FP seen in the parameter-space at that combination of
    for i=1:num_x
        for j=1:num_y
            perc_match = 0;
            for k=1:num_istim
                if target(k) == FP_domain(k+4*(t-1),j,i)
                    perc_match = perc_match + 1;
            perc_match = perc_match/num_istim;
            if perc_match > perc_match_best && perc_match >= match_thres
                perc_match_best = perc_match;
                num_best = 1;
                x_best = i;
                y_best = j;
                best_matches = [];
                best_matches = zeros(num_x*num_y,2);
                best_matches(num_best,1) = x_best;
                best_matches(num_best,2) = y_best;
            elseif perc_match == perc_match_best && perc_match >= match_thres
                x_best = i;
                y_best = j;
                num_best = num_best + 1;
                best_matches(num_best,1) = x_best;
                best_matches(num_best,2) = y_best;
            match_level(j,i) = perc_match;
    if t==1
        num_best1 = num_best;
        best_matches1 = best_matches;
        perc_match_best1 = perc_match_best;
        num_best2 = num_best;
        best_matches2 = best_matches;
        perc_match_best2 = perc_match_best;


%areas where the areas of best fit in both i_stim sections match
best_matches_final = intersect(best_matches1,best_matches2,'rows');
best_matches_final1 = zeros(num_best1,2);
best_matches_final2 = zeros(num_best2,2);

perc_match_best = (perc_match_best1*4 + perc_match_best2*3)/7.0;

num_most = max(num_best1,num_best2);
for i=1:num_best1
    best_matches_final1(i,1) = xmin + (best_matches1(i,1)-1)*dx;
    best_matches_final1(i,2) = ymin + (best_matches1(i,2)-1)*dy;
x_size = num_x;
y_size = num_y;
gA_axis = linspace(ymin,ymax,num_y);
gsub_axis = linspace(xmin,xmax,num_x);

%Display the areas of best fit in section 1, section 2, and both section 1
%and 2 (i.e. the real areas of best fit)

figure('name','Area of best fit1');
scat = scatter(best_matches_final1(:,2),best_matches_final1(:,1),10,'filled');

xlim([xmin xmax]);
ylim([ymin ymax]);

for i=1:num_best2
    best_matches_final2(i,1) = xmin + (best_matches2(i,1)-1)*dx;
    best_matches_final2(i,2) = ymin + (best_matches2(i,2)-1)*dy;

figure('name','Area of best fit2');
scat = scatter(best_matches_final2(:,2),best_matches_final2(:,1),10,'filled');

xlim([xmin xmax]);
ylim([ymin ymax]);

if size(best_matches_final) > 0
    perc_match_best = (perc_match_best1*4 + perc_match_best2*3)/7.0;
    perc_match_best = 0.0

for i=1:size(best_matches_final)
    best_matches_final(i,1) = xmin + (best_matches_final(i,1)-1)*dx;
    best_matches_final(i,2) = ymin + (best_matches_final(i,2)-1)*dy;

figure('name','Area of best fit - Total');
scat = scatter(best_matches_final(:,2),best_matches_final(:,1),10,'filled');

xlim([xmin xmax]);
ylim([ymin ymax]);

Loading data, please wait...