Problem on Implementing wang buzsaki model

Support
,
1

Description of problem

I’m going to implement Gutnisky et al (2017) ‘Mechanisms underlying a thalamocortical transformation during active tactile sensation’ paper which uses populational network based on modified Wang Buzsaki model. But unfortunately, I dosen’t work as I expected.

Main problem is variables on HH type neurons get diverges during simulation. I tried some variable changes since I think divergent was originated from original ODEs but It dosen’t get the problem worse.

if you have any tips on a better way to trobleshoot, it would be really appreciated!

Thanks

Minimal code to reproduce problem

from brian2 import *
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable
''' Inhibitory neuron '''

Cm = 1 * uF  # /cm**2
gL = 0.1 * msiemens
EL = -65 * mV
ENa = 55 * mV
EK = -90 * mV
gNa = 100 * msiemens
gK = 40 * msiemens
gKZ = 0 * msiemens

t_ampa = 2 * ms
t_gaba = 3 * ms

'''synapse decay time constant'''

t_ampa = 2 * ms
t_gaba = 3 * ms

eqs = '''
       dv/dt = (-gNa*m**3*h*(v-ENa) -gK*n**4*(v-EK) -gL*(v-EL) - gKZ*z*(v-EK)- I_ext) / Cm : volt 

       I_ext = I_ampa + I_gaba : amp

       dI_ampa/dt = -I_ampa/ t_ampa : amp
       dI_gaba/dt = -I_gaba/ t_gaba : amp


       m = alpham/(alpham+betam) : 1
       alpham = -0.1/mV*(v+35*mV)/(exp(-0.1/mV*(v+35*mV))-1)/ms : Hz
       betam = 4*exp(-(v+60*mV)/(18*mV))/ms : Hz

       dh/dt = 5*(alphah*(1-h)-betah*h) : 1
       alphah = 0.07*exp(-(v+58*mV)/(20*mV))/ms : Hz
       betah = 1./(exp(-0.1/mV*(v+28*mV))+1)/ms : Hz

       dn/dt = 5*(alphan*(1-n)-betan*n) : 1
       alphan = -0.01/mV*(v+34*mV)/(exp(-0.1/mV*(v+34*mV))-1)/ms : Hz
       betan = 0.125*exp(-(v+44*mV)/(80*mV))/ms : Hz

       dz/dt = (z_inf - z)/(60*ms) : 1
       z_inf = 1/(1+exp(-0.7*((v/(1 * mV) + 30)))) : 1
       '''


L_5_i = NeuronGroup(150, eqs, threshold='v>-48.8*mV', refractory='v>-48.8*mV', method='rk4')
L_5_i.v = EL
L_5_i.h = 1


''' excitatory neuron '''
gL = 0.05 * msiemens
gKZ = 0.5 * msiemens  # /cm**2

L_5_e = NeuronGroup(1600, eqs, threshold='v>-39.7*mV', refractory='v>-39.7*mV', method='rk4')
L_5_e.v = EL
L_5_e.h = 1


''' L4 connection '''

g_et = 0.15 * msiemens
g_it = 0.2 * msiemens
g_ee = 0.2 * msiemens
g_ie = 0.6 * msiemens
g_ei = 0.7 * msiemens
g_ii = 0.55 * msiemens

K_ET = 50
K_IT = 75
K_EE = 200
K_IE = 400
K_EI = 25
K_II = 25

# reversal potential
V_AMPA = 0 * mV
V_GABA = - 85 * mV

tau_syn_all = 1

Synapse_model_ampa = '''
           w: siemens (constant)
           I = w * (v_pre - v_post) : amp
       '''
Synapse_model_gaba = '''
           w: siemens (constant)
           I= w * (v_pre - v_post) : amp
       '''

C_i_1 = Synapses(L_5_e, L_5_e, Synapse_model_ampa, on_pre='I_ampa+=I', delay=0.001 * second)

C_i_1.connect(p=0.125)
C_i_1.w = tau_syn_all * g_ee / (sqrt(K_EE))

C_i_2 = Synapses(L_5_e, L_5_i, Synapse_model_ampa, on_pre='I_ampa+=I', delay=0.001 * second)

C_i_2.connect(p=0.25)
C_i_2.w = tau_syn_all * g_ie / (sqrt(K_IE))

C_i_3 = Synapses(L_5_i, L_5_e, Synapse_model_gaba, on_pre='I_gaba+=I', delay=0.0005 * second)

C_i_3.connect(p=0.1666)
C_i_3.w = tau_syn_all * g_ei / (sqrt(K_EI))

