Spatial discretization for radial diffusion of Ca

[Darshan]

Hi,

I am using a calcium accumulation mechanism with shells for modelling the calcium dynamics in a cell. The diameter of the soma is different from that of the dendrites and the axon. This implies that the shell widths are different for all the sections. If there are calcium activated mechanisms like BK channels on the cell membrane then the calcium seen by them will be much higher in case of axons and dendrites and lower for soma because of different shell sizes. Smaller the width of the shell less is the volume and thus, more calcium change/accumulation in the shell.
Is it then correct to use the calcium at the outermost shell to activate these mechanisms or one should use a fixed depth (e.g. 0.1 micrometre) and find the corresponding shell for each segment and use the calcium accumulated in that shell as the activating calcium and as the calcium transient generated in the cell? Changing BK channel conductance is another option for decreasing BK currents which are generated by very high calcium accumulation.


Regards,
Darshan

[ted]

Good question.

The d_lambda rule is a practical way to design spatial discretization of the cable equation for the longitudinal spread of electrical signals in neurons, but there isn't a similar rule to guide discretizing a model of radial diffusion. The issue can only be settled empirically, i.e. through computational experimentation. Here's how I'd go about it.

I'd start with a very simple model cell--call it "model 0"--that is just a single section with these properties:
L your choice, but 10 um is sufficient
diam will be an independent variable
Ra unimportant
cm 1 uf/cm2
insert the v-gated ca channel of your choice, with the density of your choice
insert the BK channel mechanism of your choice, with the density of your choice
attach an SEClamp with rs = 1e-3 megohms, dur1 = 1 ms, amp1 = your resting potential, dur2 = 1 ms, amp2 = your resting potential + 100 mV (to emulate a spike), dur3 = 1e9 ms, amp3 = your resting potential

To this I would then insert a ca accumulation mechanism. I'd start with the crude 4 compartment version, in which the thickness of each shell is a fixed fraction of section radius. This should have buffering, radial diffusion, and maybe active transmembrane extrusion of ca. No SERCA or other complications. I'd call the resulting model cell "model 1".

I'd record and plot the following vs. time:
v
cai, which should be the calcium concentration in the outermost shell
the conductance of the BK channel mechanism
the current generated by the BK mechanism

I'd do this for cell diameters of 0.2, 1, 5, and 25 um.

These simulations don't have to go on forever--just long enough that the decaying phases of their cai transients have converged to within, say, about 5-10% of each other.

Then I'd make a new calcium accumulation mechanism, based on the crude 4 compartment version, but with only 2 shells (NANN = 2). I'd give this mechanism and its NMODL source code a distinctive name to avoid confusing it with the original. If the original file and SUFFIX are cadif.mod and cadifus, the new names would be cadif2.mod and cadifus2.

I'd then make a new model cell based on model 0, but change diam to 0.2 um, insert cadifus2, call the result model 2, and use it to generate plots of v, cai, the BK conductance, and the BK-generated current, then compare those with the results from the original 4 compartment model.

I'd repeat this in a model cell called model 10 that has diam = 1 um, using a ca accumulation mechanism that has 10 shells. Then model 50 with diam = 5 um and NANN = 50, and finally model 250 with diam = 25 um and NANN = 250. Each of these models will have shells that are 0.1 um thick.

Of course there will be differences between the results produced by the 4 compartment model and these "0.1 um shell models"--especially during the early part of the cai transient. The question is whether those differences are significant. It's the time course of the BK conductance that matters, much more than the time course of cai, because the BK conductance produces a current that affects the time course of v in an unclamped model cell. So what's "significant"? Depends on the function your BK current is serving in the particular cell of interest to you, and how big it is compared to the other currents that contribute to that function. That's hard, if not impossible, to judge in a model as simple as the one proposed above. However, BK tends to participate in phenomena that occur over the time course of many ms, so the integral of BK current over many ms (10? 50? 100?) should serve as a rough indication of its contribution to membrane potential over that time interval. This integral will probaby be much less dependent on radial discretization than the detailed time course of BK current.

