Setting the basal level for Intracellular Calcium Concentrat

[auatai]

Dear Sir,

I am working on a intracellular calcium diffusion, buffering and pump mechanism and am using the cdp.mod file (from the link provided at Errata for the Enhanced Preprint of "Expanding NEURON's Repertoire of Mechanisms with NMODL")

I have added this mechanism in the model "Pyramidal Neuron Deep, Superficial; Aspiny, Stellate (Mainen and Sejnowski 1996)"
Accession Number: 2488

which has other mod files likes ca.mod and cad.mod
Now the cad.mod file has a statement wherein cai is initialized to cainf which itself was initialized to 100nM
I
I now want to set the basal cai at 116 nano Molars. I tried doing it by changing the cainf parameter value in cad.mod file from 100nM to 116nM, however when i ran the simulation and checked the cai plot in the soma it was still reflecting the basal level at 50 nM. I don't know how the basal level is 50nM. And this does not change when i change the cainf value.

Can you please guide me about how to set the basal cai to 116 nM in this context.

thanks in advance

[ted]

The first thing to do is to make sure your model does not contain a conflict that prevents accurate results. cad.mod and cdp.mod should not be inserted into the same section because both contain
USEION ca . . . WRITE cai . . .
statements. This means that each of these mechanisms is using a different set of equations to compute a value for cai. The result will be nonsense. Decide which of these calcium accumulation mechanisms you want to use, and eliminate the other one.

The answer to your original question depends on which mechanism you decide to use.

[auatai]

Sir,
I have decided to use the cdp mechanism. Now, how to set the basal level of cai.

regards

[ted]

In principle there are several ways to do this. Most intellectually satisfying is to have an analytic solution that relates one parameter "x" of the model to all the others that might affect resting cai. Then it is a simple matter of plugging in numbers and grinding out the particular value for "x" that will ensure the desired resting cai.

But with any model cell that involves one or more voltage-gated calcium currents and this cdp mechanism, one must choose which parameter to adjust, and the parameter may not have the same value in all compartments (because other paramters that affect cai may not have the same values in all compartments). A compartment with high membrane calcium conductance/permeability will need some adjustment of the pump's forward rate constants to bring cai back to the desired resting level, but that will affect they dynamics of cai. The alternative is to reduce the membrane calcium conductance or permeability of that compartment, but that may interfere with the response of the cell to changes of membrane potential.

So what is to be done? Just leave all parameters as they are, and preface each simulation with a long "warmup" run that allows cai to asymptotically settle toward whatever level it may finally reach?

In such a case, a better approach may be to create a density mechanism that does just one thing: acts as a constant sink or source of calcium. One would then adjust, during initialization, the intensity of this sink or source so that it is equal but opposite in sign to the sum of all other calcium sources. The result is a model that has the desired resting calcium concentration, without having to perturb parameters that would disturb model dynamics.

This strategy falls into the broad category of what might be called "initialization to desired steady state." A similar approach, presented in chapter 8 of The NEURON Book, uses a distributed current source to achieve a desired resting potential.

Here's the NMODL source code for such a mechanism. The purpose of the NONSPECIFIC_CURRENT ix is to ensure that this "calcium source mechanism" does not affect charge balance--we only want it to affect the local calcium concentration.

Code: Select all

NEURON {
  SUFFIX casrc
  USEION ca WRITE ica
  NONSPECIFIC_CURRENT ix
  RANGE amp
}

UNITS {
  (mA)    = (milliamp)
}

PARAMETER {
  amp     (mA/cm2)
}

ASSIGNED {
  ica     (mA/cm2)
}

BREAKPOINT {
  ica = amp
  ix = -ica
}

Here's how you would use this mechanism.

At the end of model setup, after you have inserted all other mechanisms that WRITE ica into your model cell, execute this statement:

Code: Select all

forall if (ismembrane("ca_ion")) insert casrc // add casrc only to those sections that need it

After that, and certainly after any load_file("nrngui.hoc") statement, but before you execute any run() statement, insert this custom init() procedure:

Code: Select all

