Conversion olfactory model from brian1 to brian2

Support
1

Description of problem

I’m trying to convert model from the article Computing with Neural Synchrony. I’dont understand how is used with brian2 connectivity which is used at this model.

Minimal code to reproduce problem

from brian2 import *

bmin,bmax=-7,-1

def odor(N):
    # Returns a random vector of binding constants
    return 10**(rand(N)*(bmax-bmin)+bmin)

def hill_function(c,K=1.,n=3.):
    '''
    Hill function:
    * c = concentration
    * K = half activation constant (choose K=1 for relative concentrations)
    * n = Hill coefficient
    '''
    return (c**n)/(c**n+K**n)

N=5000 # number of receptors

# Odors
seed(31415) # Get the same neurons every time
intensity=3000.
c1=odor(N)
c2=odor(N)
c0=c1
I1,I2=intensity,intensity

# Odor plumes (fluctuating concentrations)
tau_plume=75*ms
eq_plumes='''
dx/dt=-x/tau_plume+(2./tau_plume)**.5*xi : 1
y=clip(x,0,inf) : 1
'''
plume=NeuronGroup(2,model=eq_plumes) # 2 odors

# Receptor neurons
Fmax=40*Hz # maximum firing rate
tau=20*ms
Imax=1/(1-exp(-1/(Fmax*tau))) # maximum input current

eq_receptors='''
dv/dt=(Imax*hill_function(c)-v)/tau : 1
c : 1  # concentrations (relative to activation constant)
'''

receptors=NeuronGroup(N,model=eq_receptors,threshold='v > 1',reset='v = 0')
receptors.c=c1

@network_operation
def odor_to_nose():
    # Send odor plume to the receptors
    receptors.c=I1*c1*clip(plume.x[0],0,Inf)+I2*c2*clip(plume.x[1],0,Inf)

# Decoder neurons
M=200
taud=8*ms
sigma=.15
eq_decoders='''
dv/dt=-v/taud + sigma*(2/taud)**.5*xi : 1
'''
decoders=NeuronGroup(2*M,model=eq_decoders,threshold='v > 1',reset='v = 0')
# First M neurons encode odor A, next M neurons encode odor B

# Synapses
syn = Synapses(receptors,decoders)
#syn=Connection(receptors,decoders,'v')
syn.connect()

# Connectivity according to synchrony partitions
bhalf=.5*(bmin+bmax) # select only those that are well activated
u=2*(log(c1)/log(10)-bhalf)/(bmax-bmin) # normalized binding constants for odor A
for i in range(M):
    which=((u>=i*1./M) & (u<(i+1)*1./M)) # we divide in M groups with similar values
    if sum(which)>0:
        w=1./sum(which) # total synaptic weight for a postsynaptic neuron is 1
    syn[:,i]=w*which

u=2*(log(c2)/log(10)-bhalf)/(bmax-bmin)
for i in range(M): # normalized binding constants for odor B
    which=((u>=i*1./M) & (u<(i+1)*1./M))
    if sum(which)>0:
        w=1./sum(which)
    syn[:,2*M-1-i]=w*which

# Record odor concentration and output spikes
O=StateMonitor(plume,'y',record=True)
S=SpikeMonitor(receptors)
S2=SpikeMonitor(decoders)

print("Odor A")
I1,I2=intensity,0
run(2*second)
print("Odor B")
I1,I2=0,intensity
run(2*second)
print("Odor B x2")
I1,I2=0,2*intensity
run(2*second)
print("Odor A + odor C")
I1,I2=intensity,intensity
old_c2=c2
c2=odor(N) # different odor
run(2*second)
print("Odor A + odor B")
I1,I2=intensity,intensity
c2=old_c2
run(2*second)

t=O.times/ms

# Figure (9B)
subplot(311) # odor fluctuations
plot(t[t<2000],O[0][t<2000],'b')
plot(t[(t>=2000) & (t<4000)],O[1][(t>=2000) & (t<4000)],'r')
plot(t[(t>=4000) & (t<6000)],2*O[1][(t>=4000) & (t<6000)],'r')
plot(t[(t>=6000) & (t<8000)],O[0][(t>=6000) & (t<8000)],'b')
plot(t[(t>=6000) & (t<8000)],O[1][(t>=6000) & (t<8000)],'k')
plot(t[(t>=8000) & (t<10000)],O[1][(t>=8000) & (t<10000)],'r')
plot(t[(t>=8000) & (t<10000)],O[0][(t>=8000) & (t<10000)],'b')
xlim(0,10000)
xticks([])
subplot(312)
raster_plot(S)
xlim(0,10000)
ylim(2500,2600) # 100 random neurons
xticks([])
subplot(313)
raster_plot(S2)
ylim(100,300)
xlim(0,10000)

