Integrating RxD Ca reactions with MOD mechanisms

[Tuoma]

I'm trying to have both MOD mechanisms and RxD reactions affect my calcium concentrations. Based on previous examples this seems to be possible but I haven't got it to work properly. Here's a MWE:

Code: Select all

from pylab import *
from neuron import h,rxd

h("""                                                                                                                                                                                                                                                                     
load_file("stdlib.hoc")                                                                                                                                                                                                                                                   
load_file("stdrun.hoc")                                                                                                                                                                                                                                                   
                                                                                                                                                                                                                                                                          
create soma                                                                                                                                                                                                                                                               
soma nseg=1                                                                                                                                                                                                                                                               
soma insert CaDynamics_E2                                                                                                                                                                                                                                                 
soma gamma_CaDynamics_E2 = 1.0                                                                                                                                                                                                                                            
soma insert can
soma insert pas                                                                                                                                                                                                                                                           
objref cvode                                                                                                                                                                                                                                                              
cvode = new CVode()                                                                                                                                                                                                                                                       
cvode.active(1)                                                                                                                                                                                                                                                           
cvode.atol(0.00005)                                                                                                                                                                                                                                                       
""")
sec = h.soma
cyt = rxd.Region([sec], name='cyt'+str(sec), nrn_region='i')
ca = rxd.Species(cyt, name='ca', charge=2, initial=0.0)
cab = rxd.Species(cyt, name='cab', charge=0, initial=1.0)
cabca = rxd.Species(cyt, name='cabca', charge=0, initial=0.0)

cabuf_reaction = rxd.Reaction(ca + cab, cabca, 1.0, 0.1)

ca_flux_rate = rxd.Parameter(cyt, initial=0)
reaction_flux = rxd.Rate(ca, ca_flux_rate)

h.tstop = 400

vecs = []; 
vecs.append(h.Vector()); vecs[-1].record(ca.nodes(sec)(0.5)[0]._ref_concentration)
vecs.append(h.Vector()); vecs[-1].record(cabca.nodes(sec)(0.5)[0]._ref_concentration)
h("""                                                                                                                                                                                                                                                                     
objref tvec, carec, icarec, Vrec                                                                                                                                                                                                                                          
tvec = new Vector()                                                                                                                                                                                                                                                       
carec = new Vector()                                                                                                                                                                                                                                                      
icarec = new Vector()                                                                                                                                                                                                                                                     
Vrec = new Vector()                                                                                                                                                                                                                                                       
cvode.record(&v(0.5),Vrec,tvec)                                                                                                                                                                                                                                           
cvode.record(&cai(0.5),carec,tvec)                                                                                                                                                                                                                                        
cvode.record(&ica(0.5),icarec,tvec)                                                                                                                                                                                                                                       
""")

h("""objref stim                                                                                                                                                                                                                                                          
soma stim = new IClamp(0.5)                                                                                                                                                                                                                                               
stim.amp = 200                                                                                                                                                                                                                                                            
stim.dur = 5                                                                                                                                                                                                                                                              
stim.del = 50                                                                                                                                                                                                                                                             
""")

h.init()

def set_param(param, val):
  param.nodes.value = val
  h.cvode.re_init()

h.cvode.event(250, lambda: set_param(ca_flux_rate, 0.005))
h.cvode.event(350, lambda: set_param(ca_flux_rate, 0))

h.continuerun(h.tstop)

f,axarr = subplots(2,2)
ax = [axarr[0,0],axarr[0,1],axarr[1,0],axarr[1,1]]
ax[0].plot(array(h.tvec),array(vecs[0]),label='Ca (mM)')
ax[0].plot(array(h.tvec),array(vecs[1]),label='CaBCa (mM)')
ax[1].plot(array(h.tvec),array(h.Vrec),label='V (mV)')
ax[2].plot(array(h.tvec),array(h.carec),label='[Cai] (mM)')
ax[3].plot(array(h.tvec),array(h.icarec),label='ica (mA/cm2)')
for a in ax: a.legend()
ax[2].set_xlabel('t (ms)'); ax[3].set_xlabel('t (ms)')
f.savefig("mwe_ca.eps")

Here, the N-type Ca2+ current is as in https://modeldb.science/195615?tab=2&fi ... an_mig.mod and CaDynamics_E2 as in https://senselab.med.yale.edu/ModelDB/s ... mod#tabs-2

The script introduces two inputs: an IClamp that elevates the membrane potential at t=50 ms that induces Ca2+ current through the N-type Ca2+ channels and a flux of RxD-species 'ca' from 250 to 350 ms. However, the Ca2+ current through the N-type Ca2+ channels does not contribute to the intracellular Ca2+ (recorded in carec), as seen in the lower-left panel. If I rename the species as

