Nan with LONGITUDINAL_DIFFUSION
Posted: Thu Dec 11, 2008 5:59 pm
I've been working on a calcium diffusion model that is very much like the one in Chapter 9 (cadifusl.mod). I have modified it by:
1) Changing all the variable names :) (sorry this happened along the way so that i could internalize it better...)
2) Adding a second calcium buffer (useful when I am modeling calcium buffering before and after a calcium indicator like Fluo-4 has been added)
3) Adding an initialization block to help find equilibrium
4) Making most of the parameters accessible from hoc
5) Adding diffusion mechanisms for the buffers and having their diffusion constant be a range variable
6) Chaning the extrusion mechanism so that the gamma parameter is the same as the gamma parameter in the single compartment model specified in Yasuda et al. (2004, Science STKE)
You'd think that just by doing these things I wouldn't break the model, but I have done it. I have it narrowed down to the line implementing longitudinal diffusion:
If the mechanisms is inserted in a model with a single section with a single segment, then it works well and passes all of my tests that show all the buffers are diffusing radially and binding beautifully. However, if i either add a section or add a segment then all the values in caShells (the calcium state variable array) are Nan except the first one. The same is true of the voltage and anything else I can "record". If I leave the multicompartment model intact and simply comment out the LONGITUDINAL_DIFFUSION line in the mechanism then the model again works wonderfully, except there's no longitudinal diffusion of course.
I've pasted in the code for the entire mechanism below. I wish I didn't have to ask someone to debug this for me, but I don't really have anything to go on. Any ideas? What can cause a model to spit out Nan's? Your help is appreciated!
1) Changing all the variable names :) (sorry this happened along the way so that i could internalize it better...)
2) Adding a second calcium buffer (useful when I am modeling calcium buffering before and after a calcium indicator like Fluo-4 has been added)
3) Adding an initialization block to help find equilibrium
4) Making most of the parameters accessible from hoc
5) Adding diffusion mechanisms for the buffers and having their diffusion constant be a range variable
6) Chaning the extrusion mechanism so that the gamma parameter is the same as the gamma parameter in the single compartment model specified in Yasuda et al. (2004, Science STKE)
You'd think that just by doing these things I wouldn't break the model, but I have done it. I have it narrowed down to the line implementing longitudinal diffusion:
Code: Select all
LONGITUDINAL_DIFFUSION j, DCa*diam*diam*vrat[j] {caShells caB1 caB2 B1free B2free}I've pasted in the code for the entire mechanism below. I wish I didn't have to ask someone to debug this for me, but I don't really have anything to go on. Any ideas? What can cause a model to spit out Nan's? Your help is appreciated!
Code: Select all
NEURON {
SUFFIX cad
USEION ca READ cai, ica WRITE cai
RANGE vrat, cainf
RANGE B1tot, B2tot, gamma
RANGE B1kd, B1kon, B2kd, B2kon
RANGE DCa, DB1, DB2
}
DEFINE nshells 4
UNITS {
(molar) = (1/liter)
(uM) = (micromolar)
(mM) = (millimolar)
(um) = (micron)
(mA) = (milliamp)
FARADAY = (faraday) (10000 coulomb)
PI = (pi) (1)
}
PARAMETER {
DCa = 0.22 (um2/ms) :0.1 is same as 1e-6 cm^2/s
DB1 = 0.0 (um2/ms) :
DB2 = 0.22 (um2/ms) :
B1kon = 1000 (/mM-ms) :same as 1e9 M-1 * s-1
B1kd = 0.003 (mM) :same as 2.3 uM
B1tot = 0.13 (mM) :0.1
B2kon = 1000 (/mM-ms) :same as 1e9 M-1 * s-1
B2kd = 0.0013 (mM) :same as 2.3 uM
B2tot = 0.0 (mM)
gamma = 300 (1/sec)
cainf = 0.0001 (mM)
}
ASSIGNED {
diam (um)
ica (mA/cm2)
cai (mM)
vrat[nshells] (1) :dimensionless volume ratio. volume of shell i of a 1 um diameter cylinder. multiply by diam^2 to get volume per um length
B1koff (/ms) :should always be B1kon * B1kd
B2koff (/ms) :should always be B2kon * B2kd
}
STATE {
caShells[nshells] (mM) <1e-10>
caB1[nshells] (mM)
B1free[nshells] (mM)
caB2[nshells] (mM)
B2free[nshells] (mM)
}
BREAKPOINT {
SOLVE state METHOD sparse
}
LOCAL factors_done
INITIAL {
if (factors_done == 0) {
factors_done = 1
factors()
}
B1koff = B1kon * B1kd
B2koff = B2kon * B2kd
FROM i=0 TO nshells-1 {
caShells[i] = cai
B1free[i] = B1tot/(1+(1/B1kd)*cai)
B2free[i] = B2tot/(1+(1/B2kd)*cai)
caB1[i] = B1tot-B1free[i]
caB2[i] = B2tot-B2free[i]
}
ica=0
SOLVE state STEADYSTATE sparse
}
LOCAL frat[nshells]
PROCEDURE factors() {
LOCAL r, dr2
r = 1/2
dr2 = r/(nshells-1)/2
vrat[0] = 0
frat[0] = 2*r
FROM i=0 TO nshells-2 {
vrat[i] = vrat[i] + PI*(r-dr2/2)*2*dr2
r = r - dr2
frat[i+1] = 2*PI*r/(2*dr2)
r = r - dr2
vrat[i+1] = PI*(r+dr2/2)*2*dr2
}
}
LOCAL dsq, dsqvol
KINETIC state {
COMPARTMENT i, diam*diam*vrat[i] {caShells caB1 caB2 B1free B2free}
:this model works as long as this line is commented out, but there is no longitudinal diffusion!
:LONGITUDINAL_DIFFUSION j, DCa*diam*diam*vrat[j] {caShells caB1 caB2 B1free B2free}
dsq = diam*diam
~ caShells[0] <-> cainf (gamma*dsq*vrat[0],gamma*dsq*vrat[0])
~ caShells[0] << (-ica*PI*diam*frat[0]/(2*FARADAY)) :ica is calcium efflux
FROM i=0 TO nshells-2 {
~ caShells[i] <-> caShells[i+1] (DCa*frat[i+1], DCa*frat[i+1])
~ B1free[i] <-> B1free[i+1] (DB1*frat[i+1], DB1*frat[i+1])
~ caB1[i] <-> caB1[i+1] (DB1*frat[i+1], DB1*frat[i+1])
~ B2free[i] <-> B2free[i+1] (DB2*frat[i+1], DB2*frat[i+1])
~ caB2[i] <-> caB2[i+1] (DB2*frat[i+1], DB2*frat[i+1])
}
FROM i=0 TO nshells-1 {
dsqvol = dsq*vrat[i]
~ caShells[i] + B1free[i] <-> caB1[i] (B1kon*dsqvol, B1koff*dsqvol)
~ caShells[i] + B2free[i] <-> caB2[i] (B2kon*dsqvol, B2koff*dsqvol)
}
cai=caShells[0]
}