How to model diffusible protein in nmodl file with kinetic?

[mazid]

Hi,
I'm currently working on a model of the mapk kinase pathways in a CA1 neuron type. I have a problem with the diffusion of a protein in the kinectics block. It's clearly doesn't work so i'm doing something wrong or i don't understand something. So here are my questions:
1- How to model diffusible protein in nmodl file with kinetic?
2- How to declare a variable that we want to diffuse but that is not an ion?

THe physical structure is a spine + a neck + a dendrite (simplified)

Here are my codes:
FIle into the spine

Code: Select all

: Calcium ion accumulation with radial and longitudinal diffusion
NEURON {
   SUFFIX cadifus
   USEION ca READ cai, ica WRITE cai
	USEION RasGTP READ RasGTPi WRITE RasGTPi VALENCE 0
   	USEION RasGDP READ RasGDPi WRITE RasGDPi VALENCE 0
   RANGE  phi, beta
}

UNITS {
   (molar) = (1/liter)
   (mM)    = (millimolar)
   (um)    = (micron)
   (mA)    = (milliamp)
   FARADAY = (faraday)  (10000 coulomb)
   PI      = (pi)       (1)
}

PARAMETER {
	:DRas = 0.65 (um2/ms)
	DRas = 1000 (um2/ms)
   phi 	= 0.25(/ms) :(1/4)
	:phi	= 13.33 (ms)
	beta = 17.402 
	ceiling	= 2	(mM)
	caiBase = 50e-6 (mM)
	kf_cam1 = 1 (/mM-ms)
	kr_cam1 =  8.4853e-3 (/ms)
	kf_cam2 =  8.4846 (/mM-ms)
	kr_cam2 =  8.4853e-3 (/ms)
	kf_cam3 =  3.6001 (/mM-ms)
	kr_cam3 =  10e-3 (/ms)
	kf_cam4 =  0.465 (/mM-ms)
	kr_cam4 =  10e-3 (/ms)
	kf_gef =  2.75 (/mM-ms)
	kr_gef =  0.008333333e-3(/ms)
	kf_rasGDP = 0.1 (/mM-ms)
	kr_rasGDP = 0.003703703703703704e-3 (/ms)
	:kf_rasGTP = 0.003703703703703704e-3 (/ms)
	:kr_rasGTP = 0
}

ASSIGNED {
	diam (um)
	ica (mA/cm2)
	cai (mM)
	 
}

STATE {
   : ca[0] is equivalent to cai
   : ca[] are very small, so specify absolute tolerance
   ca       (mM) <1e-10>
   ca_modif (mM)
  Cam (mM)
  Cam_Ca (mM) 
  Cam_Ca2 (mM)
  Cam_Ca3 (mM)
  Cam_Ca4 (mM)
  Gef  (mM)
  Gef_activated  (mM)
  RasGDPi (mM)
  RasGTPi (mM)

}

BREAKPOINT { SOLVE state METHOD sparse }


INITIAL {
	ca = caiBase
	ca_modif = 0
	Cam = 0.02
	Cam_Ca = 0
	Cam_Ca2 = 0
	Cam_Ca3 = 0
	Cam_Ca4 = 0
	Gef = 0.0006
	Gef_activated = 0
	RasGDPi = 0.0006
	RasGTPi = 10
	
	}
	
KINETIC state {
   	COMPARTMENT PI*diam {RasGDPi RasGTPi}
	LONGITUDINAL_DIFFUSION DRas*diam {RasGDPi}
	LONGITUDINAL_DIFFUSION DRas*diam {RasGTPi}
   ~ ca << ( ((- beta * ica) /diam)- (phi * (cai - caiBase)) ) 
    cai = ca 
    ca_modif = ca - 0.00045
  ~ ca_modif + Cam <-> Cam_Ca (kf_cam1, kr_cam1)
  ~ ca_modif + Cam_Ca <-> Cam_Ca2 (kf_cam2, kr_cam2)
  ~ ca_modif + Cam_Ca2 <-> Cam_Ca3 (kf_cam3, kr_cam3)
  ~ ca_modif + Cam_Ca3 <-> Cam_Ca4 (kf_cam4, kr_cam4)
  ~ Cam_Ca4 + Gef <-> Gef_activated (kf_gef, kr_gef)
  ~ RasGDPi <-> RasGTPi (kf_rasGDP*Gef_activated, kr_rasGDP)
  
  :~ RasGTP <-> RasGDP (kf_rasGTP, kr_rasGTP)
  

}

