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;
%Sim_Analysis.m
%Arjun Balachandar 2016
%Display and analyze data generated by running AutoSim.m, which simulates
%neuron(s) with various input parameters (see AutoSim.m). Used to visualize
%simulation data and create plots/figures

clc;
clear all;
close all;

%Load files generated after running AutoSim.m simulations
%'_ioff' refers to pre-stimulus current used to change resting Vm. If 0,
%then no pre-stimulus current used.

%Various files to load
%load('AutoSim_distim1_ioff0.mat');
%load('AutoSim_istim085_distim5_ioff0_dgA0');
%load('AutoSim_istim070_distim5_ioff0_dgA0_maxgsub15');
%load('AutoSim_distim5_ioff0_dgA0'); % ****------ CONTAINS iSTIM70 -------****
%load('AutoSim_distim5_ioff0.mat');
load('AutoSim_istim060_distim10_ioff20_dgA0_maxgsub15'); %i=60, i_off=0-20, 15x15
planeBF = @PlaneBF;

%Variable to keep constant
ind = 1; %take slice through i_stim = set-value
if num_istim > 1 %if many i_stim used by AutoSim, then will create a slice through data for each i_stim (see below)
    value = 0;
else
    value = max_istim; %if only one level of stimulus i_stim, then use this as the 'value' for creating slice (see below)
end

dx = d_gA;
dy = d_gsub;

grid_size = 15; %Width of grid (i.e. gA-gsub plot square length)

%Dimensions of grids for plots produced later
xmax = grid_size;
xmin = 0;
ymax = grid_size;
ymin = 0;