proc init() {
  ca0_ca_ion = whatever_you_want_it_to_be
  forall if (ismembrane("casrc")) amp_casrc = 0 // so it doesn't contribute to total ica(x)
  finitialize(v_init)
  forall if (ismembrane("ca_ion")) {
    // iterate over all internal nodes (segments) and adjust amp_casrc accordingly
    for (x,0) amp_casrc(x) = -ica(x)
  }
  if (cvode.active()) {
    cvode.re_init()
  } else {
    fcurrent()
  }
}

Now, every time you execute run(), NEURON will first set ca0_ca_ion (the default initial value of cai) to whatever value you want (this is important because the pump in cdp generates a current whose magnitude depends on cai), and then it will adjust amp_casrc in all segments of all sections with mechanisms that WRITE ica so that the net calcium flux into or out of each segment is 0.

This should all work, but it is always possible that I mistyped something or even (gasp) made a mistake, so be sure to test it on a simple model with just one or two sections first, then try it on your whole model and test again to make sure it works properly.

Thanks for using NEURON, and thanks also for asking a question that raises some very important issues about initializing models that contain ion accumulation mechanisms.

[auatai]

Dear Sir,

I am delighted to receive such a enlightening and elaborate reply. I can't thank you enough for showing me the way forward.

With warm regards

[auatai]

Dear Sir,

I tried initializing cai as suggested by you. Firstly the NMODL source code threw up an error on compilation "Translating casrc.mod into casrc.c
Warning: ix undefined. (declared within VERBATIM?) ....." . To resolve this I added the statement "ix (mA/cm2)" within the assigned block. (I hope it was the correct thing to do) . Anyway, now the mod file compiled properly.

Thereafter I added the "insert casrc" statement and then the init() procedure(in which i had set the parameter value of ca0_ca_ion = 116.0e-6). The init procedure was added into my hoc code just prior to the run() statement. I then ran the simulation but the result was still the same. The basal cai as seen in the cai plot was stuck at 50 nM.

Can you please help me out. I have copied the hoc code that i am using down below

Code: Select all

objref st,dendritic,dendritic_only
load_file("nrngui.hoc")
objectvar gr
create soma
access soma
tstop = 3000
steps_per_ms = 40
dt = 0.025
// ----------------passive & active membrane ---------------------
ra        = 150
global_ra = ra
rm        = 30000
c_m       = 0.75
cm_myelin = 0.04
g_pas_node = 0.02
v_init    = -70
celsius   = 37
Ek = -90
Ena = 60
gna_dend = 20
gna_node = 30000
gna_soma = gna_dend
gkv_axon = 2000
gkv_soma = 200
gca = .4
gkm = .1
gkca = 3
gca_soma = .7
gkm_soma = .1
gkca_soma = .3
// -----------Graph to observe cai-------------------------------
proc addgraph(){
	gr=new Graph()
	gr.size(0,tstop,-0.00001,0.0004)
	gr.xaxis()
	gr.yaxis()
	gr.addvar($s1,3,0)
	gr.addvar($s2,2,0)
	gr.save_name("graphList.")
	graphList.append(gr)
	}
// ------------Axon Geometry-------------------------------------
n_axon_seg = 5
create soma,iseg,hill,myelin[2],node[2]
proc create_axon() {
	create iseg,hill,myelin[n_axon_seg],node[n_axon_seg]
	soma{
  		equiv_diam = sqrt(area(.5)/(4*PI)) 	// area = equiv_diam^2*4*PI
  		}
	if (numarg()) equiv_diam = $1
	iseg{                					// initial segment between hillock + myelin
		L = 15
		nseg = 5
		diam = equiv_diam/10        		// see Sloper and Powell 1982, Fig.71
		}
	hill{                
		L = 10
        	nseg = 5
     	       diam(0:1) = 4*iseg.diam:iseg.diam
  		}
	// construct myelinated axon with nodes of ranvier
	for i=0,n_axon_seg-1 {
		myelin[i]{	         				// myelin element
			nseg = 5
			L = 100
			diam = iseg.diam         
			}
		node[i]{           					// nodes of Ranvier
			nseg = 1
			L = 1.0           
			diam = iseg.diam*.75    		// nodes are thinner than axon
         		}
		}
	soma connect hill(0), 0.5
	hill connect iseg(0), 1
	iseg connect myelin[0](0), 1
	myelin[0] connect node[0](0), 1
	for i=0,n_axon_seg-2  { 
		node[i] connect myelin[i+1](0), 1 
		myelin[i+1] connect node[i+1](0), 1
		}
	}
