Units balance in rate-based Jansen and Rit model

[rth]

I’m playing with Jansen and Rit model. I’m trying to reproduce Fig.3 from the original paper, Jansen and Rit 1995.

attempt #1

In the paper, all parameters have units, and units (seems) aren’t balanced, at least I couldn’t figure out how to convince Brian that equations are OK.

defaultclock.dt = .1*ms
a, b  = 100./second,    50./second
e0    = 5./second
v0    = 6.*mV
r0    = 0.56/mV
A,B,C = 3.25*mV, 22.*mV, 135

nstim = TimedArray(rnd.randn(70000),.1*ms)
equs = """
dy0/dt = y3                            : volt/second
dy3/dt = (A        * Sp -2*y3 -y0*a)*a : volt
dy1/dt = y4                            : volt/second
dy4/dt = (A*(p+ C2 * Se)-2*y4 -y1*a)*a : volt
dy2/dt = y5                            : volt/second
dy5/dt = (B   * C4 * Si -2*y5 -y2*b)*b : volt

p  = P0+nstim(t)*(300-P0) : 1

Sp = e0/(1+exp(r0*(v0 - (y1-y2)    ))) : 1
Se = e0/(1+exp(r0*(v0 -  C1*y0     ))) : 1
Si = e0/(1+exp(r0*(v0 -  C3*y0     ))) : 1
C1            : 1
C2 = 0.8 *C1  : 1
C3 = 0.25*C1  : 1
C4 = 0.25*C1  : 1
P0            : 1
"""

n1 = NeuronGroup(3,equs,method='euler')

n1.C1[0] = 135
n1.C1[1] = 270
n1.C1[2] = 675
n1.P0    = 120

sm1 = StateMonitor(n1,['y4','y1','y3','y0','y5','y2'],record=True)


run(7*second)

gv='epi'
figure(1,figsize=(16,16))
idx1 = where(sm1.t/second>2.)[0]
ax = None
for o,p in enumerate([0,1,2]):
    if o == 0: ax = subplot(311)
    else     :subplot(311+o,sharex=ax)
    if 'e' in gv :plot(sm1.t[idx1]/second, sm1[p].y1[idx1],'g-')
    if 'p' in gv :plot(sm1.t[idx1]/second, sm1[p].y0[idx1],'b-')
    if 'i' in gv :plot(sm1.t[idx1]/second, sm1[p].y2[idx1],'r-')
show()

This code returns an error:

brian2.units.fundamentalunits.DimensionMismatchError: Expression "v0 - C1 * y0" uses inconsistent units ("v0" has unit V; "C1 * y0" has unit V/s)

attempt #2

Well, we can assume that units are OK and set everything in unit (1)

defaultclock.dt = .1*ms
a, b  = 100.,    50.
e0    = 5.
v0    = 6.
r0    = 0.56
A,B,C = 3.25, 22., 135

nstim = TimedArray(rnd.randn(70000),.1*ms)
equs = """
dy0/dt = y3                           /second : 1
dy3/dt = (A        * Sp -2*y3 -y0*a)*a/second : 1
dy1/dt = y4                           /second : 1
dy4/dt = (A*(p+ C2 * Se)-2*y4 -y1*a)*a/second : 1
dy2/dt = y5                           /second : 1
dy5/dt = (B   * C4 * Si -2*y5 -y2*b)*b/second : 1

p  = P0+nstim(t)*(300-P0) : 1

Sp = e0/(1+exp(r0*(v0 - (y1-y2)    ))) : 1
Se = e0/(1+exp(r0*(v0 -  C1*y0     ))) : 1
Si = e0/(1+exp(r0*(v0 -  C3*y0     ))) : 1
C1            : 1
C2 = 0.8 *C1  : 1
C3 = 0.25*C1  : 1
C4 = 0.25*C1  : 1
P0            : 1
"""

n1 = NeuronGroup(3,equs,method='euler')

n1.C1[0] = 135
n1.C1[1] = 270
n1.C1[2] = 675
n1.P0    = 120

sm1 = StateMonitor(n1,['y4','y1','y3','y0','y5','y2'],record=True)


run(7*second)

gv='epi'
figure(1,figsize=(16,16))
idx1 = where(sm1.t/second>2.)[0]
ax = None
for o,p in enumerate([0,1,2]):
    if o == 0: ax = subplot(311)
    else     :subplot(311+o,sharex=ax)
    if 'e' in gv :plot(sm1.t[idx1]/second, sm1[p].y1[idx1],'g-')
    if 'p' in gv :plot(sm1.t[idx1]/second, sm1[p].y0[idx1],'b-')
    if 'i' in gv :plot(sm1.t[idx1]/second, sm1[p].y2[idx1],'r-')
show()

producing something like that

pretty different from the original paper, even when looking at P-population.

attempt #3

We can try to use a bit different formulation from Thomas Knosche review, Touboul et al. 2011, or David & Friston 2003 and ask Brian balance units.

defaultclock.dt = .1*ms
a, b  = 100./second,    50./second
te,ti = 1./a, 1./b
e0    = 5
v0    = 6
r0    = 0.56
e1    = e0
r1    = r0

