.. currentmodule:: brian2

.. Jansen_Rit_1995_single_column:

Example: Jansen_Rit_1995_single_column
======================================


        .. only:: html

            .. |launchbinder| image:: http://mybinder.org/badge.svg
            .. _launchbinder: https://mybinder.org/v2/gh/brian-team/brian2-binder/master?filepath=examples/frompapers/Jansen_Rit_1995_single_column.ipynb

            .. note::
               You can launch an interactive, editable version of this
               example without installing any local files
               using the Binder service (although note that at some times this
               may be slow or fail to open): |launchbinder|_

        

[Jansen and Rit 1995 model](https://link.springer.com/content/pdf/10.1007/BF00199471.pdf) (Figure 3) in Brian2.

Equations are the system of differential equations number (6) in the original paper. 
The rate parameters $a=100 s^{-1}$ and $b=200 s^{-1}$ were changed to excitatory $\tau_e = 1000ms/a =10ms$  and inhibitory $\tau_i = 1000ms/b =20ms$ time constants as in 
[Thomas Knosche review](https://link.springer.com/referenceworkentry/10.1007%2F978-1-4614-6675-8_65), 
[Touboul et al. 2011](https://direct.mit.edu/neco/article-abstract/23/12/3232/7717/Neural-Mass-Activity-Bifurcations-and-Epilepsy?redirectedFrom=fulltext), or 
[David & Friston 2003](https://www.sciencedirect.com/science/article/pii/S1053811903004579). 
Units were removerd from parameters $e_0$, $v_0$, $r_0$, $A$, $B$, and $p$ to stop Brian's confusion.

Ruben Tikidji-Hamburyan 2021 (rth@r-a-r.org)

::

    
    
    
    from numpy import *
    from numpy import random as rnd
    from matplotlib.pyplot import *
    from brian2 import *
    
    defaultclock.dt = .1*ms   #default time step
    
    te,ti    = 10.*ms, 20.*ms #taus for excitatory and inhibitory populations
    e0       = 5.             #max firing rate
    v0       = 6.             #(max FR)/2 input
    r0       = 0.56           #gain rate
    A,B,C    = 3.25, 22., 135 #standard parameters as in the set (7) of the original paper
    P,deltaP = 120, 320.-120  #random input uniformly distributed between 120 and
                              #320 pulses per second
    
    # Random noise
    nstim = TimedArray(rnd.rand(70000),2*ms)
    
    # Equations as in the system (6) of the original paper
    equs = """
    dy0/dt = y3                             /second : 1
    dy3/dt = (A        * Sp -2*y3 -y0/te*second)/te : 1
    dy1/dt = y4                             /second : 1
    dy4/dt = (A*(p+ C2 * Se)-2*y4 -y1/te*second)/te : 1
    dy2/dt = y5                             /second : 1
    dy5/dt = (B   * C4 * Si -2*y5 -y2/ti*second)/ti : 1
    p  = P0+nstim(t)*dlP                   : 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
    dlP           : 1
    """
    
    n = NeuronGroup(6,equs,method='euler') #creates 6 JR models for different connectivity parameters
    
    #set parameters as for different traces on figure 3 of the original paper
    n.C1[0]  = 68
    n.C1[1]  = 128
    n.C1[2]  = C
    n.C1[3]  = 270
    n.C1[4]  = 675
    n.C1[5]  = 1350
    #set stimulus offset and noise magnitude
    n.P0     = P
    n.dlP    = deltaP
    
    #just record everything
    sm = StateMonitor(n,['y4','y1','y3','y0','y5','y2'],record=True)
    
    #Runs for 5 second
    run(5*second,report='text')
    
    
    #This code goes over all models with different parameters and plot activity of each population.
    
    figure(1,figsize=(22,16))
    idx1 = where(sm.t/second>2.)[0]
    
    o = 0
    for p in [0,1,2,3,4,5]:
        if o == 0: ax = subplot(6,3,1)
        else     :subplot(6,3,1+o,sharex=ax)
        if o == 0: title("E")
        plot(sm.t[idx1]/second, sm[p].y1[idx1],'g-')
        ylabel(f"C={n[p].C1[0]}")
        if o == 15: xlabel("Time (seconds)")
        subplot(6,3,2+o,sharex=ax)
        if o == 0: title("P")
        plot(sm.t[idx1]/second, sm[p].y0[idx1],'b-')
        if o == 15: xlabel("Time (seconds)")
        subplot(6,3,3+o,sharex=ax)
        if o == 0: title("I")
        plot(sm.t[idx1]/second, sm[p].y2[idx1],'r-')
        if o == 15: xlabel("Time (seconds)")
        o += 3
    
    show()
    
    

.. image:: ../resources/examples_images/frompapers.Jansen_Rit_1995_single_column.1.png