two compartment Traub model, share membrane voltage

[jcmocte]

Hi,

I create a two compartment model for Traub neuron (soma and dendrite). So for SOMA i have mytraub.mod file and for DENDRITE i have CaDen.mod, K_C.mod and KAHP.mod files
My problem is when i connect this two sections (SOMA + DENDRITE) the output for the soma is the same than if i dont make the connection, so it seems that soma is not taking into account the voltage from dendrite. So why the SOMA output is the same if i connect the dendrite or not??. I simulate by separate each section and they work fine.

The SOMA depends of Na and K ions
The DENDRITE depends of Ca and K ions

The voltage in dendrite indeed changed, so it seem to take into account the soma, but in the other way around does not
How is the mechanism for share the voltage between sections??, is it automatically when i use the "connect" command??

Thanks for any help



This is my HOC file

Code: Select all

load_file("nrngui.hoc")		

create somaT,denT
access somaT
somaT {
  nseg = 1
  diam = 43.3	// (um)
  L = 43.3		// (um)
  Ra = 100e3	// (ohm*cm)
  cm = 3		// (uF/cm2) , capacitance
  
  insert traub
}

denT {
  nseg = 1
  diam = 2.5	// (um)
  L = 11		// (um)
  Ra = 100e3	// (ohm*cm)
  cm = 3		// (uF/cm2) , capacitance
 
  insert caden
  insert kahp
  insert kc
}

objectvar stim
somaT stim = new IClamp(0.5)
    stim.del = 10	// (mseg)
    stim.dur = 50	// (mseg)
    stim.amp = 2.3 //0.1767   //2.3 	// (nA)

// +++++++++++++++++++++
// CONNECT COMPARTMENTS
// +++++++++++++++++++++
connect somaT(1),denT(0)

This is mytraub.MOD file

Code: Select all

TITLE mytraub5.mod
UNITS {
	(mA) = (milliamp)
	(mV) = (millivolt)
	(S) = (mho)
	(umho) = (micromho)
} : end UNITS

NEURON {
	SUFFIX traub
	USEION na WRITE ina 
	USEION k READ ik WRITE ik	
    NONSPECIFIC_CURRENT itl	
	RANGE gnabar,gkbar,gna,gk,ena,ek 		: this variables are added the SUFFIX
	RANGE gtl,etl
	GLOBAL 	alphaM,betaM,alphaH,betaH,alphaP,betaP  
} : end NEURON

: Constants OR variables which can be changed from HOC program
PARAMETER {
	ena	= 50.0	(mV)
	ek		= -85.0	(mV)	
	etl 	= -70.0 (mV)	
        
	gnabar 	= 0.030   (S/cm2)
	gkbar		= 0.015	  (S/cm2)	
	gtl 		= 0.0003  (S/cm2)		

	maflag 		= 3
	malphaA 	= -0.32
	malphaB		= -4.0
	malphaV0	= 13.1
	mbflag 		= 3
	mbetaA 		= 0.28
	mbetaB		= 5.0
	mbetaV0		= 40.1
	
	haflag 		= 1
	halphaA 	= 0.128
	halphaB		= -18
	halphaV0	= 17.
	hbflag 		= 2
	hbetaA 		= 4.
	hbetaB		= -5.
	hbetaV0		= 40.
	
	: MODIFICADO channel "p" added
	paflag		= 4
	palphaA 	= -0.016
	palphaB		= -5.0
	palphaV0	= 35.1
	pbflag		= 4
	pbetaA 		= 0.25
	pbetaB		= -40.0
	pbetaV0		= 20.0
} : end PARAMETER

: Variables calculated by mechanism itself OR 
: computed by NEURON
ASSIGNED {
	v (mV)
	ina (mA/cm2)		
	ik (mA/cm2)	
	itl (mA/cm2)		
	gna (S/cm2)
	gk	(S/cm2)
	mtraub_inf
	mtraub_tau (ms)
	htraub_inf 
	htraub_tau (ms)
	ptraub_inf 	
	ptraub_tau (ms)
	alphaM betaM
	alphaH betaH
	alphaP betaP
} : end ASSIGNED 

