a hemi-sopherical model

[along82]

What I am trying to model is a single cell with simple semi-sopherical 3d morphology.
It is said that NEURON treats section as straight cylindrial colum, so I image that I have to implement many cylindrial colum shaped section units to approximate a single semi-sopherical shaped neuron. Does it make sense or is there any other good alternatives?

[ted]

You mean "hemispherical"?
Do you need that particular shape to achieve a particular surface/volume
ratio (i.e. is an ion accumulation mechanism involved?), or this just for
aesthetics? If the latter, use pt3d specification of geometry. If the former,
then forget about pt3d--just set L = diam/2, but be sure to base the surface/
volume ratios of adjacent compartments on spherical geometry, not cylindrical.

[along82]

Thanks for reply.

You mean "hemispherical"?
++++++++++++++++++++++++++++++++
Exactly!

Do you need that particular shape to achieve a particular surface/volume ratio (i.e. is an ion accumulation mechanism involved?),
++++++++++++++++++++++++++++++++
Exactly, what I am doing is a single neuronal model for detailed analysis of intracellular calcium dynamics

If the former, then forget about pt3d--just set L = diam/2, but be sure to base the surface/volume ratios of adjacent compartments on spherical geometry, not cylindrical.
++++++++++++++++++++++++++++++++
I don't understand "base the surface/volume ratios of adjacent compartments on spherical geometry". Are there any items in the reference manul or I have to define some new variable for this shape by myself?
Actually in my plan, I would want to firstly divide the hemispherical neuron into 100 segments in the longitudial direction, then each segment will be further divided into 101 shells in the radial direction. But a little overwhelmed on how to define this neuron geometry at this moment.

[ted]

I don't understand "base the surface/volume ratios of adjacent compartments on spherical geometry". Are there any items in the reference manul or I have to define some new variable for this shape by myself?

My statement proposed a diffusional geometry that assumed a particular set
of boundary conditions, but apparently not what you have in mind. What
boundary conditions are you considering, that would require the geometry
that you propose?

[along82]

Hi, ted
You konw, it is only a single model for intracellular calcium dynamics description. I use hemi-sopherical geometry only because my cell mostly looks like that. So it seems a simple model and I really wish it is a simple one. I never think about any boundary conditions things. Actually, the only question is how to define this hemi-sopherical geometry, and then I can divide it into segments and shells, followed by incorporating calcium signal mechanisms.
Another problem:
I create a NMODL mechanism kdr and included it into my hoc file, however, it can't run through and the error messages like below:
at line 32 in file kdr.mod:
SOLVE states
Error at section location soma(0.5)
The increment in the independent variable is less than machine roundoff error
./i386/special: scopmath library error
It is strange since it is just fine when I include an IClamp stimuli, but never succeed for VClamp or SEClamp stimuli
my hoc file
load_file("nrngui.hoc")
2
3 // create sections
4 create soma
5 access soma
6
7 // set properties of sections
8 soma {
9 // Neuron Geometry
10 nseg = 1
11 diam = 20
12 L = 22.5
13 Ra = 123.0
14
15 // Delayed Rectifying Potassium Channel
16 insert kdr
24 }
25
26 // voltage-clamp stimuli for validation of voltage-gated channels
27 objectvar stim
28 soma stim = new SEClamp(0.5)
29
30 stim.dur1 = 100
31 stim.amp1 = -80
32 stim.dur2 = 100
33 stim.amp2 = 40
34 stim.dur3 = 100
35 stim.amp3 = -80
36 stim.rs = 1e-5 // actually 10 ohm
37
38 tstop = 800

my mod file:
TITLE Delayed Rectifying Potassium Channel

UNITS {
(mV) = (millivolt)
(mA) = (milliamp)
(S) = (siemens)
}

NEURON {
SUFFIX kdr
USEION k READ ek WRITE ik
RANGE g
}

PARAMETER {
g = 0.002 (S/cm2)
ek = -60
}

ASSIGNED {
v (mV)
ik (mA/cm2)
malpha (/ms)
mbeta (/ms)
}

STATE {
m
}

BREAKPOINT {
SOLVE states
ik = g*m*m*m*m*(v-ek)
}

