Inhibition of bAPs and Ca2+ spikes in a multi-compartment pyramidal neuron model (Wilmes et al 2016)

 Download zip file 
Help downloading and running models
Accession:187603
"Synaptic plasticity is thought to induce memory traces in the brain that are the foundation of learning. To ensure the stability of these traces in the presence of further learning, however, a regulation of plasticity appears beneficial. Here, we take up the recent suggestion that dendritic inhibition can switch plasticity of excitatory synapses on and off by gating backpropagating action potentials (bAPs) and calcium spikes, i.e., by gating the coincidence signals required for Hebbian forms of plasticity. We analyze temporal and spatial constraints of such a gating and investigate whether it is possible to suppress bAPs without a simultaneous annihilation of the forward-directed information flow via excitatory postsynaptic potentials (EPSPs). In a computational analysis of conductance-based multi-compartmental models, we demonstrate that a robust control of bAPs and calcium spikes is possible in an all-or-none manner, enabling a binary switch of coincidence signals and plasticity. ..."
Reference:
1 . Wilmes KA, Sprekeler H, Schreiber S (2016) Inhibition as a Binary Switch for Excitatory Plasticity in Pyramidal Neurons. PLoS Comput Biol 12:e1004768 [PubMed]
Citations  Citation Browser
Model Information (Click on a link to find other models with that property)
Model Type: Neuron or other electrically excitable cell;
Brain Region(s)/Organism: Neocortex; Hippocampus;
Cell Type(s): Hippocampus CA1 pyramidal GLU cell; Neocortex V1 L6 pyramidal corticothalamic GLU cell;
Channel(s):
Gap Junctions:
Receptor(s):
Gene(s):
Transmitter(s):
Simulation Environment: NEURON; Python;
Model Concept(s): Dendritic Action Potentials; Synaptic Plasticity; Synaptic Integration;
Implementer(s): Wilmes, Katharina A. [katharina.wilmes at googlemail.com];
Search NeuronDB for information about:  Hippocampus CA1 pyramidal GLU cell; Neocortex V1 L6 pyramidal corticothalamic GLU cell;
#!/usr/bin/env python

"""simulation"""

# _title_     : sim.py
# _author_     : Katharina Anna Wilmes
# _mail_     : katharina.anna.wilmes __at__ cms.hu-berlin.de


# --Imports--
import sys
import neuron
import matplotlib.pyplot as plt
import numpy as np
import matplotlib.mlab as mlab
import random

import nrn
import datetime
import random
import os
from neuron import h



# --load Neuron graphical user interface--
if not ( h('load_file("nrngui.hoc")')):
    print "Error, cannot open NEURON gui"


# general functions
def truncated_gauss(N, mu, sigma=0.05, a=0, b=2):
    """Return a N random numbers from a truncated (a,b) Gaussian distribution."""

    pos = np.zeros(N)
    i = 0
    while i<N:
        x = random.gauss(mu,sigma)
        if a <= x <= b:
            pos[i] = round(x,2)
            i += 1
    x = np.linspace(0,1,100)
    return pos

def gauss(N, mu=0, sigma=0.05):
    """Return a N random numbers from a truncated (a,b) Gaussian distribution."""

    pos = np.zeros(N)
    i = 0
    while i<N:
        x = random.gauss(mu,sigma)
        pos[i] = round(x,2)
        i += 1
    x = np.linspace(0,1,100)
    return pos

