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;
/
StanleyEtAl2013
matlab_plotting
matlab_libdave1
Wavelet
struct_masc
stats_dave.m
                            
function s = stats_dave (s)

format compact;
global sig;         %This will be used later in the curve fitting
global baseline_compare
clean_memory = 0;
clean_filtered = 0;
plotting = 0;       %Turn plotting on/off
baseline_compare = 0;    % Controls whether the baseline is plotted against 
                                % the current frequency spectrum


if isfield (s,'statsdata') == 1
    if s.statsdata.mean == mean(s.data)
        already_calculated = 1;
    else
        already_calculated = 0;
    end
end

already_calculated = 0;       %Hard coding to force recalculation of values
                              %Essential for debugging!

if ~already_calculated
    bin_duration = 5; %seconds
    s.statsdata.mean = mean(s.data);
    [s.statsdata.std s.statsdata.stderr] = dave_binoverlap_stats(s.datatimes, s.datafilt2,bin_duration,'std(x)');
    [s.statsdata.var s.statsdata.varerr] = dave_binoverlap_stats(s.datatimes, s.datafilt2,bin_duration,'var(x)');
    [s.statsdata.skew s.statsdata.skewerr] = dave_binoverlap_stats(s.datatimes, s.datafilt2,bin_duration,'skewness(x,0)');       %Setting flag to zero makes an unbiased estimator
    [s.statsdata.kurt s.statsdata.kurterr] = dave_binoverlap_stats(s.datatimes, s.datafilt2,bin_duration,'kurtosis(x,0)');


    % Histogram
    IQR = iqr(s.datafilt2);
    len = length(s.datafilt2);
    spacing = 2*IQR*len^(-1/3);   % Estimate the appropriate number of bins using Freedman-Draconis ruls
    nbins = ceil((max(s.datafilt2) - min(s.datafilt2))/spacing);
    if nbins == Inf
        nbins = 3;
    end
    [s.statsdata.nhist s.statsdata.binloc] = hist(s.datafilt2, nbins);

    % Wavelet & FFT
    [wvf wvf_val] = wvSpect(s.datatimes, s.data);
    wvf_val = sqrt(wvf_val);
    s.fft.wvf = fliplr(wvf);
    s.fft.wvfft_val = fliplr(wvf_val);
    [s.fft.f s.fft.fft_val] = dave_binoverlap_FFT(s.datatimes, s.data, 5.0);   % Taken from.... Me!!!
%     [s.fft.f s.fft.fft_val] = daveFFT(s.datatimes, s.data, 1);   % Taken from Steinmitz/Koch
    s.fft.fft_h = s.fft.fft_val(1:round(end/2));
    s.fft.f_h = s.fft.f(1:round(end/2));
    temp = round(length(s.fft.f)/2); s.fft.f = s.fft.f(1:temp); s.fft.fft_val = s.fft.fft_val(1:temp);

    
    %   Power statistics
        lowfreq_min = 0.2;
        lowfreq_max = 5;
            lowrange = find( (lowfreq_min <= s.fft.f) .* (s.fft.f < lowfreq_max) );
        midfreq_min = lowfreq_max;
        midfreq_max = 20;
            midrange = find( (midfreq_min <= s.fft.f) .* (s.fft.f < midfreq_max) );
        highfreq_min = 50;
        highfreq_max = 100;
            highrange = find( (highfreq_min <= s.fft.f) .* (s.fft.f < highfreq_max) );
        s.specstd.low = sum( abs(s.fft.fft_val(lowrange)).^2 ).^(1/2);
        s.specstd.mid = sum( abs(s.fft.fft_val(midrange)).^2 ).^(1/2);
        s.specstd.high = sum( abs(s.fft.fft_val(highrange)).^2 ).^(1/2);
        s.specstd.lowerr = 0;
        s.specstd.miderr = 0;
        s.specstd.higherr = 0;
        