show()

from brian2 import *
import numpy

N=5000
# Coincidence detectors
sigma=.15
taud=8*ms
# Connections
Nsynapses=50
w0=150./(0.02*N)
# STDP
factor=0.05
a_pre=0.06*factor
b_post=-1.*factor
tau_pre=3*ms
# Intrinsic plasticity: non-specific weight increase
IP_period=10*ms
IP_rate=-b_post*5*Hz # target firing rate = 5 Hz
# Simulation control
record_period=1*second
duration=100*second

bmin,bmax=-7,-1

def odor(N):
    # Returns a random vector of binding constants
    return 10**(rand(N)*(bmax-bmin)+bmin)

def hill_function(c,K=1.,n=3.):
    '''
    Hill function:
    * c = concentration
    * K = half activation constant (choose K=1 for relative concentrations)
    * n = Hill coefficient
    '''
    return (c**n)/(c**n+K**n)

N=5000 # number of receptors

seed(31415) # Get the same neurons every time
intensity=3000.

# Odor plumes
tau_plume=75*ms
eq_plumes='''
dx/dt=-x/tau_plume+(2./tau_plume)**.5*xi : 1
y=clip(x,0,inf) : 1
'''
plume=NeuronGroup(1,model=eq_plumes) # 1 odor

# Receptor neurons
Fmax=40*Hz # maximum firing rate
tau=20*ms
Imax=1/(1-exp(-1/(Fmax*tau))) # maximum input current

eq_receptors='''
dv/dt=(Imax*hill_function(c)-v)/tau : 1
c : 1  # concentrations (relative to activation constant)
'''

receptors=NeuronGroup(N,model=eq_receptors,threshold='v>1',reset='v=0')

@network_operation
def odor_to_nose():
    # Send odor plume to the receptors
    receptors.c=I1*c1*clip(plume.x[0],0,Inf)

odors=[odor(N),odor(N)] # two odors
c1=odors[0]
stimuli=[]
# A random odor is presented every 200 ms
@network_operation(clock=EventClock(dt=200*ms))
def change_odor():
    global c1
    nodor=randint(len(odors))
    c1=odors[nodor]
    stimuli.append((float(defaultclock.t),float(nodor)))

# Decoder neurons
M=30
eq_decoders='''
dv/dt=-v/taud + sigma*(2/taud)**.5*xi : 1
'''
decoders=NeuronGroup(M,model=eq_decoders,threshold='v>1',reset='v=0')
S2=SpikeMonitor(decoders)

# Random synapses
syn=Connection(receptors,decoders,'v',sparseness=Nsynapses*1./N,weight=w0)

# STDP
eqs_stdp='''
dApre/dt=-Apre/tau_pre : 1
Apost : 1
'''
pre='''
Apre+=a_pre
#w+=0
'''
post='''
Apost+=0
w+=Apre+b_post*w
'''
stdp=STDP(syn,eqs_stdp,pre=pre,post=post,wmax=Inf)
MC=ConnectionMonitor(syn,store=True,clock=EventClock(dt=record_period))

@network_operation(EventClock(dt=IP_period))
def intrinsic_plasticity(): # synaptic scaling
    # Increases weights of all synapses
    syn.W.alldata+=syn.W.alldata*IP_rate*IP_period

# Record the evolution of weights
weights=[]
@network_operation(EventClock(dt=record_period))
def recordW():
    Z=syn.W[:,0].copy()
    weights.append(Z)

I1=intensity
print "Started"
run(duration,report="text")

# Save data
wsave=[(t,M.todense()) for (t,M) in MC.values]
numpy.save("weights.npy",array(zip(*wsave)[1])) # 3D array (t,i,j)
numpy.save("spikesout.npy",array(S2.spikes))
numpy.save("stimuli.npy",array(stimuli))

What you have aready tried

Expected output (if relevant)

Actual output (if relevant)

Full traceback of error (if relevant)

2

Hi @serge25. Please have a look at the documentation for the Brian1 → Brian2 conversion, in particular at Synapses (Brian 1 –> 2 conversion) — Brian 2 2.5.1 documentation. You will also most likely not need any of the network_operations anymore, since Brian’s TimedArray and run_regularly are now more powerful to change things during the simulation (Input stimuli — Brian 2 2.5.1 documentation). Also, the network_operation storing the synaptic weights can be replaced by a simple StateMonitor. Let us know if anything in the documentation is unclear and/or if you run into any issues.

