using diam and L params inside MODL file

NMODL and the Channel Builder.
Post Reply
ziemek
Posts: 45
Joined: Thu May 23, 2019 8:02 am
Location: Warsaw, Poland
Contact:

using diam and L params inside MODL file

Post by ziemek »

Hi,

I want to use compartment's params (diam, L) inside my MODL mechanism "mymodl":

Code: Select all

NEURON {
    SUFFIX mymodl
}

UNITS {
    (um) = (micron)
}

ASSIGNED {
    diam (um)
    L (um)
}

STATE {
    my_diam (um)
    my_L (um)
}

BREAKPOINT {
    my_diam = diam
    my_L = L
}
hoc file looks like that:

Code: Select all

load_file("nrngui.hoc")
create soma
access soma
{
    nseg = 1
    L = 2.5
    diam = 0.5
    insert mymodl
}
tstop = 3
run()
Hovewer after run hoc file, when I check my_diam and my_l params inside the console:

Code: Select all

oc>my_diam_mymodl
    0.5 
oc>my_L_mymodl
    0 
So clearly diam works but L not. Do you know why? Is it forbidden to use compartment's params other than diam?
ted
Site Admin
Posts: 6287
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: using diam and L params inside MODL file

Post by ted »

NMODL code can access properties of the segment to which it is attached. Examples include range variables, such as v, diam, and transmembrane currents. Another is segment surface area, which is known to NMODL as area. To access any of these from inside a mod file, declare them in the ASSIGNED block. L is a section parameter and is not accessible in this way. I have not seen a situation in which it is necessary to access L from inside a mod file. If you describe what you are trying to do, maybe I can tell you how to accomplish it without explicitly referring to L.
ziemek
Posts: 45
Joined: Thu May 23, 2019 8:02 am
Location: Warsaw, Poland
Contact:

Re: using diam and L params inside MODL file

Post by ziemek »

OK, so I would like to model Ca2+ radial diffusion from the top of the head of the dendritic spine -> through the neck -> to the dendrite (and in the dendrite radially and axially - as well, so Ca2+ may diffuse to the next spine). Basically I would like to model a whole-cell Ca2+ diffusion.

EDIT:
I realized that I had used cylindric morphology too extensively - if I tried to compute multi-shell diffusion through cylindric compartment as it went through 2D rectangular compartment. But is there any other way to compute multi-shell diffusion from the top to the bottom of the dendritic spine?

Nevertheless - below is the description of what I done, even though it's quite incorrect.

Example image: https://drive.google.com/file/d/1MO5Aoh ... sp=sharing

Detailed explanation of the image:
  • head and neck are separated rectangles. Each has diam, L and inner rectangular shells (diffusion is rectangular, meaning - there are no inner annulus-like shells, but horizontal, inner rectangles - stacked from the top to the bottom)
  • dendrite is a circle of annulus inner shells (just like from the Neuron Book example)

Example MODL file of a radial Ca2+ diffusion in the spine head (with stacked, rectangular inner shells):

Code: Select all


NEURON {
	SUFFIX cadifusrect
	USEION ca READ cai, ica WRITE cai
	GLOBAL TotalBuffer
	RANGE cai0
	THREADSAFE
}

DEFINE NANN  4

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

PARAMETER {
        L (um)
	DCa = 0.6	(um2/ms)
	k1buf	= 100 (/mM-ms)
	k2buf	= 0.1	 (/ms)
	TotalBuffer = 0.003 (mM)
	cai0 = 50e-6 (mM)
}

ASSIGNED {
        diam (um)
	ica		(mA/cm2)
	cai		(mM)
	Kd		(/mM)
	B0		(mM)
}

STATE {
	ca[NANN]		(mM) <1e-6>	: ca[0] is equivalent to cai
	CaBuffer[NANN]	(mM)
	Buffer[NANN]	(mM)
}

BREAKPOINT {
	SOLVE state METHOD sparse
}