: State variables to be solved
STATE { m h p }

: Equations which calculates the ions currents
: In Yiwei thesis, pag 98, equation 4.20(a)

BREAKPOINT {
  SOLVE states METHOD cnexp
  mh(v)
  gna = gnabar * mtraub_inf*mtraub_inf*h
  gk = gkbar * p	
  ina = gna*(v-ena) 
  ik = gk*(v-ek) 
  itl = gtl*(v-etl) 
} : end BREAKPOINT

: Initial values
INITIAL { 
 	mh(v)
	m = mtraub_inf
	h = htraub_inf
	p = ptraub_inf
}

DERIVATIVE states {
	mh(v)	
	m' = alphaM * ( 1 - m ) - betaM * m 
	h' = alphaH * ( 1 - h ) - betaH * h 
	p' = alphaP * ( 1 - p ) - betaP * p
 }

: +++++++++++++++++++++++++++++++++++++++++++++++++++
: +++++++++++++++++++++++++++++++++++++++++++++++++++
: ASSIGNMENT PROCEDURES
: +++++++++++++++++++++++++++++++++++++++++++++++++++
: +++++++++++++++++++++++++++++++++++++++++++++++++++

UNITSOFF
:-------------------------------------------------------------------
: NOTE : 0 = m and 1 = h
PROCEDURE mh (v(mV)) {
	
	alphaM = -0.32*(56.9+v) /(exp(-(56.9+v)/4.0)-1)
	betaM = 0.28*(29.9+v) / (exp((29.9+v)/5.0)-1)
	mtraub_inf = alphaM / (alphaM + betaM)

	alphaH = 0.128*exp(-(v+53.0)/18.0)
	betaH = 4.0 /(exp(-(30.0+v)/5.0)+1.0)
	htraub_inf = alphaH / (alphaH + betaH)
	
	alphaP = -0.016*(34.9+v)/(exp(-(34.9+v)/5.0)-1.0)
	betaP = 0.25*exp(-(v+50.0)/40.0)
	ptraub_inf = alphaP / (alphaP + betaP)
	
} : end PROCEDURE mh (v)
UNITSON

This is CaDen.MOD file

Code: Select all

UNITS { 
	(mV) = (millivolt) 
	(mA) = (milliamp) 
}
 
NEURON { 
	SUFFIX caden
	USEION ca WRITE ica,cai
	RANGE  gcamax, is, gs
    GLOBAL alphaS,betaS
}

PARAMETER { 
  gcamax = 0.010 	(mho/cm2)
  eca    = 70.0	(mV) 
  scale_ca   = 0.13 
  time_ca    = 0.075  
} 

ASSIGNED { 
  v 	(mV)
  gs	(mho/cm2)
  myica	(mA/cm2) 
  ica	(mA/cm2) 
  alphaS betaS 
}
 
STATE { s cai }

INITIAL { 
	s = .009  :alphaS / ( alphaS + betaS ) 	
	cai = 0.2
}

BREAKPOINT { 
  SOLVE states METHOD cnexp
  gs = gcamax * s * s
  ica = gs * (v-eca)
  myica = 10.0 * s * s * (v-eca)
}
 
DERIVATIVE states { 
	settables(v)
	s' = alphaS * ( 1 - s ) - betaS * s 
	cai' = -(scale_ca * myica) - (time_ca * cai)
}

UNITSOFF 

PROCEDURE settables(v) { 
	alphaS = 1.6 / (exp( -(v + 5.0)/13.89 )+1)
	betaS = 0.02*(v+18.9)/(exp( (v + 18.9)/5.0 )-1)
}
UNITSON

This is KAHP.MOD file

Code: Select all

UNITS { 
	(mV) = (millivolt) 
	(mA) = (milliamp) 
}
 
NEURON { 
	SUFFIX kahp
	USEION k READ ik WRITE ik 
	USEION ca READ cai
	RANGE  gKAHPmax,gkahp
    GLOBAL alphaQ,betaQ
}

PARAMETER { 
  gKAHPmax = 0.0008 (mho/cm2)
  ek = -85.0 		(mV) 
} 