But as a general remark: I would strongly suggest not to start with the full model. Instead, do it step by step, e.g. try to first recreate the receptor neurons, and if they are working convert the decoder neurons, then add the synapses, then add the plasticity, etc. Otherwise you would have to work with a broken model until you fix all the issues, and cannot run anything.

3

Thank you for your answer. My first question (i don’t find at documentation): how implement user function in equation at NeuronGroup ?

def hill_function(c,K=1.,n=3.):
    '''
    Hill function:
    * c = concentration
    * K = half activation constant (choose K=1 for relative concentrations)
    * n = Hill coefficient
    '''
    return (c**n)/(c**n+K**n)

eq_receptors='''
dv/dt=(Imax*hill_function(c)-v)/tau : 1
c : 1  # concentrations (relative to activation constant)
'''
4

You don’t need to implement this as a user function (not quite sure why/whether it was necessary in Brian 1), you can directly write: dv/dt = (Imax*(c**n)/(c**n+K**n) - v)/tau : 1 in the equations.

5

Thank you. But what about this part:

# Connectivity according to synchrony partitions
bhalf=.5*(bmin+bmax) # select only those that are well activated
u=2*(log(c1)/log(10)-bhalf)/(bmax-bmin) # normalized binding constants for odor A
for i in range(M):
    which=((u>=i*1./M) & (u<(i+1)*1./M)) # we divide in M groups with similar values
    if sum(which)>0:
        w=1./sum(which) # total synaptic weight for a postsynaptic neuron is 1
    syn[:,i]=w*which

u=2*(log(c2)/log(10)-bhalf)/(bmax-bmin)
for i in range(M): # normalized binding constants for odor B
    which=((u>=i*1./M) & (u<(i+1)*1./M))
    if sum(which)>0:
        w=1./sum(which)
    syn[:,2*M-1-i]=w*which
6

There’s probably a way to write this more compactly, but as explained in the page I linked earlier, in Brian 2 you have to assign to a weight variable that you declare yourself in the equations. Then, you’d replace e.g. syn[:, i] = w*which by syn.w[:, i] = w*which (assuming you called the weight variable w).

7

Thank you @mstimberg . The first part of the model now works well. But there are questions on the second part:
how can I translate this piece of code into Brian2

@network_operation(EventClock(dt=IP_period))
def intrinsic_plasticity(): # synaptic scaling
    # Increases weights of all synapses
    syn.w+=syn.w*IP_rate*IP_period

and saving data:

wsave=[(t,M.todense()) for (t,M) in Syn.w]
numpy.save("weights.npy",array(zip(*wsave)[1])) # 3D array (t,i,j)
numpy.save("spikesout.npy",array(S2.spikes))
numpy.save("stimuli.npy",array(stimuli))

spikes=array(S2.spikes) # i,t
n,t=spikes[:,0],spikes[:,1]
9

Please tell me @mstimberg if I translated the code from Brian 1 to Brian 2 correctly:

Brian1:

MC=ConnectionMonitor(syn,store=True,clock=EventClock(dt=record_period))
wsave=[(t,M.todense()) for (t,M) in MC.values]
numpy.save("weights.npy",array(zip(*wsave)[1])) # 3D array (t,i,j)
numpy.save("spikesout.npy",array(S2.spikes))
numpy.save("stimuli.npy",array(stimuli))
..................................................................................
t,odor=numpy.load("stimuli.npy").T
W=numpy.load("weights.npy")
spikes_out=numpy.load("spikesout.npy")
weights=W[-1,:,:] # final weights

# Analyze selectivity at the end of the STDP simulation
ispikes=spikes_out[:,0] # indexes of neurons that spiked
tspikes=spikes_out[:,1] # spike timings

Brian2:

# Create a matrix to store the weights and fill it with NaN
Wout = np.full((len(receptors), len(decoders)), np.nan)
# Insert the values from the Synapses object
Wout[syn.i[:], syn.j[:]] = syn.w[:]

np.save("weights.npy",Wout)
np.save("spikesout.npy",array((S2.i[:],S2.t[:]/ms)))
np.save("stimuli.npy",array(stimuli))
..............................................
# Loads information from the STDP simulation
t,odor=np.load("stimuli.npy").T
Win=np.load("weights.npy")
spikes_out=np.load("spikesout.npy")
weights=Win[:,:] # final weights