%         original frequencies
%         lowfreq_min = 0.2;
%         lowfreq_max = 10;
%             lowrange = find( (lowfreq_min <= s.fft.f) .* (s.fft.f < lowfreq_max) );
%         midfreq_min = lowfreq_max;
%         midfreq_max = 100;
%             midrange = find( (midfreq_min <= s.fft.f) .* (s.fft.f < midfreq_max) );
%         highfreq_min = 1000;
%         highfreq_max = 5000;
%             highrange = find( (highfreq_min <= s.fft.f) .* (s.fft.f < highfreq_max) );
%         s.specstd.low = sum( abs(s.fft.fft_val(lowrange)).^2 ).^(1/2);
%         s.specstd.mid = sum( abs(s.fft.fft_val(midrange)).^2 ).^(1/2);
%         s.specstd.high = sum( abs(s.fft.fft_val(highrange)).^2 ).^(1/2);
%         s.specstd.lowerr = 0;
%         s.specstd.miderr = 0;
%         s.specstd.higherr = 0;



%         Jacobson Power stats
        lowfreq_min = 0.2;
        lowfreq_max = 2;
            jlowrange = find( (lowfreq_min <= s.fft.f) .* (s.fft.f < lowfreq_max) );
        midfreq_min = 5;
        midfreq_max = 100;
%         midfreq_min = 15;
%         midfreq_max = 35;
            jmidrange = find( (midfreq_min <= s.fft.f) .* (s.fft.f < midfreq_max) );
        s.specstd.jlow = sum( abs(s.fft.fft_val(jlowrange)).^2 ).^(1/2);
        s.specstd.jmid = sum( abs(s.fft.fft_val(jmidrange)).^2 ).^(1/2);
        s.specstd.jlowerr = 0;
        s.specstd.jmiderr = 0;
            

    %   Wavelt Power statistics
        lowfreq_min = 0.2;
        lowfreq_max = 10;
            lowrange = find( (lowfreq_min <= s.fft.wvf) .* (s.fft.wvf < lowfreq_max) );
        midfreq_min = lowfreq_max;
        midfreq_max = 100;
            midrange = find( (midfreq_min <= s.fft.wvf) .* (s.fft.wvf < midfreq_max) );
        highfreq_min = 1000;
        highfreq_max = 5000;
            highrange = find( (highfreq_min <= s.fft.wvf) .* (s.fft.wvf < highfreq_max) );
        s.wvstd.low = sum( abs(s.fft.wvfft_val(lowrange)) );
        s.wvstd.mid = sum( abs(s.fft.wvfft_val(midrange)) );
        s.wvstd.high = sum( abs(s.fft.wvfft_val(highrange)) );
        s.wvstd.lowerr = 0;
        s.wvstd.miderr = 0;
        s.wvstd.higherr = 0;


    %//////////////
    %Normalize the power/amplitude of the white noise to that of the
    %original signal (optional)
    % noise = noise / std(noise) * std(s.statsdata.data); %normalize power
    %noise = (noise / sum(abs(noise))) * sum(abs(s.statsdata.data));    %normalize amplitude
    %///////////////////// (end optional)

    if isfield (s,'noise')
        s.statsnoise.mean = mean(s.noise);
        s.statsnoise.std = std(s.noise);
        s.statsnoise.var = var(s.noise);
        s.statsnoise.skew = skewness (s.noise, 0);        %Setting flag to zero makes an unbiased estimator
        s.statsnoise.kurt = kurtosis (s.noise, 0);
        [s.statsnoise.nhist s.statsnoise.binloc] = hist(s.noise, nbins);
        [s.fftnoise.f s.fftnoise.fft_val] = dave_binoverlap_FFT(s.noisetimes, s.noise, 10);

        %Normalize histogram so that it reports approx the same # of data points
        s.statsnoise.nhist = (s.statsnoise.nhist / sum(s.statsnoise.nhist)) * sum(s.statsdata.nhist);
    end



    % Curve fitting to log-log
    intlist = sortrows (s.interval);
