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
Hi @czx. The exponential Euler integrator assumes that each variable defined by a differential equation only linearly depends on itself. It then solves the equation for each variable using the analytical solution over a time step, assuming that all other variables are constant over the time step. This is in general a good method for HH-type equations, and it should work here as well (see e.g. this paper for details). I have to have a closer look, but I think the problem here has to do with the fact that the solution for n depends on t (via ema). I think one cannot directly apply the exponential Euler algorithm in its usual form here. The solution would be to either use an explicit method that doesn’t have any requirements – I would not use euler here, since it will lead to instabilities quickly, but rather rk2 or rk4. The other option would be to mark the alpha and beta functions as (constant over dt), i.e.
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 (contant over dt)
betan = .5*exp((10.*mV-(v+(em1/eT)**2*EK0)+VT)/(40.*mV))/ms : second**-1 (constant over dt)
This means that they will be calculated once at the beginning of the time step, instead of incorporating their changes over a time step into the solution for n. Please note that this would not change anything in many cases. In particular, an euler integration already considers everything constant over a time step and for “exponential euler”, alpha and beta are usually considered constant over the time step as well, since all variables that they depend on (usually only v) are considered constant. In your case, you have a dependency on t which changes things. But again, your dependence on t is quite “slow”, since it has a time constant of roughly 100ms. Considering it as constant over a 0.1ms time step therefore seems to be very reasonable to me.
Hope that makes sense, let us know whether this fixes your issue!
PS: Note that posts/comments here are markdown-formatted, so for better readability you should enclose Python code in triple backticks like this:
Much gratitude to your reply! Now I know the definition of Exponential Euler integrator, and as suggested, mark the alpha and beta functions as (constant over dt), and the negative n value issue has been solved. However, according to the differential equation of n, there should not be negative values, no matter using numerical or analytical methods. Using (constant over dt) cannot explain this phenomenon. What’s more, except euler, explicit methods such as rk2 and rk4 are not able to calculate the spikes of that neuron, since the leaky reversal potential El here is written as
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
to accommodate the changes in concentrations due to alterations of Nav and Kv. Besides, when substituding em1 with -0.5 ema in the equations, the negative n value issue was also solved, and the results were highly close to that following your suggestion. That makes it seem not like the problem with t-dependence (via ema).
Mathematically/analytically, the variable n should certainly be between 0 and 1 and never negative. But the numerical method and the time step definitely matter. Take this very simple example with constant \alpha and \beta values:
\frac{dn}{dt} = \alpha(1-n)-\beta n
Let’s use a somewhat extreme example, say that \alpha=100/s and \beta=100000/s. Let’s assume our value for n is at 0.01, and we use the Euler method with a time step of 1ms. The update step will then be
Applying this to the current value of n, will lead to a negative value. In general, for explicit methods, the time scale of the variables need to be the same order of magnitude (or slower) than the integration time step. This is not the case above, where the \beta induces a time scale of 0.01ms, much faster than the 1ms of the integration time step. This means that euler, rk2, rk4, etc. will all fail to correctly integrate the equation above (in different ways), while the exact and exponential Euler methods would give the same, correct, solution.
I’m afraid I don’t understand this comment – rk2 and rk4 can “solve” anything that “euler” can, all three are explicit methods, just of different order.
I still need to look more closely what (if anything) goes wrong in the exponential Euler method. The time-dependence seemed to be the most obvious possibility. This does not exclude that changing something else in the equations can make things work correctly.
Thanks for your detailed illustration! I’ll try it later.
As for your example, I think it is too extreme to happen, since I extracted the alpha and beta values from the same case where negative n happened, and found both of them had a value of below 800/ s, like the figure shows:
These alpha and beta values should not make n negative, if exponential Euler is used, assuming the value for n is 0, with a time step of 0.1 ms. But it did happened. That’s why I was confused about this issue .
Hi @czx. I agree that exponential Euler should not make it go below zero in this case. I was just reponding to your statement
However, according to the differential equation of n, there should not be negative values, no matter using numerical or analytical methods.
The numerical method (and the parameters, etc.) do certainly matter – that there should be no negative values “according to the differential equation” is not enough.