# Analyze selectivity at the end of the STDP simulation
ispikes=spikes_out[0,:] # indexes of neurons that spiked
tspikes=spikes_out[1,:] # spike timings

I ask because I do not fully understand how the previous version of the simulator was stored in these objects and the simulation of the model is not yet obtained.

10

Thank you for your attention. Everything worked out.
The only question is how to write the following piece of code more elegantly and efficiently, using the capabilities of the simulator:

# Connectivity according to synchrony partitions
bhalf=.5*(bmin+bmax) # select only those that are well activated
u=2*(log(c1)/log(10)-bhalf)/(bmax-bmin) # normalized binding constants for odor A
for i in range(M):
    which=((u>=i*1./M) & (u<(i+1)*1./M)) # we divide in M groups with similar values
    if sum(which)>0:
        w=1./sum(which) # total synaptic weight for a postsynaptic neuron is 1
    syn.w[:, i] = w*which

u=2*(log(c2)/log(10)-bhalf)/(bmax-bmin)
for i in range(M): # normalized binding constants for odor B
    which=((u>=i*1./M) & (u<(i+1)*1./M))
    if sum(which)>0:
        w=1./sum(which)
    syn.w[:,2*M-1-i]=w*which
11

Hi @serge25. Sorry for not having followed the forum closely lately, but glad that you figured things out on your own! Regarding the matrix that you store in weights.npy, note that in the Brian1 version, it has values for several time points, whereas in the Brian2 version, you are only storing the final weight matrix. The later analysis code only seems to use the final weight matrix, so that should not be an issue.