%     if length(s.data) > 200000      %Includes only RO cases
%         max_freq = 5000;
%     else                            %Includes baseline and PI cases
%         max_freq = max(s.fft.f);
%     end

    min_freq = 1;
    max_freq = 100;
    
    %Fit FFT
    intlist = [0 min_freq; s.interval; max_freq Inf];
    %[const_est beta_est fitlist] = fit_betastandard_scaling (s.fft.f, s.fft.fft_val, intlist);
    [const_est beta_est fitlist] = fit_betalog (s.fft.f, s.fft.fft_val, intlist);
    
    
    s.general_beta_est.beta_est = beta_est;
    s.general_beta_est.beta_esterr = 0;
    s.general_beta_est.const_est = const_est;
    s.fft.fitlist = fitlist;
    
    %Fit Wavelet
    intlist = [0 min_freq; max_freq Inf];
    [const_est beta_est fitlist] = fit_betalog (s.fft.wvf, s.fft.wvfft_val, intlist);
    %[const_est beta_est fitlist] = fit_betastandard (s.fft.wvf, s.fft.wvfft_val, intlist);
    
    s.general_beta_est.wvbeta_est = beta_est;
    s.general_beta_est.wvbeta_esterr = 0;
    s.general_beta_est.wvconst_est = const_est;
    s.fft.wvfitlist = fitlist;



    % Fitting PDF
    [coefs_out resnorm_out] = fit_hist (s.datafilt2,  s.statsdata.binloc, s.statsdata.nhist);
%     [mu sig] = normfit (s.datafilt2);
%     [coefs3] = gamfit (abs(s.datafilt2));
%     [coefs4] = wblfit (abs(s.datafilt2));
    [coefs_out5 resnorm_out5] = fit_gendist (s.datafilt2,  s.statsdata.binloc, s.statsdata.nhist);
%     [coefs_out6 resnorm_out6] = fit_cauchy (s.datafilt2,  s.statsdata.binloc, s.statsdata.nhist);    
%     [coefs_out7 resnorm_out7] = fit_cauchy_gengaus (s.datafilt2,  s.statsdata.binloc, s.statsdata.nhist);    
    coefs_out7 = [6.6472e+03 0.1564 1.7534 0.5327];        % fit_cauchy_gengaus crashes sometimes, so i hardcode this instead

    
    s.statsdata.pdfcoefs = coefs_out;
%     s.statsdata.pdfcoefs2 = [mu sig];
%     s.statsdata.pdfcoefs3 = coefs3;
%     s.statsdata.pdfcoefs4 = coefs4;
    s.statsdata.pdfcoefs5 = coefs_out5;
%     s.statsdata.pdfcoefs6 = coefs_out6;
    s.statsdata.pdfcoefs7 = coefs_out7;
    
    s.statsdata.gauspdf = coefs_out(2);
    s.statsdata.gauspdferr = 0;
    
    coefs_out(2);
    resnorm_out;

    if isfield (s,'noise')
            [coefs_in resnorm_in] = fit_hist (s.noise,  s.statsnoise.binloc, s.statsnoise.nhist);            
            s.statsnoise.pdfcoefs = coefs_in;
            coefs_in(2);
            resnorm_in;
    end
    
%     echo on
%     data = s.data(:);
%     Demitre_psd_scaling(data)
%     echo off

    s.statsdata.stdspread = 0;
    s.statsdata.skewspread = 0;
    s.statsdata.kurtspread = 0;
    s.statsdata.pdfcoefsspread = 0;
    s.statsdata.general_beta_est.beta_estspread = 0;
    s.statsdata.general_beta_est.wvbeta_estspread = 0;
    s.specstd.lowspread = 0;
    s.specstd.midspread = 0;
    s.specstd.higherr = 0;
    s.specstd.jlowspread = 0;
    s.specstd.jmidspread = 0;    

end


if plotting == 1
    stats_plot (s);
	%stats_plot_loglog (s);    
end

clear global baseline_compare


if clean_memory == 1
   s.fft = rmfield (s.fft, 'f');
   s.fft = rmfield (s.fft, 'fft_val');
   s.fft = rmfield (s.fft, 'fitlist'); 
   s.fft = rmfield (s.fft, 'f_h');
   s.fft = rmfield (s.fft, 'fft_h');
end


if clean_filtered == 1
   s = rmfield(s, 'datafilt');
   s = rmfield(s, 'datafilt2');
end

end





Loading data, please wait...