Structure-dynamics relationships in bursting neuronal networks revealed (Mäki-Marttunen et al. 2013)

 Download zip file 
Help downloading and running models
This entry includes tools for generating and analyzing network structure, and for running the neuronal network simulations on them.
1 . Mäki-Marttunen T, Acimovic J, Ruohonen K, Linne ML (2013) Structure-dynamics relationships in bursting neuronal networks revealed using a prediction framework. PLoS One 8:e69373 [PubMed]
Citations  Citation Browser
Model Information (Click on a link to find other models with that property)
Model Type: Realistic Network;
Brain Region(s)/Organism:
Cell Type(s): Neocortex L5/6 pyramidal GLU cell; Abstract Wang-Buzsaki neuron; Abstract integrate-and-fire leaky neuron;
Channel(s): I Na,p; I Na,t; I K; I K,leak; I M;
Gap Junctions:
Receptor(s): GabaA; AMPA; NMDA;
Transmitter(s): Gaba; Glutamate;
Simulation Environment: NEST; MATLAB; Python;
Model Concept(s): Bursting; Synchronization; Short-term Synaptic Plasticity; Methods; Connectivity matrix;
Implementer(s): Maki-Marttunen, Tuomo [tuomomm at];
Search NeuronDB for information about:  Neocortex L5/6 pyramidal GLU cell; GabaA; AMPA; NMDA; I Na,p; I Na,t; I K; I K,leak; I M; Gaba; Glutamate;
%function [BL,BS,start_t,BL_abs,mFr,Rs,Fs] = burstanalysis(firings,dt,maxISI,minSpikes,minElectrodes,ignoreBegin,sigma,divisionpars)
%  Analyses the bursting properties using the population spike train "firings"
%  Input:
%    firings - a Nx2 matrix; first column has the spike times (in ms) and the
%              second the corresponding neuron indices.
%    binwidth - the binsize for calculating the burst profile
%    maxISI - the maximum distance between spikes that can belong to the same
%             burst
%    minSpikes - the minimum number of spikes that a burst must have
%    minUnique - the minimum number of individual neurons that must spike
%                in a burst
%    ignoreBegin - ignore the spikes whose spike time is less than this value
%                  (in ms)
%    sigma - the st. deviation of Gaussian with which the profile is smoothed
%  Output:
%    BL - burst lengths
%    BS - burst sizes (how many spikes in a burst)
%    start_t - times of the first spike of the bursts
%    BL_abs - absolute burst lengths (duration from first to last spike)
%    mFr - maximum firing rates of the bursts
%    Rs - lengths of the rising slopes
%    Fs - lengths of the falling slopes
% Tuomo Mäki-Marttunen
% Last modified 8.1.2013

function [BL,BS,start_t,BL_abs,mFr,Rs,Fs] = burstanalysis(firings,dt,maxISI,minSpikes,minUnique,ignoreBegin,sigma)

if nargin < 2 || isempty(dt)
    dt = 1;
if nargin < 3 || isempty(maxISI)
    maxISI = 25;
if nargin < 4 || isempty(minSpikes)
    minSpikes = 40;
if nargin < 5 || isempty(minUnique)
    minUnique = 30; %the minimum number of spikers in a burst
if nargin < 6 || isempty(ignoreBegin)
    ignoreBegin = 0; %ignore x ms from the beginning
if nargin < 7 || isempty(sigma)
    sigma = 2.5; %ms

BS = [];
start_t = [];
BL_abs = [];
mFr = [];
Rs = [];
Fs = [];

firings(firings(:,1) < ignoreBegin,:) = [];
[~,ord] = sort(firings(:,1));
firings = firings(ord,:);
times = firings(:,1); %firings times (in ms)
dtimes = times(2:end)-times(1:end-1); %ISIs (in ms)
chind = [find(dtimes > maxISI); length(times)]; %the indices ahead of which lies and inter-burst period.
                                                %The last index is included so that the last burst is also included
sd = sigma/dt; %how many bins is the SD of the convoluting Gaussian
NemptySteps = round(sd*5); %consider the Gaussian up to +-5*sigma

for j=1:length(chind)-1
    chind_j = 0; %this makes it possible to study the first possible burst too. the first chind is thus 0
    if j > 1
        chind_j = chind(j-1);
    SCnew = chind(j)-chind_j; %number of spikes in the considered burst
    nEl = length(unique(firings(chind_j+1:chind(j+1),2))); %number of unique spikers in the considered burst
    if SCnew < minSpikes || nEl < minUnique

    BL_absnew = times(chind(j)) - times(chind_j+1); % and burst length (in ms)

    mappi = zeros(1,floor(BL_absnew/dt)+1+2*NemptySteps); %the length of the burst profile is the number of
                                                          %bins required for the duration of the burst +-5*sigma
    mappiprofile = mappi; %mappi will have the raw global firing rate, and mappiprofile the smoothened
    %Save the number of spikes in bins of length dt to vector mappi:
    for i=1:length(mappi)
        mappi(i) = sum(1+floor((times(chind_j+1:chind(j)) - times(chind_j+1))/dt) == i-NemptySteps);
    %Smooth the vector mappi using Gaussian of SD sigma and save to mappiprofile
    gw = 1/sqrt(2*pi*sd^2)*exp(-1/2/(sd^2)*(-length(mappi)+1:length(mappi)-1).^2);
    for i=1:length(mappi)
      mappiprofile(i) = mappi*gw((1:length(mappi))+length(mappi)-i)';
    BS = [BS, SCnew]; %detect this as a burst, save spike count
    BL_abs = [BL_abs, BL_absnew]; % and absolute burst length (in ms)

    [a,peakind] = max(mappiprofile); %the height and index (in mappiprofile) of the peak
    mappiprofile = [0 mappiprofile 0]; %put zeros into the ends to make sure that a half-width exists
    % Calculate the widths of the rising slope and falling slope (both
    % using half (default) and quarter (application) mFr
    crossind1 = find(mappiprofile(1:end-1) < a/2 & mappiprofile(2:end) >= a/2);
    crossind2 = find(mappiprofile(1:end-1) >= a/2 & mappiprofile(2:end) < a/2);

    mFr = [mFr a/dt];
    Rs = [Rs dt*(peakind - (min(crossind1) + (abs(mappiprofile(min(crossind1))-a/2) > abs(mappiprofile(min(crossind1)+1)-a/2))))];
    Fs = [Fs dt*((max(crossind2) + (abs(mappiprofile(max(crossind2))-a/2) > abs(mappiprofile(max(crossind2)+1)-a/2)))-peakind)];
    %Fs = [Fs dt*((min(crossind2) + (abs(mappiprofile(min(crossind2))-a/2) > abs(mappiprofile(min(crossind2)+1)-a/2)))-peakind)];
    start_t = [start_t times(chind_j+1)]; %the start time of the burst
    if length(BL_abs)>50
BL = Rs + Fs;