ASSIGNED { 
  v 	(mV)
  gkahp	(mho/cm2)
  ik	(mA/cm2)
  cai  
  alphaQ betaQ 
  alphaS
}
 
STATE { q }

INITIAL { 
	q = 0.01
}

BREAKPOINT { 
  SOLVE states METHOD cnexp
  gkahp = gKAHPmax * q
  ik = gkahp*(v-ek)
}
 
DERIVATIVE states { 
	rates( cai )
	q' = alphaQ * ( 1 - q ) - betaQ * q 
}

UNITSOFF 

PROCEDURE rates(cai) { 
  if ((0.2e-4)*cai < 0.01) {
    alphaQ = (0.2e-4)*cai
  } else {
    alphaQ = 0.01
  }
  betaQ = 0.001
}

UNITSON

This is K_C.MOD file

Code: Select all

UNITS { 
	(mV) = (millivolt) 
	(mA) = (milliamp) 
}
 
NEURON { 
	SUFFIX kc
	USEION k READ ik WRITE ik 
	USEION ca READ cai
	RANGE  gkcmax, gk
    GLOBAL alphaC,betaC
}

PARAMETER { 
  gkcmax = 0.015 	(mho/cm2)
  gcamax = 0.010 	(mho/cm2)
  glden	 = 0.0003  (mho/cm2)
  gKAHPmax = 0.0008 	(mho/cm2)

  Vthc = -20.0  (mV)
  ek = -85.0	(mV) 
  Ca_th = 1.0
} 

ASSIGNED { 
  v 	(mV)
  gk 	(mho/cm2)
  ik	(mA/cm2) 
  cai
  alphaC betaC	
}
 
STATE { c }

INITIAL { 
	c = .007  :alphaC / ( alphaC + betaC )
}

BREAKPOINT { 
  SOLVE states METHOD cnexp

  if( 0.004 * cai < Ca_th ) {	: this comparison solves equation 5.15j
    gk = gkcmax * c * 0.004 * cai
  } else {
    gk = gkcmax * c
  }
  ik = gk * (v-ek)
}
 
DERIVATIVE states { 
	settables(v) 
	c' = alphaC * ( 1 - c ) - betaC * c 
}

UNITSOFF 

PROCEDURE settables(v) { 

	if( v < Vthc ) {
		alphaC = (2 / 37.95) * ( exp( ( v + 60 ) / 11 - ( v + 63.5 ) / 27 ) )
		betaC  = 2 * exp( -(v + 63.5 ) / 27 ) - alphaC 
		: Note that there is typo in the paper - missing minus sign in the front of 'v'
	}else{
		alphaC = 2 * exp(-(v + 63.5 ) / 27 )
		betaC  = 0
	}
}

UNITSON

[ted]

why the SOMA output is the same if i connect the dendrite or not??

Is it really the same? You have compared the results and found that there is absolutely no effect on soma.v, at all? the results are numerically identical?

The effect of soma on den.v is very small; I suspect that is because cm and Ra are so large. cm is 3 times larger than the usual value that one sees in models (1 uf/cm2), and Ra is about 1000 times larger than the usual value (100-160 ohm cm). With such large values, there is a strong low-pass filter for the spread of electrical signals in your model.

How is the mechanism for share the voltage between sections??, is it automatically when i use the "connect" command??

Voltage isn't "shared." Coupling between adjacent compartments is by the spread of current between them. connect tells NEURON that there is an electrical connection between the two sections.

A couple of comments about your model's code:
access somaT
makes somaT the "default" compartment. That's OK if there's something special about somaT, that is, if somaT is the most important compartment in the model.

connect somaT(1),denT(0)
has two problems. First, it makes the soma the child of denT, which is not what you'd expect to happen if somaT is the most important compartment in the model. It's probably not a good idea, because it makes denT the "root" of the model, that is, the section that has no parent. Second, denT's 0 end is the root node of the model, and anatomical distance from the root node increases as you move from the 0 end of denT toward its 1 end, but it decreases as you move from the 0 end of somaT toward somaT's 1 end. This can cause problems with space plots (plots that show how a variable depends on distance) and with code that sets up inhomogeneous biophysical properties (i.e. properties that depend on anatomical distance from a reference point). It would be better to
connect denT(0),somaT(1)

