Hippocampal CA3 network and circadian regulation (Stanley et al. 2013)

 Download zip file 
Help downloading and running models
Accession:142104
This model produces the hippocampal CA3 neural network model used in the paper below. It has two modes of operation, a default mode and a circadian mode. In the circadian mode, parameters are swept through a range of values. This model can be quite easily adapted to produce theta and gamma oscillations, as certain parameter sweeps will reveal (see Figures). BASH scripts interact with GENESIS 2.3 to implement parameter sweeps. The model contains four cell types derived from prior papers. CA3 pyramidal are derived from Traub et al (1991); Basket, stratum oriens (O-LM), and Medial Septal GABAergic (MSG) interneurons are taken from Hajos et al (2004).
Reference:
1 . Stanley DA, Talathi SS, Parekh MB, Cordiner DJ, Zhou J, Mareci TH, Ditto WL, Carney PR (2013) Phase shift in the 24-hour rhythm of hippocampal EEG spiking activity in a rat model of temporal lobe epilepsy. J Neurophysiol 110:1070-86 [PubMed]
Model Information (Click on a link to find other models with that property)
Model Type: Realistic Network; Neuron or other electrically excitable cell;
Brain Region(s)/Organism: Hippocampus; Medial Septum;
Cell Type(s): Hippocampus CA3 pyramidal GLU cell; Hippocampus CA3 interneuron basket GABA cell; Hippocampus CA3 stratum oriens lacunosum-moleculare interneuron; Hippocampus septum medial GABAergic neuron;
Channel(s): I Na,t; I A; I K; I h; I K,Ca; I Calcium;
Gap Junctions:
Receptor(s): GabaA; AMPA;
Gene(s):
Transmitter(s): Gaba; Glutamate;
Simulation Environment: GENESIS; MATLAB;
Model Concept(s): Epilepsy; Brain Rhythms; Circadian Rhythms;
Implementer(s): Stanley, David A ;
Search NeuronDB for information about:  Hippocampus CA3 pyramidal GLU cell; Hippocampus CA3 interneuron basket GABA cell; GabaA; AMPA; I Na,t; I A; I K; I h; I K,Ca; I Calcium; Gaba; Glutamate;
function s = zero_cross_ISI (s)

format compact;
dt = s.datatimes(2) - s.datatimes(1);


datafilt_lowpass = qif (s.datatimes, s.datafilt, [0 1]); datatimes_lowpass = s.datatimes;
[datatimes_lowpass datafilt_lowpass] = lowpass_avg (s.datatimes, datafilt_lowpass, 100);

% s.data_diff = diff(datafilt_lowpass);
% s.datatimes_diff = s.datatimes (1:length(s.data)-1);














% tr = 0.0:0.02:0.22; %threshold range
% mir = 0.01:0.01:0.02; %minimium interval range
lpstd = std(datafilt_lowpass);
% tr = 1.1 * lpstd;
tr = 0.65;                          % This value seems to work great for ic7gap and ic7gapsyn. Can decrease lowpass filter freq
mir = 0.000001;


figure; hold on;
ds = 5;
num = 1:length(datafilt_lowpass);
plotds (s.datatimes, s.data - mean(s.data), 'b', ds); hold on;
% plot (dt * wkeep(num, length(s.datafilt_nobase), 'c'), s.datafilt_nobase, 'g:');
plotds ( datatimes_lowpass, datafilt_lowpass - mean(datafilt_lowpass), 'r', ds);
plotds (dt * num, tr + zeros(1, length(num)), 'k', ds);
legend('unfiltered', 'filtered', ['stdev=' num2str(lpstd)]);

% return;



if length(mir) > 1
    dmir = mir(2) - mir(1);
    dtr = tr(2) - tr(1);
    amat = zeros(length(tr), length(mir));
end