To count as significant, differences would have to be larger than the experimental uncertainty in BK channel density. Most ion channel densities vary from location to location ‌in a given cell, not to mention from cell to cell even within the same class of cell, let alone from report to report in the experimental literature. Uncertainty on the order of +-20% is common if not the rule.

So I wouldn't be too concerned if two models that use different discretization strategies produce BK currents whose integrals differ by 30 - 40%.

[Darshan]

Thank you Ted for the quick and detailed explanation.

I was trying to implement equal shell widths for the models. I am using the cabpump.mod file which has a membrane pump, radial diffusion and buffering.

I have used this formula for calculating the volume of a shell: PI*(r^2-rc^2)*1 where rc is the radius of the core which remains below the shell being considered. The length of cylinder was taken as 1 um. The code is

Code: Select all

PROCEDURE coord() {
	LOCAL r, dr2
	r = 1/2				:starts at edge (half diam)
	dr2 = r/(NANN)			:full thickness of annulus
	vol[0] = 0
	frat[0] = 2*r
	FROM i=0 TO NANN-1 {
		vol[i] = PI*(r^2-(r-dr2)^2)*1	:volume of the shell = PI*(outer_radius^2-inner_radius^2)*1(um)
		frat[i] = 2*PI*r/(dr2)			:exterior edge of annulus divided by distance between centers
		r = r - dr2
	}
}

Is it correct? In the original cabpump file, a slightly different approach is used for calculating the volume of shells. The volume of the interior half and exterior half of a shell are calculated separately. And it uses a different formula for calculation which is not clear to me. The volume is calculated for a unit length cylinder for a given shell width according to the NEURON book.

Original code for shell vol. calculation:

Code: Select all

FROM i=0 TO NANN-2 {
		vol[i] = vol[i] + PI*(r-dr2/2)*2*dr2	:interior half
		r = r - dr2
		frat[i+1] = 2*PI*r/(2*dr2)	:exterior edge of annulus
					               : divided by distance between centers
		r = r - dr2
		vol[i+1] = PI*(r+dr2/2)*2*dr2	:outer half of annulus
	}

I simulated the models with the modified shell widths and I can see the differences in the cai for models 1 (4 compartment model), 2, 10, 50 and 250 (models 2,5,10 and 250 are 0.1 um models). The difference in peak cai is evident in all the 0.1 um models. But the decay phase of models 2 and 10 0.1 um models show some similarity when compared with original model 1 (which I think is mostly because of the extrusion mechanism).