// ----------------------------Cell Initialization--------------------------------
proc init_cell(){
	// passive
	forall{
		insert pas
		Ra = ra 
		cm = c_m 
		g_pas = 1/rm
		e_pas = v_init
		}
	// exceptions along the axon
	forsec "myelin" cm = cm_myelin
	forsec "node" g_pas = g_pas_node
	// na+ channels
	forall insert na
	forsec dendritic gbar_na = gna_dend
	forsec "myelin" gbar_na = gna_dend
	hill.gbar_na = gna_node
	iseg.gbar_na = gna_node
	forsec "node" gbar_na = gna_node
	// kv delayed rectifier channels
	iseg { insert kv  gbar_kv = gkv_axon }
	hill { insert kv  gbar_kv = gkv_axon }
	soma { insert kv  gbar_kv = gkv_soma }
	// dendritic channels
	forsec dendritic_only{
	    insert km    gbar_km  = gkm
    	insert kca   gbar_kca = gkca
    	insert ca    gbar_ca = gca
    	insert cdp
		}
  	soma{
		insert km  gbar_km = gkm_soma
		insert kca gbar_kca = gkca_soma
		insert ca  gbar_ca = gca_soma
    	gbar_na = gna_soma
		insert cdp2
		}
	forall if(ismembrane("k_ion")) ek = Ek
	forall if(ismembrane("na_ion")){
		ena = Ena
		// seems to be necessary for 3d cells to shift Na kinetics -5 mV
		vshift_na = -5
		}
	forall if(ismembrane("ca_ion")){
    	eca = 140
	    ion_style("ca_ion",0,1,0,0,0)
		vshift_ca = 0
		}
	forall if(ismembrane("ca_ion")) insert casrc // add casrc only to those sections that need it
	}
// ----------------Init------------------------------------------
proc init(){
	ca0_ca_ion = 116.0
    forall if (ismembrane("casrc")) amp_casrc = 0 // so it doesn't contribute to total ica(x)
    finitialize(v_init)
    forall if (ismembrane("ca_ion")) {
    	// iterate over all internal nodes (segments) and adjust amp_casrc accordingly
        for (x,0) amp_casrc(x) = -ica(x)
    	}
    if (cvode.active()) {
    	cvode.re_init()
    	} else {
        fcurrent()
    	}
	}
// -------------------------------Load cell structure---------------------------
proc load_3dcell(){
	// $s1 filename
	forall delete_section()
	xopen($s1)
	access soma
	dendritic = new SectionList()
	// make sure no compartments exceed 50 uM length
	forall {
    	        diam_save = diam
    	        n = L/50
    	        nseg = n + 1
    	        if (n3d() == 0) diam = diam_save
    	        dendritic.append()
		}    
	dendritic_only = new SectionList()
	forsec dendritic dendritic_only.append()
	soma  dendritic_only.remove()
	create_axon()
	init_cell()
	init()
	}
// -------------------------Run Simulation------------------------
load_3dcell("cells/lcAS3.hoc")
st=new IClamp(.5)
st.dur = 1000
st.del = 100
st.amp = 0.082
nrnmainmenu()
addgraph("apath[4].cai(0.5)","soma.cai(0.5)")
nrncontrolmenu()

[ted]

