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
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.