Adding a BAPTA buffer...

[gartland]

hi
i realize this may be a simple question but i have been really banging my head over it. i am using a single compartment, five-channel model with ca diffusion and pump dynamics (capump.mod) that i found on ModelDB. i would like to be able to add a specific buffer to the model and keep track of the ammount of bound buffer and a new cai value. the kinetic eqn is this:

cai + free <--> bound (k1, k2)

i have included my best attempt at a MOD file below. should this work conceptually. i thought the easiest way to do this would be with a KINETICS block but i couldn't figrue out how to do this either. this code doesn't work. does anybody know why? thanks for the help!

-Andrew

Code: Select all

TITLE model of fluorescence changes for oregon green BAPTA-1

UNITS {
	(mv) = (millivolt)
	(mA) = (milliamp)
	(molar) = (1/liter)			
	(mM)	= (millimolar)
	(nM)	= (nanomolar)
	(um)	= (micron)
}

NEURON {
	SUFFIX dye
	USEION ca READ cai WRITE cai
	RANGE conc, kd, tau, bound
}

PARAMETER {
	conc = 0.0001	(mM)
	kd = 170		(nM)
	tau = 1			(ms)
}


ASSIGNED {
	:kd = k2/k1
	k2 (/ms)
	k1 (/nM /ms)
	ref
	
}

INITIAL {
	bound = kd
	free = conc-bound
}

STATE {
	cai		(mM)
	free	(mM)
	bound	(mM)
}

BREAKPOINT {
	SOLVE rxn METHOD euler
}

DERIVATIVE rxn {
	k1 = 1/(kd*tau)
	k2 = kd*k1

	cai'=-k1*cai*free + k2*bound
	free'=-k1*cai*free + k2*bound
	bound'=k1*cai*free - k2*bound
}

[Raj]

I have no context available to test your code in.

However I have a few remarks, which hopefully bring you one step further:
1. Running modlunit on your mod-file reveals errors in your units, which you should fix.
2. In your init procedure you use a conservation law: free = conc-bound. This conservation law will hold elsewhere as well and allows you to eliminate one of your differential equations, you can chose free' or bound'.
3. I only see kinetics, not a trace of diffusion. If you want include diffusion you might want to have a look at the nmodl documentation (P 34 and further): http://www.neuron.yale.edu/neuron/stati ... /nmodl.htm. The model described there includes longitudinal and radial diffusion for free calcium. It has an endogeneous buffer instead of a dye, which from a modeling perspective is pretty much the same except it does not light up and it seems to be immobile.

[gartland]

thank you for replying and offering to help me out.

1) i fixed the units errors. thanks for spotting that.

2) will a conservation law like the one in INITIAL hold throughout the simulation or does it only get called at the begining?

3) the diffusion and calcium pump is in a different mod file authored by Alain Destexhe. it seems to be the one that most people use for simple diffusion with a time constant.

4) when i run my hoc file there is a warning that two mechanisms are attempting to write to cai. it happens at the line where i "insert dye". how do multiple k-channels get away with both altering ik?

5) my idea is to not only model calcium concentration in the compartment (which i have done successfully with other peoples ca-channels and the diffusion above), but also to model the way in which i measure [ca] with a fluorescent indicator dye. also i would like to see how the exogenous buffer (the dye) affects the [ca]. do you think i have written a mod file to do this? how would you do it?

6) now my code runs and produces non sensical values. any ideas for improvement?

btw its very exciting to me that i can correspond with a scientist in amsterdam. forums like this are fantastic!

thank you,

-Andrew

Code: Select all

TITLE model of fluorescence changes for oregon green BAPTA-1

UNITS {
	(mv) = (millivolt)
	(mA) = (milliamp)
	(molar) = (1/liter)			
	(mM)	= (millimolar)
	(nM)	= (nanomolar)
	(um)	= (micron)
}

NEURON {
	SUFFIX dye
	USEION ca READ cai WRITE cai
	RANGE conc, kd, tau, bound, free
}

PARAMETER {
	conc = 0.1		(mM)
	kd = 0.000170	(mM)
	tau = 1			(ms)
}


ASSIGNED {
	:kd = k2/k1
	k2		(/ms)
	k1		(/mM /ms)
	free	(mM)
}

INITIAL {
	bound = kd
	free = conc - bound
}

STATE {
	cai		(mM)
	bound	(mM)
}

BREAKPOINT {
	SOLVE rxn METHOD euler
}

DERIVATIVE rxn {
	k1 = 1/(kd*tau)
	k2 = kd*k1

	cai'=-k1*cai*free + k2*bound
	bound'=k1*cai*free - k2*bound
}

[/code]

[Raj]

2) will a conservation law like the one in INITIAL hold throughout the simulation or does it only get called at the begining?