Regarding the connection code, one major disadvantage of the code is that it creates all the synapses, and only sets the weights of a small number of them to a non-zero value. This is very inefficient in Brian2, both for the construction phase, and for the simulation itself. Instead, you can only create the necessary synapses, and then set the weights to a normalized value (this is what w./sum(which) is doing in the above code). Here’s how I’d write this in Brian 2, storing the binding constants for the two odors directly in the neuron model for the receptors (note that this is not the full code, it is meant to be added to the existing definitions of receptors, syn, etc.:

receptors = NeuronGroup(N, '''
                           ...
                           binding_A : 1 (constant) # normalized binding constants for odor A
                           binding_B : 1 (constant) # normalized binding constants for odor B
                           ''', ...)
bhalf=.5*(bmin+bmax)
receptors.binding_A = '2*(log(10**(rand()*(bmax - bmin) + bmin))/log(10) - bhalf)/(bmax - bmin)'
receptors.binding_B = '2*(log(10**(rand()*(bmax - bmin) + bmin))/log(10) - bhalf)/(bmax - bmin)'
decoders = NeuronGroup(2*M, ...)
syn = Synapses(receptors, decoders, 'w : 1', ...)

# Partition the decoder neurons into groups
syn.connect('j < M  and binding_A_pre >= j * 1.0/M         and binding_A_pre < (j + 1) * 1.0/M')
syn.connect('j >= M and binding_B_pre >= (2*M-1-j) * 1.0/M and binding_B_pre < (2*M-j) * 1.0/M')
# Normalize weights by number of incoming synapses
syn.w = '1.0 / N_incoming'

Hope that makes some sense! Note that the above code already includes what was previously defined in the odor function, i.e. that function is no longer necessary.

12

Thank you @mstimberg . I added this piece of code to the project, but the response of the model does not match the original one. What could be wrong?

from brian2 import *

bmin,bmax=-7,-1

def odor(N):
    # Returns a random vector of binding constants
    return 10**(rand(N)*(bmax-bmin)+bmin)


N=5000 # number of receptors

# Odors
seed(31415) # Get the same neurons every time
intensity=3000.
c1=odor(N)
c2=odor(N)
c0=c1
I1,I2=intensity,intensity

K=1.
n=3.

# Odor plumes (fluctuating concentrations)
tau_plume=75*ms
eq_plumes='''
dx/dt=-x/tau_plume+(2./tau_plume)**.5*xi : 1
y=clip(x,0,inf) : 1
'''
plume=NeuronGroup(2,model=eq_plumes) # 2 odors

# Receptor neurons
Fmax=40*Hz # maximum firing rate
tau=20*ms
Imax=1/(1-exp(-1/(Fmax*tau))) # maximum input current

eq_receptors='''
dv/dt=(Imax*(c**n)/(c**n+K**n)-v)/tau : 1
c : 1  # concentrations (relative to activation constant)
binding_A : 1 (constant) # normalized binding constants for odor A
binding_B : 1 (constant) # normalized binding constants for odor B
'''

receptors=NeuronGroup(N,model=eq_receptors,threshold='v > 1',reset='v = 0')
receptors.c=c1

@network_operation
def odor_to_nose():
    # Send odor plume to the receptors
    receptors.c=I1*c1*clip(plume.x[0],0,Inf)+I2*c2*clip(plume.x[1],0,Inf)

# Decoder neurons
M=200
taud=8*ms
sigma=.15

bhalf=.5*(bmin+bmax)

receptors.binding_A = '2*(log(10**(rand()*(bmax - bmin) + bmin))/log(10) - bhalf)/(bmax - bmin)'

receptors.binding_B = '2*(log(10**(rand()*(bmax - bmin) + bmin))/log(10) - bhalf)/(bmax - bmin)'

eq_decoders='''
dv/dt=-v/taud + sigma*(2/taud)**.5*xi : 1
'''
decoders=NeuronGroup(2*M,model=eq_decoders,threshold='v > 1',reset='v = 0')
# First M neurons encode odor A, next M neurons encode odor B

# Synapses
syn = Synapses(receptors,decoders, model='w : 1', on_pre='v +=w')


# Partition the decoder neurons into groups

syn.connect('j < M  and binding_A_pre >= j * 1.0/M         and binding_A_pre < (j + 1) * 1.0/M')

syn.connect('j >= M and binding_B_pre >= (2*M-1-j) * 1.0/M and binding_B_pre < (2*M-j) * 1.0/M')

# Normalize weights by number of incoming synapses

syn.w = '1.0 / N_incoming'

# Record odor concentration and output spikes
O=StateMonitor(plume,'y',record=True)
S=SpikeMonitor(receptors)
S2=SpikeMonitor(decoders)

print("Odor A")
I1,I2=intensity,0
run(2*second)
print("Odor B")
I1,I2=0,intensity
run(2*second)
print("Odor B x2")
I1,I2=0,2*intensity
run(2*second)
print("Odor A + odor C")
I1,I2=intensity,intensity
old_c2=c2
c2=odor(N) # different odor
run(2*second)
print("Odor A + odor B")
I1,I2=intensity,intensity
c2=old_c2
run(2*second)

t=O.t/ms

subplot(311) # odor fluctuations
plot(t[t<2000],O.y[0][t<2000],'b')
plot(t[(t>=2000) & (t<4000)],O.y[1][(t>=2000) & (t<4000)],'r')
plot(t[(t>=4000) & (t<6000)],2*O.y[1][(t>=4000) & (t<6000)],'r')
plot(t[(t>=6000) & (t<8000)],O.y[0][(t>=6000) & (t<8000)],'b')
plot(t[(t>=6000) & (t<8000)],O.y[1][(t>=6000) & (t<8000)],'k')
plot(t[(t>=8000) & (t<10000)],O.y[1][(t>=8000) & (t<10000)],'r')
plot(t[(t>=8000) & (t<10000)],O.y[0][(t>=8000) & (t<10000)],'b')
xlim(0,10000)
xticks([])
subplot(312)
plot(S.t/ms, S.i, '.')
xlim(0,10000)
ylim(2500,2600) # 100 random neurons
xticks([])
subplot(313)
plot(S2.t/ms, S2.i, '.')
ylim(100,300)
xlim(0,10000)

show()

This is the output:
index2

In the original model:
index3

13

Dear @mstimberg , can you please say what is wrong at this code?

14

Hi @serge25. Sorry, it took me longer than expected to come back to the code.

The problem with your version was that as I mentioned in my comment, the code that sets binding_A and binding_B already contains the random number generation done by the odor function. In your code, you used it alongside the previous odor function for c1 and c2, which meant that the connection pattern and the response to odors weren’t matched to each other anymore. The easiest quick fix would be to set the binding_A and binding_B variables based on the c1 and c2 variables like this:

receptors.binding_A = 2*(np.log(c1)/np.log(10) - bhalf)/(bmax - bmin)
receptors.binding_B = 2*(np.log(c2)/np.log(10) - bhalf)/(bmax - bmin)

I had another look at the full code, though, and I think the most elegant solution would be to get rid of the network_operation and odor function completely, and handle everything via Brian2 mechanisms. If you are interested in that (I’ll probably clean up the code a bit and put it into the Brian2 documentation as an example), here’s the full code using this approach:

Code 
from brian2 import *

bmin,bmax=-7,-1

N=5000 # number of receptors

# Odors
seed(31415) # Get the same neurons every time
intensity=3000.
I1,I2=intensity,intensity

K=1.
n=3.

# Odor plumes (fluctuating concentrations)
tau_plume=75*ms
eq_plumes='''
dx/dt=-x/tau_plume+(2./tau_plume)**.5*xi : 1
y=clip(x,0,inf) : 1
'''
plume=NeuronGroup(2,model=eq_plumes) # 2 odors

# Receptor neurons
Fmax=40*Hz # maximum firing rate
tau=20*ms
Imax=1/(1-exp(-1/(Fmax*tau))) # maximum input current

eq_receptors='''
dv/dt=(Imax*(c**n)/(c**n+K**n)-v)/tau : 1
c = I1*c1*y1 + I2*c2*y2 : 1  # concentrations (relative to activation constant)
y1 : 1 (linked)  # random plume for component 1
y2 : 1 (linked)  # random plume for component 2
c1 : 1 (constant)
c2 : 1 (constant)
binding_A = 2*(log(c1)/log(10) - bhalf)/(bmax - bmin) : 1  # normalized binding constants for odor A
binding_B = 2*(log(c2)/log(10) - bhalf)/(bmax - bmin) : 1  # normalized binding constants for odor B
'''
receptors=NeuronGroup(N,model=eq_receptors,threshold='v > 1',reset='v = 0')
receptors.y1=linked_var(plume, 'y', index=np.zeros(N, dtype=int))
receptors.y2=linked_var(plume, 'y', index=np.ones(N, dtype=int))
receptors.c1 = '10**(rand()*(bmax-bmin)+bmin)'
receptors.c2 = '10**(rand()*(bmax-bmin)+bmin)'

# Decoder neurons
M=200
taud=8*ms
sigma=.15

bhalf=.5*(bmin+bmax)

eq_decoders='''
dv/dt=-v/taud + sigma*(2/taud)**.5*xi : 1
'''
decoders=NeuronGroup(2*M,model=eq_decoders,threshold='v > 1',reset='v = 0')
# First M neurons encode odor A, next M neurons encode odor B

# Synapses
syn = Synapses(receptors,decoders, model='w : 1', on_pre='v +=w')

# Partition the decoder neurons into groups

syn.connect('j < M  and binding_A_pre >= j * 1.0/M         and binding_A_pre < (j + 1) * 1.0/M')
syn.connect('j >= M and binding_B_pre >= (2*M-1-j) * 1.0/M and binding_B_pre < (2*M-j) * 1.0/M')

# Normalize weights by number of incoming synapses
syn.w = '1.0 / N_incoming'

# Record odor concentration and output spikes
O=StateMonitor(plume,'y',record=True)
S=SpikeMonitor(receptors)
S2=SpikeMonitor(decoders)

print("Odor A")
I1,I2=intensity,0
run(2*second)
print(receptors.c[:5])
print("Odor B")
I1,I2=0,intensity
run(2*second)
print("Odor B x2")
I1,I2=0,2*intensity
run(2*second)
print("Odor A + odor C")
I1,I2=intensity,intensity
old_c2 = receptors.c2[:]
receptors.c2 = '10**(rand()*(bmax-bmin)+bmin)'
run(2*second)
print("Odor A + odor B")
I1,I2=intensity,intensity
receptors.c2 = old_c2
run(2*second)

t=O.t/ms

subplot(311) # odor fluctuations
plot(t[t<2000],O.y[0][t<2000],'b')
plot(t[(t>=2000) & (t<4000)],O.y[1][(t>=2000) & (t<4000)],'r')
plot(t[(t>=4000) & (t<6000)],2*O.y[1][(t>=4000) & (t<6000)],'r')
plot(t[(t>=6000) & (t<8000)],O.y[0][(t>=6000) & (t<8000)],'b')
plot(t[(t>=6000) & (t<8000)],O.y[1][(t>=6000) & (t<8000)],'k')
plot(t[(t>=8000) & (t<10000)],O.y[1][(t>=8000) & (t<10000)],'r')
plot(t[(t>=8000) & (t<10000)],O.y[0][(t>=8000) & (t<10000)],'b')
xlim(0,10000)
xticks([])
subplot(312)
plot(S.t/ms, S.i, '.')
xlim(0,10000)
ylim(2500,2600) # 100 random neurons
xticks([])
subplot(313)
plot(S2.t/ms, S2.i, '.')
ylim(100,300)
xlim(0,10000)

show()
15

@mstimberg Thank you.