HH model with Passive dendrite??

[jcmocte]

Hi ,

I have just started with NEURON, because i want to use it in order to compare with the results of my own code written in MATLAB and some results from a paper. The thing is i made some test and i got the same results in MATLAB and in the paper, but i can not get the same result in NEURON.

So i just want to modelate a HH section connect with Passive dendrite, so i have 2 sections
When i simulate soma voltage and dendrite voltage, the dendrite voltage appears attenuated with respect to soma voltage
BUT when i simulate with NEURON, i've got the same soma and dendrite voltage, the voltage in dendrite is not attenuated, why is that??

Thanks for any help you can provide

Here is the result in MATLAB:



Here is the result in NEURON:




the code i made in NEURON is very simple:

Code: Select all

load_file("nrngui.hoc")		//Load the simulator

create somaHH,denpas
access somaHH
somaHH {
  nseg = 1
  diam = 43.3	// (um)
  L = 43.3		// (um)
  Ra = 100.0	// (ohm*cm)
  cm = 1		// (uF/cm2) , capacitance
  
  insert hh
  gnabar_hh = 0.120 		// (S/cm2)
  gkbar_hh =  0.036
  gl_hh = .0003
  ena =	45.0				// (mV)
  ek =  -82.0
  el_hh = -59.0
}

denpas {
    nseg = 1
    diam = 2.5
    L = 11
    Ra = 100
	cm = 3		// Dendrite capacitance
    insert pas
    g_pas = .00001
    e_pas = -70.0
}

objectvar stim
somaHH stim = new IClamp(0.5)
    stim.del = 30	// (mseg)
    stim.dur = 30	// (mseg)
    stim.amp = 0.68  	// (nA)

// +++++++++++++++++++++
// CONNECT COMPARTMENTS
// +++++++++++++++++++++
connect denpas(0), somaHH(1)

And here is the code in MATLAB:

Code: Select all

close all;
clear;

%% STEP 1: Parameters values
% +++++++++++++++++++++++++++++++
MODE = 0;   % 0 = HH
Iamp = 0.68e-3;   % uA/cm2 (amplitude) or uA (multicompartment case)
Vrest = -70; % mV
GNA = 120;   % mS/cm2
GK = 36;     % mS/cm2
GLEAK = 0.3;  % mS/cm2
ENA = 45;     % mV
EK = -82;     % mV
ELEAK = -59;   % mV
Csoma = 1;        % uF/cm2
Cden = 3;        % uF/cm2

ti_I = 30;      %time initial current applied (mseg)
tf_I = 60;      %time final current applied (mseg)
Tinit = 0;      %mseg
Tfinal= 100;    %mseg
hstep = 0.1;   %mseg

NUM_COMP = 2;     % number of compartments
PASS_DEN = 1;     % 1=passive dendrite, 0=active dendrite
Eden = -70;        % Vrest of dendrite
GCA = 10;
GK_AHP = 0.8;
GK_C = 15;
ECA = 70;

a_soma = 43.3e-4;     % Radius of soma (cm)
a_den = 2*2.5e-4;       % Radius of cilinder ('a' in thesis) (cm)
dx = 11e-4;           % delta x, length of cilinder (cm)
Ra = 100;           % citoplasm resistivity (KOhm*cm)
Rm = 100;           % Passive resistance dendrite (KOhms/cm2)
Aden = 2*pi*a_den*dx;  %Area dendrite cilinder
Asoma = pi*a_soma^2;

t=Tinit:hstep:Tfinal;
V=zeros(length(t),NUM_COMP);
m=zeros(length(t),NUM_COMP);
h=zeros(length(t),NUM_COMP);
n=zeros(length(t),NUM_COMP);
C=zeros(1,NUM_COMP);
Iinj=zeros(size(t))+(t>=ti_I & t<=tf_I)*Iamp;



%% STEP 2 and 3: Initialize V. And initiliaze m,h & n to steady state
% +++++++++++++++++++++++++++++++
V(1,1) = Vrest;						% initial V for soma
V(1,2:end) = Eden;                    % initial V for dendrite
m(1,:) = alpha_m(Vrest,MODE)/(alpha_m(Vrest,MODE)+beta_m(Vrest,MODE));	% initial m (steady-equation)
h(1,:) = alpha_h(Vrest,MODE)/(alpha_h(Vrest,MODE)+beta_h(Vrest,MODE));	% initial h (steady-equation)
n(1,:) = alpha_n(Vrest,MODE)/(alpha_n(Vrest,MODE)+beta_n(Vrest,MODE));	% initial n (steady-equation)
C(1) = Csoma;					% Capacitance for soma
C(2:end) = Cden;                % Capacitance for dendrite


%% STEP 4: LOOP
% +++++++++++++++++++++++++++++++
for i=1:length(t)-1
    for l=1:NUM_COMP
    
        gNa = GNA * (m(i,l)^3) * h(i,l);  
        gK  = GK  * (n(i,l)^4);           
        
        if(l==1)
            Ome = (gNa*ENA + gK*EK + GLEAK*ELEAK + (Iinj(i)/Asoma)); %HH soma
            Phi = (gNa + gK + GLEAK);
            Tao = a_soma/(Ra*a_soma*a_soma);
            A(i,l) = (Tao)*V(i,l+1) + Ome;
            B(i,l) = (Tao)+Phi;
        elseif(l==NUM_COMP)
            Ome = Eden/(Rm);    %Passive dendrite
            Phi = 1/(Rm);
            Tao = a_den/(Ra*dx*dx);
            A(i,l) = (Tao)*V(i,l-1) + Ome;
            B(i,l) = (Tao)+Phi;
        end        
    
    %Exponential Euler's rule
        m(i+1,l) = first_kinetics(alpha_m(V(i,l),MODE),beta_m(V(i,l),MODE),m(i,l),hstep);
        h(i+1,l) = first_kinetics(alpha_h(V(i,l),MODE),beta_h(V(i,l),MODE),h(i,l),hstep);
        n(i+1,l) = first_kinetics(alpha_n(V(i,l),MODE),beta_n(V(i,l),MODE),n(i,l),hstep);
        V(i+1,l) = (A(i,l)/B(i,l)) + exp(-B(i,l)*hstep/C(l))*(V(i,l)-(A(i,l)/B(i,l)));        
                
    end
