Basis for temporal filters in the cerebellar granular layer (Roessert et al. 2015)

 Download zip file 
Help downloading and running models
Accession:168950
This contains the models, functions and resulting data as used in: Roessert C, Dean P, Porrill J. At the Edge of Chaos: How Cerebellar Granular Layer Network Dynamics Can Provide the Basis for Temporal Filters. It is based on code used for Yamazaki T, Tanaka S (2005) Neural modeling of an internal clock. Neural Comput 17:1032-58
Reference:
1 . Rössert C, Dean P, Porrill J (2015) At the Edge of Chaos: How Cerebellar Granular Layer Network Dynamics Can Provide the Basis for Temporal Filters. PLoS Comput Biol 11:e1004515 [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: Cerebellum;
Cell Type(s):
Channel(s):
Gap Junctions:
Receptor(s):
Gene(s):
Transmitter(s):
Simulation Environment: C or C++ program;
Model Concept(s): Rate-coding model neurons; Reservoir Computing;
Implementer(s): Roessert, Christian [christian.a at roessert.de];
# -*- coding: utf-8 -*-
"""
Created on Tue Nov  8 19:12:05 2011

@author: chris

"""

from __future__ import division

#import sys
#sys.path.insert(0, "/home/chris/lib/python")  # for new matplotlib!!!

from pylab import *
from numpy import round, random, any
from units import *
import time
from NeuroTools import stgen
import h5py
import os
import shutil

         
def make_colormap(seq):
    import matplotlib.colors as mcolors
    """Return a LinearSegmentedColormap
    seq: a sequence of floats and RGB-tuples. The floats should be increasing
    and in the interval (0,1).
    """
    seq = [(None,) * 3, 0.0] + list(seq) + [1.0, (None,) * 3]
    cdict = {'red': [], 'green': [], 'blue': []}
    for i, item in enumerate(seq):
        if isinstance(item, float):
            r1, g1, b1 = seq[i - 1]
            r2, g2, b2 = seq[i + 1]
            cdict['red'].append([item, r1, r2])
            cdict['green'].append([item, g1, g2])
            cdict['blue'].append([item, b1, b2])
    return mcolors.LinearSegmentedColormap('CustomMap', cdict)
    
def create_colnoise(t, sexp, cutf = None, seed = None, onf = None):
    '''x = create_colnoise(t, sexp, cutf)
    Make coloured noise signal
    t    = vector of times
    sexp = spectral exponent - Power ~ 1 / f^sexp
    cutf = frequency cutoff  - Power flat (white) for f <~ cutf
    
    cutf == None OR sexp == 0: white noise
    
    output:  mean = 0, std of signal = 1/2 => max 95% is 1
    '''
    
    from scipy import signal
    
    nt = len(t)
    dt = (t[-1] - t[0]) / (nt - 1)
    random.seed(seed)
    
    #if sexp == 0:  # no smooth cutoff
    #    N = 10001   # number of filter taps
    #    xp = int((N-1)/2)
    #    x = random.standard_normal(size = len(t)+xp)   
    #else:
    #    x = random.standard_normal(size = shape(t)) 
        
    x = random.standard_normal(size = shape(t))/2
    #print std(x)
    if cutf == None:
        pass
    else:
        
        if sexp == 0:  # no cutoff
            pass
        

        elif sexp == -1:  # sharp cutoff
            x, freq, freq_wp, freq_used = create_multisines(t, freq_used=None, cutf=cutf, onf=onf)
            #print std(x)
            x = x / std(x) / 2         # std of signal = 1, mean = 0
                        
        else:  # smooth cutoff   
            
            x = fft(x)
            f = fftfreq(nt, dt)
            x[nonzero(f == 0)] = 0               # remove zero frequency contribution
        
            i = nonzero(f != 0)              # find indices i of non-zero frequencies
            x[i]=x[i] / (cutf ** 2 + f[i] ** 2) ** (sexp / 4)  # using i allows cutf = 0
        
            x = real(ifft(x))      # ignore imaginary part (numerical error)
            x = x / std(x) / 2         # std of signal = 1/2, mean = 0
            
    return x


def create_multisines(t, freq_used=array([1]), cutf = None, onf = None):
    """
    This function will produce a colored noise signal using the time points 
    as defined in array t. The signal will be constructed using sinosoids with 
    frequencies as defined in the array freq_used and randomized phases.

    The output consists of signal x        
    """

    tstop = t[-1]
    df = 1 / tstop # frequency stepsize
    dt = t[2] - t[1]
    data_points = len(t) # number of time or frequency steps  

    vector_freq = zeros(data_points) 
    vector_phase = zeros(data_points)

    if cutf != None:
        f = arange(0,data_points)*df
        if onf != None:
            freq_used = f[nonzero((f <= cutf) & (f >= onf))]
            #print freq_used
        else:            
            freq_used = f[nonzero((f <= cutf) & (f > 0))]
        
    index_f_used = round(freq_used / df).astype('int') # get indices of used frequencies in frequency vector
    
    index_fneg_used = (data_points - index_f_used).astype('int') # indices of negative frequencies 

    index_fall_used = concatenate((index_f_used, index_fneg_used)) # indices of pos+neg frequencies

    vector_freq[index_fall_used] = data_points / 2 # each frequency used ???
    
    phase = 2*pi*(random.rand(len(freq_used),1)-0.5) # pick phases randomly shifted by +-pi (sould there by another 2* to shift +-2pi???)

    vector_phase[index_f_used] = phase # assign positive phases to full vector
    vector_phase[index_fneg_used] = -phase # assign negative phases to full vector

    freqtemp = vector_freq * exp(1j * vector_phase) # generate frequency domain response
    x = real(ifft(freqtemp)) # convert into time domain

    #print "- Number of msine frequencies: " + str(2 * std(x) ** 2)
    
    noise_data = x/max(abs(x)) # scale so that signal amplitude is 1
    
    freq = fftfreq(data_points, dt)[ 0:round(data_points / 2) ] # positive frequency vector
   
    noise_power = abs(fft(noise_data))[ 0:round(data_points / 2) ] # compute noise power
    freq_wp = find(noise_power > 2 * std(noise_power)) # threshold to discriminate the indexes of peak frequencies
    freq_used = freq[freq_wp] # vector of used frequencies [Hz]
    
    return noise_data, freq, freq_wp, freq_used 
    

def create_singlesine(fu = 5, amp = 0.5, ihold = 1, dt = 0.025*ms, periods = 10, minlength = 1*s, t_prestim = 2*s, l_zeros = 2):
    """
    This function will produce a single sine signal of frequency fu with holding current ihold
    Signal has at least the length periods*T (s) or minlength (s).
    Use stimulate with pre stimulus of length t_prestim (s)        
    """
    
    fs = 1 / dt  # sampling rate 
    
    tnext = 0
    # delay for no noise input
    start_zeros =  zeros(l_zeros * fs)   
    t_zeros = tnext + arange(0, l_zeros, dt)
    
    tnext = t_zeros[-1] + dt
    l_pre_signal = ceil(t_prestim / (1. / fu)) * 1. / fu # length of pre stimulus should to be at least t_prestim seconds but with length of full periods
    t_pre_signal = arange(0, l_pre_signal, dt) # create pre time vector
    
    pre_signal = amp * sin(2 * pi * t_pre_signal * fu) # create pre signal vector
    t_pre_signal = t_pre_signal + tnext
    
    tnext = t_pre_signal[-1] + dt
    l_t = max(minlength, periods * 1 / fu) # length of input_signal: stimulate for at least periods*T or minlength
    t_input_signal = arange(0, l_t, dt) # create stimulus time vector
    
    #window = sin(2 * pi * t_input_signal * 1/l_t/2)  # not really good if nonlinear membrane function!!!!
    input_signal = amp * sin(2 * pi * t_input_signal * fu)
    t_input_signal = t_input_signal + tnext
    
    i_start = len(start_zeros) + len(pre_signal)
    i_stop = len(start_zeros) + len(pre_signal) + len(input_signal) 
    
    tnext = t_input_signal[-1] + dt
    l_post_signal = 1 # length of post stimulus should only be 1 s, equivalent to 1 Hz lower bound for spike cutoff 
    t_post_signal = arange(0, l_post_signal, dt) # create pre time vector
    post_signal = amp * sin(2 * pi * t_post_signal * fu) # create pre signal vector
    t_post_signal = t_post_signal + tnext

    all_data = concatenate((start_zeros, pre_signal, input_signal, post_signal)) # combine all
    t = concatenate((t_zeros, t_pre_signal, t_input_signal, t_post_signal)) # combine all
    t1 = arange(0, size(all_data) * dt, dt) # time vector of stimulus [s]  
    
    i_startstop = array([i_start, i_stop])
    t_startstop = array([t[i_start], t[i_stop]])
    
    iholdvec = concatenate((zeros(1 * fs), ones(len(all_data) - 1 * fs) * ihold)) # create holding current vector
    #iholdvec = concatenate((zeros(1 * fs), ones(len(all_data) - 2 * fs) * ihold, zeros(1 * fs))) # create holding current vector
    
    stimulus = all_data + iholdvec # create full stimulus vector
    
    return t, stimulus, i_startstop, t_startstop

        
    

def syn_kernel(t, tau1 = 5, tau2 = 10):    

    if tau1 == 0:
        
        G = exp(-t/tau2)    

    else:
        
        if (tau1/tau2 > .9999):
            tau1 = .9999*tau2
    
        tp = (tau1*tau2)/(tau2 - tau1) * log(tau2/tau1)
        factor = -exp(-tp/tau1) + exp(-tp/tau2)
        factor = 1/factor                        
        G = factor * (exp(-t/tau2) - exp(-t/tau1))
    
    return G
    

def construct_Stimulus(noise_data, fs, amp = 1, ihold = 0, tail_points = 2, delay_baseline = 4):
    """
    Construct Stimulus from cnoise/msine input and other parameters.
    """

    #inin = 8 # stimulate before with 10 s of signal
    inin = np.array((len(noise_data)/fs)*0.1).clip(max=8)

    stim_data = concatenate((noise_data[-inin*fs:], noise_data))  # increase length of stimulus # no normalization here: / max(abs(noise_data))
    
    stimulus = concatenate((concatenate((zeros(round(delay_baseline*fs)), amp * stim_data)), zeros(round(tail_points*fs))))  # construct stimulus
            
    iholdvec = concatenate((zeros(round(fs)), ones(round(len(stimulus) - 1 * fs)) * ihold))
    stimulus = stimulus + iholdvec
    
    dt = 1 / fs
    t = arange(0, len(stimulus) * dt,dt)  # time vector of stimulus [s]
            
    t_startstop = np.array([inin+delay_baseline, inin+delay_baseline+len(noise_data)/fs])
    
    return stimulus, t, t_startstop
    

def construct_Pulsestim(dt = 0.025e-3, pulses = 1, latency = 10e-3, stim_start = 0.02, stim_end = 0.02, len_pulse = 0.5e-3, amp_init = 1, amp_next = None):
    """
    Construct a pulse stimulus in the form of |---stim_start---|-len_pulse-|--(latency-len_pulse)--|...|-len_pulse-|--(latency-len_pulse)--|--stim_end--|
    For stim_end shorter than pulse: stim_end = stim_end + len_pulse
    """
    #print dt
    
    fs = 1 / dt
    
    if len(np.shape(amp_next)) > 0:
        pulses = len(amp_next)
        amp_vec = amp_next
    else:
        if amp_next == None:
            amp_vec = np.ones(pulses)*amp_init
        else:
            amp_vec = np.ones(pulses)*amp_next
            amp_vec[0] = amp_init

    if len(np.shape(latency)) > 0: 
        pass
    else:
        latency = np.ones(pulses)*latency
        
    if len(amp_vec) != len(latency):    
        raise ValueError('amp_vec and latency vectors do not have the same size!!!')
        
     
    if stim_end < len_pulse:
        print "From construct_Pulsestim: stim_end shorter than pulse, setting stim_end = stim_end + len_pulse"
        stim_end = stim_end + len_pulse
           
    ivec = zeros(round((stim_start + sum(latency) + stim_end)*fs))  # construct zero vector to begin with 
    t = arange(0, len(ivec))*dt  # construct time vector
    
    ivec[round(stim_start*fs):round((stim_start+len_pulse)*fs)] = amp_vec[0]
                    
    for i in range(1, pulses):
        
        ivec[round((stim_start+sum(latency[0:i]))*fs):round((stim_start+sum(latency[0:i])+len_pulse)*fs)] = amp_vec[i]
                               
    return t, ivec

    
def slugify(value):
    """
    Normalizes string, converts to lowercase, removes non-alpha characters,
    and converts spaces to hyphens.
    
    From Django's "django/template/defaultfilters.py".
    """
    
    for c in r'[]/\;,><&*:%=+@!#^\'()$| ?^':
        value = value.replace(c,'')
        
    return value

    
def nans(shape, dtype=float):
    a = np.empty(shape, dtype)
    a.fill(np.nan)
    return a

    
def nanmean(a, dim=0):
    return np.mean(np.ma.masked_array(a,np.isnan(a)),dim).filled(np.nan)  


def rw_hdf5(filepath = "data_dict.hdf5", data_dict = None, export = False):
    
    if data_dict == None: # load
        
        data_dict = {}
        
        print "rw_hdf5 load:", filepath
        f = h5py.File(filepath, 'r')

        for i in f.items():
            data_dict[i[0]]  = np.array(i[1])    
            
        f.close() 
            
        if export:
            if not os.path.exists(export):
                os.makedirs(export)
            shutil.copy(filepath, export)

    else:
        
        f = h5py.File(filepath, 'w')
        
        for name  in data_dict:

            data0 = data_dict[name]
            
            if (type(data0) == np.ndarray):
                pass
            else:
                data0 = np.array(data0)
            
            #print data0
            if len(np.shape(data0)) < 1:
                f.create_dataset(name, data=data0)
            else:
                f.create_dataset(name, data=data0, compression='lzf')
                

        f.close() 
    
    return data_dict
    

# test code
if __name__ == '__main__': 
    
    pass