[jcmocte]

Is it really the same? You have compared the results and found that there is absolutely no effect on soma.v, at all?

Yes they are the same

Voltage isn't "shared." Coupling between adjacent compartments is by the spread of current between them. connect tells NEURON that there is an electrical connection between the two sections

Well i mean, when you are solving the cable equation for many sections, you should now the voltage of your "neighbor" sections right?, yes, you spread the current but in order to calculate the current you need to know the voltage
So thats i want to know, if this "neighbor voltages" are taken into consideration between sections?? because in my MOD files i only see declarated the "v" variable, but i dont know if this "v" is calculated by every section independently or if both "v" for every section are automatically combined and take into consideration when i make the connection between them

Another questions:

what does mean the "i" and "o" suffix for ions??, e.g. cai, cao, ki, ko, etc.

When i run the utility "modlunit" it gives me a warning about one state variable, that it is 1/1000 when it should be 1, should this affect my results??
I just want to know how important is to be successful (without errors) in the unit test??

thanks again for your help

[ted]

I tried reproducing your results, and here's what I found.

First, it looks like your model's initial condition is not v = -65 mV. Was it -70 mV?

Is it really the same? You have compared the results and found that there is absolutely no effect on soma.v, at all?

Yes they are the same

Not so. I ran two simulations--a "control" simulation with the model as is, and a "test" simulation in which Ra was set to 1e9 in both sections so that coupling between the two sections would be much smaller--and compared the "control" and "test" somaT.v(0.5). Result: within the first 0.5 ms, there were noticeable differences in the third and fourth decimal places between the somaT.v(0.5) recorded under "control" and "test" conditions. The differences grew progressively larger throughout the course of the simulation, so that they became quite obvious for the last half of the simulation when the traces were superimposed in the same plot. That could not have happened if denT had no effect on somaT. The effect occurs because of conservation of charge. somaT can't affect denT unless charge spreads from somaT to denT. If that happens, the loss of charge from somaT will affect the time course of somaT.v(0.5).

By the way, there are some things in the model's source code that may be errors. Here is a partial list:
Why is Ra = 1e5? This is 3 orders of magnitude larger than the values commonly used for cytoplasmic resistivity.
Specific membrane capacitance cm is 3 uf/cm2. This seems large.
In caden.mod this formula
cai' = -(scale_ca * myica) - (time_ca * cai)
omits diameter (which is needed to take surface/volume ratio into account) and Faraday's constant--was this developed from a model implemented in some other simulation environment?
For sections that have the caden mechanism, eca will not be 70 mV. Instead, it will be calculated from the Nernst equation, using default values for extracellular Ca concentration and temperature (2 mM and 6.3 deg C, respectively) since there are no statements that would assign other values to these variables.

As I mentioned in a different thread, it would be "better" if the connect statement
connect somaT(1),denT(0)
were written
connect denT(0), somaT(1)
so that
--the soma is the root section
and
--the child's 0 end is connected to the parent's 1 end
simply because this would avoid possible future problems if you wanted to use space plots (plots of a range variable vs. distance from a reference point) or assign spatially varying parameters whose values are governed by distance from a reference point.

So thats i want to know, if this "neighbor voltages" are taken into consideration between sections??

The answer is yes. NEURON sets up the equations for the spatially discretized cable equation, according to the specification of length, diameter, branched architecture, and properties of membrane and cytoplasm that are specified by the user in hoc or Python. It automatically takes care of boundary conditions.

because in my MOD files i only see declarated the "v" variable

and this is the local v, i.e. membrane potential in a compartment (segment).

but i dont know if this "v" is calculated by every section independently or if both "v" for every section are automatically combined and take into consideration when i make the connection between them

The system equations for a model consist of a collection of algebraic and ordinary differential equations that are solved simultaneously.

