.. currentmodule:: brian2 .. example_6_COBA_with_astro: Example: example_6_COBA_with_astro ================================== .. only:: html .. |launchbinder| image:: http://mybinder.org/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() .. image:: ../resources/examples_images/frompapers.Stimberg_et_al_2018.example_6_COBA_with_astro.1.png