initialization of calcium concentration in a segment with attached point process

Anything that doesn't fit elsewhere.
Post Reply
rth
Posts: 41
Joined: Thu Jun 21, 2012 4:47 pm

initialization of calcium concentration in a segment with attached point process

Post by rth »

I have a model of a V1 pyramidal neuron, with simplified morphology. Somatic, basal and apical compartments have calcium channels and calcium dynamics.
Zip archive with all mod files and test-model example can be downloaded here https://storage.r-a-r.org/index.php/s/jLCz6m4k3TXFgmK

Only a few channels use calcium

Code: Select all

CaDynamics_E2.mod:      USEION ca READ ica WRITE cai
Ca_HVA.mod:             USEION ca READ eca WRITE ica
Ca_LVAst.mod:           USEION ca READ eca WRITE ica
SK_E2.mod:              USEION ca READ cai
In this model external calcium concentration may be different from the standard, as well as initial internal calcium concentration, therefor, ino_style is set for all segments in all sections:

Code: Select all

        for secname in self.sections:
            if secname == "axon": continue
            h.ion_style("ca_ion",3,2,1,1,0,sec=self.sections[secname])
            self.sections[secname].cao = 1.2
            for seg in self.sections[secname]:
                seg.cai = 0.001
The basal and apical dendrites bear synapses with NMDA and AMPA glutamate receptors. These synapses also have synaptic depression and they are plastic. I combine Tsodyks et al 1998 model https://senselab.med.yale.edu/ModelDB/s ... mod#tabs-2 with my implementation of Graupner and Brunel 2007 https://www.pnas.org/content/109/10/3991.short calcium-dependent learning. The synapse is a point processes which reads cai, cao and write ica, as well as NONSPECIFIC_CURRENT

Code: Select all

COMMENT

AMPA, NMDA, tonic glutamate, and GAMAb
synaptic depression is taken from T&M model
https://senselab.med.yale.edu/ModelDB/showmodel?model=3815&file=/tsodyks/tmgsyn.mod#tabs-2

ENDCOMMENT

NEURON {
	POINT_PROCESS AMPAandNMDAplastic
	: -- SYNAPSE -- :
	:parameters
	RANGE AMPAt1, AMPAt2, AMPAnorm, NMDAt1, NMDAt2, NMDAnorm, e, AMPAfrac, NMDAfrac, AMPAfracca, NMDAfracca
	:synaptic depression
	RANGE tau_1, tau_rec, u0
	NONSPECIFIC_CURRENT i
	USEION ca READ cai,cao WRITE ica
	: -- PLASTICITY -- :
    : STDP
    RANGE gammad, gammap, thetap, thetad, taurho, rho12, rhoinit
    : STDP FLAGs
    RANGE bistable, plastic
	:readouts
	RANGE i, ampai, nmdai, ampag, nmdag, xrho
    
}

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

PARAMETER {
	AMPAt1       = 0.1 (ms) <1e-9,1e9>
	AMPAt2       = 10 (ms) <1e-9,1e9>
	AMPAnorm     = 0
	AMPAfrac     = 0.5 () <0,1>
	NMDAt1       =  9 (ms) <1e-9,1e9>
	NMDAt2       = 50 (ms) <1e-9,1e9>
	NMDAnorm     = 0
	NMDAfrac     = 0.5 () <0,1> 
	e            = 0 (mV)  
	AMPAfracca   = 0.13                 : fraction of current that is ca ions; Srupuston &al 95
	NMDAfracca   = 0.13                 : fraction of current that is ca ions; Srupuston &al 95
	: tau_1 was the same for inhibitory and excitatory synapses
	: in the models used by T et al.
	tau_1        = 3 (ms) < 1e-9, 1e9 >
	: tau_rec = 800 ms for excitatory
	tau_rec      = 800 (ms) < 1e-9, 1e9 >
	: u0 amount of x->y
	u0           = 0.04 <0,1>
    : STDP parameters 
    gammad = 200 <0,50000>              : gain of depression
    gammap = 321.808 <0,25000>          : gain of potentiation
    thetad = 1 (mM) <1e-9,1e9>          : threshold of depression
    thetap = 1.3 (mM) <1e-9,1e9>        : threshold of potentiation
    taurho = 150e3 (ms) <25e3,2500e3>   : time constant of STDP 
    rho12  = 0.5                        : middle point of cubic polynomial term
    rhoinit = 0.1                       : initial value for rho
    : STDP FLAGs
    bistable = 0 <0,1>                  : FLAG if 0 - STDP is not bistable, if 1 - STDP has bistability.
    plastic = 1 <0,1>                   : FLAG if 1 - synapse is plastic, i.e. STDP is active.

}