C_i_4 = Synapses(L_5_i, L_5_i, Synapse_model_gaba, on_pre='I_gaba+=I', delay=0.0005 * second)

C_i_4.connect(p=0.1666)
C_i_4.w = tau_syn_all * g_ii / (sqrt(K_II))

''' external input '''

''''VPM input'''
Synapse_model_thalamic = '''
           w: siemens (constant)
           I= w * (V_AMPA - v_post) : amp
       '''

P_1 = PoissonGroup(200, rates=14 * Hz)

S_1 = Synapses(P_1, L_5_e, Synapse_model_thalamic, on_pre='I_ampa+=I', delay=0.001 * second)

S_1.connect(p=0.25)
S_1.w = tau_syn_all * g_et / (sqrt(K_ET))

S_2 = Synapses(P_1, L_5_i, Synapse_model_thalamic, on_pre='I_ampa+=I', delay=0.001 * second)

S_2.connect(p=0.5)
S_2.w = tau_syn_all * g_it / (sqrt(K_IT))

print('variable completed')

'''run simulation'''

A = SpikeMonitor(L_5_e)
B = PopulationRateMonitor(L_5_e, name='L5_Excitatory_neurons')

C = SpikeMonitor(L_5_i)
D = PopulationRateMonitor(L_5_i, name='L5_inhibitatory_neurons')

trace = StateMonitor(L_5_e, 'v', record=[0, 1, 2])
trace_1 = StateMonitor(L_5_i, 'v', record=[0, 1, 2])

run(3 * second)

''' spike statistics '''

print('1 simulation complete')

plot(trace.t / ms, trace[0].v)
plot(trace.t / ms, trace_1[0].v)

xlabel('Time (ms)')
ylabel('v')
show()

Actual output (if relevant)

WARNING    neurongroup's variable 'z' has NaN, very large values, or encountered an error in numerical integration. This is usually a sign that an unstable or invalid integration method was chosen. [brian2.groups.group.invalid_values]
WARNING    neurongroup's variable 'n' has NaN, very large values, or encountered an error in numerical integration. This is usually a sign that an unstable or invalid integration method was chosen. [brian2.groups.group.invalid_values]
WARNING    neurongroup's variable 'h' has NaN, very large values, or encountered an error in numerical integration. This is usually a sign that an unstable or invalid integration method was chosen. [brian2.groups.group.invalid_values]
WARNING    neurongroup's variable 'v' has NaN, very large values, or encountered an error in numerical integration. This is usually a sign that an unstable or invalid integration method was chosen. [brian2.groups.group.invalid_values]
WARNING    neurongroup's variable 'I_ampa' has NaN, very large values, or encountered an error in numerical integration. This is usually a sign that an unstable or invalid integration method was chosen. [brian2.groups.group.invalid_values]
WARNING    neurongroup's variable 'I_gaba' has NaN, very large values, or encountered an error in numerical integration. This is usually a sign that an unstable or invalid integration method was chosen. [brian2.groups.group.invalid_values]
WARNING    neurongroup_1's variable 'z' has NaN, very large values, or encountered an error in numerical integration. This is usually a sign that an unstable or invalid integration method was chosen. [brian2.groups.group.invalid_values]
WARNING    neurongroup_1's variable 'n' has NaN, very large values, or encountered an error in numerical integration. This is usually a sign that an unstable or invalid integration method was chosen. [brian2.groups.group.invalid_values]
WARNING    neurongroup_1's variable 'h' has NaN, very large values, or encountered an error in numerical integration. This is usually a sign that an unstable or invalid integration method was chosen. [brian2.groups.group.invalid_values]
WARNING    neurongroup_1's variable 'v' has NaN, very large values, or encountered an error in numerical integration. This is usually a sign that an unstable or invalid integration method was chosen. [brian2.groups.group.invalid_values]
WARNING    neurongroup_1's variable 'I_ampa' has NaN, very large values, or encountered an error in numerical integration. This is usually a sign that an unstable or invalid integration method was chosen. [brian2.groups.group.invalid_values]
WARNING    neurongroup_1's variable 'I_gaba' has NaN, very large values, or encountered an error in numerical integration. This is usually a sign that an unstable or invalid integration method was chosen. [brian2.groups.group.invalid_values]
1 simulation complete

Figure_1

2

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.

3

I didn’t know that constants in equations get same values in this code.

Its quite simple though but if you didn’t point it out, I will spend some extra weeks to deal with it.

Thank you for giving me a full solution in this implementation.

1 Like