for min_int = mir
    a = [];
    for thresh = tr
        thresh;
        [ints EPSPamps indicies] = down_up_ints (datatimes_lowpass, datafilt_lowpass, thresh);
    %     ints = down_up_ints (s.datatimes_diff, s.data_diff, thresh);

    %     ints = gamrnd(2,2,1,50000);
        keeping = find(ints>min_int);
        ints = ints(keeping);
        indicies = indicies(keeping+1);
        EPSPamps = EPSPamps(keeping);
        
        
        hold on; plot(s.datatimes(indicies),tr*ones(length(indicies),1),'r*');
        
        % EPSP Amplitude Histogram
        IQR = iqr(EPSPamps);
        len = length(EPSPamps);

        %  Friedman Diaconis Rule
        spacing = 2*IQR*len^(-1/3);                       % Estimate the appropriate number of bins
        nbins = ceil ((max(EPSPamps) - min(EPSPamps))/spacing); % using Freedman-Draconis ruls

        sp = spacing;

        [epsp_hist epsp_bins] = hist(EPSPamps, min(EPSPamps):sp:(max(EPSPamps)));
        figure; h1 = bar (epsp_bins, epsp_hist,'b');


        % ISI Times Histogram
        IQR = iqr(ints);
        len = length(ints);

    %     % Sturges' Formula
    %     nbins = log2(len) + 1;
    %     spacing = (max(ints) - min(ints)) / nbins;

    %      % Scott Rule
    %      spacing = 3.5*std(ints)*len^(-1/3);   % Estimate the appropriate number of bins
    %      nbins = ceil ((max(ints) - min(ints))/spacing); % using Freedman-Draconis ruls    

        %  Friedman Diaconis Rule
        spacing = 2*IQR*len^(-1/3);                       % Estimate the appropriate number of bins
        nbins = ceil ((max(ints) - min(ints))/spacing); % using Freedman-Draconis ruls

        sp = max(dt, spacing);

        [nhist binloc] = hist(ints, min(ints):sp:max(ints));

    
        fprintf ('Thresh = %g; min_int = %g\n',thresh,min_int);

        [coefs resnorm] = fit_poisson (ints, binloc, nhist);
        [area_under_curve resnorm2] = fit_poisson_scale (ints, binloc, nhist);
%         area_under_curve = sum(nhist)*sp;

        figure; h1 = bar (binloc, nhist,'b'); hold on;
        plot (binloc, area_under_curve*exppdf(binloc, mean(ints)),'g');
        h2 = plot (binloc, poisson_exp ([coefs(1) coefs(6)], binloc), 'r');    
        legend (['thresh=' num2str(thresh) ' mint=' num2str(min_int) 'mean=' num2str(mean(ints))] ,['a=' num2str(coefs(2)) ' b=' num2str(coefs(3)) ' max=' num2str(coefs(4)) ' err=' num2str(resnorm)]);

        % Make a array (old code, 1d)
        a = [a coefs(2)];

        % Make a matrix (2d)
        if length(mir) > 1
            miindex = round((min_int-min(mir))/dmir + 1);
            trindex = round((thresh-min(tr))/dtr + 1);
            amat(trindex, miindex) = coefs(2);
        end
    end


    
    if length(a) > 1
        figure; bar(tr, a);
        xdiff = max(tr) - min(tr);
        ydiff = max(a) - min(a);
        axis ([(min(tr)-0.25*xdiff) (max(tr)+0.25*xdiff) min(0.8, min(a)-0.25*ydiff) max(a)]);
    end
end

if length(mir)>1
    s.amat = amat;
    imagesc (mir, tr, amat);
end

end



function [times lowpass] = lowpass_avg (datatimes, data, filt_freq)

    % Set constants
    filt_time = 1/filt_freq;
    dt = datatimes(2) - datatimes(1);
    len = length (data);
    
    % Design filter
    filt_size = round(filt_time / dt);
    filt = ones(1, filt_size) / filt_size;
    
    % Apply filter
    fout = conv (data, filt);
    lkeep = len - filt_size;
    lowpass = wkeep (fout, lkeep, 'c');
    times = (0:lkeep-1)*dt;

end