Clearly you have a lot of ideas about what should go into your model (that's good), and you want to do it all at once. But the hoc code contains several actual or incipient bugs. Among other things, there are multiple create statements, scattered at different lines, such that if all of them were executed, the same sections would be created more than once, destroying previous effort. It is also too disorganized for efficient development and debugging. Throw it away and start from scratch, following the guidelines in this post
http://www.neuron.yale.edu/phpBB/viewto ... f=15&t=950
especially items 5, 7, and 10.

Do not bother adding a calcium accumulation mechanism such as cdp, or adding a custom initialization, until you have a well-structured program that creates a model whose other properties (anatomy, geometry, and biophysics) are "correct" (i.e. a close match to your conceptual model).

Once your complex model is at the point where it seems ready for insertion of a ca accumulation mechanism, put it aside and set up a toy model that has just one section and only two or three biophysical mechanisms--maybe pas, cdp, and a voltage-gated calcium current. Use this toy model as the testbed for your custom initialization strategy. After you have implemented and debugged your initialization strategy with the toy model, apply what you have learned to your complex model--and test again to make sure it is done correctly.

This may seem a difficult path to follow, but it is the way to master your own code, and after some initial delay your progress will be much smoother and faster.

[auatai]

Thanks Ted. I'll try and do what you have suggested.

with warm regards

[ted]

The "code cleanup" part of your task may be somewhat tedious but at least it should be straightforward. Proper initialization of ionic concentrations is one of the more difficult aspects of model specification and usage, which is why it should be tackled only after everything else is working correctly.

[auatai]

Dear Ted,
As suggested by you, i cleaned up the code after starting from scratch. Also as you had advised, i made a toy model with just the soma, into which i inserted the cdp,ca and pas mechanisms. I also inserted the casrc mechanism that you had provided.
This is the code for the toy model.

Code: Select all

load_file("nrngui.hoc")
create soma
soma      {nseg=2 diam=13.2 L= 8.2 }
//---------------- passive & active membrane properties--------------------
ra        = 150
rm        = 30000
c_m       = 0.75
v_init    = -70
celsius   = 37
Ek = -90
Ena = 60
gca_soma = .3
//---------------------Mechanisms-----------------------------------
proc init_cell() {
	soma{
    	insert pas
    	Ra = ra 
    	cm = c_m 
    	g_pas = 1/rm
    	e_pas = v_init
		insert ca
		insert cdp
		}
	forall if(ismembrane("ca_ion")) {
    	eca = 140
    	ion_style("ca_ion",0,1,0,0,0)
    	vshift_ca = 0
  		}
	forall if(ismembrane("ca_ion")) insert casrc // add casrc only to those sections that need it
	}
//------------------Init-------------------------------------------
proc init(){
	ca0_ca_ion = 116.0e-6
    forall if (ismembrane("casrc")) amp_casrc = 0 // so it doesn't contribute to total ica(x)
    finitialize(v_init)
    forall if (ismembrane("ca_ion")) {
    	// iterate over all internal nodes (segments) and adjust amp_casrc accordingly
        for (x,0) amp_casrc(x) = -ica(x)
    	}
    if (cvode.active()){
    	cvode.re_init()
    	} else {
        fcurrent()
    	}
	}
//---------------SIMULATION FLOW CONTROL---------------------------------
access soma
objref stim,shape,cgr,vgr
tstop = 3000
steps_per_ms = 40
dt = 0.025
stim=new IClamp(.5)
stim.dur = 1000
stim.del = 500
stim.amp = 0.082

init()
init_cell()

// show cell
shape = new PlotShape()
shape.size(-300,300,-300,300)
shape.show(2)

// show cai plot
cgr=new Graph()
cgr.size(0,tstop,-0.0001,0.0005)
cgr.xaxis() cgr.yaxis()
cgr.addvar("soma.cai(0.5)",2,0) 
cgr.save_name("graphList.")
graphList.append(cgr)

// show volt spikes
vgr=new Graph()
vgr.size(0,tstop,-90,50)
vgr.xaxis() vgr.yaxis()
vgr.addvar("soma.v(0.5)",1,0)
vgr.save_name("graphList.")
graphList.append(vgr)

nrncontrolmenu()

However, on running the code, I got the following error message

Code: Select all

oc>at line 69 in file cdp.mod:
  SOLVE state METHOD sparse
Error at section location soma(0.25)
Convergence not achieved in maximum number of iterations
/cad/neuron/nrn/x86_64/bin/nrniv: scopmath library error
 near line 0
 {run()}
        ^
        fadvance()
      advance()
    step()
  continuerun(3000)
and others

I guessed that this had something to do with the small size of the soma, so I resized it to the following dimensions (Diam=15, L=40). This time the code worked and I got the following output. (Was I right in my approach of resizing the soma or should i do something else?)

http://www.flickr.com/photos/48556469@N07/4447265946/

And this was the plot for cai

http://www.flickr.com/photos/48556469@N07/4447265950/
As seen above, the basal level for cai is still stuck at 50 nM. I am totally clueless as to where this figure of 50 nM is coming from. Is it part of some internal initialization within the neuron simulator?

I would highly appreciate it, if you could guide me on this issue.

For ready refernce, I have pasted the code for the cdp.mod and ca.mod mechanism that I have used
cdp.mod

Code: Select all

: Calcium ion accumulation with radial and longitudinal diffusion and pump

NEURON {
  SUFFIX cdp
  USEION ca READ cao, cai, ica WRITE cai, ica
  RANGE ica_pmp
  GLOBAL vrat, TotalBuffer, TotalPump
    : vrat must be GLOBAL--see INITIAL block
    : however TotalBuffer and TotalPump may be RANGE
}

DEFINE Nannuli 4

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

PARAMETER {
  DCa   = .6 (um2/ms)
  k1buf = 9000 (/mM-ms) : Yamada et al. 1989
  k2buf = 10 (/ms)
  TotalBuffer = .09  (mM)

  k1    = 1   (/mM-ms)
  k2    = 0.005   (/ms)
  k3    = 1       (/ms)
  k4    = 0.005   (/mM-ms)
  : to eliminate pump, set TotalPump to 0 in hoc
  TotalPump = 0.95e-11  (mol/cm2)
  fact = 0.75
}

ASSIGNED {
  diam      (um)
  ica       (mA/cm2)
  ica_pmp   (mA/cm2)
  ica_pmp_last   (mA/cm2)
  parea     (um)     : pump area per unit length
  cai       (mM)
  cao       (mM)
  vrat[Nannuli]  (1) : dimensionless
                     : 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
  Kd        (/mM)
  B0        (mM)
}

CONSTANT { volo = 1e10 (um2) }

STATE {
  : ca[0] is equivalent to cai
  : ca[] are very small, so specify absolute tolerance
  : let it be ~1.5 - 2 orders of magnitude smaller than baseline level
  ca[Nannuli]       (mM) <1e-7>
  CaBuffer[Nannuli] (mM) <1e-5>
  Buffer[Nannuli]   (mM) <1e-5>
  pump              (mol/cm2) <1e-15>
  pumpca            (mol/cm2) <1e-15>
}

BREAKPOINT {
  SOLVE state METHOD sparse
  ica_pmp_last = ica_pmp
  ica = ica_pmp
}

LOCAL factors_done

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
   }

  Kd = k1buf/k2buf
  B0 = TotalBuffer/(1 + Kd*cai)

  FROM i=0 TO Nannuli-1 {
    ca[i] = cai
    Buffer[i] = B0
    CaBuffer[i] = TotalBuffer - B0
  }

  parea = PI*diam