Maybe it would help to read about how NEURON works. Have you read
Hines, M.L. and Carnevale, N.T.
The NEURON simulation environment.
Neural Computation 9:1179-1209, 1997
(available for download or online browsing from links at http://www.neuron.yale.edu/neuron/nrnpubs)
?

what does mean the "i" and "o" suffix for ions??, e.g. cai, cao, ki, ko, etc.

"i"ntra- and "e"xtracellular concentrations.

When i run the utility "modlunit" it gives me a warning about one state variable, that it is 1/1000 when it should be 1, should this affect my results??

Good question. What is the exact error message, and what is the source code?

[jcmocte]

Hi Ted, thanks for your time

I already put the connection order as you mention before : connect denT(0), somaT(1)
However it does not make change in the result.

The main goal is that i want to implement the model in "Intrinsic and Network Rhythmogenesis in a Reduced Traub Model for CA3 Neurons" by Pinsky and Rinzel, so i want to implement their formulas in MATLAB and compare it with NEURON, however i am not sure if my MATLAB program is wrong or i am misunderstanding the parameters in NEURON, so that is my problem.

I want to program a soma connected to a dendrite (2-compartment model)
So the way i "coupled" my two sections in MATLAB is using the cable equation as follow. For the soma i solved this equation:

Ome = (gNa*ENA + gK*EK + GLEAK*ELEAK + (Iinj(i)/Asoma));
Phi = (gNa + gK + GLEAK);
Tao = a_soma/(Ra*a_soma*a_soma); %a_soma is the radius of my soma
A(i,l) = (Tao)*V(i,l+1) + Ome; %V(i,l+1) is the voltage for right compartment (dendrite)
B(i,l) = (Tao)+Phi;

And for dendrite i solved this equation:

Ome = (Gca*ECA + Gk_ahp*EK + Gk_c*EK + (Iinj_den(i)/Aden) ); % in this case current in dendrite is zero
Phi = (Gca + Gk_ahp + Gk_c);
Tao = a_den/(Ra*dx*dx); %a_den is the radius of my dendrite and dx is the lenght of dendrite
A(i,l) = (Tao)*V(i,l-1) + Ome; %V(i,l-1) is the voltage for left compartment (soma)
B(i,l) = (Tao)+Phi;

And for each compartment i solve the equations using exponential Euler method:

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))));

However i dont know how NEURON works "under the hood" with the "coupling" compartments, i have doubts in the "Tao" variable and i dont know i am manage well the geometric parameters of the sections

Below i answer your questions and explain more my situation:

First, it looks like your model's initial condition is not v = -65 mV. Was it -70 mV?

Yes my initial condition is -70 mV

By the way, there are some things in the model's source code that may be errors. Here is a partial list:
Why is Ra = 1e5? This is 3 orders of magnitude larger than the values commonly used for cytoplasmic resistivity.
Specific membrane capacitance cm is 3 uf/cm2.

Yes i think Ra is just 100 but capacitance is 3 because the paper i am working with, they manage that capacitance

In caden.mod this formula
cai' = -(scale_ca * myica) - (time_ca * cai)
omits diameter (which is needed to take surface/volume ratio into account) and Faraday's constant--was this developed from a model implemented in some other simulation environment?

In the paper i am working they just manage the formula as i put before, at the final i put all the equations and parameters which i want to implement and these equations were extracted for the paper mentioned

When i run the utility "modlunit" it gives me a warning about one state variable, that it is 1/1000 when it should be 1, should this affect my results??

Good question. What is the exact error message, and what is the source code?[/quote]

This is the message when i check "mytraub.mod" file

Checking units of ./mytraub4.mod

units: 1
units: 1000 /sec
The units of the previous two expressions are not conformable
at line 110 in file ./mytraub.mod
m' = alphaM * ( 1 - m ) - betaM * m<<ERROR>>


And for CaDen.mod, KAHP.mod and K_C.mod i got the message:

cai must have the units, milli/liter, instead of .
at line 17 in file ./KAHP.mod
NEURON<<ERROR>> {

Here is the link for the paper and the summary for the equations:

http://ccc.inaoep.mx/~jcmoctezuma/docs/Traub.pdf

[jcmocte]

By the way these are my files for MATLAB implementation

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

[lneisenman]

FWIW, Bill Lytton implemented a version of this model in NEURON.

https://senselab.med.yale.edu/modeldb/S ... odel=35358

Larry

[jcmocte]

thanks Larry,

I already have a look to this model, but it is to huge, the main HOC file is about 4,000 lines, and the way how they implement the mechanism is a little strange for me, for example in one MOD files they just declared some parameters and include the INC file, in some files they have repetitive equations and so on. So it was difficult to me follow it

[jcmocte]

In fact i based my MOD files in theirs, i also put my parameters but i get a totally different output, so it was cumbersome try to fix their files to my parameters
so thats why i want to implement by my self this model and make sure i am implementing the cable equation in the right way

[ted]

Bill Lytton is pretty careful about what he does and how he does it. The comments in his README file
https://senselab.med.yale.edu/modeldb/S ... b12\README
suggest that, for various technical reasons, it may not be possible to replicate the Pinsky-Rinzel results perfectly with NEURON. If that's the case, how can such a NEURON implementation be useful as a standard for comparison with a Matlab implementation of their model? unless, of course, the Matlab implementation differs from the original Pinsky-Rinzel implementation in the same way as the NEURON implementation does.

So unless you have a deep and abiding interest in the Pinsky-Rinzel model itself, it might be a more fruitful investment of your time and effort to try replicating some other model.

[ted]

i dont know how NEURON works "under the hood" with the "coupling" compartments

There's no great mystery about it. NEURON uses the central difference approximation for discretizing space. Given the fact that the Pinsky-Rinzel model has only two compartments, that simplifies the problem. Use two sections, each with nseg = 1. All of the membrane properties of each section will be treated as if they are lumped into the middle of the section. The resistance between the two compartments will be half of the axial resistance of one compartment plus half of the axial resistance of the other compartment. Done.

In caden.mod this formula
cai' = -(scale_ca * myica) - (time_ca * cai)
omits diameter (which is needed to take surface/volume ratio into account) and Faraday's constant--was this developed from a model implemented in some other simulation environment?

In the paper i am working they just manage the formula as i put before, at the final i put all the equations and parameters which i want to implement and these equations were extracted for the paper mentioned

NEURON has a "policy" or "standard way" of dealing with density mechanisms: their currents are reported in density units, but total current is surface area times current density. XPP doesn't have such a policy. If the equations used by Pinsky and Rinzel assume that currents and conductances and capacitances are in "absolute" units (not mA/cm2 or mho/cm2 or uf/cm2), you can't just plop them into NMODL code for a density mechanism and expect things to magically turn out OK. You have to make sure that each compartment in your NEURON implementation has the same total membrane current, total ionic conductance, and total membrane capacitance as in P&R's implementation. You have to figure out what surface areas your sections must have, and what conductance and capacitance densities they should have, and you have to choose the correct section diameter and length and Ra so that the series resistance between your compartments will be the same as P&R used. It's a pain but it's the way things are.

When i run the utility "modlunit" it gives me a warning about one state variable, that it is 1/1000 when it should be 1, should this affect my results??

Good question. What is the exact error message, and what is the source code?

This is the message when i check "mytraub.mod" file

Checking units of ./mytraub4.mod

units: 1
units: 1000 /sec
The units of the previous two expressions are not conformable
at line 110 in file ./mytraub.mod
m' = alphaM * ( 1 - m ) - betaM * m<<ERROR>>

This isn't a fatal error, but you need to eliminate the error message so that modlunit can check the rest of the mod file (it stops at the first error it finds). To eliminate this error, you have to specify that alphaM and betaM have units (/ms). Do this in the ASSIGNED block with statements of the form
alphaM (/ms)
Then run modlunit again and see if it finds anything else to complain about.

And for CaDen.mod, KAHP.mod and K_C.mod i got the message:

cai must have the units, milli/liter, instead of .
at line 17 in file ./KAHP.mod
NEURON<<ERROR>> {

This might be a big problem, or it might not. In whatever block that declares cai (ASSIGNED? PARAMETER?) you need to make sure that it is declared like this:
cai (milli/liter)
Then check again with modlunit and see if it finds any other errors.

Good luck. This can be tedious. But (if and) when everything works, ah, success!

[jcmocte]

Hi Ted,

Well even if running the modlunit to the Bill Lytton original files, i still have the same warnings and errors, so i supposed that it was ok, but i will try the changes you mentioned.

The reason why is i want to implement this in NEURON is that i want to compare my results with some well know and stable tool as NEURON, other way, how do i know if my program is right??

Also i try to simulate the HH model with a passive membrane, but still i have different results from NEURON and my program, so that why i am interested in to know how represent my equations in a 2-compartment model, because maybe i am misunderstanding something. Mainly in the "coupling" process when i take into consideration Ra, soma radius and dendrite length.

Sorry for insist but i think this is important for my issues. For example i always do for the soma

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_soma is soma radius
A(i,l) = (Tao)*V(i,l+1) + Ome;
B(i,l) = (Tao)+Phi;

And for dendrite i do this:

Ome = Eden/(Rm); %Passive dendrite
Phi = 1/(Rm);
Tao = a_den/(Ra*dx*dx); %a_den is the radius of the dendrite and dx is the length
A(i,l) = (Tao)*V(i,l-1) + Ome;
B(i,l) = (Tao)+Phi;

And then for each one i solved by exponential Euler like this:

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)));