INITIAL {
settables(v)
m = malpha/(malpha+mbeta)
}

DERIVATIVE states {
settables(v)
m' = ((malpha*(1-m)) - (mbeta*m))
}

UNITSOFF
PROCEDURE settables(v (mV)) {
TABLE malpha, mbeta FROM -100 TO 100 WITH 200

malpha = 10*(60-v)/(exp((60-v)/20)-1)
mbeta = 5*exp(-(v+60)/20)
}

UNITSON

[ted]

Lots of issues here, some deep, some not.


If you are serious about modeling ion accumulation and diffusion, it will serve
you well to learn how to approach such problems in a principled manner. Two
things you absolutely must know: how to take advantage of the geometry in
which accumulation and diffusion are happening, and how to deal with
boundary conditions. Chopping space into hundreds of little bits without
regard to geometry or boundary conditions is not the way to go. For a
hemispherical geometry in which solute exchange between cell and bath can
occur only through the curved surface, and there are no intracellular
complications that destroy radial symmetry, you may be able to represent
diffusion as exchange between nested hemispherical shells.

But don't just take my word for it. Try to find someone in your near
vicinity who can collaborate with you on the problem, or at least help
guide your readings. A good book is Crank's "The Mathematics of
Diffusion."


The mod file contains two items that should be revised. First, the units of
v should be specified as (mV) so modlunit won't complain. Second, the
formula for malpha becomes 0/0 at v = 60. This indeterminate condition
can be eliminated by using Lhospital's rule, i.e. replacing
malpha = 10*(60-v)/(exp((60-v)/20)-1)
with
malpha = 10*vtrap(60-v, 20)
and adding this FUNCTION

Code: Select all

FUNCTION vtrap(x,y) {  :Traps for 0 in denominator of rate eqns.
        if (fabs(x/y) < 1e-6) {
                vtrap = y*(1 - x/y/2)
        }else{
                vtrap = x/(exp(x/y) - 1)
        }
}

However, neither of these changes will eliminate the error message you
saw. To get rid of the error message, just don't make the SEClamp's
series resistance (rs) so small. Voltage control is perfectly adequate with
rs = 1e-3. Alternatively, use adaptive integration (CVode).


A caveat about your kdr mechanism: it's really a fast rectifier. Activates in
a millisecond or less.


Next, an initialization issue: default initialization with the standard run
system is a "voltage clamp initialization" to the value specified by the
variable v_init, i.e. as if all compartments were held for a long time at
v_init. The default value of v_init is -65 mV. But this doesn't make much
sense if your SEClamp's amp1 is -80 mV--i.e. why initialize the model to
one potential if, as soon as the simulation starts to execute, the SEClamp
is immediately going to force v to a different level? So decide if you want
to change v_init to -80, or if you would rather change amp1 to -65.


Minor stylistic items about postings that include code:
There's no need to include line numbers in code listings--destroys the
ability to cut and paste, so interferes with convenient reproducibility of
results.
If you encapsulate code in

Code: Select all

 . . . 

markup (select code, then
click on the "Code" button at the top of the "Message body" area), then
it will be displayed with indents preserved, which improves readability.

[along82]

Hi, ted
Thank you so much for your suggestions, I will learn some preliminary knowledge on the mathematics of diffusion before going on this model...
I think there is a mistake in your LHospital function

Code: Select all

FUNCTION vtrap(x,y) {  :Traps for 0 in denominator of rate eqns. 
        if (fabs(x/y) < 1e-6) { 
                vtrap = y*(1 - x/y/2) 
             }else{ 
                vtrap = x/(exp(x/y) - 1) 
        } 
}

vtrap should be equal to y rather than y*(...) when x/y=0

The rs parameter of SEClamp didn't fix my problem, but it finally work when I use cnexp Method to solve my states equations, just like below:

Code: Select all

BREAKPOINT {
          SOLVE states METHOD cnexp 
          ik  = g*m*m*m*m*(v-ek)
   }

[ted]

You're absolutely right--SOLVE states METHOD cnexp
makes it work properly. How did I overlook the absence of a METHOD
specifier in the SOLVE statement? How careless of me to miss that. And how
gracious of you to overlook what happens to y*(1 - x/y/2) as x/y -> 0.