LOCAL factors_done
LOCAL frat
LOCAL da

INITIAL {
	MUTEXLOCK
	if (factors_done == 0) {
		factors_done = 1
		da = diam/(NANN-1)
                frat = L/da  : area_of_diffusion/thickness_of_the_shell
	}
	MUTEXUNLOCK

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

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

KINETIC state {
	COMPARTMENT i, da*L {ca CaBuffer Buffer} 
	LONGITUDINAL_DIFFUSION i, DCa*da*L {ca}

	~ ca[0] << ((-ica*da*L)/(2*FARADAY))

	FROM i=0 TO NANN-2 { : radial diffusion
		~ ca[i] <-> ca[i+1] (DCa*frat, DCa*frat)
	}

	FROM i=0 TO NANN-1 { : calcium buffering
		~ ca[i] + Buffer[i] <-> CaBuffer[i] (k1buf*da*L, k2buf*da*L)
	}

	cai = ca[0]
}

Key features of the code above:
  • da = diam/(NANN-1)
  • ~ da is a single shell thickness
  • shell_volume = da*L
  • ~ shell_volume is used in: COMPARTMENT, LONGITUDINAL_DIFFUSION and CaBuffer computations
  • frat = L/da
  • ~ I understand that frat is forward and backward rates
  • ~ as I understand the equation (form the Neuron Book as well) it's just a ratio of area_of_diffusion/thickness_of_the_shell
  • ~ since the diffusion starts from the top of the spine head, nominator is not the circumference of the head, but rather its length
RxD:
Yesterday I had a reflection that my approach could be wrong at all, because maybe I should use RxD for the whole cell Ca2+ diffusion, rather then simplified MODL-like modelling?

But I've never used RxD and I have no experience about its features and computational capabilities. So I would be very glad for your comments and help - what's the best approach for my scenario here?
ted
Site Admin
Posts: 6287
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: using diam and L params inside MODL file

Post by ted »

Models of diffusion can be implemented with NMODL, but NMODL is best suited for nearly cylindrical geometries--structures in which diameter changes very slowly along the length of a neurite--and under conditions in which concentration gradients and solute fluxes are either predominantly longitudinal or predominantly radial.

Assuming spine head diameter is about 1 um and effective Ca diffusion constant in cytoplasm is about 20 um2/s, if calcium enters a spine head through a cluster of channels over the course of about 1 ms it will take on the order of (x^2)/(2*D) = 10 ms for [Ca] to equilibrate throughout the spine head (yes, I know the formula assumes 3 D diffusional geometry, but all I need here is a rough estimate), after which d[Ca]/dt will be governed by three processes:
flux through the spine neck (1 dimensional diffusion) to the dendritic shaft
efflux through the spine head membrane (pump, Na-Ca exchange)
uptake by organelles (primarily SER)
These three processes can be handled by NMODL reasonably easily and accurately.

What's not so easily handled are the details of diffusion in the spine head during the first ~10 ms, especially at the junction between the spine head and neck. Invoking LONGITUDINAL_DIFFUSION in NMODL couples corresponding "concentric shells" between adjacent sections (incidentally, L is not needed in NMODL code that includes LONGITUDINAL_DIFFUSION--see Modeling diffusion with kinetic schemes in chapter 9 of the NEURON Book). If the spine head's nseg is 1, the outermost shell of the spine neck will "see" the increase of [Ca] in the outermost shell in the spine head far too early. I suppose one could make the spine head's nseg equal to the number of radial shells, and restrict ligand- and/or v-gated Ca channels to a distal compartment of the spine head. A 3x3 spatial grid would improve temporal resolution to about 1 ms, and a 5x5 grid (5 = # shells = nseg) would improve it to about 0.4 ms. And of course you'd want to have a similar grid in the spine neck (actually you'd need to use the same number of concentric shells in the neck, or you won't get longitudinal diffusion from head to neck). And we haven't even begun to consider what to do about the dendritic shaft . . .

These headaches go away if you implement calcium accumulation and diffusion with rxd. And in the latest version of NEURON, rxd's performance with 3 D diffusional geometry is much improved.
ziemek
Posts: 45
Joined: Thu May 23, 2019 8:02 am
Location: Warsaw, Poland
Contact:

Re: using diam and L params inside MODL file

Post by ziemek »

Thanks for your help and comments, they are very valuable to me.

I started with RxD and it works nice but I have 4 follow-up questions and one comment.

Questions:
  • if you don't specify geometry for the Region - it assumes the whole compartment's segment as a "single shell"?
  • If I do specify geometry - it takes the volume of the compartment's segment and then compute geometry with RxD for the same volume but the new geometry?
  • Why shells diffusion through MultiCompartmentReaction are defined in molecules instead of moles? Is it because volumes of different compartments may differ?
  • Why avogadro's number is scalled by 1e-18? Is it just a liter scale factor?

Code: Select all

mM_to_mol_per_um = 6.0221409e+23 * 1e-18
The comment:
if I have the same cylindric geometry for the head and the neck as earlier, it seems to me that the diffusion from the head to the neck is still a problem.
  • if I left cylindric geometry as earlier (head's and neck's L param is horizontal and diam is assumed to be vertical) - then how to make a diffusion from the head to the neck from the bottom shell of the head?
  • For the other hand - if I rotate head and neck 90 degree (L is vertical, diam horizontal) the nseg actually defines "shells number" (as you mentioned earlier) so any shell diffusion with RxD will go horizontally instead of vertically.
So I still don't know how to deal with that. Should I use 3D geometry of the spine? But how in that case define shells in RxD?
adamjhn
Posts: 54
Joined: Tue Apr 18, 2017 10:05 am

Re: using diam and L params inside MODL file

Post by adamjhn »

1] The default geometry is rxd.geometry.inside. This is effectively a single shell for the segment.

