.. currentmodule:: brian2 .. example_5_astro_ring: Example: example_5_astro_ring ============================= .. 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_5_astro_ring.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 5: Astrocytes connected in a network. Intercellular calcium wave propagation in a ring of 50 astrocytes connected by bidirectional gap junctions (see Goldberg et al., 2010) :: from brian2 import * import plot_utils as pu set_device('cpp_standalone', directory=None) # Use fast "C++ standalone mode" ################################################################################ # Model parameters ################################################################################ ### General parameters duration = 4000*second # Total simulation time sim_dt = 50*ms # Integrator/sampling step ### 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 O_delta = 0.6*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_ex = 0.09*umolar/second # Maximal exogenous IP3 flow 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 ################################################################################ # Model definition ################################################################################ defaultclock.dt = sim_dt # Set the integration time ### Astrocytes astro_eqs = ''' dI/dt = J_delta - J_3K - J_5P + J_ex + J_coupling : mmolar 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 # Exogenous stimulation (rectangular wave with period of 50s and duty factor 0.4) stimulus = int((t % (50*second))<20*second) : 1 delta_I_bias = I - I_bias*stimulus : mmolar J_ex = -F_ex/2*(1 + tanh((abs(delta_I_bias) - I_Theta)/omega_I)) * sign(delta_I_bias) : 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 # External IP_3 drive I_bias : mmolar (constant) ''' N_astro = 50 # Total number of astrocytes in the network astrocytes = NeuronGroup(N_astro, astro_eqs, method='rk4') # Asymmetric stimulation on the 50th cell to get some nice chaotic patterns astrocytes.I_bias[N_astro//2] = 1.0*umolar astrocytes.h = 0.9 # Diffusion between astrocytes astro_to_astro_eqs = ''' delta_I = I_post - I_pre : mmolar J_coupling_post = -F/2 * (1 + tanh((abs(delta_I) - I_Theta)/omega_I)) * sign(delta_I) : mmolar/second (summed) ''' astro_to_astro = Synapses(astrocytes, astrocytes, model=astro_to_astro_eqs) # Couple neighboring astrocytes: two connections per astrocyte pair, as # the above formulation will only update the I_coupling term of one of the # astrocytes astro_to_astro.connect('j == (i + 1) % N_pre or ' 'j == (i + N_pre - 1) % N_pre') ################################################################################ # Monitors ################################################################################ astro_mon = StateMonitor(astrocytes, variables=['C'], record=True) ################################################################################ # Simulation run ################################################################################ run(duration, report='text') ################################################################################ # Analysis and plotting ################################################################################ plt.style.use('figures.mplstyle') fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(6.26894, 6.26894 * 0.66), gridspec_kw={'left': 0.1, 'bottom': 0.12}) scaling = 1.2 step = 10 ax.plot(astro_mon.t/second, (astro_mon.C[0:N_astro//2-1].T/astro_mon.C.max() + np.arange(N_astro//2-1)*scaling), color='black') ax.plot(astro_mon.t/second, (astro_mon.C[N_astro//2:].T/astro_mon.C.max() + np.arange(N_astro//2, N_astro)*scaling), color='black') ax.plot(astro_mon.t/second, (astro_mon.C[N_astro//2-1].T/astro_mon.C.max() + np.arange(N_astro//2-1, N_astro//2)*scaling), color='C0') ax.set(xlim=(0., duration/second), ylim=(0, (N_astro+1.5)*scaling), xticks=np.arange(0., duration/second, 500), xlabel='time (s)', yticks=np.arange(0.5*scaling, (N_astro + 1.5)*scaling, step*scaling), yticklabels=[str(yt) for yt in np.arange(0, N_astro + 1, step)], ylabel='$C/C_{max}$ (cell index)') pu.adjust_spines(ax, ['left', 'bottom']) pu.adjust_ylabels([ax], x_offset=-0.08) plt.show() .. image:: ../resources/examples_images/frompapers.Stimberg_et_al_2018.example_5_astro_ring.1.png