I used nseg = 5 (odd) instead of 4 (Though, it doesn't play any role in the changing shell cai (?)).

Regards,
Darshan

[ted]

I was trying to implement equal shell widths for the models.

A natural first thought, and that would be fine if concentration in the shells was calculated with a first-order correct method. However, the spatial discretization scheme implemented in nrn/demo/release/cabpump.mod and in 9.10.1 Modeling diffusion with kinetic schemes uses the central difference approximation for the second spatial derivative of calcium concentration in the radial direction. Consequently,

"The outermost shell is half as thick as the others so that [Ca2+] will be second-order correct with respect to space at the surface of the segment."

(page 246 of The NEURON Book, paragraph 2, sentence 3). That's what we want because the BK channels are in the membrane, not in the cytoplasm. Making all shells have equal thickness will introduce an error because the value assigned to cai will have been calculated at the "wrong" location in space.

I simulated the models with the modified shell widths and I can see the differences in the cai for models 1 (4 compartment model), 2, 10, 50 and 250 (models 2,5,10 and 250 are 0.1 um models). The difference in peak cai is evident in all the 0.1 um models.

That makes sense. At what time does it occur?

the decay phase of models 2 and 10 0.1 um models show some similarity when compared with original model 1 (which I think is mostly because of the extrusion mechanism).

Yes, the decay time course of cai is partly radial diffusion but the slowest component of decay will be due to the pump.

For a semiquantitative indicator of similarity, the most relevant results to compare between different mechanisms would be time course and integral of BK current over a phenomenologically-relevant time interval--integral would probably be best if the primary phenomenon that BK current affects has a slow time course, i.e. many ms.

But the ultimate test is what happens when one of these ca accumulation mechanisms is used in a more "complete" model, that is, one that has all of the other active and passive currents. If BK current is always small compared to the other currents, then it won't matter which of these ca accumulation mechanisms is used. If the model fires spikes with long interspike intervals, or bursts with long interburst intervals, the ISI or interburst interval may be sensitive to which model is used. In any case, the ultimate question is whether using one mechanism or another causes a significant qualitative change in model behavior. If the answer is no, then my inclination would be to use the simplest mechanism because it requires the fewest assumptions and is computationally efficient.

I used nseg = 5 (odd) instead of 4

nseg should be 1 for a 1 compartment model. Perhaps you mean the parameter that specifies the number of shells, which is NANN in cabpump.mod and Nannuli in the NEURON book example.

[Darshan]

Hi Ted,

Thanks for reply. I will get back to you once I have implemented and tested in the 'complete model' along with integral of BK current and their time course.

That makes sense. At what time does it occur?

The peak times with the models I used above fall within 1 ms of each other. I will test the models after doing the corrections.

nseg should be 1 for a 1 compartment model.

I misunderstood your statement in the first reply:

I'd start with the crude 4 compartment version.

I thought you were referring to nseg and not NANN.

Regards,
Darshan

[Darshan]

Hi,
I have a 1-D shell calcium dynamics model with intracellular release and uptake mechanisms. I wanted to implement the volume distribution for cytoplasm and ER similar to what has been done in rxd.This can give me changes in the cytoplasm and ER calcium for different simulations.
My model assumes that each shell’s 83% volume is occupied by cytoplasm and 12% by ER. I have implemented this by using the following code for the compartment statement of KINETIC block for my calcium dynamics mod file:

Code: Select all

KINETIC state {
	COMPARTMENT ii, (1+bbr)*diam*diam*vol[ii]*0.83 {ca} :Cytoplasmic volume for each shell, bbr is buffer binding ratio, vol[ii] are volumes of shells :(as implemented in cabpump.mod)
	COMPARTMENT jj, diam*diam*vol[jj]*0.17 {caer}         :ER volume for each shell 

caer is the calcium concentration of ER initialized to a resting value.

For one of the mechanisms, SERCA, I have the following code (in the KINETIC block):

Code: Select all

dsq = diam*diam      
	 
FROM i=0 TO NANN-1{
	 dsqvol = dsq*vol[i]*0.83
	 dsqvoler = dsq*vol[i]*0.17            

	: SERCA pump
	jserca[i] = (-vmaxsr*ca[i]^2 / (ca[i]^2 + kpsr^2)) : jserca is the flux, vmaxsr is the maximum rate of release, Kpsr is a constant
	~ ca[i] << (dsqvol*jserca[i])                      : used to make changes in cytoplasmic calcium
	~ caer[i] << (-dsqvoler*jserca[i])                 : used to make changes in ER calcium
}

I have used similar codes for other mechanisms. The code doesn't give errors with modlunit and mknrndll.
Is this approach correct and will give me correct values for calcium changes in both ER and cytoplasm?

[ted]

83% volume is occupied by cytoplasm and 12% by ER

12% or 17%?

The COMPARTMENT statements look correct.

caer is the calcium concentration of ER initialized to a resting value.

Not sure I understand why it isn't sufficient to call caer the free calcium concentration in the ER; it's a state variable, not an initialization constant, right?

For one of the mechanisms, SERCA, I have the following code (in the KINETIC block)

This also looks quite reasonable, and the fact that modlunit passes it is reassuring. The proof that it works properly would be to perform a very simple test run and analyze the results, even if only for a few short time steps, to verify that they make sense, then perturb one of the parameters (buffering for example) and run a new simulation to make sure that simulation results are appropriately altered. I have seen a case in which everything looked correct but the simulation results were bizarre, all because of a simple typographical error that was easily overlooked by multiple expert programmers/NEURON users (the cause was a - sign that should have been a +).

[Darshan]

Thank you, Ted for the quick reply!

Sorry for the typo. Yes, it is 17%.

Yes, caer is a state variable and not a constant. I have initialized it with a resting value in the INITIAL block.

I will analyze the results by changing parameters and check to see if I get the desired output. I will also look thoroughly for any sign and typographical errors. You are correct, they are sometimes very hard to find.