'''
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()