ASSIGNED {
	v (mV)
	i (nA)
	ik (nA)
	ampag (uS)
	nmdag (uS)
	ampai (nA)
	nmdai (nA)
	AMPAfactor
	NMDAfactor
	mgblock (1)
	ica	  (nA)
	x
	cai     (mM)
	cao     (mM)
	celsius
}

STATE {
	ampaA (uS)
	ampaB (uS)
	nmdaA (uS)
	nmdaB (uS)
	xrho (1)
}

INITIAL {

	LOCAL tp, factor

	if (AMPAt1/AMPAt2 > 0.9999) { AMPAt1 = 0.9999*AMPAt2 }
	if (AMPAt1/AMPAt2 < 1e-9  ) { AMPAt1 = AMPAt2*1e-9   }
	ampaA = 0
	ampaB = 0
	if (AMPAnorm > 0.5){
		AMPAfactor = 1./(AMPAt2 - AMPAt1)
	} else {
		tp = (AMPAt1*AMPAt2)/(AMPAt2 - AMPAt1)*  log(AMPAt2/AMPAt1)
		factor = -exp(-tp/AMPAt1) + exp(-tp/AMPAt2)
		AMPAfactor = 1./factor
	}

	if (NMDAt1/NMDAt2 > 0.9999) { NMDAt1 = 0.9999*NMDAt2 }
	if (NMDAt1/NMDAt2 < 1e-9  ) { NMDAt1 = NMDAt2*1e-9   }
	nmdaA = 0
	nmdaB = 0
	if (NMDAnorm > 0.5 ){
		NMDAfactor = 1./(NMDAt2 - NMDAt1)
	} else {
		tp = (NMDAt1*NMDAt2)/(NMDAt2 - NMDAt1)*  log(NMDAt2/NMDAt1)
		factor = -exp(-tp/NMDAt1) + exp(-tp/NMDAt2)
		NMDAfactor = 1./factor
	}
	xrho = rhoinit
}

BREAKPOINT {

	LOCAL i_ca2p

	SOLVE state METHOD cnexp
	if (plastic){
		ampag = xrho*(ampaB - ampaA)
		nmdag = xrho*(nmdaB - nmdaA)
	} else {
		ampag = ampaB - ampaA
		nmdag = nmdaB - nmdaA
	}
	: Jahr Stevens 1990 J. Neurosci
	mgblock = 1.0 / (1.0 + 0.28 * exp(-0.062(/mV) * v) )
	ampai =  ampag * (v - e) *           (1-AMPAfracca)
	nmdai =  nmdag * (v - e) * mgblock * (1-NMDAfracca)
	i_ca2p = 0
	if(AMPAfracca>0.0){ : assuming only NMDA pumps calcium
		i_ca2p = i_ca2p +  ampag * ghkg(v,cai,cao,2) *           AMPAfracca
	}
	if(NMDAfracca>0.0){ : assuming only NMDA pumps calcium
		i_ca2p = i_ca2p +  nmdag * ghkg(v,cai,cao,2) * mgblock * NMDAfracca
	}
	ica = i_ca2p
	i = ampai + nmdai
}

