HH mechanism implemented in Matlab

[brocke]

Dear group,

I've tried to replicate voltage and current traces in Matlab with HH mechanism, but the traces don't look similar. Probably, there is a trick, and it was discussed already. If it's so, excuse me for this, I really searched and didn't find.
I have a simple one isopotentional section and HH mechanism, smth like this:

Code: Select all

class hh_cell(object) :
....
    def geom(self) :
        soma = h.Section()
      #  all values are in microns
        soma.L= 30
        soma.diam = 30 
        soma.nseg = 1 # -> isopotential soma
        self.soma = soma
        self.soma_area = h.area(0.5)*1e-8 #in cm2
    def biophys(self) :
        self.soma.insert('hh')
        self.soma.gnabar_hh = (5.4*1e-7)/self.soma_area #0.0191 default .12 S/cm2
        self.soma.gkbar_hh = (5.4*1e-8)/self.soma_area #0.00191 default .036 S/cm2
        self.soma.gl_hh = 1.0/(8.333*1e8*self.soma_area) #0.0000424 default .0003 S/cm2
....
cell = hh_cell()
h.finitialize(-70)
neuron.run(200)

In Matlab the main code looks like:

Code: Select all

initials =  [-0.07; 0.028905534475191907; 0.7540796658225246; 0.24458654944007166];

options  = odeset('Reltol',3e-5,'Abstol',1e-20,'MaxOrder', 1 ); 

% --------- Simulation
[t y]    = ode15s(@eq_HH,T_span,initials,options);

...
function out = eq_HH(t,input) 
...
CM = 1.0*1e-6*soma_A;    %F                 
gm1 = 1.0/(8.333*1e8); %S

ENa = 50.0*1e-3; %V
gNa  = 5.4e-7;%S      
EK  = -77.0*1e-3; %V
gK   = 5.4e-8;    %S   
Em  = -54.3*1e-3; %V

V1 = input(1);  
m = input(2);
h = input(3);
n = input(4);

alpha(1) = 100000*vtrap((-0.040-V1),0.01);%0.1*vtrap(-(V1+40.0),10.0);
beta(1)  = 4000/exp((0.065+V1)/0.0180);%4.0*exp(-(V1+65.0)/18.0);
%------ h
alpha(2) = 70/exp((V1+0.065)*50);%0.07*exp(-(V1+65.0)/20.0);
beta(2)  = 1000/(exp(1000*(-V1-0.035))+1);%1.0 /(exp(-(V1+35.0)/10.0)+1);
%------ n
alpha(3) = 10000*vtrap(-0.055-V1,0.01);%0.01*vtrap(-(55.0+V1 ),10.0);
beta(3)  = 125/exp((0.065+V1)/0.080);%0.125*exp(-(65.0+V1 )/80.0); 

for i=1:
    %probabilities of the chanel gates being opened
    dP_dt(i) = alpha(i)*(1-P(i))-beta(i)*P(i);
end

dV1_dt = -((V1-ENa)*gNa*h*m^3+(V1-EK)*gK*n^4+(V1-Em)*gm1)/CM;
out=[dV1_dt; dP_dt'];

I used ode15s ans SI units to solve the system in matlab. Initial values were taken from NEURON after finitialize.

By "don't look similar" I meant, that there are more spikes in matlab and values are higher.
I will be thankful for any help/suggestions.

[brocke]

Interestingly,
when I do this trick in Matlab:

Code: Select all

V1 = V1*1e3;
%------ m
alpha(1) = 1e3*0.1*vtrap((-V1-40.0),10.0);
beta(1)  = 1e3*4.0*exp((-V1-65.0)/18.0);
%------ h
alpha(2) = 1e3*0.07*exp((-V1-65.0)/20.0);
beta(2)  = 1e3*1.0 /(exp((-V1-35.0)/10.0)+1);
%------ n
alpha(3) = 1e3*0.01*vtrap((-55.0-V1 ),10.0);
beta(3)  = 1e3*0.125*exp((-65.0-V1 )/80.0); 
V1 = V1*1e-3;

then traces are exactly the same as in NEURON.
Could it be possible that the roundoff error in Matlab lead to such a different behaviour?

[ted]

The first question to ask is:
are you sure that what's in the computer is the same as what you think is in the computer? Specifically, are the geometries and biophysical properties identical to what you think they should be?

The second question, which is related to the first, is:
do your Python statements for NEURON produce a computational model whose properties are identical to the model specified by your Matlab statements?

The third question is whether the numerical integration methods have comparable accuracies so that the should be expected to generate results that are quantitatively similar to each other. However, this can be addressed only after the first two questions have been answered "yes".

[brocke]

The first question to ask is:
are you sure that what's in the computer is the same as what you think is in the computer? Specifically, are the geometries and biophysical properties identical to what you think they should be?

I believe yes.

The second question, which is related to the first, is:
do your Python statements for NEURON produce a computational model whose properties are identical to the model specified by your Matlab statements?

I believe yes.

The third question is whether the numerical integration methods have comparable accuracies so that the should be expected to generate results that are quantitatively similar to each other. However, this can be addressed only after the first two questions have been answered "yes".

I believe again yes, since I tried also cvode with the same values or order, atol, and rtol.
The only difference I can see is that NEURON, in case of HH mechanism, operates with values of the order of milli, and Matlab (for my personal simplicity) I used SI units.

[ted]

The null hypothesis is that every new observation is artifact. Have you proven that what you believe is true?

[brocke]

Finally success. After comparing for the 5th time all the digits, I've found an extra 0 in exp term in Matlab code:

Code: Select all

beta(2)  = 1000/(exp(1000*(-V1-0.035))+1);%1.0 /(exp(-(V1+35.0)/10.0)+1);

Sorry for this post. I don't think it can be of an help for anyone. Shall I remove it?

thank you Ted.

[ted]

It's good to see that you found and fixed the problem.

I'd rather this thread remained intact for several reasons. For one thing, it illustrates how easy it is to make a small mistake, and how difficult it can be to find such an error. Also it is a reminder that such errors are far from rare.

More particularly, even though this thread doesn't contain a lot of code, the code that is present might be useful to some other Matlab user who is getting started with NEURON, or vice versa.

And it illustrates the idea of validating simulation results by replication with different source code, something that doesn't get near the attention or effort it deserves.