File into the dendrite

Code: Select all

: Calcium ion accumulation with radial and longitudinal diffusion
NEURON {
   SUFFIX cadifus_dend
   USEION ca READ cai, ica WRITE cai
   USEION RasGTP READ RasGTPi WRITE RasGTPi VALENCE 0
   USEION RasGDP READ RasGDPi WRITE RasGDPi VALENCE 0
   RANGE  phi, beta
}

UNITS {
   (molar) = (1/liter)
   (mM)    = (millimolar)
   (um)    = (micron)
   (mA)    = (milliamp)
   FARADAY = (faraday)  (10000 coulomb)
   PI      = (pi)       (1)
}

PARAMETER {
	:DRas = 0.65 (um2/ms)
	DRas = 1000 (um2/ms)
   phi 	= 0.075(/ms) :(1/4)
	:phi	= 13.33 (ms)
	beta = 17.402 
	ceiling	= 2	(mM)
	caiBase = 50e-6 (mM)
	kf_cam1 = 1 (/mM-ms)
	kr_cam1 =  8.4853e-3 (/ms)
	kf_cam2 =  8.4846 (/mM-ms)
	kr_cam2 =  8.4853e-3 (/ms)
	kf_cam3 =  3.6001 (/mM-ms)
	kr_cam3 =  10e-3 (/ms)
	kf_cam4 =  0.465 (/mM-ms)
	kr_cam4 =  10e-3 (/ms)
	kf_gef =  2.75 (/mM-ms)
	kr_gef =  0.008333333e-3(/ms)
	kf_rasGDP = 0 (/mM-ms)
	kr_rasGDP = 0.003703703703703704e-3 (/ms)
	:kf_rasGTP = 0.003703703703703704e-3 (/ms)
	:kr_rasGTP = 0
}

ASSIGNED {
	diam (um)
	ica (mA/cm2)
	cai (mM)
	 
}

STATE {
   : ca[0] is equivalent to cai
   : ca[] are very small, so specify absolute tolerance
   ca       (mM) <1e-10>
   ca_modif (mM)
  Cam (mM)
  Cam_Ca (mM) 
  Cam_Ca2 (mM)
  Cam_Ca3 (mM)
  Cam_Ca4 (mM)
  Gef  (mM)
  Gef_activated  (mM)
  RasGDPi (mM)
  RasGTPi (mM)
 
}

BREAKPOINT { SOLVE state METHOD sparse }


INITIAL {
	ca = caiBase
	ca_modif = 0
	Cam = 0.02
	Cam_Ca = 0
	Cam_Ca2 = 0
	Cam_Ca3 = 0
	Cam_Ca4 = 0
	Gef = 0.0006
	Gef_activated = 0
	RasGDPi = 0.0006
	RasGTPi = 0
	
	}
	
KINETIC state {
   	COMPARTMENT PI*diam {RasGDPi RasGTPi}
	LONGITUDINAL_DIFFUSION DRas*diam {RasGDPi}
	LONGITUDINAL_DIFFUSION DRas*diam {RasGTPi}
   ~ ca << ( ((- beta * ica) /diam)- (phi * (cai - caiBase)) ) 
    cai = ca 
    ca_modif = ca - 0.00045
  ~ ca_modif + Cam <-> Cam_Ca (kf_cam1, kr_cam1)
  ~ ca_modif + Cam_Ca <-> Cam_Ca2 (kf_cam2, kr_cam2)
  ~ ca_modif + Cam_Ca2 <-> Cam_Ca3 (kf_cam3, kr_cam3)
  ~ ca_modif + Cam_Ca3 <-> Cam_Ca4 (kf_cam4, kr_cam4)
  ~ Cam_Ca4 + Gef <-> Gef_activated (kf_gef, kr_gef)
  ~ RasGDPi <-> RasGTPi (kf_rasGDP*Gef_activated, kr_rasGDP)

  :~ RasGTP <-> RasGDP (kf_rasGTP, kr_rasGTP)
  

}