DERIVATIVE state {
	LOCAL rhoinf
	rhoinf =  -bistable*xrho*(1.-xrho)*(rho12-xrho)+gammap*(1.-xrho)*xheav(cai-thetap)-gammad*xrho*xheav(cai-thetad)
:             |      STDP bistability              |            Potentiation          |          Depression        |
	ampaA' = -ampaA/AMPAt1
	ampaB' = -ampaB/AMPAt2
	nmdaA' = -nmdaA/NMDAt1
	nmdaB' = -nmdaB/NMDAt2
	xrho'  =  rhoinf/taurho
	:DB>>
	printf("DB> cai=%g, cao=%g\n",cai,cao)
	:<<DB
}

NET_RECEIVE(weight (uS), y, z, tsyn (ms) ) {
INITIAL {
	y = 0
	z = 0
	tsyn = t
}
	: first calculate z at event-
	:   based on prior y and z
	z = z*exp(-(t - tsyn)/tau_rec)
	z = z + ( y*(exp(-(t - tsyn)/tau_1) - exp(-(t - tsyn)/tau_rec)) / ((tau_1/tau_rec)-1) )
	: now calc y at event-
	y = y*exp(-(t - tsyn)/tau_1)


	x = 1-y-z
	y = y + x*u0
	tsyn = t
	ampaA = ampaA + weight*x*u0*AMPAfactor*AMPAfrac
	ampaB = ampaB + weight*x*u0*AMPAfactor*AMPAfrac
	nmdaA = nmdaA + weight*x*u0*NMDAfactor*NMDAfrac
	nmdaB = nmdaB + weight*x*u0*NMDAfactor*NMDAfrac
}

: from
:  http://senselab.med.yale.edu/modeldb/ShowModel.asp?model=144490&file=\bpap\ghk.inc

COMMENT
    GHK function that returns effective driving force
    Slope at low voltages is 1
    z needs to be set as a PARAMETER
ENDCOMMENT

FUNCTION ghkg(v(mV), ci(mM), co(mM), Z) (mV) {
    LOCAL xi, f, exi, fxi
    f = R*(34.+273.15)/(Z*(1e-3)*FARADAY)
    xi = v/f
    exi = exp(xi)
    if (fabs(xi) < 1e-4) {
        fxi = 1 - xi/2
    }else{
        fxi = xi/(exi - 1)
    }
    ghkg = f*((ci/co)*exi - 1)*fxi
}

FUNCTION ghk(v(mV), ci(mM), co(mM), Z) (.001 coul/cm3) {
    LOCAL xi, f, exi, fxi
    f = R*(34.+273.15)/(Z*(1e-3)*FARADAY)
    xi = v/f
    exi = exp(xi)
    if (fabs(xi) < 1e-4) {
        fxi = 1 - xi/2
    }else{
        fxi = xi/(exi - 1)
    }
    ghk = (.001)*Z*FARADAY*(ci*exi - co)*fxi
}

FUNCTION xheav( x (1) ) (1) {
    if ( x >= 0) {
        xheav = 1
    } else {
        xheav = 0
    }
}

The problem appears during the initialization of this synapse. If this synapse connected to a segment, inner and outer calcium concentrations are set into some default numbers, even though cinit in ion_style is set to zero.

Let's connect it a middle of apical dendrite:

Code: Select all

xsyn = h.AMPAandNMDAplastic(0.5,sec=xn.sections["apic"])
The loop through all segments just before and just after init, shows what happens.

Code: Select all

    for secname in xn.sections:
        if secname == "axon": continue
        logging.debug("Section {} cao={}".format(secname,xn.sections[secname].cao) )
        for seg in xn.sections[secname]:
            logging.debug(" > {}, {}".format(seg.cai, seg.cao) )

    
    logging.debug("--- NEURON INIT ---")
    h.finitialize()
    h.fcurrent()
    h.frecord_init()
    for secname in xn.sections:
        if secname == "axon": continue
        logging.debug("Section {} cao={}".format(secname,xn.sections[secname].cao) )
        for seg in xn.sections[secname]:
            logging.debug(" > {}, {}".format(seg.cai, seg.cao) )