: Manually computed initalization of pump
: assumes that cai has been equal to cai0_ca_ion for a long time
:  pump = TotalPump/(1 + (cai*k1/k2))
:  pumpca = TotalPump - pump
: If possible, instead of using formulas to calculate pump and pumpca,
: let NEURON figure them out--just uncomment the following four statements
  ica=0
  ica_pmp = 0
  ica_pmp_last = 0
  SOLVE state STEADYSTATE sparse
: This requires that pump and pumpca be constrained by the CONSERVE
: statement in the STATE block.
: If there is a voltage-gated calcium current, 
: this is almost certainly the wrong initialization. 
: In such a case, first do an initialization run, then use SaveState
: On subsequent runs, restore the initial condition from the saved states.
}

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] {ca CaBuffer Buffer}
  COMPARTMENT (1e10)*parea {pump pumpca}
  COMPARTMENT volo {cao}
  LONGITUDINAL_DIFFUSION i, DCa*diam*diam*vrat[i] {ca}

  :pump
  ~ ca[0] + pump <-> pumpca  (k1*parea*(1e10), k2*parea*(1e10))
  ~ pumpca <-> pump + cao    (k3*parea*(1e10), k4*parea*(1e10))
  CONSERVE pump + pumpca = TotalPump * parea * (1e10)
  ica_pmp = 2*FARADAY*(f_flux - b_flux)/parea

  : all currents except pump
  : ica is Ca efflux
  ~ ca[0] << (-(ica - ica_pmp_last)*PI*fact*diam/(2*FARADAY))
  FROM i=0 TO Nannuli-2 {
    ~ ca[i] <-> ca[i+1]  (DCa*frat[i+1], DCa*frat[i+1])
  }
  dsq = diam*diam
  FROM i=0 TO Nannuli-1 {
    dsqvol = dsq*vrat[i]
    ~ ca[i] + Buffer[i] <-> CaBuffer[i]  (k1buf*dsqvol, k2buf*dsqvol)
  }
  cai = ca[0]
}