The molecules in diffusion are RasGTP.

Your help will be greatly appreciated.
Thank you in advance.

[ted]

First suggestion: there is no substantive difference between these two models--they differ only in the values of a very few parameters--so use a single mod file for both of them. This will simplify development and debugging. Make the default parameters suitable for either spine or dendrite, and preserve the other parameter values as comments, e.g.
kf_rasGDP = 0.1 (/mM-ms) : for spine, 0 for dendrite
Insert a comment at the very top that indicates which location the default parameter values are suitable for, e.g.
: diffusion of RasGTP and RasGDP
: default parameter values are for spine
: dendritic values are indicated in comments

Second suggestion: the last time I implemented a "second messenger accumulation/diffusion" mechanism, the USEION statement required specification of a nonzero value for VALENCE, e.g.
USEION foo READ fooi WRITE fooi VALENCE 1
It's not going to mess up charge balance equations, as long as you don't also WRITE a 2nd messenger current--that is, don't do
USEION foo READ fooi WRITE fooi, ifoo VALENCE 1
Well, maybe even that would work as long as you don't calculate a nonzero value for ifoo in the INITIAL or BREAKPOINT block.

Third suggestion: check units for consistency with modlunit. If you find an inconsistency, figure out how to fix it. modlunit's first whack at the mod file produces this message

Code: Select all

Checking units . . .
1-15 m3/sec : relevant area * diffusion constant must
   be micron2 micron2/ms (1-21 m4/s) at line 87 in file xd.mod
	LONGITUDINAL_DIFFUSION DRas*diam<<ERROR>> {RasGDPi}

and if you look close, and compare with example 9.8 in the NEURON book, you'll see what needs to be done to the LONGITUDINAL_DIFFUSION statements. If you don't have the book, pick up this pdf
https://www.neuron.yale.edu/ftp/ted/boo ... xedref.pdf

After you fix that, check with modlunit again and see if the error goes away. If modlunit reports another error,

Code: Select all

REPEAT
  fix the error
  check with modlunit
UNTIL modlunit doesn't complain

By the way, just for future reference, you may want to download a recent gzipped tar file of NEURON's source code e.g.
http://www.neuron.yale.edu/ftp/neuron/v ... 232.tar.gz
then expand it and look at the examples in nrn/share/examples/nrniv/nmodl

[mazid]

Thank you for the quick reply! I'll do as you say and will let you know the results.

[mazid]

Hi Mr Ted,
I've done some debugging and come out with something that's working i.e i've a diffusion of RasGTPi.

Code: Select all

: Calcium ion accumulation with longitudinal diffusion
: diffusion of RasGTP and RasGDP
: default parameter values are for spine
: dendritic values are indicated in comments
: USEION foo READ fooi WRITE fooi VALENCE 1
NEURON {
   SUFFIX cadifus
   USEION ca READ cai, ica WRITE cai
   USEION RasGTP READ RasGTPi WRITE RasGTPi VALENCE 1
   USEION RasGDP READ RasGDPi WRITE RasGDPi VALENCE 1
   RANGE  phi, beta, kf_rasGDP
}

UNITS {
   (molar) = (1/liter)
   (mM)    = (millimolar)
   (um)    = (micron)
   (mA)    = (milliamp)
   FARADAY = (faraday)  (10000 coulomb)
   PI      = (pi)       (1)
}

PARAMETER {
	DRas = 0.65 (um2/ms)
	:DRas = 1000 (um2/ms)
   phi 	= 0.25(/ms) :(1/4)
	:phi	= 13.33 (ms)
	beta = 17.402  
	ceiling	= 2	(mM)
	caiBase = 50e-6 (mM)
	kf_cam1 = 1 (/mM-ms)
	kr_cam1 =  8.4853e-3 (/ms)
	kf_cam2 =  8.4846 (/mM-ms)
	kr_cam2 =  8.4853e-3 (/ms)
	kf_cam3 =  3.6001 (/mM-ms)
	kr_cam3 =  10e-3 (/ms)
	kf_cam4 =  0.465 (/mM-ms)
	kr_cam4 =  10e-3 (/ms)
	kf_gef =  2.75 (/mM-ms)
	kr_gef =  0.008333333e-3(/ms)
	kf_rasGDP = 0.1 (/mM-ms) : for spine, 0 for dendrite
	kr_rasGDP = 0.003703703703703704e-3 (/ms)
}