Here what this code prints in the log.

Code: Select all

2020-07-27 00:19:15,301:178   DEBUG   :Section apic cao=1.2
2020-07-27 00:19:15,301:180   DEBUG   : > 0.001, 1.2
2020-07-27 00:19:15,302:180   DEBUG   : > 0.001, 1.2
2020-07-27 00:19:15,302:180   DEBUG   : > 0.001, 1.2
2020-07-27 00:19:15,302:180   DEBUG   : > 0.001, 1.2
2020-07-27 00:19:15,302:180   DEBUG   : > 0.001, 1.2
2020-07-27 00:19:15,302:180   DEBUG   : > 0.001, 1.2
2020-07-27 00:19:15,302:180   DEBUG   : > 0.001, 1.2
2020-07-27 00:19:15,302:180   DEBUG   : > 0.001, 1.2
2020-07-27 00:19:15,302:180   DEBUG   : > 0.001, 1.2
2020-07-27 00:19:15,302:180   DEBUG   : > 0.001, 1.2
2020-07-27 00:19:15,302:180   DEBUG   : > 0.001, 1.2
2020-07-27 00:19:15,302:180   DEBUG   : > 0.001, 1.2
2020-07-27 00:19:15,302:180   DEBUG   : > 0.001, 1.2
2020-07-27 00:19:15,302:180   DEBUG   : > 0.001, 1.2
2020-07-27 00:19:15,302:180   DEBUG   : > 0.001, 1.2
2020-07-27 00:19:15,302:180   DEBUG   : > 0.001, 1.2
2020-07-27 00:19:15,302:180   DEBUG   : > 0.001, 1.2
2020-07-27 00:19:15,302:180   DEBUG   : > 0.001, 1.2
2020-07-27 00:19:15,302:180   DEBUG   : > 0.001, 1.2
2020-07-27 00:19:15,302:180   DEBUG   : > 0.001, 1.2
2020-07-27 00:19:15,302:180   DEBUG   : > 0.001, 1.2
2020-07-27 00:19:15,302:180   DEBUG   : > 0.001, 1.2
2020-07-27 00:19:15,302:180   DEBUG   : > 0.001, 1.2
2020-07-27 00:19:15,302:180   DEBUG   : > 0.001, 1.2
2020-07-27 00:19:15,302:180   DEBUG   : > 0.001, 1.2
2020-07-27 00:19:15,302:180   DEBUG   : > 0.001, 1.2
2020-07-27 00:19:15,302:180   DEBUG   : > 0.001, 1.2
2020-07-27 00:19:15,302:180   DEBUG   : > 0.001, 1.2
2020-07-27 00:19:15,303:180   DEBUG   : > 0.001, 1.2
2020-07-27 00:19:15,303:180   DEBUG   : > 0.001, 1.2
2020-07-27 00:19:15,303:180   DEBUG   : > 0.001, 1.2
2020-07-27 00:19:15,303:180   DEBUG   : > 0.001, 1.2
2020-07-27 00:19:15,303:180   DEBUG   : > 0.001, 1.2
2020-07-27 00:19:15,303:180   DEBUG   : > 0.001, 1.2
2020-07-27 00:19:15,303:180   DEBUG   : > 0.001, 1.2
2020-07-27 00:19:15,303:180   DEBUG   : > 0.001, 1.2
2020-07-27 00:19:15,303:180   DEBUG   : > 0.001, 1.2
2020-07-27 00:19:15,303:180   DEBUG   : > 0.001, 1.2
2020-07-27 00:19:15,303:180   DEBUG   : > 0.001, 1.2
2020-07-27 00:19:15,303:180   DEBUG   : > 0.001, 1.2
2020-07-27 00:19:15,303:180   DEBUG   : > 0.001, 1.2
2020-07-27 00:19:15,303:180   DEBUG   : > 0.001, 1.2
2020-07-27 00:19:15,303:180   DEBUG   : > 0.001, 1.2
2020-07-27 00:19:15,303:180   DEBUG   : > 0.001, 1.2
2020-07-27 00:19:15,303:180   DEBUG   : > 0.001, 1.2
2020-07-27 00:19:15,303:180   DEBUG   : > 0.001, 1.2
2020-07-27 00:19:15,303:180   DEBUG   : > 0.001, 1.2
2020-07-27 00:19:15,303:180   DEBUG   : > 0.001, 1.2
2020-07-27 00:19:15,303:180   DEBUG   : > 0.001, 1.2
2020-07-27 00:19:15,304:180   DEBUG   : > 0.001, 1.2
2020-07-27 00:19:15,304:180   DEBUG   : > 0.001, 1.2
2020-07-27 00:19:15,304:180   DEBUG   : > 0.001, 1.2
2020-07-27 00:19:15,304:180   DEBUG   : > 0.001, 1.2
2020-07-27 00:19:15,304:178   DEBUG   :Section bas0 cao=1.2
.....
2020-07-27 00:19:15,305:183   DEBUG   :--- NEURON INIT ---
2020-07-27 00:19:15,305:189   DEBUG   :Section soma cao=1.2
2020-07-27 00:19:15,305:191   DEBUG   : > 0.001, 1.2
2020-07-27 00:19:15,306:189   DEBUG   :Section apic cao=2.0
2020-07-27 00:19:15,306:191   DEBUG   : > 0.001, 1.2
...
2020-07-27 00:19:15,307:191   DEBUG   : > 0.001, 1.2
2020-07-27 00:19:15,307:191   DEBUG   : > 5e-05, 2.0
2020-07-27 00:19:15,307:191   DEBUG   : > 0.001, 1.2
.....
I use NEURON -- VERSION 7.6.7 HEAD (603da174) 2019-04-19 on Ubuntu (like) Linux OS.
ted
Site Admin
Posts: 5793
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: initialization of calcium concentration in a segment with attached point process

