Description of problem
I have a model running perfectly fine on brian 2.5 (python 3.8) and I would like it to run on brian & python latests versions. However, the model does not reproduce at all the results from version 2.5. The miniEPSPs (mini excitatory post synaptic potentials) do not work. Those are defined as custom events, appearing following a poisson process. The rate of arrival is defined in my axosomatic compartment, as a function dependent on the last spike (here, pyramidal spike but this is not important). Depending on that rate, the custom_event (defining the time of a poisson process) is triggered and this triggers in return changes in the synapses (releasing small potentials).
Minimal code to reproduce problem
The custom_event (defined when instantiating the axosomatic compartment of our model):
PY_soma = NeuronGroup(N_PY,Soma_eqs,method='rk4',threshold='v>40*mV',refractory=3*ms,events={'custom_poisson_PY':'rand()<mean_rate_PY*dt','custom_poisson_IN':'rand()<mean_rate_IN*dt'})
Equation for the axosomatic compartment of the model where the mean_rate is defined:
Soma_eqs = '''
v = -68*mV*init_timestep + truev*(1-init_timestep) : volt
init_timestep = init_timedarray(t) : 1
truev = (vgap + (R * s_soma) * g2) / (1 + (R * s_soma) * g1) : volt
g1 = g_na *(m_na ** 3) * h_na * Phi_m + g_k * n_k * Tad + g_nap * m_nap : siemens * meter**-2
g2 = g_na *(m_na ** 3) * h_na * Phi_m * E_na + g_k * n_k * Tad * E_k + g_nap * m_nap * E_na + 6.74172*uA*cm**-2 + Iext : amp*meter**-2
I_na = Phi_m * g_na * (m_na ** 3) * h_na * (v - E_na) : amp * meter**-2
dm_na/dt = -(m_na - m_nainf) / tau_m_na : 1
dh_na/dt = -(h_na - h_nainf) / tau_h_na : 1
alpham_na = (0.182/ms * (v + 25*mV)/mV / (1 - exp(-(v + 25*mV)/9/mV))) * int(abs((v-10*mV)/(-35*mV)) > 1e-6) + (0.182/ms * 9) * int(abs((v-10*mV)/(-35*mV)) < 1e-6) : Hz
betam_na = (-0.124/ms * (v + 25*mV)/mV / (1 - exp((v + 25*mV)/9/mV))) * int(abs((-v+10*mV)/(35*mV)) > 1e-6) + (0.124/ms * 9) * int(abs((-v+10*mV)/(35*mV)) < 1e-6) : Hz
alphah_na = (0.024/ms * (v + 40*mV)/mV / (1 - exp(-(v + 40*mV)/5/mV))) * int(abs((v-10*mV)/(-50*mV)) > 1e-6) + (0.024/ms * 5) * int(abs((v-10*mV)/(-50*mV)) < 1e-6) : Hz
betah_na = (-0.0091/ms * (v + 65*mV)/mV / (1 - exp((v + 65*mV)/5/mV))) * int(abs((-v+10*mV)/(75*mV)) > 1e-6) + (0.0091/ms * 5) * int(abs((-v+10*mV)/(75*mV)) < 1e-6) : Hz
h_nainf = 1 / (1 + exp((v + 55*mV)/6.2/mV)) : 1
tau_h_na = (1 / (alphah_na + betah_na)) / Qt : second
m_nainf = alpham_na / (alpham_na + betam_na) : 1
tau_m_na = (1 / (alpham_na + betam_na)) / Qt : second
I_nap = g_nap * m_nap * (v - E_na) : amp * meter**-2
dm_nap/dt = -(m_nap - m_napinf)/tau_m_nap : 1
m_napinf = 0.02 / (1 + exp(-(v + 42*mV)/5/mV)) : 1
I_k = Tad * g_k * n_k * (v - E_k) : amp * meter**-2
dn_k/dt = -(n_k - n_kinf) / tau_n_k : 1
alphan_k = (0.02/mV) * (v - 25*mV) / (1 - exp(-(v - 25*mV)/(9*mV))) : 1
betan_k = (-0.002/mV) * (v - 25*mV) / (1 - exp((v - 25*mV)/(9*mV))) : 1
n_kinf = alphan_k / (alphan_k + betan_k) : 1
tau_n_k = (1*msecond / (alphan_k + betan_k)) / Qt : second
Iext : amp * meter**-2
vgap : volt
g_na : siemens * meter**-2
g_k : siemens * meter**-2
g_nap : siemens * meter**-2
x : meter
t_last_spike_PY : second
t_last_spike_IN : second
mean_rate_PY = (log((t-t_last_spike_PY + 50*ms)/(50*ms))/400*kHz)*int((t-t_last_spike_PY)>70*ms) : Hz
mean_rate_IN = (log((t-t_last_spike_PY + 50*ms)/(50*ms))/400*kHz)*int((t-t_last_spike_PY)>70*ms) : Hz
mean_rate_GABAA = (log((t-t_last_spike_IN + 50*ms)/(50*ms))/400*kHz)*int((t-t_last_spike_IN)>70*ms) : Hz
'''
The synapse from where the issue is coming from:
def syn_ampa(source,target,syntype,surfacetarget,connection_pattern,g_syn,A_mEPSP,poissontype,connection_pattern_mini,synapses_received):
eq_syn_AMPA = '''_post = (D * g_syn * (1/surfacetarget) * (W/mM) * (v_post - E_syn)) / N_incoming : amp * meter**-2 (summed)
dW/dt = alpha * (1*mM - W) * T - beta * W : mM (clock-driven)
T = A_syn * int((t_last_spike + tmax - t) > 0*ms)* int((t - t_last_spike) > 0*ms) : mM
D : 1
q1 : second
g_syn : siemens
E_syn = 0*mV : volt
t_last_spike : second
alpha = 0.94*ms**-1*mM**-1 : second**-1*mM**-1
beta = 0.18*kHz : Hz
U = 0.073 : 1
tau = 700*ms : second
tmax = 0.3*ms : second
A_syn = 0.5*mM : mM
connection_pattern_mini : 1
surfacetarget : meter**2
synapses_received : 1
'''
eq_syn_AMPA2 = '''_post = (D * A_mEPSP * (1/surfacetarget) * (W_Poisson/mM) * (v_post - E_syn)) / N_incoming : amp * meter**-2 (summed)
dW_Poisson/dt = alpha * (1*mM - W_Poisson) * T_Poisson - beta * W_Poisson : mM (clock-driven)
T_Poisson = (A_syn * int((t_last_spike_Poisson_PY + tmax - t) > 0*ms)* int((t - t_last_spike_Poisson_PY) > 0*ms))*int(connection_pattern_mini == 0) + (A_syn * int((t_last_spike_Poisson_IN + tmax - t) > 0*ms)* int((t - t_last_spike_Poisson_IN) > 0*ms))*int(connection_pattern_mini == 1) : mM
t_last_spike_Poisson_PY : second
t_last_spike_Poisson_IN : second
A_mEPSP : siemens
'''
pre_code_ampa = '''q1 = ((t - t_last_spike) - tmax)
D = 1 - (1 - D*(1 - U)) * exp(-q1/tau)*int((-q1/tau)>-10)*int((-q1/tau)<10)
t_last_spike = t
t_last_spike_PY_pre = t
'''
poisson_PY = '''
q1 = ((t - t_last_spike) - tmax)
t_last_spike_Poisson_PY = t
D = (1 - (1 - D*(1 - 0)) * exp(-q1/tau)*int((-q1/tau)>-10)*int((-q1/tau)<10))*int(connection_pattern_mini == 0) + D*int(connection_pattern_mini == 1)
W0 = W
'''
poisson_IN = '''
q1 = ((t - t_last_spike) - tmax)
t_last_spike_Poisson_IN = t
D = (1 - (1 - D*(1 - 0)) * exp(-q1/tau)*int((-q1/tau)>-10)*int((-q1/tau)<10))*int(connection_pattern_mini == 1) + D*int(connection_pattern_mini == 0)
W0 = W
'''
S=Synapses(source,target,model=(syntype+eq_syn_AMPA)+(poissontype+eq_syn_AMPA2),method='rk4',on_pre={'pre': pre_code_ampa,'mini_PY': poisson_PY,'mini_IN': poisson_IN},on_event={'mini_PY':'custom_poisson_PY','mini_IN':'custom_poisson_IN'})
if connection_pattern=='':
S.connect()
else :
S.connect(condition=connection_pattern, skip_if_invalid=True)
S.surfacetarget = surfacetarget
S.g_syn = g_syn
S.A_mEPSP = A_mEPSP
S.connection_pattern_mini = connection_pattern_mini
S.synapses_received = synapses_received
# print("AMPA synapses from "+str(source)+" to "+str(target))
# print(S.N_outgoing_pre)
# print(S.N_incoming_post)
return S
What you have aready tried
I have considered it being a ran() error, but the custom_event actually happen correctly (can be easily tracked using a monitor like M0=EventMonitor(PY_soma,'custom_poisson_PY',record=True)). The issue is coming from the âintegrationâ between the custom_event and the synapse itself: the variables that should be updated when this event happen are simply not (all of them, not just a single one). I have gone through the release notes of both 2.6 and 2.7, and canât find anything that changed that would impact the synapse.
Expected output (if relevant)
Update of variables in âpoisson_PYâ of the synapse, as the event âcustom_poisson_PYâ happen.
Actual output (if relevant)
No update of variables
Thank you in advance for your help.