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