How to call my ionic current variable implemented in NMODL in HOC

[MarvinM]

Hi all,

This is my first post, and I am looking forward to engaging in fruitful discussions in this excellent curated forum.

I have a problem.

Right now, I am trying to calculate the ATP consumption in the NEURON CA1 model I am using (https://www.sciencedirect.com/science/a ... 8021002975) by recording the current for different channel mechanisms that contribute to calcium and sodium flux across the cell membrane. For that, I am using open-source code from a Renaud Jolivet paper (https://journals.plos.org/ploscompbiol/ ... bi.1007226) to read out currents of interest. So far, I succeeded in tracking the sodium current of diverse mechanisms with specified ina currents in my mechanisms and disentangling the sodium part from my ih and leak currents. However, I am having a hard time to track the sodium part of my synaptic mechanisms (NMDAR & AMPAR).

This is the code I am using in HOC:

Code: Select all

//**************************************************************//
// New advance procedure to read out currents
//**************************************************************//
objref vec_sodium, vec_sodium_leak, vec_sodium_hd, vec_calcium, vec_sodium_comp, vec_sodium_ampa, vec_sodium_nmda

itot_na = 0
ina_leak = 0
ina_h = 0
ina_ampa = 0
ina_nmda = 0
ina_comp = 0
itot_ca = 0


vec_sodium = new Vector()
vec_sodium_leak = new Vector()
vec_sodium_hd = new Vector()
vec_sodium_comp = new Vector()
vec_sodium_ampa = new Vector()
vec_sodium_nmda = new Vector()
vec_calcium = new Vector()
vec_sodium.record(&itot_na)
vec_sodium_leak.record(&ina_leak)
vec_sodium_hd.record(&ina_h)
vec_sodium_comp.record(&ina_comp)
vec_sodium_ampa.record(&ina_ampa)
vec_sodium_nmda.record(&ina_nmda)
vec_calcium.record(&itot_ca)
proc read_compartment_current(){
    if (ismembrane("na_ion")){
        ina_comp = ina($1) * 1e-2 * area($1)
        ina_leak = g_pas * (Ek - v_init) / (Ek - Ena) * ( v - Ena) * area($1)
        ina_h = 0
        if (ismembrane("hd")){
            ina_h = ghdbar_hd * (Ek - Eh) / (Ek - Ena) * (v - Ena) * area($1) 
        }
        ina_syn = 0
        //if (ismembrane("ghkampa")){
            //ina_syn = AMPAsyn.Pmax * Ek / (Ek - Ena) * (v - Ena) * area($1)
            //AMPAsyn[j].Pmax * Ek / (Ek - Ena) * (v - Ena) * area($1)
            //}
                    // INa, syns = gsyns *VK / (VK – VNa) * (V – VNa) with gsyns = gAMPA for AMPA type synapses
                    // or INa, nmda = gmax *VK / (VK – VNa) * (V – VNa) * 1/(1 + eta*mg*exp(-gamma*V))
                    // With eta=0.05 mg=2 gamma=0.06 from the NMDA kinetics or fitting for NMDA type synapses
        itot_na = itot_na + ina_comp + ina_leak + ina_h + ina_syn
        }
    if (ismembrane("ca_ion")){
        ica_comp = ica($1) * 1e-2 * area($1)
        itot_ca = itot_ca + ica_comp
    }
}

proc advance(){
    fadvance()
    itot_na = 0
    itot_ca = 0
    forall {
        for (x) {
            if (x==0 || x==1) { continue }
            read_compartment_current(x)
        }
    }
}

My approach can be summarized as, creating empty vectors for the currents, iterating over the cell segments, defining a function to read out currents, having a look if an ion or channel type exists in the respective segment, and using the math to derive the current.

Now to the core of my problem, as I am randomly distributing synapses on the dendritic tree, enumerated from 1 to 80, I don't know in which order they show up. Thus, I can't iterate through my list of point processes (example, AMPAsyn[j].Pmax) in numerical order during my iteration of the segments to look for the conductance to calculate the current. At least, this didn't work out for me. I forgot to mention, another caveat is that I also assign a range of conductances. Hence, I can't specify the current without looking it up in the respective segment. Instead, I would like to implement the sodium part calculation in the NMODL file itself, so I can call the current from the NMODL file directly without the need to iterate through the point processes, which seems not feasible to me with my configuration.

Unfortunately, I am a little bit overwhelmed with calling variables specified in NMODL in HOC.

For example, this is the code for my AMPAR mechanism (not much different from the NMDAR one):

Code: Select all

COMMENT

The kinetics part is obtained from Exp2Syn of NEURON.

Two state kinetic scheme synapse described by rise time taur, and
decay time constant taud. The normalized peak conductance is 1.
Decay time MUST be greater than rise time.

The solution of A->G->bath with rate constants 1/taur and 1/taud is

 A = a*exp(-t/taur) and
 G = a*taud/(taud-taur)*(-exp(-t/taur) + exp(-t/taud))
	where taur < taud

If taud-taur -> 0 then we have a alphasynapse.
and if taur -> 0 then we have just single exponential decay.

The factor is evaluated in the initial block such that an event of
weight 1 generates a peak conductance of 1.

Because the solution is a sum of exponentials, the coupled equations
can be solved as a pair of independent equations by the more efficient
cnexp method.

GHK based ionic currents for AMPA current added by Rishikesh Narayanan.

ENDCOMMENT

NEURON {
	POINT_PROCESS ghkampa
	USEION na WRITE ina
	USEION k WRITE ik

	RANGE taur, taud
	RANGE iampa
	RANGE P, Pmax
}

UNITS {
	(nA) = (nanoamp)
	(mV) = (millivolt)
	(uS) = (microsiemens)
	(molar) = (1/liter)
	(mM) = (millimolar)
	FARADAY = (faraday) (coulomb)
	R = (k-mole) (joule/degC)
}

PARAMETER {
	taur=2 (ms) <1e-9,1e9>
	taud = 10 (ms) <1e-9,1e9>
	nai = 18	(mM)	: Set for a reversal pot of +55mV
	nao = 140	(mM)
	ki = 140	(mM)	: Set for a reversal pot of -90mV
	ko = 5		(mM)
	celsius		(degC)
	Pmax=1e-6   (cm/s)	
}

ASSIGNED {
	ina     (nA)
	ik      (nA)
	v (mV)
	P (cm/s)
	factor
	iampa	(nA)

	Area (cm2)
}

STATE {
	A (cm/s)
	B (cm/s)
}

INITIAL {
	LOCAL tp
	if (taur/taud > .9999) {
		taur = .9999*taud
	}
	A = 0
	B = 0
	tp = (taur*taud)/(taud - taur) * log(taud/taur)
	factor = -exp(-tp/taur) + exp(-tp/taud)
	factor = 1/factor
	Area=1
}

BREAKPOINT {
	SOLVE state METHOD cnexp
	P=B-A

: Area is just for unit conversion of ghk output

	ina = P*ghk(v, nai, nao,1)*Area	
	ik = P*ghk(v, ki, ko,1)*Area
	iampa = ik + ina
}

DERIVATIVE state {
	A' = -A/taur
	B' = -B/taud
}

FUNCTION ghk(v(mV), ci(mM), co(mM),z) (0.001 coul/cm3) {
	LOCAL arg, eci, eco
	arg = (0.001)*z*FARADAY*v/(R*(celsius+273.15))
	eco = co*efun(arg)
	eci = ci*efun(-arg)
	ghk = (0.001)*z*FARADAY*(eci - eco)
    ina_syn = ina
}

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

NET_RECEIVE(weight (uS)) { 	: No use to weight, can be used instead of Pmax,
							: if you want NetCon access to the synaptic
							: conductance.
	state_discontinuity(A, A + Pmax*factor)
	state_discontinuity(B, B + Pmax*factor)
}

Notice, that there is already the calculation of the sodium part of the current implemented: ina = P*ghk(v, nai, nao,1)*Area under BREAKPOINT and previously defined as an ASSIGNED variable. Now, I would like to know how I can call the sodium part from the mechanism in HOC when I am iterating over the segments and encounter a segment with the existence of the channel mechanism. I am uncertain how to specify the variable, i.e., in RANGE, GLOBAL, or a VERBATIM block. And, how to call it thereafter.

I hope you can help me.

All the best,

Marvin

[MarvinM]

Also, this is my try to iterate over my synaptic point processes nested in my segment for-loop (latter not shown here, but see above).

Code: Select all

ina_syn = 0
        if (ismembrane("ghkampa")){
            for(j=0;j<80;j+=1) {
                ina_syn = AMPAsyn[j].Pmax * Ek / (Ek - Ena) * (v - Ena) * area($1)
                }
            }

What I receive is: NEURON: AMPAsyn not an array variable

[ted]

There may be some ways to reduce complexity, but I'll have to ask some questions.

Here's a big suggestion (I don't know why Lavzin et al. didn't do this):
The mod files for several mechanisms generate nonspecific currents. Presumably you want to attribute those currents to corresponding ion fluxes, yes or no? If yes, why not replace the NONSPECIFIC_CURRENT declarations with appropriate USEION statements? If you do that, then every segment that uses ionic species foo will automatically keep track of the net current attributable to foo, called ifoo, and its units will be mA/cm2. For any given segment, reduces the number of vectors to the number of ionic species that move through the membrane, and entirely eliminates the need to iterate over point processes.

Of course, if you really want to focus on just the currents generated by point processes, all you have to do is invent some new ionic species--e.g. inap and ikp (p for point process), or maybe inas and iks (s for synapse). I like the latter better because it focuses on the anatomy of the cell and the network in which it is embedded, NOT on the totally nonbiological and irrelevant fact that you happened to implement your model with NEURON. Then you can record inas and iks to keep track of synaptic na and k fluxes, and ina and ik to discover the na and k fluxes through all other na and k channels. And units will be mA/cm2. Again, you don't have to iterate over the point processes--just record inas and iks in those sections that have inas and iks--and you can set up those recordings AFTER synapses have been attached to your model.

Regarding those USEION statements--

Usually a NONSPECIFIC_CURRENT has a fixed reversal potential. If you're reusing a published model, replacing the NONSPECIFIC_CURRENT statement with
USEION na READ ena WRITE ina
USEION k READ ek WRITE ik
is very likely to make your revised model generate results that differ from those produced by the published original. Why? The gna/gk ratio (which should be based on experimental reports or at least on justifiable assumptions) and the values of ena and ek (which are specified outside of the mod file) will probably make the reversal potential of the revised mechanism different from the original mechanism. Especially if ena and/or ek actually vary in the course of a simulation. Reviewers will not like this at all. It would be better to declare
USEION na WRITE ina
USEION k WRITE ik
and in the PARAMETER block declare
gnagkratio = 10 : or whatever is reasonable to assume
ena = somevalue
ek = othervalue
where ena and ek are assigned values that are resonable and that, combined with the gnagkratio, produce the same reversal potential that was used in the original mod file.

[MarvinM]

Hi Ted,

Thank you very much for your reply. I like your proposal a lot and think that, as you mentioned, it should reduce the complexity.

It is right, that I want to keep track of the currents of nonspecific current and their corresponding ion fluxes. For that, I calculated the sodium part of the nonspecific I-h current from Magee (https://www.jneurosci.org/content/18/19/7613.short). Besides, the synaptic mechanisms (AMPAR&NMDAR) from Rishi Narayanan generate specific currents.

After reading your reply, I would like to go over what I've done so far.
Firstly, I modified the two synaptic NMODL files by declaring
USEION na WRITE inas
USEION k WRITE iks
Moreover, I adapted the subsequent notations of ina and ik in the RANGE, ASSIGNED, and BREAKPOINT blocks.
Then, I tried to recompile.

Code: Select all

NEURON {
	POINT_PROCESS ghkampa
	USEION na WRITE inas
	USEION k WRITE iks

	RANGE taur, taud
	RANGE iampa, inas, iks
	RANGE P, Pmax
}

UNITS {
	(nA) = (nanoamp)
	(mV) = (millivolt)
	(uS) = (microsiemens)
	(molar) = (1/liter)
	(mM) = (millimolar)
	FARADAY = (faraday) (coulomb)
	R = (k-mole) (joule/degC)
}

PARAMETER {
	taur=2 (ms) <1e-9,1e9>
	taud = 10 (ms) <1e-9,1e9>
	nai = 18	(mM)	: Set for a reversal pot of +55mV
	nao = 140	(mM)
	ki = 140	(mM)	: Set for a reversal pot of -90mV
	ko = 5		(mM)
	celsius		(degC)
	Pmax=1e-6   (cm/s)	
}

ASSIGNED {
	inas     (nA)
	iks      (nA)
	v (mV)
	P (cm/s)
	factor
	iampa	(nA)

	Area (cm2)
}

STATE {
	A (cm/s)
	B (cm/s)
}

INITIAL {
	LOCAL tp
	if (taur/taud > .9999) {
		taur = .9999*taud
	}
	A = 0
	B = 0
	tp = (taur*taud)/(taud - taur) * log(taud/taur)
	factor = -exp(-tp/taur) + exp(-tp/taud)
	factor = 1/factor
	Area=1
}

BREAKPOINT {
	SOLVE state METHOD cnexp
	P=B-A

: Area is just for unit conversion of ghk output

	inas = P*ghk(v, nai, nao,1)*Area
	iks = P*ghk(v, ki, ko,1)*Area
	iampa = ik + ina
}

DERIVATIVE state {
	A' = -A/taur
	B' = -B/taud
}

FUNCTION ghk(v(mV), ci(mM), co(mM),z) (0.001 coul/cm3) {
	LOCAL arg, eci, eco
	arg = (0.001)*z*FARADAY*v/(R*(celsius+273.15))
	eco = co*efun(arg)
	eci = ci*efun(-arg)
	ghk = (0.001)*z*FARADAY*(eci - eco)
}

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

NET_RECEIVE(weight (uS)) { 	: No use to weight, can be used instead of Pmax,
							: if you want NetCon access to the synaptic
							: conductance.
	state_discontinuity(A, A + Pmax*factor)
	state_discontinuity(B, B + Pmax*factor)
}

Previously, I called the ina current in the respective segments in HOC by ina($1) (see below for example) to record the ionic current. However, a modification such as inas($1) is not functioning. When I now try to recompile the synaptic NMODL files with mknrndll it will give me the error: "inas is not a valid ionic variable for na at line 32 in file ghkampa.mod", and terminates the compilation. Furthermore, without compiling it does not recognize the variable: "NEURON: inas undefined function".

How could I modify the mechanism to invoke the changes?

Code: Select all

proc read_compartment_current(){
    if (ismembrane("na_ion")){
        ina_comp = ina($1) * 1e-2 * area($1)
        ina_leak = g_pas * (Ek - v_init) / (Ek - Ena) * ( v - Ena) * area($1)
        ina_h = 0
        if (ismembrane("hd")){
            ina_h = ghdbar_hd * (Ek - Eh) / (Ek - Ena) * (v - Ena) * area($1) 
        }
        ina_syn = 0
        if (ismembrane("ghkampa")||ismembrane("ghknmda")){
            ina_syn = inas($1) * 1e-2 * area($1)
        }

[MarvinM]

In addition, I can declare inas a NONSPECIFIC_CURRENT, and implement the equation to calculate it in NMODL - similar to how I showed that I calculated the inas and ikas function before. By that, I can compile the NMODL file. Unfortunately, I still receive the error that inas is undefined, when I try to invoke it in HOC.

[ted]

I modified the two synaptic NMODL files by declaring
USEION na WRITE inas
USEION k WRITE iks

In a USEION statement, the variable name declared in the USEION clause dictates the form of the names declared in the READ and WRITE clauses. Also, NEURON only "knows" about na, k, and ca, so any declaration of a user-created species must include a VALENCE clause. These statements will work:
USEION nas WRITE inas VALENCE 1
USEION ks WRITE iks VALENCE 1

You'll find a lot more about NMODL in chapter 9 of The NEURON Book. If you don't have the book, read https://www.neuron.yale.edu/ftp/neuron/ ... odl400.pdf

I can declare inas a NONSPECIFIC_CURRENT, and implement the equation to calculate it in NMODL

Right, and then you'll have keep track of all synapses that are attached to all segments, and iterate over them to discover how much current they deliver to each segment, etc. etc. etc.. The whole point of
USEION nas WRITE inas VALENCE 1
is that it eliminates the need to iterate over synaptic instances. Each segment that has any mechanisms that use nas will automatically have its own inas variable whose value is the sum of all inas currents generated those mechanisms. The units of that variable will be mA/cm2, so (in pseudocode)

Code: Select all

for each section in the model cell
  for each segment in this section
    net inas in nA = inas in mA/cm2 * segment area in um2 / 100

I should mention that NEURON has a function called area() that returns a segment's surface area in um2. Here's a hoc implementation of the pseudocode in the form of a function that returns the total inas (in nA) for a model cell

Code: Select all

func inas_total() { local tmp
  tmp = 0
  forall for (x,0) tmp += inas(x) * area(x) / 100
  return tmp
}

and python users are welcome to add a python implementation to this discussion thread.

Got an appointment now, but I'll be back soon with more comments about your most recent posts.

< edited after initial post--was in too much of a hurry, had to come back and revise the function slightly >

[ted]

Before digging into implementational details, it is a good idea to verify that what we know, or think we know, makes sense. This is especially important when it comes to mathematical formulas and units.

A synapse generates a current that is specified in "absolute units" (e.g. nA in NEURON) rather than density units (e.g mA/cm2 in NEURON). The GHK formula for synaptic current carried by ionic species X with charge z has the form

Ix = Px v (z^2 F^2/RT) (Xi - Xo exp(-z v F/RT)) / (1 - exp(-z v F/RT))

(e.g. equations 4a-c on page 513 of
Mayer ML and Westbrook GL. Permeation and block of N-methyl-D-aspartic acid receptor channels by divalent cations in mouse cultured central neurones. J Physiol. 1987 Dec;394:501-27. doi: 10.1113/jphysiol.1987.sp016883. PMID: 2451020)

Let's start by verifying that the right and left hand sides of this formula are "dimensionally consistent at the conceptual level". The left hand side is in units of current, which is the same as charge/time, so the right hand side must also be in units of charge/time--but is it?

Start by eliminating the RHS terms that are dimensionless (z and the exponentials). This leaves

Px v (F^2/RT) (Xi - Xo)

v is in units of electrical potential, and F/RT is in 1/(electrical potential), so these cancel each other, and we have

Px F (Xi - Xo)

F is in units of charge/mole and (Xi - Xo) is in units of mole/volume, so we are left with

Px charge/volume

which tells us that Px must be in units of volume/time.

the synaptic mechanisms (AMPAR&NMDAR) from Rishi Narayanan generate specific currents

So you got the original NMDA code from modeldb.science/244848 which is associated with this paper
Basak R and Narayanan R. Active dendrites regulate the spatiotemporal spread of signaling microdomains. PLoS Comput Biol. 2018 Nov 1;14(11):e1006485. doi: 10.1371/journal.pcbi.1006485. eCollection 2018 Nov. PMID: 3038374

They represent the NMDA-dependent current as the sum of sodium, potassium, and calcium components. These are described by equations 8-10, which correspond to the formula for Ix, but instead of a single Px term they use the product of three terms:

• MgB(v), a dimensionless term which represents voltage-dependent Mg channel block
• Pbar_NMDA which they call "maximum permeability of the NMDAR". They don't specify its units, but the corresponding term in their AMPA-dependent current has units of nm/s (apparently nanometers/second; elsewhere in the paper they use nm to mean nanometers)
• a dimensionless scale factor (which they call PNa, PK, or PCa) that is related to permselectivity of the NMDAR-gated channels

They also use a fourth term s(t), which is dimensionless, to specify the time course of channel opening and closing.

So here are several inconsistencies that the authors and reviewers missed. First, the left hand sides of equations 8-10 are in units of charge/time, but the right hand sides are in charge/area. That's a problem. Second, if Pbar_NMDA is truly the "maximum permeability of the NMDAR", then how can the NMDAR's PNa and PK both be 1? Not to mention that its PCa is 10.6!

The AMPAR-dependent current is represented as the sum of sodium and potassium components (equations 14 and 15). In these equations, net permeability is the product of two terms:

• Pbar_AMPA which they call "maximum permeability of the AMPAR . . . whose default value was set at 20 nm/s" (nanometers/second)
• a dimensionless scale factor (PNa or PK) related to permselectivity of the AMPAR-gated channels

There is also a third term, which is dimensionless, that specifies the time course of channel opening and closing.

And of course the same inconsistencies that affect the NMDAR equations also affect the AMPAR equations. It is possible that "20 nm/s" was a typographical error--did the authors mean that to be "20/s" ?. But that wouldn't resolve the inconsistencies surrounding channel permeabilities, especially the issue of how can it make sense for a channel's permeability to a single ionic species far exceed that channel's "maximum permeability"? The authors' NMODL source code will reveal what they actually did.

< revised slightly after initial post for consistency of units >

[MarvinM]

Thank you very much for your extensive elaboration on the mathematical formulas and units, Ted!

This is most welcome and encourages us to think through the mechanisms and theory behind it again.

To sketch out what I've done since your last post, let me elucidate.

I used
USEION nas WRITE inas VALENCE 1
USEION ks WRITE iks VALENCE 1
to declare the new ionic species. Because of that, I could record the current in HOC via inas(x) with:

Code: Select all

    if (ismembrane("nas_ion")) {
        ina_syn = inas($1) * area($1) / 100
    }

Which is already a success in my book.

Thereafter, I used:

Code: Select all

inas = P*ghk(v, nai, nao,1)*Area

to calculate the ionic current, based on the ionic concentrations inside and outside the cell using the GHK equation. However, I am worried about the compatibility with my new ionic species as the GHK formula is using the variables nai and nao which are default notations in NEURON if I am not misdoing.

Moreover, I also used the conventional:

Code: Select all

inas = P*ek/(ek - ena)*(v - ena)*Area

with newly defined equilibrium potentials, similar to the results obtained by GHK in the original mod file.

I also paid attention that I set up the recordings after I attached the synapses to my model.

Nonetheless, the ina current I record is empty. Do you have any more ideas about what I've done wrong and where the problem might be hidden?

Furthermore, I wanted to inquire: Do you have an idea, why:

Code: Select all

ina_syn = 0
        if (ismembrane("ghkampa")){
            for(j=0;j<80;j+=1) {
                ina_syn = AMPAsyn[j].Pmax * Ek / (Ek - Ena) * (v - Ena) * area($1)
                }
            }

is not working in the first place and gives me the error: NEURON: AMPAsyn not an array variable.
Even though I created the enumerated point processes before by:

Code: Select all

forsec "apical" {
        if (raddist(0)<=300){
            for(x) {
                if(x!=0 && x!=1) {
                	// print "x ", x
                    for(j=0;j<e_synloc.x[comp_count];j+=1) {
                        
                        AMPAsyn[e_syncount]= new ghkampa(x)
                        AMPAsyn[e_syncount].taur = 2
                        AMPAsyn[e_syncount].taud = 10
                        AMPAsyn[e_syncount].Pmax = e_gvalue.x[comp_count]
                        
                        NMDAsyn[e_syncount]= new ghknmda(x)
                        NMDAsyn[e_syncount].taur = 5
                        NMDAsyn[e_syncount].taud = 50
                        NMDAsyn[e_syncount].Pmax = e_gvalue.x[comp_count] * NAR
                        
                        conne_a[e_syncount] = new NetCon(nile, AMPAsyn[e_syncount], 0,0,0)
                        conne_n[e_syncount] = new NetCon(nile, NMDAsyn[e_syncount], 0,0,0)
                        
                        e_syncount+=1
                        //print x
                    }
                    comp_count+=1
                    // print "n comps", comp_count
                }    
            }
            print raddist(x)
        }
    }

< corrected some errors after the initial post >

[MarvinM]

Update

I found the error regarding the synaptic iteration in hoc.

Code: Select all

ina_syn = 0
        if (ismembrane("ghkampa")){
            for(j=0;j<80;j+=1) {
                ina_syn = AMPAsyn[j].Pmax * Ek / (Ek - Ena) * (v - Ena) * area($1)
                }
            }

I declared the objref statement later and overlooked the inconsistency.

[ted]

< revised from initial post >

The chief problem with the models of NMDAergic and AMPAergic synaptic transmission used and published by Basak and Narayanan (2018) is they are implemented as point processes but their parameters and variables have units that are appropriate not for point processes but for density (AKA distributed) mechanisms. It may be understandable why they did this--AFAIK all previous examples of NMODL code for currents whose IV relationships are governed by the GHK constant field equation are density mechanisms, and use units appropriate for density mechanisms. For me this was so confusing that I kept falling back into thinking of code for a density mechanism, and my initial writeup of this post reflected that impression.

But no, ghknmda.mod and ghkampa.mod really do define point processes.

The problems with these files are
1. The use of inappropriate units means that the numerical values of key parameters, variables, and scale factors are questionable.
and
2. This is confounded by numerous inconsistencies of units not only between the left- and right-hand sides of many assignment statements, but also between the terms on the right hand side of many assignment statements.

Simulations that use the mechanisms described by ghknmda.mod and ghkampa.mod might produce perfectly correct numerical results, but that remains to be demonstrated by someone else. Before those files are reused, they need to be revised to resolve the units inconsistencies, and that may require changing the numerical values of key parameters.

You might well ask "what units should a point process use?" For currents, the answer is easy: absolute units (e.g. nA), not density currents (mA/cm2). But what units should permeability have? For density mechanisms, which describe ion fluxes through surfaces, permeability is in velocity units (distance/time, e.g. cm/s), but what should it be for fluxes through one or more discrete channels, in other words point processes? All textbook and almost all journal article treatments of membrane currents described by the constant field equation that I have seen treat the current as if it were a flux through a membrane, not as a current through one or more discrete channels, and for the most part they omit any discussion of units. From the GHK current equation one can show that it must be in units of bulk flow (volume/time, e.g. m3/s), but it is easier to read about it in this article
The concept of channel permeability.
W.F. Pickard
Journal of Mathematical Biology 22:11-19, 1985.
https://link.springer.com/content/pdf/1 ... 276543.pdf

I now have a draft of ghknmda.mod that passes testing with modlunit with the exception of one parameter: Pmax, whose units and numerical value are both incorrect. I know what the units must be (m3/sec), but not what its numerical value should be. I expect to resolve this and post a corrected version of ghknmda.mod soon, and a corrected version of ghkampa.mod should follow quickly.

[ted]

Well, maybe not so quickly, allowing for intrusions like a grant application. Here at last are properly written implementations of AMPAergic and NMDAergic synaptic transmission https://www.neuron.yale.edu/ftp/ted/neuron/ampanmda.zip along with demonstration programs that illustrate usage via hoc and Python. Be sure to read usage.txt. Let me know if you have any questions.

[MarvinM]

Thank you, Ted! I truly appreciate the generous spending of your time. 

Your comments helped me understand better the difference between point processes and density mechanisms. Also, it helped me to reflect on my current readout procedure and how to disentangle the currents.

I will try to incorporate your implementation of the AMPAergic and NMDAergic synaptic transmission and observe if it will make any difference in my simulation results.

[ted]

Here is the primary consideration: reusing code that is known to have serious problems is inadvisable for many reasons. It wastes the time and effort of whoever reuses it, because it jeopardizes publication of results, not to mention support of current or future work. And it's not good for computational neuroscience.

Some might ask "Well, is this a serious problem?" Hmm, let's see a show of hands: who in the audience thinks it's a serious problem if a key parameter uses the wrong units and has a numerical value that is multiple orders of magnitude too big or too small? "Well, for the moment, let's try to ignore that. Is there some other problem?" (the sound of jaws dropping open)

"Who's going to know?" whispers someone in a small, dark room with the water running and a radio playing.

Anybody who bothers to check, thanks to the widespread adoption of the FAIR principles by scientific journals and institutional (governmental and non-governmental) funding sources. The number of journals that require model source code to be made publicly accessible is large and growing larger.