2] Specifying a geometry changes the way rxd handle the same space, e.g.
A region with rxd.Shell(lo, hi) is treated as a hollow cylinder extending from lo multiplied by the segments diameter, to hi multiplied by the segments diameter.

3] MultiCompartmentReaction has to deal with a flux between two compartments depending on the surface area between them. As the compartments may have different volumes, it would complicate things for users to specify the flux in the usually units of (amount/volume/time) as they would have to indicate which volume they intended to be used.
Instead the rates are given in amount (molecules) per surface area (μm²) per time (ms) then rxd uses the relative volumes to compute the flux.

The desired flux for radial diffusion of calcium from Fick's law is (D/dr)Ca
With the diffusion coefficient D (μm²/ms) and calcium Ca (mM or mol/m³) and shell thickness dr (μm)

Multiplying by Avogadro's number and by 1e-18 (note in future releases of NEURON these constants will be available as rxd.constants.NA and neuron.units.m**-3) gives the concentration Ca in (number of molecules per μm³) and the rate has the desired units i.e;

(D/dr)*Ca * rxd.constants.NA * neuron.units.m**-3 has units of number molecules per μm² per ms

4]
A common approach for radial diffusion is to define the same shells (as a proportion of the diameter) in each section, diffusion occurs between matched shells, i.e. outer to outer, inner to inner, etc. This does not make a lot of sense when there is an abrupt change in diameter or when one section is connected at an angle to another.

The simplest way to implement such a solution is to define the same regions (shells and borders) for both the head and neck, then let rxd use longitudinal diffusion (vertical in your schematic), for flux between the segments. Longitudinal diffusion is provided by rxd just by specifying a diffusion coefficient in the rxd.Species declaration. This is how diffusion between the dendrite and soma occurs in the radial diffusion example.

Alternatively if you suspect the abrupt change in diameter or the angle between the neck and the dendrite makes mapping of shells an inappropriate solution, you can solve it in 3D, by calling rxd.set_solve_type. This will be slower and require more memory.
Here is a 3D and 1D/3D hybrid example.
Post Reply