A,B,C = 3.25, 22., 135

nstim = TimedArray(rnd.randn(70000)*250,.1*ms)

equs_v1 = """
dy0/dt = y3                           /ms : 1
dy3/dt = (A        * Sp -2*y3 -y0*a*ms)*a : 1
dy1/dt = y4                           /ms : 1
dy4/dt = (A*(p+ C2 * Se)-2*y4 -y1*a*ms)*a : 1
dy2/dt = y5                           /ms : 1
dy5/dt = (B   * C4 * Si -2*y5 -y2*b*ms)*b : 1
p  = P0+nstim(t) *(300-P0) : 1
Sp = e0/(1+exp(r0*(v0 - (y1-y2)    ))) : 1
Se = e0/(1+exp(r0*(v0 -  C1*y0     ))) : 1
Si = e0/(1+exp(r0*(v0 -  C3*y0     ))) : 1
C1            : 1
C2 = 0.8 *C1  : 1
C3 = 0.25*C1  : 1
C4 = 0.25*C1  : 1
P0            : 1
"""

equs_v2 = """
dy0/dt = y3                             /ms : 1
dy3/dt = (A        * Sp -2*y3 -y0/te*ms)/te : 1
dy1/dt = y4                             /ms : 1
dy4/dt = (A*(p+ C2 * Se)-2*y4 -y1/te*ms)/te : 1
dy2/dt = y5                             /ms : 1
dy5/dt = (B   * C4 * Si -2*y5 -y2/ti*ms)/ti : 1
p  = P0+nstim(t)*(300-P0): 1
Sp = e1/(1+exp(r1*(v0 - (y1-y2)    ))) : 1
Se = e1/(1+exp(r1*(v0 -  C1*y0     ))) : 1
Si = e1/(1+exp(r1*(v0 -  C3*y0     ))) : 1
C1            : 1
C2 = 0.8 *C1  : 1
C3 = 0.25*C1  : 1
C4 = 0.25*C1  : 1
P0            : 1
"""

n1 = NeuronGroup(3,equs_v1,method='euler')
n2 = NeuronGroup(3,equs_v2,method='euler')

n1.C1[0] = n2.C1[0] = 128
n1.C1[1] = n2.C1[1] = 270
n1.C1[2] = n2.C1[2] = 675
n1.P0    = n2.P0    = 120

sm1 = StateMonitor(n1,['y4','y1','y3','y0','y5','y2'],record=True)
sm2 = StateMonitor(n2,['y4','y1','y3','y0','y5','y2'],record=True)

run(7*second,report='text')
gv='epi'
figure(1,figsize=(16,16))
idx1 = where(sm1.t/second>2.)[0]
idx2 = where(sm2.t/second>2.)[0]
o = 0
for p in [0,1,2]:
    for sm,idx in [(sm1,idx1), (sm2,idx2)]:
        if o == 0: ax = subplot(321)
        else     :subplot(321+o,sharex=ax)
        if 'e' in gv :plot(sm.t[idx]/second, sm[p].y1[idx],'g-')
        if 'p' in gv :plot(sm.t[idx]/second, sm[p].y0[idx],'b-')
        if 'i' in gv :plot(sm.t[idx]/second, sm[p].y2[idx],'r-')
        o += 1
show()

To my surprise, this version generates completely different dynamics,

and for just Pyramidal cells,

which doesn’t make any sense!

Any ideas/suggestions?
I can promise to commit a working model whenever @mstimberg will tell

[rth]

After some playing with the JR model, I got it working. The graphs are not reproduced in Fig3 precisely, but closed enough to stop having fun with the code. I used \tau_e and \tau_i with ms units as in Thomas Knosche review, Touboul et al. 2011, or David & Friston 2003. Units were remover from parameters e_0, v_0, r_0, A, B, and p to stop Brian confusion.

The public gist with Jupiter notebook is available here: JR1995_SingleColumn_Figure3.ipynb

The dynamics for each population (E, P, and I) in networks with 6 different connection scales C (see labels on the left axes) are shown below:

[mstimberg]

This is looking great! Unfortunately I still did not find time to look into this in more detail, but I’m happy you figured things out by yourself It would be great to share this in a bit more discoverable way than here as a reply. As I just mentioned in the thread on Using the "Siegart formula" for predicting rate in LIF networks - #3 by adam , it could be become a blog post, or we could have a dedicated place to share examples like this (in addition to or replacing the example section in the docs), … Would be happy to hear thoughts about this question!

[rth]

a bit more discoverable way than here as a reply

@mstimberg could you please elaborate a bit more “discoverable way”? I added more comments that explain each parameter/formula in the gist. Should I assume that a reader knows the original J&R model and have access to the paper so that I can refer to equations/parameters there? Is the code clean enough to go to Brian documentation?

[mstimberg]

@rth

@mstimberg could you please elaborate a bit more “discoverable way”?

Oh I just meant that someone interested in that model would maybe not come across the link because it is a bit “hidden” a few replies deep in the support section. I did not mean that the example needs more work.

And yes, looks perfect to go into the documentation, please open a pull request (but it might take a while before the merge due to the summer break)!