# 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
```