class Simulation(object):
    """
    Objects of this class control a simulation. Example of use:
    >>> cell = Cell()
    >>> sim = Simulation(cell)
    >>> sim.go()
    """
    def __init__(self, cell, params):
        self.cell = cell
        self.sim_time = params['sim_time'].value
        self.dt = params['dt'].value
        self.celsius = params['celsius'].value
        self.v_init = params['v_init'].value
        self.high_res = params['high_res']
        self.theta = params['theta'].value
        self.go_already = False
        self.SE = False
        self.vector = False
        self.shunt = h.List()
        self.nc = h.List()
        self.ns = h.List()
        self.v_rest = dict()
        self.syn = h.List()
        self.nc_syn = h.List()
        self.ns_syn = h.List()
        self.inhibition = h.List()
        self.nc_inh = h.List()
        self.ns_inh = h.List()
        self.excitation = h.List()
        self.nc_exc = h.List()
        self.ns_exc = h.List()
        self.stimlist = h.List()
        self.vlist = h.List()
        self.vtimelist = h.List()
        self.stimnumber = 0
        self.mu = []
        self.interneuron = h.List()


    def set_colored_noise_stim(self,params):
        noise = neuron.h.IClamp(params.noise.pos)
        noise.delay = params.noise.delay
        noise.dur = 1e9
        VecStim = h.Vector()  #the ramp vector: linearly increasing from 0 to I_max
        VecStim.play(noise._ref_amp, h.dt)  #playing the ramp vector into the iramp.amp variable




    def set_current(self, vector):
        """
        gets an array with the current that should be delivered by the electrode
        that is named stim
        """
        self.vector = True
        global vList
        timepoints = vector.size    # =self.sim_time/self.dt
        v = h.Vector(timepoints)
        v_time = h.Vector(timepoints)
        for i in range(timepoints):
            v.x[i] = vector[i] # set var to value in vector
        for i in range(timepoints):
            v_time.x[i] = i * self.sim_time / float(timepoints) # =i*self.dt
        self.vlist.append(v)
        self.vtimelist.append(v_time)

    def get_Ramp(self,params):
        """
        defines a stim ramp with duration (10) ms from 0 to max (1) nA,
        concatenated to zeros for sim_time-duration (40-10=30) ms
        """
        ramp = np.arange(0, params.max.value,params.max.value / (params.duration.value / self.dt))
        vector = np.concatenate((ramp, np.zeros(((self.sim_time - params.duration.value) / self.dt))))
        return vector


    def get_Alpha(self,params):
        """
        defines Alpha-function stimulus
        """

        x = np.arange(0, params.duration.value, self.dt)
        if self.SE == True:
            alpha = (self.cell.E + params.peak.value * (x / params.tau.value) * np.exp( - (x - params.tau.value) / params.tau.value))
            vector = np.concatenate((alpha,-70 * np.ones(((self.sim_time - params.duration.value) / self.dt))))
        else:
            alpha = params.max.value * (x / params.tau.value) * np.exp( - (x - params.tau.value) / params.tau.value)
            vector = np.concatenate((alpha, np.zeros(((self.sim_time - params.duration.value) / self.dt))))
        plt.figure()
        plt.plot(np.arange(0, self.sim_time, self.dt), vector)
        plt.show()
        return vector


    def set_Poisson(self,params, rate = None):
        print "set Poisson"
        if rate == None:
            rate = params.rate.value
        else:
            print rate
        syn_pos = params.pos.value
        g = params.g.value
        E = params.E.value
        tau1 = params.tau1.value
        tau2 = params.tau2.value
        number = params.number.value

        if params.location == 'soma' or params.type == 'excitation':
            synapse = neuron.h.Exp2Syn(self.cell.soma(0.5))
        elif params.location == 'apical':
            pos = (syn_pos / 2.0) + (syn_pos-0.5)
            synapse = neuron.h.Exp2Syn(self.cell.apical(pos))
        else:
            raise ValueError

        synapse.e = E
        synapse.tau1 = tau1
        synapse.tau2 = tau2


        ns = neuron.h.NetStim(0.5)
        ns.noise = 1
        ns.start = params.delay.value
        ns.number = number
        ns.interval = 1000.0/rate
        w = float(g)
        nc = h.NetCon(ns,synapse,1,1,w)

        self.syn.append(synapse)
        self.ns_syn.append(ns)
        self.nc_syn.append(nc)



    def set_distributed_shunt(self, params, shunt_pos, shunt_weight, compartment, number, pos_sigma, pattern, timing_sigma, delay, no_reps = 1, interval = None):

        shunt = h.List()
        nc = h.List()
        ns = h.List()
        pos = truncated_gauss(number,shunt_pos,pos_sigma,a=0,b=1) # params.number = number of shunts
        timing = gauss(number, 0, timing_sigma)

        pos.sort()
        for x in pos:
            if x <= 0.2:
                x = x * 5.0
                shunt.append(neuron.h.Exp2Syn(self.cell.apical_prox(x)))
            else:
                x = (x-0.2) * (1.0+(1.0/4.0))
                shunt.append(neuron.h.Exp2Syn(self.cell.apical(x)))
                sec = self.cell.apical.name()

        shunt_count = 0
        for s in shunt:
            s.tau1 = params.tau1.value # ms rise time
            s.tau2 = params.tau2.value # ms decay time
            s.e = params.e.value # mV reversal potential
            ns_i = neuron.h.NetStim(0.5)
            ns_i.number = 1 # (average) number of spikes
            if pattern == 'none':
                # for temporal patterns add delay to shunts at later pos
                ns_i.start = delay # ms (most likely) start time of first spike
            elif pattern == 'gauss':
                ns_i.start = delay + timing[shunt_count]
            elif pattern == 'random':
                ns_i.start = delay - 1 + (2.0 / number) * randi[shunt_count]
            elif pattern == 'sequence':
                ns_i.start = delay - 1 + (2.0 / number) * shunt_count
            else:
                ns_i.start = delay - 1 + (2.0 / number) * (number - shunt_count)
            w = (1.0 / number) * float(shunt_weight)
            ns.append(ns_i)
            nc.append(h.NetCon(ns_i,s,1,0,w))
            shunt_count += 1

        self.shunt = shunt
        self.ns = ns
        self.nc = nc

        return pos, timing


    def set_Shunt(self, params, shunt_pos, shunt_weight, compartment, source = False, delay=None, no_reps = 1, interval = None):

        """
        sets a single or a number of shunts at the defined shunt_pos (e.g. 0.3)
        of the shunt_compartment (e.g. "a", "o", or "b") with a particular shunt_weight (in uS)
        and a delay as specified in params.delay.value unless the function is called
        with delay.
        """
        # apical
        if compartment == 'a':
            if shunt_pos <= 0.2:
                pos = shunt_pos * 5.0
                print pos
                shunt = neuron.h.Exp2Syn(self.cell.apical_prox(pos))
                sec = self.cell.apical_prox.name()
                print sec
            else:
                pos = (shunt_pos-0.2) * (1.0+(1.0/4.0))
                shunt = neuron.h.Exp2Syn(self.cell.apical(pos))
                sec = self.cell.apical.name()
        # oblique
        elif compartment == 'o':
            shunt = neuron.h.Exp2Syn(self.cell.oblique_branch(shunt_pos-1))
            sec = self.cell.oblique_branch.name()
        # apical branch
        elif compartment == 'b':
            shunt = neuron.h.Exp2Syn(self.cell.branches[0](shunt_pos))
            sec = self.cell.branches[0].name()
        # basal
        elif compartment == 'basal':
            shunt = neuron.h.Exp2Syn(self.cell.basal_main(shunt_pos))
            sec = self.cell.basal_main.name()
        # basal branch
        elif compartment == 'basalb':
            shunt = neuron.h.Exp2Syn(self.cell.basal_branches[0](shunt_pos-1))
            sec = self.cell.basal_branches.name()
        else:
            shunt = neuron.h.Exp2Syn(self.cell.branches[1](shunt_pos-1))
        print params
        shunt.tau1 = params['tau1'].value
        shunt.tau2 = params['tau2'].value
        shunt.e = params['reversal'].value
        self.shunt_sec = sec
        ns = neuron.h.NetStim(0.5)
        ns.number = no_reps
        if interval is not None:
            ns.interval = interval
        w = float(shunt_weight)
        if delay == None: # if no value was provided
            ns.start = params.delay.value # ms (most likely) start time of first spike
        else:
            ns.start = delay
        # if inhibition is timed to activity of post-synaptic neuron
        if source:
            nc = h.NetCon(self.cell.soma(0.5)._ref_v, shunt, 1,1,w, sec = self.cell.soma)
        else:
            nc = h.NetCon(ns,shunt,1,0,w)

        self.shunt.append(shunt)
        # list makes it possible to have multiple shunts without overriding
        self.ns.append(ns)
        self.nc.append(nc)




    def set_synaptic_recording(self, switch = True, all=True):
        self.Inh_distalRec = neuron.h.Vector()
        self.nc[0].record(self.Inh_distalRec) # shunting inhibition from set_shunt
        if len(self.nc)>1:
            self.Inh_distalRec2 = neuron.h.Vector()
            self.nc[1].record(self.Inh_distalRec2) # shunting inhibition from set_shunt

        i = 1
        if switch:
            self.Inh2_distalRec = neuron.h.Vector()
            self.nc_syn[i].record(self.Inh2_distalRec)
            i +=1
        if all:
            self.ExcRec = neuron.h.Vector()
            self.InhRec = neuron.h.Vector()
            self.nc_syn[i].record(self.ExcRec)
            self.nc_syn[i+1].record(self.InhRec)

    def set_synaptic_current_rec(self):
        self.syn_current = neuron.h.Vector()
        self.syn_current.record(self.shunt[0]._ref_i)

    def set_synaptic_cond_rec(self):
        self.syn_cond = neuron.h.Vector()
        self.syn_cond.record(self.shunt[0]._ref_g)

    def set_rec_pos(self,rec_pos):
        self.rec_pos = rec_pos


    def set_highestres_recording(self):

        '''Record Time'''
        self.rec_t = neuron.h.Vector()
        self.rec_t.record(neuron.h._ref_t)

        '''record voltage for all compartments'''
        self.rec_va = neuron.h.Vector()
        self.rec_v = neuron.h.Vector()
        self.rec_v1 = neuron.h.Vector()
        self.rec_v2 = neuron.h.Vector()
        self.rec_v3 = neuron.h.Vector()
        self.rec_v4 = neuron.h.Vector()
        self.rec_v5 = neuron.h.Vector()
        self.rec_v6 = neuron.h.Vector()
        self.rec_v7 = neuron.h.Vector()
        self.rec_v8 = neuron.h.Vector()
        self.rec_v9 = neuron.h.Vector()
        self.rec_v10 = neuron.h.Vector()
        self.rec_v11 = neuron.h.Vector()
        self.rec_v12 = neuron.h.Vector()
        self.rec_v13 = neuron.h.Vector()
        self.rec_v14 = neuron.h.Vector()
        self.rec_v15 = neuron.h.Vector()
        self.rec_v16 = neuron.h.Vector()
        self.rec_v17 = neuron.h.Vector()
        self.rec_v18 = neuron.h.Vector()
        self.rec_v19 = neuron.h.Vector()
        self.rec_v20 = neuron.h.Vector()
        self.rec_v21 = neuron.h.Vector()
        self.rec_v22 = neuron.h.Vector()
        self.rec_v23 = neuron.h.Vector()
        self.rec_v24 = neuron.h.Vector()
        self.rec_v25 = neuron.h.Vector()
        self.rec_v26 = neuron.h.Vector()
        self.rec_v27 = neuron.h.Vector()
        self.rec_v28 = neuron.h.Vector()
        self.rec_v29 = neuron.h.Vector()
        self.rec_v30 = neuron.h.Vector()
        self.rec_v31 = neuron.h.Vector()
        self.rec_v32 = neuron.h.Vector()
        self.rec_v33 = neuron.h.Vector()
        self.rec_v34 = neuron.h.Vector()
        self.rec_v35 = neuron.h.Vector()
        self.rec_v36 = neuron.h.Vector()
        self.rec_v37 = neuron.h.Vector()
        self.rec_v.record(self.cell.soma(0.5)._ref_v)
        self.rec_va.record(self.cell.ais(0.5)._ref_v)
        '''record apical and oblique if existent'''
        self.rec_v1.record(self.cell.apical_prox(0.3)._ref_v)
        self.rec_v2.record(self.cell.apical_prox(0.6)._ref_v)
        self.rec_v3.record(self.cell.apical_prox(0.9)._ref_v)
        self.rec_v4.record(self.cell.apical(0.15)._ref_v)
        self.rec_v5.record(self.cell.apical(0.3)._ref_v)
        self.rec_v6.record(self.cell.apical(0.45)._ref_v)
        self.rec_v7.record(self.cell.apical(0.6)._ref_v)
        self.rec_v8.record(self.cell.apical(0.75)._ref_v)
        self.rec_v9.record(self.cell.apical(0.9)._ref_v)


        self.rec_vo1 = neuron.h.Vector()
        self.rec_vo2 = neuron.h.Vector()
        self.rec_vo3 = neuron.h.Vector()
        self.rec_vo4 = neuron.h.Vector()
        self.rec_vo5 = neuron.h.Vector()
        self.rec_vo6 = neuron.h.Vector()
        self.rec_vo7 = neuron.h.Vector()
        self.rec_vo8 = neuron.h.Vector()
        self.rec_vo9 = neuron.h.Vector()
        self.rec_vo1.record(self.cell.oblique_branch(0.1)._ref_v)
        self.rec_vo2.record(self.cell.oblique_branch(0.2)._ref_v)
        self.rec_vo3.record(self.cell.oblique_branch(0.3)._ref_v)
        self.rec_vo4.record(self.cell.oblique_branch(0.4)._ref_v)
        self.rec_vo5.record(self.cell.oblique_branch(0.5)._ref_v)
        self.rec_vo6.record(self.cell.oblique_branch(0.6)._ref_v)
        self.rec_vo7.record(self.cell.oblique_branch(0.7)._ref_v)
        self.rec_vo8.record(self.cell.oblique_branch(0.8)._ref_v)
        self.rec_vo9.record(self.cell.oblique_branch(0.9)._ref_v)

        self.rec_vbasal1 = neuron.h.Vector()
        self.rec_vbasal2 = neuron.h.Vector()
        self.rec_vbasal3 = neuron.h.Vector()
        self.rec_vbasal4 = neuron.h.Vector()
        self.rec_vbasal5 = neuron.h.Vector()
        self.rec_vbasal6 = neuron.h.Vector()
        self.rec_vbasal7 = neuron.h.Vector()
        self.rec_vbasal8 = neuron.h.Vector()
        self.rec_vbasal9 = neuron.h.Vector()
        self.rec_vbasal10 = neuron.h.Vector()
        self.rec_vbasal11 = neuron.h.Vector()
        self.rec_vbasal12 = neuron.h.Vector()
        self.rec_vbasal13 = neuron.h.Vector()
        self.rec_vbasal14 = neuron.h.Vector()
        self.rec_vbasal1.record(self.cell.basal_main(0.1)._ref_v)
        self.rec_vbasal2.record(self.cell.basal_main(0.2)._ref_v)
        self.rec_vbasal3.record(self.cell.basal_main(0.3)._ref_v)
        self.rec_vbasal4.record(self.cell.basal_main(0.4)._ref_v)
        self.rec_vbasal5.record(self.cell.basal_main(0.5)._ref_v)
        self.rec_vbasal6.record(self.cell.basal_main(0.6)._ref_v)
        self.rec_vbasal7.record(self.cell.basal_main(0.7)._ref_v)
        self.rec_vbasal8.record(self.cell.basal_main(0.8)._ref_v)
        self.rec_vbasal9.record(self.cell.basal_main(0.9)._ref_v)
        self.rec_vbasal10.record(self.cell.basal_branches[0](0.1)._ref_v)
        self.rec_vbasal11.record(self.cell.basal_branches[0](0.3)._ref_v)
        self.rec_vbasal12.record(self.cell.basal_branches[0](0.5)._ref_v)
        self.rec_vbasal13.record(self.cell.basal_branches[0](0.7)._ref_v)
        self.rec_vbasal14.record(self.cell.basal_branches[0](0.9)._ref_v)


        '''record branches'''
        self.rec_v20.record(self.cell.branches[0](0.05)._ref_v)
        self.rec_v21.record(self.cell.branches[0](0.1)._ref_v)
        self.rec_v22.record(self.cell.branches[0](0.15)._ref_v)
        self.rec_v23.record(self.cell.branches[0](0.2)._ref_v)
        self.rec_v24.record(self.cell.branches[0](0.25)._ref_v)
        self.rec_v25.record(self.cell.branches[0](0.3)._ref_v)
        self.rec_v26.record(self.cell.branches[0](0.35)._ref_v)
        self.rec_v27.record(self.cell.branches[0](0.4)._ref_v)
        self.rec_v28.record(self.cell.branches[0](0.45)._ref_v)
        self.rec_v29.record(self.cell.branches[0](0.5)._ref_v)
        self.rec_v30.record(self.cell.branches[0](0.55)._ref_v)
        self.rec_v31.record(self.cell.branches[0](0.6)._ref_v)
        self.rec_v32.record(self.cell.branches[0](0.65)._ref_v)
        self.rec_v33.record(self.cell.branches[0](0.7)._ref_v)
        self.rec_v34.record(self.cell.branches[0](0.75)._ref_v)
        self.rec_v35.record(self.cell.branches[0](0.8)._ref_v)
        self.rec_v36.record(self.cell.branches[0](0.85)._ref_v)
        self.rec_v37.record(self.cell.branches[0](0.9)._ref_v)

        self.rec_distal = neuron.h.Vector()
        self.rec_distal.record(self.cell.second_branches[0](0.5)._ref_v)



        if self.cell.ifca:

            self.rec_icaL = neuron.h.Vector()
            self.rec_icaL1 = neuron.h.Vector()
            self.rec_icaL2 = neuron.h.Vector()
            self.rec_icaL3 = neuron.h.Vector()
            self.rec_icaL4 = neuron.h.Vector()
            self.rec_icaL5 = neuron.h.Vector()
            self.rec_icaL6 = neuron.h.Vector()
            self.rec_icaL7 = neuron.h.Vector()
            self.rec_icaL8 = neuron.h.Vector()
            self.rec_icaL9 = neuron.h.Vector()
            self.rec_icaL10 = neuron.h.Vector()
            self.rec_icaL11 = neuron.h.Vector()
            self.rec_icaL1.record(self.cell.apical_prox(0.2)._ref_ica)
            self.rec_icaL2.record(self.cell.apical_prox(0.6)._ref_ica)
            self.rec_icaL3.record(self.cell.apical_prox(1)._ref_ica)
            self.rec_icaL4.record(self.cell.apical(0.4)._ref_ica)
            self.rec_icaL5.record(self.cell.apical(0.8)._ref_ica)
            self.rec_icaL6.record(self.cell.branches[0](0.1)._ref_ica)
            self.rec_icaL7.record(self.cell.branches[0](0.5)._ref_ica)
            self.rec_icaL8.record(self.cell.branches[0](0.9)._ref_ica)
            self.rec_icaL9.record(self.cell.oblique_branch(0.1)._ref_ica)
            self.rec_icaL10.record(self.cell.oblique_branch(0.3)._ref_ica)
            self.rec_icaL11.record(self.cell.oblique_branch(0.5)._ref_ica)


            self.rec_ina0 = neuron.h.Vector()
            self.rec_ina1 = neuron.h.Vector()
            self.rec_ina2 = neuron.h.Vector()
            self.rec_ina3 = neuron.h.Vector()
            self.rec_ina4 = neuron.h.Vector()
            self.rec_ina5 = neuron.h.Vector()
            self.rec_ina6 = neuron.h.Vector()
            self.rec_ina7 = neuron.h.Vector()
            self.rec_ina8 = neuron.h.Vector()
            self.rec_ina9 = neuron.h.Vector()
            self.rec_ina10 = neuron.h.Vector()
            self.rec_ina11 = neuron.h.Vector()
            self.rec_ina0.record(self.cell.ais(0.5)._ref_ina)
            self.rec_ina1.record(self.cell.apical_prox(0.2)._ref_ina)
            self.rec_ina2.record(self.cell.apical_prox(0.6)._ref_ina)
            self.rec_ina3.record(self.cell.apical_prox(0.9)._ref_ina)
            self.rec_ina4.record(self.cell.apical(0.4)._ref_ina)
            self.rec_ina5.record(self.cell.apical(0.8)._ref_ina)
            self.rec_ina6.record(self.cell.branches[0](0.1)._ref_ina)
            self.rec_ina7.record(self.cell.branches[0](0.5)._ref_ina)
            self.rec_ina8.record(self.cell.branches[0](0.9)._ref_ina)
            self.rec_ina9.record(self.cell.oblique_branch(0.1)._ref_ina)
            self.rec_ina10.record(self.cell.oblique_branch(0.3)._ref_ina)
            self.rec_ina11.record(self.cell.oblique_branch(0.5)._ref_ina)


    def get_spike_dict(self):
        spike_dict = spiketrainutil.netconvecs_to_dict(self.t_vec, self.id_vec)
        return spike_dict

    def get_spike_times(self,voltage):
        voltage = voltage - (-70)
        above_threshold = voltage > 50 # gives a vector of zeros and ones
        threshold_crossings = np.diff(above_threshold) # gives ones when 1 changes to 0 or reverse
        strokes = np.nonzero(threshold_crossings)
        upstrokes = strokes[0][::2]
        print upstrokes
        downstrokes = strokes[0][1::2]
        spike_times = np.zeros(len(downstrokes))
        for i in range(len(downstrokes)):
            spike_time = upstrokes[i] + np.argmax(voltage[upstrokes[i]:downstrokes[i]])
            spike_times[i] = int(spike_time)
        spike_times = spike_times.astype(int)
        spike_times = spike_times * self.dt
        print "spike_times: " + str(spike_times)
        return spike_times


    def get_idx(self,time):
        "transforms time [ms] into timestep"
        idx = time*np.int(1/self.dt)
        return idx

    def get_spike_train(self,voltage):
        st = self.get_spike_times(voltage)
        spike_train = np.zeros(len(voltage))
        spike_train[st] = 1
        return spike_train


    def get_rate(self,start,end):
        idx_start = self.get_idx(start)
        idx_end = self.get_idx(end)
        voltage = np.array(self.rec_v) - (-70) # I used to measure this at the ais
        voltage = voltage[idx_start:idx_end]
        above_threshold = voltage > 50 # gives a vector of zeros and ones
        threshold_crossings = np.diff(above_threshold) # gives ones when 1 changes to 0 or reverse
        number_of_spikes = np.shape(np.nonzero(threshold_crossings))[1]*0.5
        if number_of_spikes == 0:
            return number_of_spikes, 0
        else:
            return number_of_spikes, (number_of_spikes * 1000.0) / (end-start)



    def get_ISI(self,voltage,start,end):
        idx_start = self.get_idx(start)
        idx_end = self.get_idx(end)
        spiketimes = self.get_spike_times(voltage[idx_start:idx_end])
        ISI = np.diff(spiketimes)
        return ISI


    def get_bursts(self, voltage, start, end):
        idx_start = self.get_idx(start)
        idx_end = self.get_idx(end)
        spiketimes = self.get_spike_times(voltage[idx_start:idx_end])
        spike_timings = spiketimes # spike_timings have unit mst

        ISI = np.diff(spike_timings)
        small_ISI = np.zeros(len(ISI))
        small_ISI[ISI<25] = 1
        small_ISI = np.append(np.array(0),small_ISI)

        ISI_boundaries = np.diff(small_ISI)
        idx = np.nonzero(ISI_boundaries)
        try:
            ISIstart = np.nonzero(ISI_boundaries==1)[0]
            ISIend = np.nonzero(ISI_boundaries==-1)[0]
            if not np.shape(ISIstart)[0] == np.shape(ISIend)[0]: # if there is not an endpoint for a startpoint
                ISIend = np.append(ISIend, np.array(len(spike_timings)-1))
                # add last spike as burst end, because there was no single spike after the last burst

        except IndexError:
            print "IndexError excepted in sim.get_bursts()"
        realbursts = (ISIend-ISIstart)>1
        ISIstart = ISIstart[realbursts]
        ISIend = ISIend[realbursts]
        print ISIstart
        print ISIend
        burststart = spike_timings[ISIstart]
        burstend = spike_timings[ISIend]
        print burststart
        print burstend
        print "end get_bursts"
        return ISIstart, ISIend, burststart, burstend

    def get_burst_stats(self, voltage, start, end):
        ISIstart, ISIend, burststart, burstend = self.get_bursts(voltage, start, end)
        diff = ISIend - ISIstart
        print "diff: " + str(diff)
        burstnumber = np.shape(np.nonzero(diff>1))[1]
        spikes_in_bursts = np.sum(diff[np.nonzero(diff>1)]) + burstnumber # spikes = ISI + 1
        print "burstnumber: " + str(burstnumber)
        print "spikes_in_bursts: " + str(spikes_in_bursts)
        if burstnumber == 0:
            burstrate = 0
        else:
            burstrate = (burstnumber * 1000.0) / (end-start)
        if spikes_in_bursts == 0:
            spikes_in_bursts_rate = 0
        else:
            spikes_in_bursts_rate = (spikes_in_bursts * 1000.0) / (end-start)

        return burstrate, spikes_in_bursts, spikes_in_bursts_rate


    def plot_bursts(self, voltage, start, end):
        ISIstart, ISIend, burststart, burstend = self.get_bursts(voltage, start, end)
        idx_burststart = self.get_idx(burststart)
        idx_burstend = self.get_idx(burstend)
        idx_start = self.get_idx(start)
        idx_end = self.get_idx(end)
        voltage = voltage[idx_start:idx_end]
        plot = Plot(self, params = P.plot)
        for i in np.arange(0,np.shape(burststart)[0]):
            print "voltage"
            print voltage[idx_burststart[i]-80:idx_burstend[i]+80]
            plot.voltage(voltage[idx_burststart[i]-80:idx_burstend[i]+80],'burst%d'%i)


    def plot_spikes(self,spikes):
        spike_list = list()
        for value in spikes.itervalues():
            print value
            spike_list.append(value)
        sp = spikeplot.SpikePlot(fig_name = "/extra/wilmes/figures/spikeplot.png")
        print spike_list
        sp.plot_spikes(spike_list, label=None, draw=True, savefig=True, cell_offset=0)


    def set_basal_recorders(self):
        '''record basal voltage for all compartments'''
        self.basallist = neuron.h.List()
        for seg in self.cell.basal_main:
            basal = neuron.h.Vector() # int(tstopms/self.timeres_python+1))
            basal.record(seg._ref_v)
            self.basallist.append(basal)

    def get_recording(self):
        self.recording = np.array((self.rec_v,self.rec_v1,self.rec_v2,self.rec_v3,
        self.rec_v4,self.rec_v5,self.rec_v6,self.rec_v7,self.rec_v8,self.rec_v9,self.rec_v10))
        return self.recording

    def get_otherbranch_recording(self):
        self.recording = np.array((self.rec_v,self.rec_v1,self.rec_v1a,self.rec_v2,
        self.rec_v2a,self.rec_v3,self.rec_v3a,self.rec_v4,self.rec_v4a,self.rec_v5,
        self.rec_v6o,self.rec_v7o,self.rec_v8o,self.rec_v9o,self.rec_v10o))
        return self.recording


    def get_highestres_recording(self):
        self.recording = np.array((self.rec_va, self.rec_v,self.rec_v1,self.rec_v2,self.rec_v3,
        self.rec_v4,self.rec_v5,self.rec_v6,self.rec_v7,self.rec_v8,self.rec_v9,
        self.rec_v20,self.rec_v21,
        self.rec_v22,self.rec_v23,self.rec_v24,self.rec_v25,self.rec_v26,self.rec_v27,
        self.rec_v28,self.rec_v29,self.rec_v30,self.rec_v31,self.rec_v32,self.rec_v33,
        self.rec_v34, self.rec_v35,self.rec_v36,self.rec_v37))

        return self.recording

    def get_oblique_recording(self):
        self.oblique_rec = np.array((self.rec_va, self.rec_v,self.rec_v2,
        self.rec_v4,self.rec_v6,self.rec_v8,self.rec_v9,self.rec_vo1,
        self.rec_vo2,self.rec_vo3, self.rec_vo4,self.rec_vo5,self.rec_vo6,
        self.rec_vo7,self.rec_vo8,self.rec_vo9))
        return self.oblique_rec

    def get_highres_oblique_recording(self):
        self.oblique_rec = np.array((self.rec_va, self.rec_v,self.rec_v1,self.rec_v2,self.rec_v3,self.rec_vo1,
        self.rec_vo2,self.rec_vo3, self.rec_vo4,self.rec_vo5,self.rec_vo6,
        self.rec_vo7,self.rec_vo8,self.rec_vo9))
        return self.oblique_rec


    def get_basal_recording(self):
        self.basal_rec = np.array((self.rec_va, self.rec_v, self.rec_vbasal1, self.rec_vbasal2,self.rec_vbasal3,
        self.rec_vbasal4,self.rec_vbasal5,self.rec_vbasal6,self.rec_vbasal7,self.rec_vbasal8,
        self.rec_vbasal9,self.rec_vbasal10, self.rec_vbasal11,self.rec_vbasal12,self.rec_vbasal13,
        self.rec_vbasal14))
        return self.basal_rec

    def get_basal_recording1(self):
        return self.basallist

    def get_ais_ina(self):
        return np.array((self.rec_ina0))

    def get_ina(self):
        self.ina = np.array((self.rec_ina1,self.rec_ina2,self.rec_ina3,
        self.rec_ina4,self.rec_ina5,self.rec_ina6,self.rec_ina7, self.rec_ina8,
        self.rec_ina9,self.rec_ina10,self.rec_ina11))
        return self.ina

    def get_ica(self):
        self.ical = np.array((self.rec_icaL1, self.rec_icaL2, self.rec_icaL3,
        self.rec_icaL4, self.rec_icaL5, self.rec_icaL6, self.rec_icaL7, self.rec_icaL8, self.rec_icaL9, self.rec_icaL10, self.rec_icaL11))
        return self.ical

    def get_ca_current(self,start,end):
        idx_start = self.get_idx(start)
        idx_end = self.get_idx(end)
        ca_current = np.array(self.rec_icaL7)
        return ca_current[idx_start:idx_end]


    def get_ica_simple(self):
        self.ica_simple = np.array((self.rec_icaL1, self.rec_icaL2, self.rec_icaL3,
        self.rec_icaL4, self.rec_icaL5, self.rec_icaL6, self.rec_icaL7, self.rec_icaL8,
        self.rec_icaL9, self.rec_icaL10, self.rec_icaL11))
        return self.ica_simple

    def get_time(self):
        return self.rec_t

    def get_stim_current(self):
        return self.rec_stim_i

    def get_stim_voltage(self):
        return self.rec_stim_vc

    def get_shunt_current(self):
        return self.rec_i

    def get_EPSC(self):
        return np.array((self.rec_EPSC))

    def get_shunt_conductance(self):
        return self.rec_g

    def get_distal_rec(self):
        return np.array((self.rec_distal))


    def get_synaptic_recording(self,kind):
        if kind == 'dendritic_inh':
            return np.array((self.Inh_distalRec))
        elif kind == 'dendritic_inh2':
            return np.array((self.Inh2_distalRec))
        elif kind == 'exc':
            return np.array((self.ExcRec))
        elif kind == 'inh':
            return np.array((self.InhRec))
        else:
            return ValueError

    def get_synaptic_current_rec(self):
        return np.array(self.syn_current)

    def get_synaptic_cond_rec(self):
        return np.array(self.syn_cond)

    def go(self, sim_time=None):
        #if self.high_res:
        #    self.set_highestres_recording()
        #    self.set_basal_recorders()
        #else:
        #    self.set_recording()
        neuron.h.dt = self.dt
        neuron.h.celsius = self.celsius

        if self.vector == True:
            if self.SE:
                self.v.play(self.stim._ref_amp1, self.v_time, 1)
            else:
                for i in range(self.stimnumber):
                    print "stimnumber: %d"%self.stimnumber
                    self.vlist[i].play(self.stimlist[i]._ref_amp, self.vtimelist[i], 1)
                    outime = np.array(self.vtimelist[i])
                    oustim = np.array(self.stimlist[i])
                    ouv = np.array(self.vlist[i])
                    print "ou"
                    print outime[::100]
                    print oustim[::100]
                    print ouv[::100]
                    print self.mu
        neuron.init()
        neuron.h.finitialize(self.v_init)

        h.fcurrent()


        neuron.run(self.sim_time)
        self.go_already = True