# Conversion olfactory model from brian1 to brian2

**URL:** <https://brian.discourse.group/t/conversion-olfactory-model-from-brian1-to-brian2/910>\
**Category:** Support\
**Created:** [16 February 2023 09:54 UTC](https://brian.discourse.group/t/conversion-olfactory-model-from-brian1-to-brian2/910 "2023-02-16T09:54:13Z")\
**Posts on this page:** 14\
**Page:** 1

<div class="post-metadata">

**Author:** ![serge25](https://avatars.discourse-cdn.com/v4/letter/s/4af34b/32.png) [@serge25](https://brian.discourse.group/u/serge25)\
**Post date:** [16 February 2023 09:54 UTC](https://brian.discourse.group/t/conversion-olfactory-model-from-brian1-to-brian2/910/1 "2023-02-16T09:54:13Z")

</div>

# Description of problem

I’m trying to convert model from the article [Computing with Neural Synchrony](https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1002561#s5). I’dont understand how is used with brian2 connectivity which is used at this model.

# Minimal code to reproduce problem

```auto
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()

```

```auto
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)

---

<div class="post-metadata">

**Author:** ![mstimberg](https://yyz2.discourse-cdn.com/free1/user_avatar/brian.discourse.group/mstimberg/32/11_2.png) [@mstimberg](https://brian.discourse.group/u/mstimberg)\
**Post date:** [16 February 2023 10:10 UTC](https://brian.discourse.group/t/conversion-olfactory-model-from-brian1-to-brian2/910/2 "2023-02-16T10:10:03Z")

</div>

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](https://brian2.readthedocs.io/en/stable/introduction/brian1_to_2/synapses.html). You will also most likely not need any of the `network_operation`s 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](https://brian2.readthedocs.io/en/stable/user/input.html)). 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.

---

<div class="post-metadata">

**Author:** ![serge25](https://avatars.discourse-cdn.com/v4/letter/s/4af34b/32.png) [@serge25](https://brian.discourse.group/u/serge25)\
**Post date:** [16 February 2023 11:02 UTC](https://brian.discourse.group/t/conversion-olfactory-model-from-brian1-to-brian2/910/3 "2023-02-16T11:02:25Z")

</div>

> [@mstimberg](#):
>
> run\_regularly

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

```auto
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)
'''

```

---

<div class="post-metadata">

**Author:** ![mstimberg](https://yyz2.discourse-cdn.com/free1/user_avatar/brian.discourse.group/mstimberg/32/11_2.png) [@mstimberg](https://brian.discourse.group/u/mstimberg)\
**Post date:** [16 February 2023 11:04 UTC](https://brian.discourse.group/t/conversion-olfactory-model-from-brian1-to-brian2/910/4 "2023-02-16T11:04:45Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![serge25](https://avatars.discourse-cdn.com/v4/letter/s/4af34b/32.png) [@serge25](https://brian.discourse.group/u/serge25)\
**Post date:** [16 February 2023 11:21 UTC](https://brian.discourse.group/t/conversion-olfactory-model-from-brian1-to-brian2/910/5 "2023-02-16T11:21:03Z")

</div>

Thank you. But what about this part:

```auto
# 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

```

---

<div class="post-metadata">

**Author:** ![mstimberg](https://yyz2.discourse-cdn.com/free1/user_avatar/brian.discourse.group/mstimberg/32/11_2.png) [@mstimberg](https://brian.discourse.group/u/mstimberg)\
**Post date:** [16 February 2023 13:31 UTC](https://brian.discourse.group/t/conversion-olfactory-model-from-brian1-to-brian2/910/6 "2023-02-16T13:31:24Z")

</div>

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`).

---

<div class="post-metadata">

**Author:** ![serge25](https://avatars.discourse-cdn.com/v4/letter/s/4af34b/32.png) [@serge25](https://brian.discourse.group/u/serge25)\
**Post date:** [17 February 2023 08:47 UTC](https://brian.discourse.group/t/conversion-olfactory-model-from-brian1-to-brian2/910/7 "2023-02-17T08:47:30Z")

</div>

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

```auto
@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:

```auto
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))

```

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

```

---

<div class="post-metadata">

**Author:** ![serge25](https://avatars.discourse-cdn.com/v4/letter/s/4af34b/32.png) [@serge25](https://brian.discourse.group/u/serge25)\
**Post date:** [20 February 2023 17:28 UTC](https://brian.discourse.group/t/conversion-olfactory-model-from-brian1-to-brian2/910/9 "2023-02-20T17:28:52Z")

</div>

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

Brian1:

```auto
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:

```auto
# 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.

---

<div class="post-metadata">

**Author:** ![serge25](https://avatars.discourse-cdn.com/v4/letter/s/4af34b/32.png) [@serge25](https://brian.discourse.group/u/serge25)\
**Post date:** [21 February 2023 08:56 UTC](https://brian.discourse.group/t/conversion-olfactory-model-from-brian1-to-brian2/910/10 "2023-02-21T08:56:14Z")

</div>

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:

```auto
# 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

```

---

<div class="post-metadata">

**Author:** ![mstimberg](https://yyz2.discourse-cdn.com/free1/user_avatar/brian.discourse.group/mstimberg/32/11_2.png) [@mstimberg](https://brian.discourse.group/u/mstimberg)\
**Post date:** [21 February 2023 16:14 UTC](https://brian.discourse.group/t/conversion-olfactory-model-from-brian1-to-brian2/910/11 "2023-02-21T16:14:28Z")

</div>

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.:

```python
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.

---

<div class="post-metadata">

**Author:** ![serge25](https://avatars.discourse-cdn.com/v4/letter/s/4af34b/32.png) [@serge25](https://brian.discourse.group/u/serge25)\
**Post date:** [24 February 2023 18:56 UTC](https://brian.discourse.group/t/conversion-olfactory-model-from-brian1-to-brian2/910/12 "2023-02-24T18:56:22Z")

</div>

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?

```auto
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](https://global.discourse-cdn.com/free1/uploads/brian/original/1X/5c8a2b45fc8e55814e063fdb5218d58a38866b77.png)

In the original model:  
 ![index3](https://global.discourse-cdn.com/free1/uploads/brian/original/1X/f3d5e597f4d139a714bee5285b356f77c2fd1795.png)

---

<div class="post-metadata">

**Author:** ![serge25](https://avatars.discourse-cdn.com/v4/letter/s/4af34b/32.png) [@serge25](https://brian.discourse.group/u/serge25)\
**Post date:** [6 March 2023 15:41 UTC](https://brian.discourse.group/t/conversion-olfactory-model-from-brian1-to-brian2/910/13 "2023-03-06T15:41:06Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![mstimberg](https://yyz2.discourse-cdn.com/free1/user_avatar/brian.discourse.group/mstimberg/32/11_2.png) [@mstimberg](https://brian.discourse.group/u/mstimberg)\
**Post date:** [13 March 2023 13:24 UTC](https://brian.discourse.group/t/conversion-olfactory-model-from-brian1-to-brian2/910/14 "2023-03-13T13:24:53Z")

</div>

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:

```python
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**
>
> ```python
> 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()
> 
> ```

---

<div class="post-metadata">

**Author:** ![serge25](https://avatars.discourse-cdn.com/v4/letter/s/4af34b/32.png) [@serge25](https://brian.discourse.group/u/serge25)\
**Post date:** [14 March 2023 05:45 UTC](https://brian.discourse.group/t/conversion-olfactory-model-from-brian1-to-brian2/910/15 "2023-03-14T05:45:36Z")

</div>

> [@mstimberg](#):
>
> ```auto
> 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
> 
> ```

@mstimberg Thank you.