Depending on the experimental situation, the length of the experiment your modelling for instance, the conservation law will be a good approximation of the real situation. So if we ignore the filling phase and photo degradation of the dye the conservation law applies. In fact it was in your original code because in the derivative block free'=-bound'. So all dye in the free state that is being bound reappears as bound dye, locking the movement of free and bound together like communicating vessels. Of course in your present code you will have to enforce it by hand by repeating free = conc - bound in the derivative block. Alternatively you can replace every occurence of free with conc-bound.

3) the diffusion and calcium pump is in a different mod file authored by Alain Destexhe. it seems to be the one that most people use for simple diffusion with a time constant.

Because the reaction and diffusion are so closely linked I wouldn't separate them in to two separate mechanisms. I think that if the mechanisms as you envision it in your mod file would work the reaction would only take place in the outer shell. You will have to make sure you have the reaction taking place in all compartments.

4) when i run my hoc file there is a warning that two mechanisms are attempting to write to cai. it happens at the line where i "insert dye". how do multiple k-channels get away with both altering ik?

My best answer is by design, but admittedly I was not involved there I am only a user like you. Anyway the mechanism you want to write directly controls a state variable, if two mechanisms control the same variable which mechanism is right? With currents its slightly different, the currents in a mechanism are contributions to an overall current and the 'state' variables involved are the membrane potential and maybe an ion-concentration.

5) my idea is to not only model calcium concentration in the compartment (which i have done successfully with other peoples ca-channels and the diffusion above), but also to model the way in which i measure [ca] with a fluorescent indicator dye. also i would like to see how the exogenous buffer (the dye) affects the [ca]. do you think i have written a mod file to do this? how would you do it?

6) now my code runs and produces non sensical values. any ideas for improvement?

I think my comment to your third point covers this.

[Raj]

By the way you might want to go to the library and leaf through Crank's book `The Mathematics of Diffusion' to get a feeling for the numerical issues surrounding diffusion probblems. I would advice not to read it from cover to cover unless you have large quantities of time at hand.

[ted]

Crank's book is indeed excellent.

I hope these somewhat belated comments add more light than heat.

Remarks about accumulation mechanisms

Be sure to read the article about NMODL that Raj mentioned. With regard
to ion accumulation mechanisms, this is quite similar to the material in
chapter 9 of The NEURON Book.

Initialization of accumulation mechanisms presents a special challenge.
Incorrect initialization can lead to undesired results, such as slow shifts
of concentrations that require prolonged simulations for transients to die
away. Chapter 8 of The NEURON Book contains much useful information
about initialization of accumulation mechanisms that can make such
models much easier to implement and debug.

Kinetic scheme representation is most natural, i.e. easiest to conceptualize
and implement. DEs should be used only for the very simplest
accumulation mechanisms. Anything that involves radial diffusion,
multiple buffers, or organelles, will be much easier to express in terms
of kinetic schemes rather than DEs.

Multiple mechanisms can assert
USEION x . . . WRITE ix
(this ensures conservation of mass when there are multiple sources of
transmembrane ionic flux). However, only one mechanism is allowed to
WRITE the corresponding concentration xi or xo.

Integration methods

NEVER euler. NEVER. You will find legacy code with BREAKPOINT blocks
that assert
SOLVE foo METHOD euler
and they may seem to work, but this is asking for trouble. Just because
it's in ModelDB doesn't mean that it's the right thing to do. Doesn't matter
who it was who wrote it, they shouldn't have, and neither should you.
No new code should be written that uses euler.
Linear ODEs use cnexp.
Nonlinear ODEs and simple ion accumulation (a la cacum.mod in
examples) should use derivimplicit.
Kinetic schemes use sparse.

Where to look for examples of working NMODL code

1. Chapter 9 of The NEURON Book.

2. Hines, M.L. and Carnevale, N.T. Expanding NEURON's repertoire of
mechanisms with NMODL. Neural Computation 12:995-1007, 2000.
Downloadable from http://www.neuron.yale.edu/neuron/nrnpubs

3. c:\nrnxx\examples\nrniv\nmodl (MSWin) or the corresponding
examples/nrniv/nmodl that UNIX/Linux/OS X users can find in the
gzipped tar sources for NEURON

4. ModelDB, with the following caveat.

General caveat about ModelDB entries

This should not be taken as a criticism of any particular model mentioned
in this thread.

The presence of a model entry in ModelDB means two things: the model
was described in a published paper, and either (a) the authors provided
the source code they used or (b) somebody else developed a
computational implementation of their own.

It should always be kept in mind that the presence of a model entry in
ModelDB should not be regarded as an endorsement of the source code's
"correctness" or programming style. There are many examples of
excellent programming in ModelDB, but also many examples that are
mediocre or worse. Also, many that are otherwise excellent fall into the
category of "legacy code" because of the use of obsolete or deprecated
features of their respective programming languages or simulation
environments.