And i just can not get the same behaviour than NERUON, and this is a simple HH model + passive dendrite. So if anyone see something wrong in my program please let me know. The program is very easy and the complete code is here

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.025;   %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 (Ohm*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

alpha_h = 0.07*exp(-(V(1,1)+70)/20);   % HH model
alpha_m = -0.1*(45+V(1,1))./(exp(-(45+V(1,1))/10)-1); 
alpha_n =  -0.01*(60+V(1,1)) ./ (exp(-(60+V(1,1))/10)-1);
beta_h = 1 ./ (exp(-(40+V(1,1))/10)+1);
beta_m = 4*exp(-(V(1,1)+70)/18);
beta_n = 0.125*exp(-(V(1,1)+70)/80);

m(1,:) = alpha_m/(alpha_m +beta_m);	% initial m
h(1,:) = alpha_h/(alpha_h +beta_h);	% initial h
n(1,:) = alpha_n/(alpha_n +beta_n);	% initial n
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
        
        alpha_h = 0.07*exp(-(V(i,l)+70)/20);   % HH model
        alpha_m = -0.1*(45+V(i,l))./(exp(-(45+V(i,l))/10)-1); 
        alpha_n =  -0.01*(60+V(i,l)) ./ (exp(-(60+V(i,l))/10)-1);
        beta_h = 1 ./ (exp(-(40+V(i,l))/10)+1);
        beta_m = 4*exp(-(V(i,l)+70)/18);
        beta_n = 0.125*exp(-(V(i,l)+70)/80);
        
        %First kinetics
        gamma = alpha_h + beta_h;
        h(i+1,l) = (alpha_h/gamma) + (h(i,l) - (alpha_h/gamma))*exp(-hstep*gamma); 
        gamma = alpha_m + beta_m;
        m(i+1,l) = (alpha_m/gamma) + (m(i,l) - (alpha_m/gamma))*exp(-hstep*gamma); 
        gamma = alpha_n + beta_n;
        n(i+1,l) = (alpha_n/gamma) + (n(i,l) - (alpha_n/gamma))*exp(-hstep*gamma); 
                
        %Exponential Euler's rule
        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
    %figure
    plot(t,V(:,l),'linewidth',1.7,'color',colorgraph(l));   
    cad=sprintf('Compartment %d',l);
    title(cad);
    hold on;
    %plot(t,Iinj/Asoma,'c--');
    grid on;
end

[ted]

Suggest you start with the very simplest model possible--single compartment, totally passive--and make sure it works properly. Then add another compartment and make sure it works. If necessary, print out intermediate results at each time step and make sure they are correct. Only after everything works for a totally passive model, is it time to add a voltage-gated current. Start by adding the delayed rectifier to a single compartment model and make sure it works. Then add a second compartment etc.