end

colorgraph='brgyk';
for l=1:NUM_COMP
    plot(t,V(:,l),'linewidth',1.7,'color',colorgraph(l));   
    cad=sprintf('Compartment %d',l);
    title(cad);
    hold on;
    grid on;
end

The complete code in MATLAB can be downloaded from here

http://ccc.inaoep.mx/~jcmoctezuma/docs/filesMATLAB.zip

[ted]

When simulation results differ from the modeler's expectations, the explanation is usually that what is in the computer is different from what is in the modeler's head--a "mismatch between the computational model and the conceptual model." That's what's going on here.

The NEURON model is working properly--it is doing what it should do, based on the properties assigned by the hoc statements. The diameters, lengths, and membrane properties of the soma and dentrite are such that these sections are electrically very short, so it would be very surprising to see a significant difference of membrane potential between them. If your MATLAB model does something different, the reason is probably that its properties differ from those of the NEURON model.

So what are the properties of the NEURON model?

Code: Select all

oc>forall psection()
somaHH { nseg=1  L=43.3  Ra=100
        /*location 0 attached to cell 0*/
        /* First segment only */
        insert morphology { diam=43.3}
        insert capacitance { cm=1}
        insert hh { gnabar_hh=0.12 gkbar_hh=0.036 gl_hh=0.0003 el_hh=-59}
        insert na_ion { ena=45}
        insert k_ion { ek=-82}
        insert IClamp { del=30 dur=30 amp=0.68}
}
denpas { nseg=1  L=11  Ra=100
        somaHH connect denpas (0), 1
        /* First segment only */
        insert morphology { diam=2.5}
        insert capacitance { cm=3}
        insert pas { g_pas=1e-05 e_pas=-70}
}

Note that all anatomical measurements are in um, Ra is in ohm cm, cm is in uf/cm2, and ionic conductance densities are in S/cm2.

And what are the properties of the MATLAB model?

Code: Select all

a_soma = 43.3e-4;     % Radius of soma (cm)
a_den = 2*2.5e-4;       % Radius of cilinder ('a' in thesis) (cm)
 . . .
Ra = 100;           % citoplasm resistivity (KOhm*cm)

so it looks like the MATLAB model has larger diameters and much larger cytoplasmic resistivity. What other differences exist? Does the MATLAB model use the same ion channel densities as the NEURON model? Do the rate constants have the same values and voltage dependence as they do in the NEURON model? Are the ionic equilibrium potentials the same?

[jcmocte]

Hi Ted,

Thanks for your answer, yes you was right, apparently i was given wrong the Ra parameter , because in fact it should be 100,000 Ohms/cm2
But still i have different results, now the output of dendrite in NEURON is too attenuated compared with soma output, and yeah i totally trust in NEURON more than my model :D ,
I made my MATLAB model because i want to understand the mathematics behind of the calculations, so according with this:

1. the "cnexp" method used by NERUON is the exponential euler method??, or which numerical methods use NEURON by default, and what other options have for solve the differential equations??

2. When you solved the cable equations mathematically usually there is a parameter that is related with the morphology of the "cilinder" and it is something like this:

var = rad*dt / (2*Ra*dx^2)
where rad and dx are the radius and lenght of the cilinder segment, and Ra is the citoplasmic resistivity

But i for example in my code i just solve the equations using : var = rad*dt / (Ra*dx^2) without the multiplication of 2, and it works fine, but when i add the "2" multiplier my results change signifincantly
So in theory i know that it must be multiplied by 2, however in practice it seems works fine without 2

So this is my confusion, where can i know how exactly neuron is computing these two parameters (the radius and length) in order to solved the cable equation??

thanks :)

[ted]

Many of these answers are in The NEURON Book, in the NEURON Forum, or in one of the papers listed on the Documentation page at NEURON's WWW site http://www.neuron.yale.edu/neuron/docs

One thing that might account for different results: the style of spatial discretization. NEURON uses the central difference approximation to the second derivative of distance. This means that each compartment has the equivalent circuit

Code: Select all

o-- rax/2 --+-- rax/2 --o
            |
        membrane
            |
         ground

so the resistance between the center of the soma and dendrite in your model is half of the axial resistance of the soma plus half the axial resistance of the dendrite. I haven't examined your code closely enough to decide whether you are using the same approach.

1. the "cnexp" method used by NERUON is the exponential euler method??

For fixed dt the default integration method is fully implicit Euler.

When you solved the cable equations mathematically usually there is a parameter that is related with the morphology of the "cilinder"

NEURON takes care of this automatically. Each segment ("compartment") is treated as a cylinder. Axial resistance is computed from cytoplasmic resistivity Ra in exactly the same way as it would be for a cylinder of length L/nseg and cross-sectional area PI*diam^2/4. Channel conductances associated with density mechanisms are computed as the product of the conductance density (in S/cm2) and the surface area in cm2 (PI*diam*L/nseg times the conversion factor 10-8 cm2/um2).