.. currentmodule:: brian

.. index::
   pair: example usage; sqrt

.. _example-frompapers_Muller_et_al_2011:

Example: Muller_et_al_2011 (frompapers)
=======================================

Interplay of STDP and input oscillations
----------------------------------------
Figure 4 from:
Muller L, Brette R and Gutkin B (2011) Spike-timing dependent plasticity and 
feed-forward input oscillations produce precise and invariant spike phase-locking. Front. 
Comput. Neurosci. 5:45. doi: 10.3389/fncom.2011.00045

Description:
In this simulation, a group of IF neurons is given a tonic DC input and a tonic AC input.
The DC input is mediated by current injection (neurons.I, line 62), and the AC input is 
mediated by Poisson processes whose rate parameters are oscillating in time. Each neuron in 
the group is given a different DC input, ensuring a unique initial phase. After two seconds 
of simulation (to integrate out any initial transients), the STDP rule is turned on 
(ExponentialSTDP, line 68), and the population of neurons converges to the theoretically 
predicted fixed point. As there is some noise in the phase due to the random inputs, the 
simulation is averaged over trials (50 in Figure 4, though 10 trials should be fine for 
testing).

The trials run in parallel on all available processors (10 trials take about
2 minutes on a modern PC).

::

    
    ### IMPORTS
    from brian import *
    import multiprocessing
    
    ### PARAMETERS
    N=5000
    M=10
    taum=33*ms                          
    tau_pre=20*ms
    tau_post=tau_pre
    Ee=0*mV
    vt=-54*mV
    vr=-70*mV
    El=-70*mV
    taue=5*ms
    f=20*Hz
    theta_period = 1/f
    Rm=200*Mohm
    a = linspace(51,65,num=M)       
    weights = .001
    ratio=1.50
    dA_pre=.01
    dA_post=.01*ratio 
    trials=10
    
    ### SIMULATION LOOP
    def trial(n): # n is the trial number
        reinit_default_clock()
        clear(True)
    
        eqs_neurons='''
        dv/dt=((ge*(Ee-vr))+Rm*I+(El-v))/taum : volt   
        dge/dt=-ge/taue : 1
        I : amp
        '''
    
        inputs = PoissonGroup(N,rates=lambda t:((.5-.5*cos(2*pi*f*t)))*10*Hz)           
        neurons=NeuronGroup(M,model=eqs_neurons,threshold=vt,reset=vr)
        neurons.I = a*pA
        synapses=Connection(inputs,neurons,'ge',weight=weights)
        neurons.v=vr
    
        S = SpikeMonitor(neurons)
        run(2*second)
        stdp=ExponentialSTDP(synapses,tau_pre,tau_post,dA_pre,-dA_post,wmax=10*weights,interactions='all',update='additive')     
        run(5*second)
        
        phase=zeros((M,200))
        for b in range(0,M):    
            tmp_phase=(S[b]%theta_period)*(360/theta_period)
            phase[b,range(0,len(tmp_phase))] = tmp_phase
            
        return phase
    
    if __name__=='__main__': # This is very important on Windows, otherwise the machine crashes!
        phase = zeros((M,200,trials))
        
        print "This will take approximately 2 minutes."
        pool=multiprocessing.Pool() # uses all available processors b
        results=pool.map(trial,range(trials))
        for i in range(trials):
            phase[:,:,i]=results[i]
        
        ### PLOTTING
        for b in range(0,M):
            m = mean(phase[b,:,:],axis=1)
            st = std(phase[b,:,:],axis=1)/sqrt(trials)
            errorbar(range(0,135), m[range(0,135)], yerr=st[range(0,135)], xerr=None,
                 fmt='-', ecolor=None, elinewidth=None, capsize=3,
                 barsabove=False, lolims=False, uplims=False,
                 xlolims=False, xuplims=False)
        
        title('STDP + Oscillations Simulation')
        xlabel('Spike Number')
        ylabel('Spike Phase (deg)')
        xlim([0, 135])
        ylim([140, 280])
        
        show()