ASSIGNED {
	diam (um)
	L (um)
	ica (mA/cm2)
	cai (mM)
}

STATE {
   ca       (mM) <1e-10>
   ca_modif (mM)
  Cam (mM)
  Cam_Ca (mM) 
  Cam_Ca2 (mM)
  Cam_Ca3 (mM)
  Cam_Ca4 (mM)
  Gef  (mM)
  Gef_activated  (mM)
  RasGDPi (mM)
  RasGTPi (mM)
}

BREAKPOINT { SOLVE state METHOD sparse }


INITIAL {
	ca = caiBase
	ca_modif = 0
	Cam = 0.02
	Cam_Ca = 0
	Cam_Ca2 = 0
	Cam_Ca3 = 0
	Cam_Ca4 = 0
	Gef = 0.0006
	Gef_activated = 0
	RasGDPi = 0.0012
	RasGTPi = 0
	}
	
KINETIC state {
   	COMPARTMENT PI*diam {RasGDPi RasGTPi}
	LONGITUDINAL_DIFFUSION DRas*diam {RasGDPi}
	LONGITUDINAL_DIFFUSION DRas*diam {RasGTPi}
   ~ ca << ( (- beta * ica)- (phi * (cai - caiBase)) ) 
    cai = ca 
    ca_modif = ca - 0.00045
  ~ ca_modif + Cam <-> Cam_Ca (kf_cam1, kr_cam1)
  ~ ca_modif + Cam_Ca <-> Cam_Ca2 (kf_cam2, kr_cam2)
  ~ ca_modif + Cam_Ca2 <-> Cam_Ca3 (kf_cam3, kr_cam3)
  ~ ca_modif + Cam_Ca3 <-> Cam_Ca4 (kf_cam4, kr_cam4)
  ~ Cam_Ca4 + Gef <-> Gef_activated (kf_gef, kr_gef)
  ~ RasGDPi <-> RasGTPi (kf_rasGDP*Gef_activated, kr_rasGDP)
  
  :~ RasGTP <-> RasGDP (kf_rasGTP, kr_rasGTP)
  

}

But the diffusion is very very very small in the dendrites and the spine. I know the

Code: Select all

	COMPARTMENT PI*diam {RasGDPi RasGTPi}
	LONGITUDINAL_DIFFUSION DRas*diam {RasGDPi}
	LONGITUDINAL_DIFFUSION DRas*diam {RasGTPi}

isn't conform to modlunits check but i did it on purpose cause i would like to make a 2D diffusion on the membrane! How could i specify it on the COMPARTMENT and LONGITUDINAL_DIFFUSION parameters? I look at EXPANDING NEURON to catch my mistake but i didn't find anything helpful!
Thank you in advance!

[ted]

diffusion is very very very small in the dendrites and the spine

I haven't had time to work through your NMODL code. However, given that the RMS displacement for 1 dimensional diffusion* over time interval t is sqrt(2*D*t), and the following diffusion coefficients
Ca in cytosol ~ 0.6 um2/ms
proteins in membrane ~ 1e-5 um2/ms
(see http://www.life.illinois.edu/crofts/bio ... sion1.html and http://www.sbcny.org/membrane_diffusion ... cients.htm)
I'd expect that, over any time interval, cytosolic Ca would diffuse about sqrt(0.6/1e-5) ~ 250 times farther than a protein in membrane would. Calcium is no speed demon, taking about 100 ms to achieve an RMS displacement of about 10 um, so a protein in membrane must be very slow indeed--. Might 100 ms be long enough for various reactions to severely attenuate whatever temporary longitudinal gradient might be induced by a localized metabotropic event?

*--Why one dimensional diffusion? Because I'm assuming that the concentration gradient is along the length of the neurite, i.e. longitudinal, not both longitudinal and circumferential. Diffusion in membrane may be at most 2 dimensional for a while (assuming very localized concentration), for which the RMS displacement is only sqrt(4*D*t)--not much faster than 1 dimensional diffusion.