Hierarchical network model of perceptual decision making (Wimmer et al 2015)

 Download zip file 
Help downloading and running models
Accession:168867
Neuronal variability in sensory cortex predicts perceptual decisions. To investigate the interaction of bottom-up and top-down mechanisms during the decision process, we developed a hierarchical network model. The network consists of two circuits composed of leaky integrate-and-fire neurons: an integration circuit (e.g. LIP, FEF) and a sensory circuit (MT), recurrently coupled via bottom-up feedforward connections and top-down feedback connections. The integration circuit accumulates sensory evidence and produces a binary categorization due to winner-take-all competition between two decision-encoding populations (X.J. Wang, Neuron, 2002). The sensory circuit is a balanced randomly connected EI-network, that contains neural populations selective to opposite directions of motion. We have used this model to simulate a standard two-alternative forced-choice motion discrimination task.
Reference:
1 . Wimmer K, Compte A, Roxin A, Peixoto D, Renart A, de la Rocha J (2015) Sensory integration dynamics in a hierarchical network explains choice probabilities in cortical area MT. Nat Commun 6:6177 [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: Neocortex;
Cell Type(s): Abstract integrate-and-fire leaky neuron;
Channel(s):
Gap Junctions:
Receptor(s): AMPA; NMDA; Gaba;
Gene(s):
Transmitter(s):
Simulation Environment: Brian; Python;
Model Concept(s): Attractor Neural Network; Winner-take-all;
Implementer(s): Wimmer, Klaus [wimmer.klaus at gmail.com];
Search NeuronDB for information about:  AMPA; NMDA; Gaba;
'''

Hierarchical network model of perceptual decision making

@author: Klaus Wimmer

wimmer.klaus@googlemail.com

'''


from brian import *
import numpy
import random as pyrandom
from numpy.random import rand as rand
from numpy.random import randn as randn
from scipy.signal import lfilter

from integration_circuit import make_integration_circuit
from sensory_circuit import make_sensory_circuit
 

def get_OU_stim(n, tau):
# UO process in discrete time => AR(1) process
      
    a = numpy.exp(-(1.0 / tau))
    i = lfilter(numpy.ones(1),[1.0, -a], numpy.sqrt(1-a*a)*randn(n))   
         
    return i


if __name__ == '__main__':

    #  initialize  
    defaultclock.reinit()
    clear(True) 


    #------------------------------------------------------------------------------ 
    # Simulation parameters 
    #------------------------------------------------------------------------------ 
    
    connect_seed = 1284                          # seed for random number generators (set before establishing network connectivity)
    stim_seed = 123                              # seed for random number generators (set before generating common part of stimulus)
    init_seed = 8190                             # seed for random number generators (set before generating private part of stimulus)

    # Timing 
    stim_on = 500.0 * ms                         # stimulus onset
    stim_off = 2500.0 * ms                       # stimulus offset   
    stim_duration = stim_off - stim_on           # duration of stimulus interval
    runtime = 3000.0 * ms                        # total simulation time


    #------------------------------------------------------------------------------
    # Construct hierarchical network
    #------------------------------------------------------------------------------ 

    # Set the seed of the random number generator
    numpy.random.seed(connect_seed) 
    pyrandom.seed(connect_seed)

    # Integration circuit
    Dgroups, Dconnections, Dnetfunctions, Dsubgroups = make_integration_circuit()
    
    # get populations from the integrations circuit
    decisionE = Dgroups['DE']
    decisionI = Dgroups['DI']
    decisionE1 = Dsubgroups['DE1']
    decisionE2 = Dsubgroups['DE2']
    decisionE3 = Dsubgroups['DE3']

    # Sensory network  
    Sgroups, Sconnections, Ssubgroups = make_sensory_circuit()

    # get sensory populations
    sensoryE = Sgroups['SE']
    sensoryI = Sgroups['SI']
    sensoryE1 = Ssubgroups['SE1']
    sensoryE2 = Ssubgroups['SE2']
    
    # Feed-forward connections from the sensory to the integration circuit
    
    wSD = 0.0036                                 # Synaptic weight of feed-forward connections from the corresponding stimulus-encoding population 
                                                 # from the sensory circuit (E1 -> D1, E2 -> D2); 
                                                 # synaptic weight (0.09 nS) is given in multiples of the leak conductance of the excitatory neurons in the integration circuit  
                                                 
    C_SE1_DE1 = Connection(sensoryE1, decisionE1, 'gea', weight=wSD, sparseness=0.2, delay=1.0 * ms)
    C_SE2_DE2 = Connection(sensoryE2, decisionE2, 'gea', weight=wSD, sparseness=0.2, delay=1.0 * ms)
	
    # Top-down feedback connections from the integration circuit to the sensory circuit
    
    b_FB = 0.0                                   # Feedback strength
    wDS = 0.004 * b_FB                           # Synaptic weight of feedback connections from the integration circuit to the sensory circuit (D1 -> E1, D2 -> E2);
                                                 # synaptic weight (0.0668 nS) is given in multiples of the leak conductance of the excitatory neurons in the sensory circuit  
    
    C_DE1_SE1 = Connection(decisionE1, sensoryE1, 'xe', weight=wDS, sparseness=0.2, delay=1.0 * ms)
    C_DE2_SE2 = Connection(decisionE2, sensoryE2, 'xe', weight=wDS, sparseness=0.2, delay=1.0 * ms)


    #------------------------------------------------------------------------------
    # Stimulus 
    #------------------------------------------------------------------------------ 
    
    # Stimulus parameters    
    I0 = 0.08 * nA                               # Mean input current for zero-coherence stimulus 
    c = 0.0                                      # Stimulus coherence (between 0 and 1) 
    mu_E1 = +0.25                                # Average additional input current to E1 at highest coherence (c = 1)
    mu_E2 = -0.25                                # Average additional input current to E2 at highest coherence (c = 1)    
    sigma = 1.0																	 # Amplitude of temporal modulations of the stimulus
    sigma_stim = 0.212 * sigma                   # S.d. of modulations of stimulus inputs
    sigma_ind = 0.212 * sigma                    # S.d. of modulations in individual inputs
    tau_stim = 20.0 * ms                         # Correlation time constant of Ornstein-Uhlenbeck process
    
    # Generate stimulus
    # set seed of random number generator (in order to generate a specific stimulus each time)
    numpy.random.seed(stim_seed)    
    pyrandom.seed(stim_seed)

    # "common part" of the stimulus
    z1 = get_OU_stim(stim_duration/ms, tau_stim/ms)
    z1 = numpy.tile(z1,(len(sensoryE1),1))
    z2 = get_OU_stim(stim_duration/ms, tau_stim/ms)
    z2 = numpy.tile(z2,(len(sensoryE2),1))

    # set seed of random number generator (in order to generate a specific stimulus each time)
    numpy.random.seed(init_seed)    
    pyrandom.seed(init_seed)

    # "private part" - part of the stimulus for each neuron, different in each trial               
    zk1 = get_OU_stim(int(stim_duration/ms * len(sensoryE1)), tau_stim/ms)  
    zk1 = numpy.asarray(zk1).reshape(len(sensoryE1), stim_duration/ms)
    zk2 = get_OU_stim(int(stim_duration/ms * len(sensoryE2)), tau_stim/ms)  
    zk2 = numpy.asarray(zk2).reshape(len(sensoryE2), stim_duration/ms)
    
    # stimulus (time series with dt = 1ms)
    # most general case: different input to each neuron in each time step 
    i1 = I0 * (1 + c * mu_E1 + sigma_stim * z1 + sigma_ind * zk1)     
    i2 = I0 * (1 + c * mu_E2 + sigma_stim * z2 + sigma_ind * zk2)     
    ii = numpy.zeros((len(sensoryI),stim_duration/ms))     # no current input to inh population

    # Stimulus-related external current input     
    myclock=Clock(dt=1*ms)
    @network_operation(myclock)
    def update_input():  
        if myclock.t >= stim_on and myclock.t < stim_off:           
            sensoryE1.I = i1[:,int( (myclock.t - stim_on) / (1 * ms))] * amp
            sensoryE2.I = i2[:,int( (myclock.t - stim_on) / (1 * ms))] * amp
            sensoryI.I = ii[:,int( (myclock.t - stim_on) / (1 * ms))] * amp
        else:
            sensoryE1.I = 0 * nA
            sensoryE2.I = 0 * nA
            sensoryI.I = 0 * nA
        
          
    #------------------------------------------------------------------------------
    # Initial conditions and Monitors
    #------------------------------------------------------------------------------

    # --- set seed of random number generator to a different value in each run
    np_seed = int(1587.47)
    numpy.random.seed(np_seed)
    py_seed = int(4736.28) 
    pyrandom.seed(py_seed)

    # ---- set initial conditions (random)
    decisionE.gen = decisionE.gen * (1 + 0.2 * rand(decisionE.__len__()))
    decisionI.gen = decisionI.gen * (1 + 0.2 * rand(decisionI.__len__()))
    decisionE.V = decisionE.V + rand(decisionE.__len__()) * 2 * mV
    decisionI.V = decisionI.V + rand(decisionI.__len__()) * 2 * mV

    # ---- set initial conditions (random)
    sensoryE.V = -50.0 * mV - 2 * mV + rand(sensoryE.__len__()) * 2 * mV
    sensoryI.V = -50.0 * mV - 2 * mV + rand(sensoryI.__len__()) * 2 * mV
    sensoryE.gea = 0.05 * (1 + rand(sensoryE.__len__()) * 0.2)
    sensoryI.gea = 0.05 * (1 + rand(sensoryI.__len__()) * 0.2)

    # record spikes of excitatory neurons
    S_DE1 = SpikeMonitor(decisionE1, record=True)
    S_DE2 = SpikeMonitor(decisionE2, record=True)
    S_SE1 = SpikeMonitor(sensoryE1, record=True)
    S_SE2 = SpikeMonitor(sensoryE2, record=True)

    # record instantaneous populations activity
    R_DE1 = PopulationRateMonitor(decisionE1, bin=20*ms)
    R_DE2 = PopulationRateMonitor(decisionE2, bin=20*ms)
    R_SE1 = PopulationRateMonitor(sensoryE1, bin=20*ms)
    R_SE2 = PopulationRateMonitor(sensoryE2, bin=20*ms)
 
 
    #------------------------------------------------------------------------------
    # Run the simulation
    #------------------------------------------------------------------------------

    # construct network       
    net = Network(Dgroups.values(), Sgroups.values(), Dconnections.values(), Sconnections.values(), 
                  Dnetfunctions, update_input, C_SE1_DE1, C_SE2_DE2, C_DE1_SE1, C_DE2_SE2,
                  S_DE1, S_DE2, S_SE1, S_SE2, R_DE1, R_DE2, R_SE1, R_SE2)
    net.prepare()
    net.run(runtime)        
      
       
    #------------------------------------------------------------------------------
    # Show results (single simulation run similar to Fig. 1b)
    #------------------------------------------------------------------------------

    fig, axs = subplots(7,1, figsize=(4,9))

		# Integration circuit	
    fig.add_axes(axs[0])
    raster_plot(S_DE1,color=(1,0,0),title="Integration circuit",xlabel="",ylabel="")
    xlim(0,runtime/ms)   
    ylim(0,len(decisionE1))
    fig.add_axes(axs[1])
    raster_plot(S_DE2,color=(0,0,1),title="",xlabel="",ylabel="")
    xlim(0,runtime/ms)   
    ylim(0,len(decisionE2))
    fig.add_axes(axs[2])
    plot(R_DE1.times/ms,R_DE1.rate,color=(1,0,0))
    plot(R_DE2.times/ms,R_DE2.rate,color=(0,0,1))
    yticks(range(0,50,10))
    xlim(0,runtime/ms)   
    ylim(0,45)
    ylabel("Rate (sp/s)")

    # Sensory circuit
    fig.add_axes(axs[3])
    raster_plot(S_SE1,color=(1,0,0),title="Sensory circuit",xlabel="",ylabel="")
    xlim(0,runtime/ms)   
    ylim(0,len(sensoryE1))
    fig.add_axes(axs[4])
    raster_plot(S_SE2,color=(0,0,1),title="",xlabel="")
    xlim(0,runtime/ms)   
    ylim(0,len(sensoryE2))
    fig.add_axes(axs[5])
    plot(R_SE1.times/ms,R_SE1.rate,color=(1,0,0))
    plot(R_SE2.times/ms,R_SE2.rate,color=(0,0,1))
    yticks(range(0,30,5))
    xlim(0,runtime/ms)   
    ylim(0,20)
    ylabel("Rate (sp/s)")
    
    # Stimulus
    fig.add_axes(axs[6])
    t = linspace(0., runtime/second, runtime/ms)
    plot(t, numpy.r_[zeros(stim_on/ms), mean(i1,0)/I0, zeros((runtime-stim_off)/ms)], color=(1,0,0))
    plot(t, numpy.r_[zeros(stim_on/ms), mean(i2,0)/I0, zeros((runtime-stim_off)/ms)], color=(0,0,1))
    xticks(range(0,3))
    yticks(range(0,2,1))
    xlim(0,runtime/second)
    ylim(0,1.5)
    xlabel("Time (s)")
    ylabel("Stimulus")

    for i in [0,1,3,4]:
        axs[i].set_yticklabels([])
    for i in range(0,6):
        axs[i].set_xticklabels([])   
    show()