Post by ted »

I see what you mean. Python statements in the L5neuron class's definition explicitly assign values to cai (1.0e-3 mM) and cao (1.2 mM) in each segment that has a mechanism with a USEION ca statement. And when you create a new L5neuron instance, those are the calcium concentrations in those segments. But calling h.finitialize() changes cai and cao in one segment of the cell: the segment to which the AMPAandNMDAplastic synaptic mechanism is attached. The new values are controlled by NEURON's built-in "default initial calcium concentration" parameters h.cai0_ca_ion and h.cao0_ca_ion (5.0e-5 mM and 2.0 mM, respectively).

You probably think this is a bug, but I think the real bug is that cai and cao didn't change to 5.0e-5 mM and 2.0 mM in ALL segments that use ca. Why do I think this? Because the CaDynamics_E2 mechanism, which controls calcium concentration, has no INITIAL block to initialize the values of cai and cao, and the USEION ca statement WRITEs only cai but not cao. So if the NMODL code doesn't say what to do about initializing cai and cao, the built-in default constants should rule.

Before I suggest a solution, here are a couple of comments:
1. Destexhe et al. 1994 may have been written at a time when it was very cumbersome, if not impossible, to initialize ionic concentrations with NMODL.
2. Maybe what I regard as the "correct" action didn't happen because of the ion_style statement.

But this is one case that doesn't need an ion_style statement. NEURON is smart enough to know what to do if you make the following changes to CaDynamics.mod, then recompile the mod files. These recommendations assume that you want cao to be a parameter, not a state variable that might change during a simulation. If that assumption is incorrect, let me know. You'll need to make a few changes to CaDynamics_E2.mod and to your main program file L5PC-db.py

To take care of initialization of cao, you could insert the statement
h.cao0_ca_ion = 1.2
anywhere you like in L5PC-db.py. I'd put it outside of the cell class definition--maybe right after the end of the cell class definition. It changes the default value of cao to 1.2 mM. This isn't elegant, because h.cao0_ca_ion affects cao in all segments with a mechanism that has a
USEION ca
statement. I'll mention a cleaner approach later.

