# Description of problem

As described in the paper DOI 10.1007/s10237-013-0543-7, I combined the mechanical strain parameter with ionic reversal potential ENa, EK and El to see the electro-mechanical coupling effect of an HH type neuron. However, the potassium conductance variable ‘n’ calculated with the method ‘exponential_euler’ was negative in the prior stage, which is not acceptable according to its definition. When ‘exponential_euler’ was altered to ‘euler’, ‘n’ became possitive during the whole stage. Thus I think it’s something wrong with ‘exponential_euler’, but I can’t find a clear format of it.

A simplified code to reproduce the problem is as follows:

# Minimal code to reproduce problem

```
from brian2 import *
prefs.codegen.target = "numpy"
import matplotlib.pyplot as plt
start_scope()
runtime = 100*ms
# electrical parameters
area = 20000*umetre**2
Cm0 = 1*ufarad*cm**-2 * area
gl = 5e-5*siemens*cm**-2 * area
El0 = -65*mV
EK0 = -90*mV
ENa0 = 50*mV
g_na = 100*msiemens*cm**-2 * area
g_kd = 30*msiemens*cm**-2 * area
VT = -63*mV
I = 10*nA #input current
# introduce mechanical properties
eT = 0.2 # strain threshold
E = 800000*pascal # modulus of neuron， calibrated under et=0.1
sigma = 1000*kilogram*meter**-1*second**-2 # assume the stress is 1 kPa
tau = 95.74*ms # the time constant of the tensile axonal active molecular reorganization, assume to be 0.09574s
eqs_HH = '''
dv/dt = (gl*(El-v) - g_na*(m*m*m)*h*(v-ENa) - g_kd*(n*n*n*n)*(v-EK) + I)/Cm : volt
dm/dt = 0.32*(mV**-1)*(13.*mV-(v+(em1/eT)**2*ENa0)+VT)/
(exp((13.*mV-(v+(em1/eT)**2*ENa0)+VT)/(4.*mV))-1.)/ms*(1-m)-0.28*(mV**-1)*((v+(em1/eT)**2*ENa0)-VT-40.*mV)/
(exp(((v+(em1/eT)**2*ENa0)-VT-40.*mV)/(5.*mV))-1.)/ms*m : 1
dn/dt = alphan*(1.-n)-betan*n : 1
dh/dt = 0.128*exp((17.*mV-(v+(em1/eT)**2*ENa0)+VT)/(18.*mV))/ms*(1.-h)-4./(1+exp((40.*mV-(v+(em1/eT)**2*ENa0)+VT)/(5.*mV)))/ms*h : 1
Cm = Cm0*sqrt(1+ema) : farad
alphan = .032*(mV**-1)*(15.*mV-(v+(em1/eT)**2*EK0)+VT)/
(exp((15.*mV-(v+(em1/eT)**2*EK0)+VT)/(5.*mV))-1.)/ms : second**-1
betan = .5*exp((10.*mV-(v+(em1/eT)**2*EK0)+VT)/(40.*mV))/ms : second**-1
ema = 0.5/E*(1-exp(-t/tau))*sigma : 1
em = sqrt(1+ema)-1 : 1
em1 = em : 1
EK = EK0*(1-(em1/eT)**2) : volt
ENa = ENa0*(1-(em1/eT)**2) : volt
El =(1+g_na*(m*m*m)*h/gl+g_kd*(n*n*n*n)/gl)*VT-(g_na*(m*m*m)*ENa+g_kd*(n*n*n*n)*EK)/gl: volt
'''
group = NeuronGroup(1, eqs_HH,
threshold='v > -40*mV',
refractory='v > -40*mV',
method='exponential_euler')
group.v = -65*mV
statemon = StateMonitor(group, ('v', 'em1', 'm', 'h', 'n'), record=True, dt=0.1*ms)
run(runtime)
plt.figure(figsize=(18, 4)) # show conductance variables
subplot(131)
plot(statemon.t/ms, statemon.m[0], color='gray')
xlabel('t (ms)')
ylabel('m');
subplot(132)
plot(statemon.t/ms, statemon.h[0], color='gray')
xlabel('t (ms)')
ylabel('h');
subplot(133)
plot(statemon.t/ms, statemon.n[0], color='gray')
xlabel('t (ms)')
ylabel('n');
show()
```

# What you have aready tried

Altered the integration method to make ‘n’ posiitive

# Expected output (if relevant)

‘n’ should be possitive at anytime

# Actual output (if relevant)

‘n’ is negative in the prior stage