num_xTot = (max_gsub-min_gsub)/dx + 1; %for modified domain size (i.e. If grid_size is less than the maximum conductance value used (smaller plot
%than present in data)
num_x = (xmax-xmin)/dx + 1;
num_y = (ymax-ymin)/dy + 1;

%set axes
gA_axis = linspace(ymin,ymax,num_y);
gsub_axis = linspace(xmin,xmax,num_x);

%Plot Colors:
rate_c_mod = zeros((num_istim)*(num_x)*(num_y),3);
type_c_mod = zeros((num_istim)*(num_x)*(num_y),3);
rate_c = zeros((num_istim)*(num_gA)*(num_gsub),3);
type_c = zeros((num_istim)*(num_gA)*(num_gsub),3);

%For rheobase (not used here)
current_slice = zeros(num_gA,num_gsub);
rheo0_1 = zeros((num_istim)*(num_gA)*(num_gsub),3);
rheo0_2 = zeros((num_istim)*(num_gA)*(num_gsub),3);
rheo0_3 = zeros((num_istim)*(num_gA)*(num_gsub),3);
rheo0_4 = zeros((num_istim)*(num_gA)*(num_gsub),3);
rheo1_2 = zeros((num_istim)*(num_gA)*(num_gsub),3);
rheo1_3 = zeros((num_istim)*(num_gA)*(num_gsub),3);
rheo1_4 = zeros((num_istim)*(num_gA)*(num_gsub),3);
rheo2_3 = zeros((num_istim)*(num_gA)*(num_gsub),3);
rheo2_4 = zeros((num_istim)*(num_gA)*(num_gsub),3);
rheo3_4 = zeros((num_istim)*(num_gA)*(num_gsub),3);

cur_istim = -1;

cur_ind = 1;
param_array_mod = zeros((num_istim)*(num_x)*(num_y),6 + nSpikeATPs);
for i=1:num_istim*num_gA*num_gsub
    if param_array(i,2) <= grid_size && param_array(i,3) <= grid_size
        param_array_mod(cur_ind,:) = param_array(i,:);
        cur_ind = cur_ind + 1;
    end
end

%For each data-point (i.e. each neuron simulated), add pertinent data to
%comprehensive array that stores color corresponding to each neuron, based
%on firing-type or firing-rate
for i=1:(num_istim)*(num_gA)*(num_gsub)
    %Firing-Type color scheme:
    if param_array(i,4) == 0 %R = Grey
        type_c(i,1) = 0.7;
        type_c(i,2) = 0.7;
        type_c(i,3) = 0.7;
    elseif param_array(i,4) == 1 %SS = Blue
        type_c(i,1) = 0;
        type_c(i,2) = 0;
        type_c(i,3) = 1;
    elseif param_array(i,4) == 2 %DO = Green
        type_c(i,1) = 0;
        type_c(i,2) = 1;
        type_c(i,3) = 0;
    elseif param_array(i,4) == 3 %Gap = Yellow
        type_c(i,1) = 1;
        type_c(i,2) = 1;
        type_c(i,3) = 0;
    elseif param_array(i,4) == 4 %RF = Red
        type_c(i,1) = 1;
        type_c(i,2) = 0;
        type_c(i,3) = 0;
    end
    
    %Firing-rate color scheme:
    if param_array(i,5) > 0
        rate_c(i,1) = (1.0/double(param_array(i,5)));
        rate_c(i,2) = 0.4;
        rate_c(i,3) = 0.4;
    else
        rate_c(i,1) = 0.7;
        rate_c(i,2) = 0.7;
        rate_c(i,3) = 0.7;
    end
end

%Modified type_c size (if grid-size is different from gA/gsub value)
for i=1:(num_istim)*(num_x)*(num_y)
    %Firing-Type color scheme:
    if param_array_mod(i,4) == 0 %R = Grey
        type_c_mod(i,1) = 0.7;
        type_c_mod(i,2) = 0.7;
        type_c_mod(i,3) = 0.7;
    elseif param_array_mod(i,4) == 1 %SS = Blue
        type_c_mod(i,1) = 0;
        type_c_mod(i,2) = 0;
        type_c_mod(i,3) = 1;
    elseif param_array_mod(i,4) == 2 %DO = Green
        type_c_mod(i,1) = 0;
        type_c_mod(i,2) = 1;
        type_c_mod(i,3) = 0;
    elseif param_array_mod(i,4) == 3 %Gap = Yellow
        type_c_mod(i,1) = 1;
        type_c_mod(i,2) = 1;
        type_c_mod(i,3) = 0;
    elseif param_array_mod(i,4) == 4 %RF = Red
        type_c_mod(i,1) = 1;
        type_c_mod(i,2) = 0;
        type_c_mod(i,3) = 0;
    end
end

%%Keeping one variable constant, geenrate slice through data such that all
%%points in slice have that constant value of the parameter meant to be
%%constant (given by 'ind')
slize_size = num_istim*num_gA*num_gsub;
if ind == 1 %if i_stim is constant (this is used primarily)
    slice_size = num_x*num_y; %num_gsub*num_gA;
    x_ind = 2;
    y_ind = 3;
    x_size = num_gA;
    y_size = num_gsub;
    param_array_slice_contour_FR = zeros(num_y,num_x);
    param_array_slice_contour_FP = zeros(num_y,num_x);
    param_array_slice_contour_ATPave = zeros(num_y,num_x);
    param_array_slice_contour_ATP = zeros(num_y,num_x);
    param_array_slice_contour_ATP2 = zeros(num_y,num_x);
elseif ind == 2 %if gsub constant (not used)
    slice_size = num_istim*num_gA;
    x_ind = 3;
    y_ind = 1;
    x_size = num_gA;
    y_size = num_istim;
    param_array_slice_contour_FR = zeros(num_gA,num_gsub);
    param_array_slice_contour_FP = zeros(num_gA,num_gsub);
    param_array_slice_contour_ATPave = zeros(num_gA,num_gsub);
    param_array_slice_contour_ATP = zeros(num_gA,num_gsub);
    param_array_slice_contour_ATP2 = zeros(num_gA,num_gsub);
elseif ind == 3 %if gA constant (not used)
    slice_size = num_istim*num_gsub;
    x_ind = 2;
    y_ind = 1;
    x_size = num_gsub;
    y_size = num_istim;
    param_array_slice_contour_FR = zeros(num_gA,num_gsub);
    param_array_slice_contour_FP = zeros(num_gA,num_gsub);
    param_array_slice_contour_ATPave = zeros(num_gA,num_gsub);
    param_array_slice_contour_ATP = zeros(num_gA,num_gsub);
    param_array_slice_contour_ATP2 = zeros(num_gA,num_gsub);
end

%note: slice_size is modified above if modified domain size used
param_array_slice = zeros(slice_size, 6 + nSpikeATPs);
temp_slice = zeros(num_y, num_x); 
prev_slice = zeros(num_y, num_x);
rheo_slice = [];
rheo_colour = [];
color_array_slice = zeros(slice_size, 3);

cur_ind = 1;
cur_ind2 = 1;
cur_istim = min_istim;
prev_fp = 0;
cur_fp = 0;


for i=1:(num_istim)*(num_gA)*(num_gsub)

    %Calculation of rheobase (not used here)
%     if param_array(i,1) == cur_istim
%         temp_slice( (param_array(i,3) - min_gA)/d_gA + 1, (param_array(i,2) - min_gsub)/d_gsub + 1) = param_array(i,4);
%     else
%         for j=0:num_gA-1
%             for k=0:num_gsub-1
%                 if temp_slice(j+1,k+1) > prev_slice(j+1,k+1)
%                     prev_fp = prev_slice(j+1,k+1);
%                     cur_fp = temp_slice(j+1,k+1);
%                     rheo_slice(cur_ind2,1) = min_gA + d_gA*j;
%                     rheo_slice(cur_ind2,2) = min_gsub + d_gsub*k;
%                     rheo_slice(cur_ind2,3) = cur_istim; %current i_stim
%                     rheo_slice(cur_ind2,4) = temp_slice(j+1,k+1);
%                     
%                     if temp_slice(j+1,k+1) == 0 %R = Grey
%                         rheo_colour(cur_ind2,1) = 0.7;
%                         rheo_colour(cur_ind2,2) = 0.7;
%                         rheo_colour(cur_ind2,3) = 0.7;
%                     elseif temp_slice(j+1,k+1) == 1 %SS = Blue
%                         rheo_colour(cur_ind2,1) = 0;
%                         rheo_colour(cur_ind2,2) = 0;
%                         rheo_colour(cur_ind2,3) = 1;
%                     elseif temp_slice(j+1,k+1) == 2 %DO = Green
%                         rheo_colour(cur_ind2,1) = 0;
%                         rheo_colour(cur_ind2,2) = 1;
%                         rheo_colour(cur_ind2,3) = 0;
%                     elseif temp_slice(j+1,k+1) == 3 %Gap = Yellow
%                         rheo_colour(cur_ind2,1) = 1;
%                         rheo_colour(cur_ind2,2) = 1;
%                         rheo_colour(cur_ind2,3) = 0;
%                     elseif temp_slice(j+1,k+1) == 4 %RF = Red
%                         rheo_colour(cur_ind2,1) = 1;
%                         rheo_colour(cur_ind2,2) = 0;
%                         rheo_colour(cur_ind2,3) = 0;
%                     end
%                     cur_ind2 = cur_ind2 + 1;
%                 end
%             end
%         end
%         prev_slice = temp_slice;
%         temp_slice = zeros(num_gA, num_gsub);
%         cur_istim = cur_istim + d_istim;
%     end
    
    %***Create 'slice' (i.e. all data such that parameter-of-interest
    %equals***
    %'value' (inputted by user) -> e.g. all points where i_stim = value)
    if param_array(i,ind) == value && param_array(i,2) <= grid_size && param_array(i,3) <= grid_size
        for j=1:(6+nSpikeATPs)
            param_array_slice(cur_ind,j) = param_array(i,j);
        end
        color_array_slice(cur_ind,1) = type_c(i,1);
        color_array_slice(cur_ind,2) = type_c(i,2);
        color_array_slice(cur_ind,3) = type_c(i,3);
        cur_ind = cur_ind + 1;
    end
end

gridX = (1:x_size);
gridY = (1:y_size);

rheo_01X = [];
rheo_01Y = [];
rheo_01 = [];

cur_ind = 1;
for i=1:size(rheo_slice)
    if rheo_slice(i,4) == 1
        rheo_01X(cur_ind) = rheo_slice(i,1);
        rheo_01Y(cur_ind) = rheo_slice(i,2);
        rheo_01(cur_ind,cur_ind) = rheo_slice(i,3);
        cur_ind = cur_ind + 1;
    end
end

%%Display figures:

%3D-scatter plot of all points/neurons (i.e. i_stim vs. gA vs. gsub)
if num_istim > 1 
    figure('name','3D Scatter Plot - Firing Pattern');
    scat3 = scatter3(param_array_mod(:,2),param_array_mod(:,3),param_array_mod(:,1),10,type_c_mod,'filled');
    ylabel('gA');
    xlabel('gsub');
    zlabel('istim');
end

%2D-scatter plot (slice) of all points/neurons with i_stim = value (set by user) -> (i.e. i_stim vs. gA vs. gsub)
figure('name','Scatter Plot Slice - Firing Pattern');
scat = scatter(param_array_slice(:,x_ind),param_array_slice(:,y_ind),10,color_array_slice,'filled');
title(['istim = ',int2str(value),', ioff = ',int2str(i_off)]);
xlabel('gsub');
ylabel('gA');

%Same as 2D scatter plot above, but using contourf (commented here since scatter plot conveys same information):
% figure('name','Contour Plot Slice - Firing Pattern');
% [aC aC] = contourf(param_array_slice_contour_FP,5);
% title(['istim = ',int2str(value),', ioff = ',int2str(i_off)]);
% xlabel('gsub');
% ylabel('gA');
% set(aC,'LineStyle','none');
% 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.

%Plot firing-rate slice using contourf
% figure('name','Contour Plot Slice - Firing Rate');
% [hC hC] = contourf(param_array_slice_contour_FR,255);
% set(hC,'LineStyle','none');
% title(['istim = ',int2str(value),', ioff = ',int2str(i_off)]);
% xlabel('gsub');
% ylabel('gA');
% 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.
% colorbar;


%plot rheo-base (not used here)
if size(rheo_slice) > 0
    figure('name','3D Scatter Plot - Rheobase');
    scat_rheo3D = scatter3(rheo_slice(:,1),rheo_slice(:,2),rheo_slice(:,3),20,rheo_colour,'filled');
    xlabel('gA');
    ylabel('gsub');
    zlabel('istim');
end






Loading data, please wait...