To initialize cai, only a few changes to CaDynamics_E2.mod are needed.
In the PARAMETER block insert
cai0 = 1e-3 (mM)

Add this INITIAL block
INITIAL {
cai = cai0
}

For the sake of stability, in the BREAKPOINT block I'd change the SOLVE method to derivimplicit (otherwise cai might fall < 0).

If in the future you'd like to have different cai in different segments, you might declare
RANGE cai0
in the NEURON block.

And that should take care of it.


A cleaner way to deal with cao would be to have CaDynamics_E2.mod control it, instead of leaving it under the control of h.cao0_ca_ion. You'd have to
--make the USEION statement WRITE cai, cao
--(optionally) make cao0 a RANGE variable
--declare cao0 and assign it a default value in the PARAMETER block
--insert cao = cao0 in the INITIAL block
--insert cao (mM) in the STATE block
--insert the ODE cao' = 0 in the DERIVATIVE block
ted
Site Admin
Posts: 5793
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: initialization of calcium concentration in a segment with attached point process

Post by ted »

One other change I'd make to L5PC-db.py: get rid of this part of the cell class definition

Code: Select all

        for secname in self.sections:
            if secname == "axon": continue
            h.ion_style("ca_ion",3,2,1,1,0,sec=self.sections[secname])
            self.sections[secname].cao = 1.2
            for seg in self.sections[secname]:
                seg.cai = 0.001
because it no longer serves a useful purpose.
rth
Posts: 41
Joined: Thu Jun 21, 2012 4:47 pm

Re: initialization of calcium concentration in a segment with attached point process

Post by rth »

Dear Ted,

Thank you so much, for your very helpful comments and suggestions. I saw the problem with the INITIAL block in the CaDynamics_E2 module but somehow ignored it. We need to control cao, because, it seems extracellular calcium concentration drifts during the development period.

I implemented all of your suggestions including control for cao, and they work great (^-^)

Thank you,
-rth
hines
Site Admin
Posts: 1601
Joined: Wed May 18, 2005 3:32 pm

Re: initialization of calcium concentration in a segment with attached point process

Post by hines »

You probably think this is a bug,
I think we need to expand the documentation about ion_style a bit in this area. ie. from
neuron.yale.edu/neuron/static/py_doc/modelspec/programmatic/ions.html?highlight=ion_style
Warning: if other mechanisms are inserted subsequent to a call of this function, the style will be “promoted” according to the rules associated with adding the used ions to the style previously in effect.
And what happened in your case is that after you called the ion_style function with the cinit slot = 0, adding the synapse that READ cai, promoted
it (for that segment only) back to 1. This was a consequence of the fact that c_style specifies cai as a STATE variable along with the promotion rules which are specified in the internal "nrn_promote" function in order as
1) c_style is the largest of its previous value and the new insertion. ie. 0,AMBIGUOUS < 1,PARAMETER < 2,ASSIGNED < 3,STATE
2) e_style is the largest of its previous value and the new insertion
3) if c_style > 0 and e_style <2 then e_style set to 2,ASSIGNED
4) if c_style is STATE then cinit=1,ON (this is what caught you)
5) if c_style is STATE and e_style is ASSIGNED then e_advance is 1,ON
6) if c_style is not AMBIGUOUS and e_style is ASSIGNED then e_init is 1,ON
ted
Site Admin
Posts: 5793
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: initialization of calcium concentration in a segment with attached point process

Post by ted »

And what happened in your case is that after you called the ion_style function with the cinit slot = 0, adding the synapse that READ cai, promoted it (for that segment only) back to 1. . . .
Which suggests that the use of ion_style is surrounded by caveats, because a later change (e.g. adding of a new mechanism that READs or WRITEs some ion concentration) can produce a mismatch between the conceptual model and what actually exists in the computer. If the mismatch is very localized--e.g. limited to just one segment, as in this particular case--the symptoms may also be very localized, so that the cause is difficult to discover.
Post Reply