ca.mod

Code: Select all

COMMENT
26 Ago 2002 Modification of original channel to allow variable time step and to correct an initialization error.
    Done by Michael Hines(michael.hines@yale.e) and Ruggero Scorcioni(rscorcio@gmu.edu) at EU Advance Course in Computational Neuroscience. Obidos, Portugal

ca.mod
Uses fixed eca instead of GHK eqn

HVA Ca current
Based on Reuveni, Friedman, Amitai and Gutnick (1993) J. Neurosci. 13:
4609-4621.

Author: Zach Mainen, Salk Institute, 1994, zach@salk.edu

ENDCOMMENT

INDEPENDENT {t FROM 0 TO 1 WITH 1 (ms)}

NEURON {
	SUFFIX ca
	USEION ca READ eca WRITE ica
	RANGE m, h, gca, gbar
	RANGE minf, hinf, mtau, htau
	GLOBAL q10, temp, tadj, vmin, vmax, vshift
}

PARAMETER {
	gbar = 0.1   	(pS/um2)	: 0.12 mho/cm2
	vshift = 0	(mV)		: voltage shift (affects all)

	cao  = 2.5	(mM)	        : external ca concentration
	cai		(mM)
						
	temp = 23	(degC)		: original temp 
	q10  = 2.3			: temperature sensitivity

	v 		(mV)
	dt		(ms)
	celsius		(degC)
	vmin = -120	(mV)
	vmax = 100	(mV)
}


UNITS {
	(mA) = (milliamp)
	(mV) = (millivolt)
	(pS) = (picosiemens)
	(um) = (micron)
	FARADAY = (faraday) (coulomb)
	R = (k-mole) (joule/degC)
	PI	= (pi) (1)
} 

ASSIGNED {
	ica 		(mA/cm2)
	gca		(pS/um2)
	eca		(mV)
	minf 		hinf
	mtau (ms)	htau (ms)
	tadj
}
 

STATE { m h }

INITIAL { 
	trates(v+vshift)
	m = minf
	h = hinf
}

BREAKPOINT {
        SOLVE states METHOD cnexp
        gca = tadj*gbar*m*m*h
	ica = (1e-4) * gca * (v - eca)
} 

LOCAL mexp, hexp

:PROCEDURE states() {
:        trates(v+vshift)      
:        m = m + mexp*(minf-m)
:        h = h + hexp*(hinf-h)
:	VERBATIM
:	return 0;
:	ENDVERBATIM
:}

DERIVATIVE states {
        trates(v+vshift)      
        m' =  (minf-m)/mtau
        h' =  (hinf-h)/htau
}

PROCEDURE trates(v) {  
                      
        
        TABLE minf, hinf, mtau, htau 
	DEPEND  celsius, temp
	
	FROM vmin TO vmax WITH 199

	rates(v): not consistently executed from here if usetable == 1

:        tinc = -dt * tadj

:        mexp = 1 - exp(tinc/mtau)
:        hexp = 1 - exp(tinc/htau)
}