function [ints zc_indices] = crossing_intervals (t, x, thresh)

    x = x - thresh;
    dt = t(2) - t(1);
    N = length(x);
    
    x_sign = ( x >= 0 ) - ( x < 0 );
    zc_list = (x_sign(1:N-1) - x_sign(2:N));
    zc_indices = find (zc_list ~= 0);
    ints = diff (zc_indices);
    ints = ints * dt;

end

function [ints ISIamps downup_indicies] = down_up_ints (t, x, thresh)

    x = x - thresh;
    dt = t(2) - t(1);
    N = length(x);
    
    x_sign = ( x >= 0 ) - ( x < 0 );
    zc_list = (x_sign(1:N-1) - x_sign(2:N));
    downup_indicies = find (zc_list == -2);
    ints = diff (downup_indicies);
    ints = ints * dt;
    
    ISIamps = [];
    x = x + thresh;
    for i = 1:length(downup_indicies)-1
        ISIamps = [ISIamps max(x(downup_indicies(i):downup_indicies(i+1)))];
    end

end




function [coefs_out resnorm_out] = fit_poisson (data, binloc, nhist)
global sig
global C_scale
global b_scale
                % sig = sqrt(alpha) * beta
                % C = Cprime * C_scale



    sig = std(data);
            % Note: variance = alpha*beta^2, so we only need to fit 2
            % parameters
    a0 = 1;
    C0 = (binloc(2) - binloc(1)) * sum(nhist);   %Set the constant multiplier to be equal to the integrated area of our data
    b0 = sig;
    Cprime0 = a0;
    bprime = a0;
    
    C_scale = C0 / Cprime0;
    b_scale = b0 / bprime;
    
    coefs_out0 = [Cprime0 bprime];
    options = optimset ('MaxFunEvals', 5000, 'TolFun', 0.000001);
    [coefs_out resnorm_out] = lsqcurvefit(@poisson_exp, coefs_out0, binloc, nhist, -Inf, Inf, options);
    
    a = 1;
    bprime = coefs_out(2);
    b = bprime * b_scale;
    Cprime = coefs_out(1);
    C = Cprime * C_scale;

    coefs_out = [Cprime a b C C_scale bprime b_scale];
    
end



function f = poisson_exp(coefs, x)
global sig
global C_scale
global b_scale
                % sig = sqrt(alpha) * beta
                % C = Cprime * C_scale

    C = coefs(1) * C_scale;
    a = 1;
    b = coefs(2) * b_scale;
%     b = sig/sqrt(a);
    

    f = C * 1/((b^a)*gamma(a)) * x.^(a-1).*exp(-x/b);

end





function [C_scale resnorm_out] = fit_poisson_scale (data, binloc, nhist)
global lambda
global C_scale
global b_scale
                % lambda = sqrt(alpha) * beta
                % C = Cprime * C_scale



    lambda = mean(data);
            % Note: variance = alpha*beta^2, so we only need to fit 2
            % parameters
    a0 = 1;
    C0 = (binloc(2) - binloc(1)) * sum(nhist);   %Set the constant multiplier to be equal to the integrated area of our data
    b0 = lambda;
    Cprime0 = a0;
    bprime = a0;
    
    C_scale = C0 / Cprime0;
    b_scale = b0 / bprime;
    
    coefs_out0 = [Cprime0];
    options = optimset ('MaxFunEvals', 5000, 'TolFun', 0.000001);
    [coefs_out resnorm_out] = lsqcurvefit(@poisson_default_sig, coefs_out0, binloc, nhist, -Inf, Inf, options);
    
    a = 1;
    b = lambda;
    Cprime = coefs_out(1);
    C = Cprime * C_scale;

%     coefs_out = [Cprime a b C C_scale bprime b_scale];
    C_scale = C;
    
end



function f = poisson_default_sig (coefs, x)
global lambda
global C_scale
global b_scale
                % lambda = sqrt(alpha) * beta
                % C = Cprime * C_scale

    C = coefs(1) * C_scale;
    a = 1;
    b = lambda;
%     b = lambda/sqrt(a);
    

    f = C * 1/((b^a)*gamma(a)) * x.^(a-1).*exp(-x/b);

end

Loading data, please wait...