.. currentmodule:: brian2

.. example_6_COBA_with_astro:

Example: example_6_COBA_with_astro
==================================


        .. only:: html

            .. |launchbinder| image:: file:///usr/share/doc/python-brian-doc/docs/badge.svg
            .. _launchbinder: https://mybinder.org/v2/gh/brian-team/brian2-binder/master?filepath=examples/frompapers/Stimberg_et_al_2018/example_6_COBA_with_astro.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|_

        

Modeling neuron-glia interactions with the Brian 2 simulator
Marcel Stimberg, Dan F. M. Goodman, Romain Brette, Maurizio De Pittà
bioRxiv 198366; doi: https://doi.org/10.1101/198366

Figure 6: Recurrent neuron-glial network.

Randomly connected COBA network (see Brunel, 2000) with excitatory synapses modulated
by release-increasing gliotransmission from a randomly connected network of astrocytes.

::

    
    from brian2 import *
    
    import plot_utils as pu
    
    set_device('cpp_standalone', directory=None)  # Use fast "C++ standalone mode"
    seed(28371)  # to get identical figures for repeated runs
    
    ################################################################################
    # Model parameters
    ################################################################################
    ### General parameters
    N_e = 3200                   # Number of excitatory neurons
    N_i = 800                    # Number of inhibitory neurons
    N_a = 3200                   # Number of astrocytes
    
    ## Some metrics parameters needed to establish proper connections
    size = 3.75*mmeter           # Length and width of the square lattice
    distance = 50*umeter         # Distance between neurons
    
    ### Neuron parameters
    E_l = -60*mV                 # Leak reversal potential
    g_l = 9.99*nS                # Leak conductance
    E_e = 0*mV                   # Excitatory synaptic reversal potential
    E_i = -80*mV                 # Inhibitory synaptic reversal potential
    C_m = 198*pF                 # Membrane capacitance
    tau_e = 5*ms                 # Excitatory synaptic time constant
    tau_i = 10*ms                # Inhibitory synaptic time constant
    tau_r = 5*ms                 # Refractory period
    I_ex = 100*pA                # External current
    V_th = -50*mV                # Firing threshold
    V_r = E_l                    # Reset potential
    
    ### Synapse parameters
    rho_c = 0.005                # Synaptic vesicle-to-extracellular space volume ratio
    Y_T = 500.*mmolar            # Total vesicular neurotransmitter concentration
    Omega_c = 40/second          # Neurotransmitter clearance rate
    U_0__star = 0.6              # Resting synaptic release probability
    Omega_f = 3.33/second        # Synaptic facilitation rate
    Omega_d = 2.0/second         # Synaptic depression rate
    w_e = 0.05*nS                # Excitatory synaptic conductance
    w_i = 1.0*nS                 # Inhibitory synaptic conductance
    # --- Presynaptic receptors
    O_G = 1.5/umolar/second      # Agonist binding (activating) rate
    Omega_G = 0.5/(60*second)    # Agonist release (deactivating) rate
    
    ### Astrocyte parameters
    # ---  Calcium fluxes
    O_P = 0.9*umolar/second      # Maximal Ca^2+ uptake rate by SERCAs
    K_P = 0.05*umolar            # Ca2+ affinity of SERCAs
    C_T = 2*umolar               # Total cell free Ca^2+ content
    rho_A = 0.18                 # ER-to-cytoplasm volume ratio
    Omega_C = 6/second           # Maximal rate of Ca^2+ release by IP_3Rs
    Omega_L = 0.1/second         # Maximal rate of Ca^2+ leak from the ER
    # --- IP_3R kinectics
    d_1 = 0.13*umolar            # IP_3 binding affinity
    d_2 = 1.05*umolar            # Ca^2+ inactivation dissociation constant
    O_2 = 0.2/umolar/second      # IP_3R binding rate for Ca^2+ inhibition
    d_3 = 0.9434*umolar          # IP_3 dissociation constant
    d_5 = 0.08*umolar            # Ca^2+ activation dissociation constant
    # --- IP_3 production
    # --- Agonist-dependent IP_3 production
    O_beta = 0.5*umolar/second   # Maximal rate of IP_3 production by PLCbeta
    O_N = 0.3/umolar/second      # Agonist binding rate
    Omega_N = 0.5/second         # Maximal inactivation rate
    K_KC = 0.5*umolar            # Ca^2+ affinity of PKC
    zeta = 10                    # Maximal reduction of receptor affinity by PKC
    # --- Endogenous IP3 production
    O_delta = 1.2*umolar/second  # Maximal rate of IP_3 production by PLCdelta
    kappa_delta = 1.5*umolar     # Inhibition constant of PLC_delta by IP_3
    K_delta = 0.1*umolar         # Ca^2+ affinity of PLCdelta
    # --- IP_3 degradation
    Omega_5P = 0.05/second       # Maximal rate of IP_3 degradation by IP-5P
    K_D = 0.7*umolar             # Ca^2+ affinity of IP3-3K
    K_3K = 1.0*umolar            # IP_3 affinity of IP_3-3K
    O_3K = 4.5*umolar/second     # Maximal rate of IP_3 degradation by IP_3-3K
    # --- IP_3 diffusion
    F = 0.09*umolar/second       # GJC IP_3 permeability
    I_Theta = 0.3*umolar         # Threshold gradient for IP_3 diffusion
    omega_I = 0.05*umolar        # Scaling factor of diffusion
    # --- Gliotransmitter release and time course
    C_Theta = 0.5*umolar         # Ca^2+ threshold for exocytosis
    Omega_A = 0.6/second         # Gliotransmitter recycling rate
    U_A = 0.6                    # Gliotransmitter release probability
    G_T = 200*mmolar             # Total vesicular gliotransmitter concentration
    rho_e = 6.5e-4               # Astrocytic vesicle-to-extracellular volume ratio
    Omega_e = 60/second          # Gliotransmitter clearance rate
    alpha = 0.0                  # Gliotransmission nature
    
    ################################################################################
    # Define HF stimulus
    ################################################################################
    stimulus = TimedArray([1.0, 1.2, 1.0, 1.0], dt=2*second)
    
    ################################################################################
    # Simulation time (based on the stimulus)
    ################################################################################
    duration = 8*second          # Total simulation time
    
    ################################################################################
    # Model definition
    ################################################################################
    ### Neurons
    neuron_eqs = '''
    dv/dt = (g_l*(E_l-v) + g_e*(E_e-v) + g_i*(E_i-v) + I_ex*stimulus(t))/C_m : volt (unless refractory)
    dg_e/dt = -g_e/tau_e : siemens  # post-synaptic excitatory conductance
    dg_i/dt = -g_i/tau_i : siemens  # post-synaptic inhibitory conductance
    # Neuron position in space
    x : meter (constant)
    y : meter (constant)
    '''
    neurons = NeuronGroup(N_e + N_i, model=neuron_eqs,
                          threshold='v>V_th', reset='v=V_r',
                          refractory='tau_r', method='euler')
    exc_neurons = neurons[:N_e]
    inh_neurons = neurons[N_e:]
    # Arrange excitatory neurons in a grid
    N_rows = int(sqrt(N_e))
    N_cols = N_e//N_rows
    grid_dist = (size / N_cols)
    exc_neurons.x = '(i // N_rows)*grid_dist - N_rows/2.0*grid_dist'
    exc_neurons.y = '(i % N_rows)*grid_dist - N_cols/2.0*grid_dist'
    # Random initial membrane potential values and conductances
    neurons.v = 'E_l + rand()*(V_th-E_l)'
    neurons.g_e = 'rand()*w_e'
    neurons.g_i = 'rand()*w_i'
    
    ### Synapses
    synapses_eqs = '''
    # Neurotransmitter
    dY_S/dt = -Omega_c * Y_S                                    : mmolar (clock-driven)
    # Fraction of activated presynaptic receptors
    dGamma_S/dt = O_G * G_A * (1 - Gamma_S) - Omega_G * Gamma_S : 1 (clock-driven)
    # Usage of releasable neurotransmitter per single action potential:
    du_S/dt = -Omega_f * u_S                                    : 1 (event-driven)
    # Fraction of synaptic neurotransmitter resources available for release:
    dx_S/dt = Omega_d *(1 - x_S)                                : 1 (event-driven)
    U_0                                                         : 1
    # released synaptic neurotransmitter resources:
    r_S                                                         : 1
    # gliotransmitter concentration in the extracellular space:
    G_A                                                         : mmolar
    # which astrocyte covers this synapse ?
    astrocyte_index : integer (constant)
    '''
    synapses_action = '''
    U_0 = (1 - Gamma_S) * U_0__star + alpha * Gamma_S
    u_S += U_0 * (1 - u_S)
    r_S = u_S * x_S
    x_S -= r_S
    Y_S += rho_c * Y_T * r_S
    '''
    exc_syn = Synapses(exc_neurons, neurons, model=synapses_eqs,
                       on_pre=synapses_action+'g_e_post += w_e*r_S',
                       method='exact')
    exc_syn.connect(True, p=0.05)
    exc_syn.x_S = 1.0
    inh_syn = Synapses(inh_neurons, neurons, model=synapses_eqs,
                       on_pre=synapses_action+'g_i_post += w_i*r_S',
                       method='exact')
    inh_syn.connect(True, p=0.2)
    inh_syn.x_S = 1.0
    # Connect excitatory synapses to an astrocyte depending on the position of the
    # post-synaptic neuron
    N_rows_a = int(sqrt(N_a))
    N_cols_a = N_a/N_rows_a
    grid_dist = size / N_rows_a
    exc_syn.astrocyte_index = ('int(x_post/grid_dist) + '
                               'N_cols_a*int(y_post/grid_dist)')
    ### Astrocytes
    # The astrocyte emits gliotransmitter when its Ca^2+ concentration crosses
    # a threshold
    astro_eqs = '''
    # Fraction of activated astrocyte receptors:
    dGamma_A/dt = O_N * Y_S * (1 - clip(Gamma_A,0,1)) -
                  Omega_N*(1 + zeta * C/(C + K_KC)) * clip(Gamma_A,0,1) : 1
    # Intracellular IP_3
    dI/dt = J_beta + J_delta - J_3K - J_5P + J_coupling              : mmolar
    J_beta = O_beta * Gamma_A                                        : mmolar/second
    J_delta = O_delta/(1 + I/kappa_delta) * C**2/(C**2 + K_delta**2) : mmolar/second
    J_3K = O_3K * C**4/(C**4 + K_D**4) * I/(I + K_3K)                : mmolar/second
    J_5P = Omega_5P*I                                                : mmolar/second
    # Diffusion between astrocytes:
    J_coupling                                                       : mmolar/second
    
    # Ca^2+-induced Ca^2+ release:
    dC/dt = J_r + J_l - J_p                                   : mmolar
    dh/dt = (h_inf - h)/tau_h                                 : 1
    J_r = (Omega_C * m_inf**3 * h**3) * (C_T - (1 + rho_A)*C) : mmolar/second
    J_l = Omega_L * (C_T - (1 + rho_A)*C)                     : mmolar/second
    J_p = O_P * C**2/(C**2 + K_P**2)                          : mmolar/second
    m_inf = I/(I + d_1) * C/(C + d_5)                         : 1
    h_inf = Q_2/(Q_2 + C)                                     : 1
    tau_h = 1/(O_2 * (Q_2 + C))                               : second
    Q_2 = d_2 * (I + d_1)/(I + d_3)                           : mmolar
    
    # Fraction of gliotransmitter resources available for release:
    dx_A/dt = Omega_A * (1 - x_A) : 1
    # gliotransmitter concentration in the extracellular space:
    dG_A/dt = -Omega_e*G_A        : mmolar
    # Neurotransmitter concentration in the extracellular space:
    Y_S                           : mmolar
    # The astrocyte position in space
    x : meter (constant)
    y : meter (constant)
    '''
    glio_release = '''
    G_A += rho_e * G_T * U_A * x_A
    x_A -= U_A *  x_A
    '''
    astrocytes = NeuronGroup(N_a, astro_eqs,
                             # The following formulation makes sure that a "spike" is
                             # only triggered at the first threshold crossing
                             threshold='C>C_Theta',
                             refractory='C>C_Theta',
                             # The gliotransmitter release happens when the threshold
                             # is crossed, in Brian terms it can therefore be
                             # considered a "reset"
                             reset=glio_release,
                             method='rk4',
                             dt=1e-2*second)
    # Arrange astrocytes in a grid
    astrocytes.x = '(i // N_rows_a)*grid_dist - N_rows_a/2.0*grid_dist'
    astrocytes.y = '(i % N_rows_a)*grid_dist - N_cols_a/2.0*grid_dist'
    # Add random initialization
    astrocytes.C = 0.01*umolar
    astrocytes.h = 0.9
    astrocytes.I = 0.01*umolar
    astrocytes.x_A = 1.0
    
    ecs_astro_to_syn = Synapses(astrocytes, exc_syn,
                                'G_A_post = G_A_pre : mmolar (summed)')
    ecs_astro_to_syn.connect('i == astrocyte_index_post')
    ecs_syn_to_astro = Synapses(exc_syn, astrocytes,
                                'Y_S_post = Y_S_pre/N_incoming : mmolar (summed)')
    ecs_syn_to_astro.connect('astrocyte_index_pre == j')
    # Diffusion between astrocytes
    astro_to_astro_eqs = '''
    delta_I = I_post - I_pre            : mmolar
    J_coupling_post = -(1 + tanh((abs(delta_I) - I_Theta)/omega_I))*
                      sign(delta_I)*F/2 : mmolar/second (summed)
    '''
    astro_to_astro = Synapses(astrocytes, astrocytes,
                              model=astro_to_astro_eqs)
    # Connect to all astrocytes less than 75um away
    # (about 4 connections per astrocyte)
    astro_to_astro.connect('i != j and '
                           'sqrt((x_pre-x_post)**2 +'
                           '     (y_pre-y_post)**2) < 75*um')
    
    ################################################################################
    # Monitors
    ################################################################################
    # Note that we could use a single monitor for all neurons instead, but this
    # way plotting is a bit easier in the end
    exc_mon = SpikeMonitor(exc_neurons)
    inh_mon = SpikeMonitor(inh_neurons)
    ast_mon = SpikeMonitor(astrocytes)
    
    ################################################################################
    # Simulation run
    ################################################################################
    run(duration, report='text')
    
    ################################################################################
    # Plot of Spiking activity
    ################################################################################
    plt.style.use('figures.mplstyle')
    
    fig, ax = plt.subplots(nrows=3, ncols=1, sharex=True, figsize=(6.26894, 6.26894*0.8),
                           gridspec_kw={'height_ratios': [1, 6, 2],
                                        'left': 0.12, 'top': 0.97})
    time_range = np.linspace(0, duration/second, int(duration/second*100))*second
    ax[0].plot(time_range, I_ex*stimulus(time_range)/pA, 'k')
    ax[0].set(xlim=(0, duration/second), ylim=(98, 122),
              yticks=[100, 120], ylabel='$I_{ex}$ (pA)')
    pu.adjust_spines(ax[0], ['left'])
    
    ## We only plot a fraction of the spikes
    fraction = 4
    ax[1].plot(exc_mon.t[exc_mon.i <= N_e//fraction]/second,
               exc_mon.i[exc_mon.i <= N_e//fraction], '|', color='C0')
    ax[1].plot(inh_mon.t[inh_mon.i <= N_i//fraction]/second,
               inh_mon.i[inh_mon.i <= N_i//fraction]+N_e//fraction, '|', color='C1')
    ax[1].plot(ast_mon.t[ast_mon.i <= N_a//fraction]/second,
               ast_mon.i[ast_mon.i <= N_a//fraction]+(N_e+N_i)//fraction,
               '|', color='C2')
    ax[1].set(xlim=(0, duration/second), ylim=[0, (N_e+N_i+N_a)//fraction],
              yticks=np.arange(0, (N_e+N_i+N_a)//fraction+1, 250),
              ylabel='cell index')
    pu.adjust_spines(ax[1], ['left'])
    
    # Generate frequencies
    bin_size = 1*ms
    spk_count, bin_edges = np.histogram(np.r_[exc_mon.t/second, inh_mon.t/second],
                                        int(duration/bin_size))
    rate = 1.0*spk_count/(N_e + N_i)/bin_size/Hz
    rate[rate<0.001] = 0.001 # Fix 0 lower bound for log scale
    ax[2].semilogy(bin_edges[:-1], rate, '-', color='k')
    pu.adjust_spines(ax[2], ['left', 'bottom'])
    ax[2].set(xlim=(0, duration/second), ylim=(0.1, 150),
              xticks=np.arange(0, 9), yticks=[0.1, 1, 10, 100],
              xlabel='time (s)', ylabel='rate (Hz)')
    ax[2].get_yaxis().set_major_formatter(ScalarFormatter())
    
    pu.adjust_ylabels(ax, x_offset=-0.11)
    
    plt.show()
    

