NEURON {
	 SUFFIX ip3dif
 	 USEION ip3 READ iip3 WRITE ip3i VALENCE 1
  	 GLOBAL vrat, DIP3, ip3i0
}

DEFINE Nannuli 4

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

PARAMETER {
  	kdegr = 0.14e-3 (/ms)  : degredation rate Finck et al & wang ey al 1995
  	DIP3 = 0.283(um2/ms)
  	ip3i0 = 0.16e-3 (mM)  : [IP3]0  initial and resting ip3i conc
}


ASSIGNED {
  	diam      (um)
 	ip3i      (mM)
	iip3      (mA/cm2)
  	vrat[Nannuli]  : numeric value of vrat[i] equals the volume 
                 : of annulus i of a 1um diameter cylinder
                 : multiply by diam^2 to get volume per um length
}

STATE {
  	: ip3[0] is equivalent to ip3i
  	: ip3[] are very small, so specify absolute tolerance
  	ip3[Nannuli]       (mM) <1e-6>
}

LOCAL factors_done

BREAKPOINT { SOLVE state METHOD sparse}


INITIAL {
   	if (factors_done == 0) {  : flag becomes 1 in the first segment
     	factors_done = 1       :   all subsequent segments will have
      	factors()              :   vrat = 0 unless vrat is GLOBAL
       }

  	ip3i = ip3i0
  	FROM i=0 TO Nannuli-1 {
    	ip3[i] = ip3i
	}
}


LOCAL frat[Nannuli]  : scales the rate constants for model geometry

PROCEDURE factors() {
  	LOCAL r, dr2
  	r = 1/2                : starts at edge (half diam)
  	dr2 = r/(Nannuli-1)/2  : full thickness of outermost annulus,
                         : half thickness of all other annuli
  	vrat[0] = 0
  	frat[0] = 2*r
 	 FROM i=0 TO Nannuli-2 {
    		vrat[i] = vrat[i] + PI*(r-dr2/2)*2*dr2  : interior half
    		r = r - dr2
    		frat[i+1] = 2*PI*r/(2*dr2)  : outer radius of annulus
                                : div by distance between centers
   		 r = r - dr2
    		vrat[i+1] = PI*(r+dr2/2)*2*dr2  : outer half of annulus
  	}
}

LOCAL dsq, dsqvol  : can't define local variable in KINETIC block
                   :   or use in COMPARTMENT statement


KINETIC state {
  	COMPARTMENT i, diam*diam*vrat[i] {ip3 ip3i0}
  	LONGITUDINAL_DIFFUSION i, DIP3*diam*diam*vrat[i] {ip3}
  	~ ip3[0] << (iip3*PI*diam*(1e4)/FARADAY)

  	FROM i=0 TO Nannuli-2 {
   		 ~ ip3[i] <-> ip3[i+1]  (DIP3*frat[i+1], DIP3*frat[i+1])
  	}

  	dsq = diam*diam
  	FROM i=0 TO Nannuli-1 {
    		dsqvol = dsq*vrat[i]
    		~ ip3[i] <-> ip3i0  (kdegr*dsqvol, kdegr*dsqvol)
  	}

  	ip3i = ip3[0]
}