Code: Select all

ca = rxd.Species(cyt, name='caFoo', charge=2, initial=0.0)

then I get the opposite: intracellular Ca2+ (recorded in carec) is affected by the N-type Ca2+ currents but not the now wrongly named RxD-species.

Any suggestions how to solve this issue?

[ted]

I don't see anything in your program that specifies the physical size of soma. The default values for a section's diam and L are 500 and 100 um, respectively. The volume of your model soma is huge--almost 2e7 um3--and the surface area is only about 1.6e5 um2, so the surface/volume ratio is 0.008. The volume of a 10 um diameter sphere is only about 520 um3 and its area about 310 um2, so the surface/volume ratio is 0.6. In the absence of buffering, for a given calcium current density, cai in a 10 um diam sphere should increase about 75 times faster than in your model cell. << revised to correct numerical errors when this post was first created >>

Simplificiation is often a useful first step in debugging. Suggest you first assign "correct" values to L and diam (pick values appropriate for whatever cell class you're thinking about), and omit the rxd stuff. Also replace the IClamp with an SEClamp so you control v directly. And instead of calling h.init(), call h.finitialize() with a numerical argument equal to whatever you want the initial membrane potential to be (-65 mV? -70? whatever you choose, be sure to use that as the SEClamp's amp1 and amp3). What happens to ica and cai during a simulation? Note that most ca entry happens after the depolarizing step ends (for the sake of others who may read this thread, why?).

[Tuoma]

Thanks, changing the compartment size and making the clamp stimulus longer indeed solved part of the problem (previously the Ca flux through the RxD rate input was very large compared to the Ca flux through the N-type channels).

So now I see that both Ca fluxes contribute to cai.

However, I'm puzzled about the observation that manipulating CaDynamics_E2 parameters doesn't seem to affect the simulations. If I set

Code: Select all

soma gamma_CaDynamics_E2 = 0.0

then ica through Ca2+ channels should no longer increase cai, right? However, I see that cai is still increased by both RxD Ca2+ flux and Ca2+ flux through N-type channels. By contrast, when I rename the RxD species 'ca' to 'caFoo' as before, then setting gamma_CaDynamics to 0 disconnects cai from Ca2+ flux through N-type channels as it should. Do you have some idea what is going on here, or am I still missing something?

[ted]

First, a side-comment: for the sake of stability, the SOLVE statement in CaDynamics_E2.mod should really use derivimplicit instead of cnexp. For further details, go to the Forum's Hot tips area and read
Integration methods for SOLVE statements

Now, returning to your latest question. The fact that setting gamma to 0 allows your simulation to show cai varying in response to nonzero ica means that something other than CaDynamics_E2 is controlling cai. Comment out every statement in your Python file that refers to CaDynamics_E2, and run a new simulation. I bet this has no effect on your simulation results. If that prediction is correct, then your problem is interference between two mechanisms that are trying to set the value of cai.

That couldn't happen if both mechanisms had been implemented with NMODL, because NEURON's computational engine would have detected it and generated an error message. The same should happen if one mechanism is implemented with NMODL and the other with rxd; failure to generate an error message under such circumstances would be a bug.

Which leads me to ask: what version of NEURON are you using?

For any given solute in a pariticular compartment, there should only be one mass balance equation in that compartment. Otherwise conservation of mass can be violated. If this is what's going on, you need to make a decision: do you want the code in CaDynamics_E2.mod to govern cai, or do you want something implemented with the rxd module to govern it? If you favoer the former, you need to revise your rxd statements so that they generate an appropriate ca flux. If you prefer the latter, you need to get rid of the CaDynamics_E2 mechanism.

[adamjhn]

Ted is right, CaDynamics_E2 writes to cai but it is then overwritten by rxd. Looking at the mod file it looks like it does a couple of things; calcium accumulation and restoration to minCa.

RxD already handles accumulation, but unlike the mod file it uses the segment's area and volume to calculate the concentration change. If you want to add restoration, a simple way to do it would be to add an rxd Rate, something like;

Code: Select all

from neuron.units import ms, mM
decay = 80 * ms     # rate of removal of calcium
minCai = 1e-4 * mM
E2 = rxd.Rate(ca, -(ca - minCai)/decay)

If you want to model calcium in an outer shell using rxd, you can specify a region with a Shell geometry,

Code: Select all

sur = rxd.Region([sec], nrn_region='i', geometry=rxd.geometry.Shell(0.9,1.0))

This will change the volumes used to calculate the concentration changes due to currents.

[Tuoma]

Thanks to both of you - yes probably better to convert the CaDynamics_E2 mechanism to RxD to avoid conflicting cai update. I'm using NEURON 7.8.2.