I plan to build this astrocyte-neuron network and keep running into some errors. There is probably something that I am overlooking. Also sidebar, I would like to meet virtually with someone to discuss brian2, my projects, and theirs and possibly create new relationships.
Hi @amanforthepeople , it’s unfortunately impossible to help you with an overly general question like that (“keep running into some errors”).
The code in the image seems to be from here: comp-glia-book/example_6_COBA_with_astro.py at master · mdepitta/comp-glia-book · GitHub (BTW, there’s a cleaner version of the same code in the Brian documentation: Example: example_6_COBA_with_astro — Brian 2 2.5.1 documentation), but that code should run correctly. The error message complains about a missing definition of astrocyte_index in the synaptic equations, but the above linked code contains that definition in synapses_eqs. So I guess you are running a different version of the code? Please be more specific and/or post some code (preferably as text surrounded by triple backtick quotes – see e.g. Basic writing and formatting syntax - GitHub Docs – and not as an image). Thanks!
Yes, I am using the concept of the “example 6” code to work on a network. Okay I will upload some code for better context, my apologies
Hi Marcel,
You are correct, I utilized the Brain documentation code. Here is the code that I was referring to; the equation is the same. I also maintained variables similar to the example because I wanted to verify if it works.
astro_eqs = '''
# ELLIPSIS BEGIN
# 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
# ELLIPSIS END
# 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, threshold='C>C_Theta', refractory='C>C_Theta',
# The gliotransmitter release happens when the threshold
# is crossed, it can be considered a "reset"
reset=glio_release,
method='rk4',
dt=1e-2*second)
# Add random initialization
astrocytes.C = 0.01*umolar #calcium released
astrocytes.h = 0.9 #gating variable
astrocytes.I = 0.01*umolar #IP3
astrocytes.x_A = 1.0 #fraction of astrocyte receptors activated
########################################################
#More Connections
#connections from synapses to astrocytes, as pathways to trigger astrocyte activation;
#connections from astrocytes to synapses as routes for gliotransmission and thereby modulation of synaptic release;
#connections between astrocytes by GJCs.
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 70um away
# (about 4 connections per astrocyte)
# INCLUDE BEGIN
astro_to_astro.connect('i != j and '
'sqrt((x_pre-x_post)**2 +'
' (y_pre-y_post)**2) < 70*um')
ast_mon = SpikeMonitor(astrocytes)
duration = 1*second
#Stimulation
run(duration, report ='text')
Hi @amanforthepeople I’m afraid it is still impossible for me to help you further. Your posted code is incomplete, in particular it misses the definition of the exc_syn synapses. In the original example, these synapses define the astrocyte_index variable in their equations, to assign each synapse to an astrocyte. As I said earlier, your error message indicates that this is not the case in your model.
We describe the approach in detail in chapter 18.2.8 of doi.org/10.1007/978-3-030-00817-8_18 (or equivalently in section 2.8 of the preprint: https://hal.sorbonne-universite.fr/hal-02325277/file/chapter18_brian.pdf).
from brian2 import *
from brian2tools import *
import matplotlib.pyplot as plt
import plot_utils as pu
start_scope()
#Model Parameters
##################################################################
#General Parameters
N_e = 20 #Number of excitatory Neurons
N_i = 5 #Number of inhibitory Neurons
N_a = 20 #Number of astrocytes
Int_dt = 0.1*ms #Integrator/sampling step
#time interval for which simulation will progress during next "step"
sigma = 0.5*mV
#Physical space metrics
#size = 3.75*mmeter #Length and width of square lattice
#distance = 50*umeter #Distance between Neurons
#Neuron specific Parameters
E_l= -70*mV #Leak Reverse potential
g_l= 5.50*nS #Leak conductance
E_e= 0*mV #Excitatory reverse potential
E_i= -80*mV #inhibitory reverse potential
I_ex= 150*pA*randn() #external current
C_m= 200*pF #membrane capacitance
tau_e= 5*ms #excitatory time constant
tau_i= 10*ms #inhibitory time constant
V_th= -55*mV #Firing threshold
V_r= E_l #Reset potential
tau_r= 5*ms #refractory potential
tau = tau_e+tau_i
### Synapse parameters
rho_c = 0.010 # 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.3/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
alpha = 0.0 # Gliotransmission nature
###############################################################
#Neuron definition
#Leaky integrate and fire model
neuron_eqs= '''
dv/dt = ((g_l*(E_l-v) + g_e*(E_e-v) + g_i*(E_i-v) + I_ex)/C_m)+ sigma*xi*tau**-0.5 : volt (unless refractory)
dg_e/dt = -g_e/tau_e : siemens #exc_post-synaptic conductance
dg_i/dt = -g_i/tau_i : siemens #inh_post-synaptic conductance
'''
#Neurongroup holding both exc and inh cells
neurons = NeuronGroup(N_e + N_i, model=neuron_eqs,
threshold = 'v>V_th', reset = 'v=V_r',
refractory = 'tau_r', method = 'euler')
#Initializing conductance and voltage
neurons.v = 'E_l + rand()*(V_th-E_l)' # Initialize the membrane potential randomly
neurons.g_e = '2*rand()*w_e' #Initialize random excitatory conductance
neurons.g_i = 'rand()*w_i' #Initialize random inhibitory conductance
#Which neurons are exc and which are inh
exc_neurons = neurons[:N_e]
inh_neurons = neurons[N_e:]
####Synapses: a product of both x_s and u_s
#x_S represents fraction available for release
#u_S represents transmitters docked for release
synapses_eqs = '''
#Releaseable neurotransmitter per single action potential(Via Ca2+ stimulation):
du_S/dt = -Omega_f *u_S : 1 (event-driven)
#fraction of synaptic neurotransmitter resources available:
dx_S/dt = Omega_d * (1 - x_S) : 1 (event-driven)
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)
U_0 : 1
# released synaptic neurotransmitter resources:
r_S : 1
# gliotransmitter concentration in the extracellular space:
G_A : mmolar
'''
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
'''
#Y_S :the concentration of synaptically-released neurotransmitter in the extracellular space
###Connection between neurons
###Establishing of source and target neurons and synapse eqs model
exc_syn = Synapses(exc_neurons, neurons, model=synapses_eqs,
on_pre=synapses_action+'g_e_post += w_e*r_S')
inh_syn = Synapses(inh_neurons, exc_neurons, model=synapses_eqs,
on_pre=synapses_action+'g_i_post += w_i*r_S')
#on_pre keyword argument denotes that the series of statements should be executed
#on arrival of a presynaptic action potential.
#In reference to g_e/g_i, excitatory (inhibitory) synapses will increase the excitatory (inhibitory) conductance
#in the postsynaptic cell whenever a presynaptic action potential arrives.
#Condition of connection
exc_syn.connect (condition='i!=j', p=0.05)
inh_syn.connect (condition='i!=j', p=0.20)
#p = probability of connection for possible neurons
###Start from "resting" condition: all synapses should have fully replenished
# neurotransmiter resources. (Reuptake?)
exc_syn.x_S = 1
inh_syn.x_S = 1
##Monitior spike of specific neurons
exc_mon = SpikeMonitor(exc_neurons)
inh_mon = SpikeMonitor(inh_neurons)
ni = 19
#Record data from a single excitatory neuron
neuron_statemon = StateMonitor(exc_neurons, ['v', 'g_e', 'g_i'], record=ni)
### Monitor synaptic variables after synapse are updated, record neuron ni
synapse_statemon = StateMonitor(exc_syn, ['u_S', 'x_S'],
record = exc_syn[ni,:], when = 'after_synapses')
defaultclock.dt = Int_dt
duration = 1*second
#Stimulation
run(duration, report ='text')
print (neuron_eqs)
print (exc_neurons)
print (inh_neurons)
print (exc_mon)
print (inh_mon)
print (neuron_statemon)
print (neurons.g_e)
print ("I_ex= " + str(I_ex))
Hi @amanforthepeople, no worries about the delay, and thanks for the code. But now I am confused: the code you posted seems to run successfully, and does not raise the error that you initially reported.
Hi @mstimberg, the code I just sent works independently, but when ran post the first set of code I sent; it falls apart
I’m sorry, but I still don’t understand: what do you mean by “when ran post the first set of code I sent”? Does the “first set of code” refer to the screenshot or the code posted in KeyError: 'The identifier "astrocyte_index_post" could not be resolved.' - #4 by amanforthepeople? Both were incomplete and cannot be run like that, so what is the actual code you are running?
Okay, so let me clarify. The first set of code is the one I just sent. with the model parameters where I initialize the neurons and their connections. The second is the one I previously sent, beginning with “astro_eqs.”
Thanks for the clarification. Can you say more about what you want to achieve? You want to first run the model for 1s without astrocytes, and then run them for 1s with astrocytes present? Regardless of the error, this will not work like this, because you cannot add new components (with some exceptions, e.g. monitors) to a network after it has been run. Instead, you can create all the objects before the first run, and then use something like
astrocytes.active = False
ecs_astro_to_syn.active = False
ecs_syn_to_astro.active = False
astro_to_astro.active = False
to deactivate them, and after the run use astrocytes.active = True, etc. to activate them for the second run.
Regarding the error about astrocyte_index_post, my previous comment in KeyError: 'The identifier "astrocyte_index_post" could not be resolved.' - #2 by mstimberg still applies: the model you adapted wants to connect synapses to astrocytes by specifying an astrocyte_index in the synaptic equations that can be used to assign specific synapses to specific astrocytes. Your code does the same in the connect statements by referring to astrocyte_index_post and astrocyte_index_pre, but your synaptic model does not define astrocyte_index, so it cannot work.