Nan with LONGITUDINAL_DIFFUSION

[gartland]

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:

Code: Select all

LONGITUDINAL_DIFFUSION j, DCa*diam*diam*vrat[j] {caShells caB1 caB2 B1free B2free}

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!

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]
}

[ted]

Start over. After each change, check with modlunit. If it reports an error, fix the error before going further. If you run into a wall at some point, let me know.

[hines]

I recently made a bug fix to the LONGITUDINAL_DIFFUSION code related to changing diameter in Version 7.0.
Could you send me your hoc,ses,mod files in a zip file that exhibits the NAN but works correctly without
LONG... I'd like to see what is going wrong. Send the zip file to michael.hines@yale.edu.

[gartland]

Unfortunately I can only offer you the mechanism above since the hoc and ses pieces are all in python. But, before I waste anyone else's time I will go through and eliminate all the other parts of the model and see if I can get *just* longitudinal diffusion to work. The LONGITUDINAL_DIFFUSION statement above looks alright though? Thank you.

[ted]

Looks OK. modlunit finds units inconsistencies elsewhere and reports the lack of a COMPARTMENT statement for
~ caShells[0] <-> cainf (gamma*dsq*vrat[0],gamma*dsq*vrat[0])
However, other errors may also be present that modlunit cannot "see" because these gagged it first. That's why it makes sense to go back to an earlier, simpler mechanism that works and passes modlunit's tests, and add further complications, one at a time, testing at each increment of new complexity.