Hi @Bolloknoon . There are a few issues with the code that I think together explain the instability. The most important problem is that the synapse model is incorrect, and therefore all your synapses are in effect inhibitory. You are currently using:
Synapse_model_gaba = '''
w: siemens (constant)
I= w * (v_pre - v_post) : amp
'''
with on_pre='I_gaba+=I'. This means that at the time step when the time step is propagated to the post-synaptic cell (i.e. after the delay), you calculate an instantaneous current based on the difference between the pre- and post-synaptic membrane potential. This will (almost) always give a positive value, since the pre-synaptic neuron is above threshold. Instead, you probably rather want to use on_pre='g_ampa += w' in the synapse and replace I_ampa and I_gaba by synaptic conductances in your neuron model:
I_ext = g_ampa*(v-V_AMPA) + g_gaba*(v-V_GABA) : amp
dg_ampa/dt = -g_ampa/ t_ampa : siemens
dg_gaba/dt = -g_gaba/ t_gaba : siemens
There’s also a small mistake concerning the parametrization of the excitatory and inhibitory neurons: both refer to the constants gL and gKZ and both will therefore use the same values (see my related answer in another post: Exploring STDP parameters through one network - #2 by mstimberg). Instead, you should use different names (e.g. gL_e and gL_i, etc.) for the different values. In order to avoid duplicating the equations, you can use the equations with the general gL/gKZ constants that you currently have, and replace them when you create the NeuronGroup, e.g.:
L_5_i = NeuronGroup(150, Equations(eqs, gL='gL_i', gKZ='gKZ_i'), ...)
Finally, you are using the rk4 algorithm, but in a model like this it might become unstable with the default simulation time step of 0.1ms – you can change it by using, e.g. defaultclock.dt = 0.025*ms (if you want to use a for many neurons StateMonitor, consider setting e.g. StateMonitor(..., dt=0.1*ms) to avoid filling up the memory with many almost identical values). The exponential_euler algorithm might be a better choice, since it is generally stable. However, you cannot directly apply it to you equations since you introduce a non-linear dependency of \frac{dv}{dt} on v via the instantaneous value for m. You could add (constant over dt) to the definition of m (which means: calculate the value of m once at the beginning at the time step, and then consider it constant) – this way you can switch to exponential_euler.