PROCEDURE rates(vm) {  
        LOCAL  a, b

        tadj = q10^((celsius - temp)/10)

	a = 0.055*(-27 - vm)/(exp((-27-vm)/3.8) - 1)
	b = 0.94*exp((-75-vm)/17)
	
	mtau = 1/tadj/(a+b)
	minf = a/(a+b)

		:"h" inactivation 

	a = 0.000457*exp((-13-vm)/50)
	b = 0.0065/(exp((-vm-15)/28) + 1)

	htau = 1/tadj/(a+b)
	hinf = a/(a+b)
}

FUNCTION efun(z) {
	if (fabs(z) < 1e-4) {
		efun = 1 - z/2
	}else{
		efun = z/(exp(z) - 1)
	}
}

[ted]

The hoc code is much cleaner. One minor problem is that proc init_cell(), which is really part of the "specification of biological properties" (it inserts mechanisms into sections!), is not executed until somewhere in the middle of a bunch of simulation parameter, simulation flow control, and user interface code. init_cell() should really be called before the custom proc init() is defined.

Another is that init_cell() calls ion_style() with a specification of how Ca concentration, reversal potential, etc. that interferes with proper initialization (it tries to force ena to be a constant and prevent Ca concentrations from changing!). This should just be omitted for now. If it turns out that you need eca to have a constant value, let's take care of that later, after Ca concentration is handled properly.

A potentially bigger problem is of my creation: a typographical error in the custom proc init(), which assigns a value to ca0_ca_ion. That should have been cai0_ca_ion.

And the use of an even value for nseg might cause you some difficulties in the future (see
Why should I use an odd value for nseg?
in the FAQ list
http://www.neuron.yale.edu/neuron/faq/general-questions
).

However, after fixing all of these issues, I see that the initial value of cai remains unaffected. Re-reading the code in cdp.mod, especially its INITIAL block, and the related discussion on page 258 of The NEURON block, I wonder if a more complex initialization process is required. It would be easy to change the pump's rate constants to force its "resting cai" to be 116 nM--for any particular resting cai and cao, the pump generates a net flux of 0 when
kf*cap = kr*cao*p
kr*cap = kf*cai*p
where numerical value of k1 = numerical value of k3 = kf
and numerical value of k2 = numerical value of k4 = kr
from which we have
kf = kr*sqrt(cao/cai)
But if one wants to leave the rate constants unchanged and use the membrane's ion channels plus a "constant ca source" to force resting cai to a desired value, the question is how to find the required intensity of the "constant ca source."

I'm going to have to think about this some more.

[auatai]

Thank you Ted for sparing the time to help me out.

Meanwhile, I'll work on the guidelines given by you.

[auatai]

Dear Sir,

Were you able to think out a way of initializing the cai? I tried doing it by adjusting the pump parameters but then the behavior of the model changed a lot from what I wanted it to be.

Regards

[auatai]

Dear Sir,

But if one wants to leave the rate constants unchanged and use the membrane's ion channels plus a "constant ca source" to force resting cai to a desired value, the question is how to find the required intensity of the "constant ca source."

I'm going to have to think about this some more.

At the outset I apologise if I am being annoying. But i need to set the basal cai and i am just not able to do so. I would be greatly obliged if you could help me find a way.

thanks and regards

[ted]

You're not being annoying. The problem has up to now stymied me. However, I finally think I see a way forward.

There are actually three difficulties to overcome. The first is to get the cdp mechanism to allow cai to be set to a value that differs from the one at which the net pump current is 0; this I have done, but only after many false starts. The second is to properly initialize the buffer concentrations, levels of pump and pumpca, and the current generated by the pump itself; this is taking more cleverness than I expected, partly because the pump's rate constants do not quite satisfy the Michaelis-Menten assumptions (the result is that pump + cai <-> pumpca doesn't replenish pumpca quite as rapidly as pumpca <-> pump + cao depletes it, and this causes a noticeable transient in the initial few ms of operation of the pump following initialization). However, I have an idea that I believe will resolve this problem, and will try it tomorrow. The third difficulty is how to deal with a model that contains other sources of transmembrane ca flux. I think this last can be handled by modifying the initialization to include a brief "initialization run," and maybe taking advantage of SaveState.