Different species parameters for different regions

Extending NEURON to handle reaction-diffusion problems.

Moderators: hines, wwlytton, ramcdougal

Post Reply
bschneiders
Posts: 34
Joined: Thu Feb 02, 2017 11:30 am

Different species parameters for different regions

Post by bschneiders »

What is the simplest way to allow a species to have different initial concentrations in different regions? Do I need to make two different "sub-species"?
For example, I am working with a buffer "buf" that has a higher buffering capacity in the dendrites than in spines. I want to try having either a larger total buffer concentration or different reaction rates in the dendrites than in spines, both yielding different initial concentrations (i.e. "b0_dend","b0_spine"). Do I need to write something like the following?:

Code: Select all

spines = rxd.Region([head,neck],nrn_region='i')
dends = rxd.Region([dend,dend1,dend2],nrn_region='i')
ca =    rxd.Species(r, initial=c0,name='ca',charge=2, d = dCa,atolscale = 10-9)
buf_dend =   rxd.Species(dends, initial=b0_dend,name='buf', d = dB,atolscale = 10-6)
buf_spine = rxd.Species(spines, initial=b0_spine,name='buf', d = dB,atolscale = 10-6)
cabuf_dend = rxd.Species(dends, initial=bc0_dend,name='cabuf', d = dB,atolscale = 10-6)
cabuf_spine = rxd.Species(spines, initial=bc0_spine,name='cabuf', d = dB,atolscale = 10-6)
I tried something like the above, but it doesn't seem to like two species having the same name "buf". However, I want it to know they are in fact the same species since the mobile buffer should be able to diffuse from dendrite to spine, and vice versa.
This question would also trickle down into defining the reactions I imagine (I am using "RxD.Rate" for all of mine). What is the best way to go about this? Thanks!
bschneiders
Posts: 34
Joined: Thu Feb 02, 2017 11:30 am

Re: Different species parameters for different regions

Post by bschneiders »

On a completely unrelated note (and I'm only posting this because Rob told me to post about it to keep the issue relevant...), are you closer to making RxD such that you don't have to restart Python every time you change an RxD-related parameter? That would be awesome!
adamjhn
Posts: 54
Joined: Tue Apr 18, 2017 10:05 am

Re: Different species parameters for different regions

Post by adamjhn »

To set inhomogeneous initial concentrations you can use a function which take an rxd.Node as an argument. Node has several properties; including x (the 1D location), sec, region and segment. If you want the buffers to diffusion between the spine and dendrite you should use the same region.

Here is a simple example with inhomogeneous initial concentrations;

Code: Select all

from neuron import h, rxd
#Parameters
b0_spine, b0_dend = 0.1, 1
bc0_spine, bc0_dend = 0.1 ,1
dB = 1.0


#Geometry
head, neck  = h.Section(), h.Section()
dend1, dend2, dend3 = h.Section(), h.Section(), h.Section()
head.connect(neck)
neck.connect(dend3)
dend3.connect(dend2)
dend2.connect(dend1)

spine_secs = [head,neck]

#Regions
cyt = rxd.Region(h.allsec(),nrn_region='i')

#Initial concentrations
init_buf = lambda nd: b0_spine if nd.sec._sec in spine_secs else b0_dend
init_cbuf = lambda nd: bc0_spine if nd.sec._sec in spine_secs else bc0_dend

#Species
buf = rxd.Species(cyt,name='buf', initial=init_buf, d=dB, atolscale=1e-6)
cabuf = rxd.Species(cyt,name='cabuf', initial=init_cbuf, d=dB, atolscale=1e-6)
Alternatively you could use FInitializeHandler to set the initial values;

Code: Select all

def init_species(species,dend0,spine0):
    for nd in species.nodes:
        nd.concentration = spine0 if nd.sec._sec in spine_secs else dend0

fih = []
fih.append(h.FInitializeHandler(1,lambda: init_species(buf,b0_dend,b0_spine)))
fih.append(h.FInitializeHandler(1,lambda: init_species(cabuf,bc0_dend,bc0_spine)))
bschneiders
Posts: 34
Joined: Thu Feb 02, 2017 11:30 am

Re: Different species parameters for different regions

Post by bschneiders »

This worked perfectly, thank you!
bschneiders
Posts: 34
Joined: Thu Feb 02, 2017 11:30 am

Re: Different species parameters for different regions

Post by bschneiders »

I have a follow up question to this. I am trying to initialize different calcium concentrations in different RxD regions (i.e. "cyt" and "er" as opposed to different neuron sections). I am following Neymotin et. al to do so ("Neuronal dendrite calcium wave model", 2015) as follows:

Code: Select all

cyt = rxd.Region(h.allsec(), nrn_region='i', geometry=rxd.FractionalVolume(1-er_cyt_ratio, surface_fraction=1))
er = rxd.Region(h.allsec(), geometry=rxd.FractionalVolume(er_cyt_ratio))
cyt_er_membrane = rxd.Region(h.allsec(), geometry=rxd.ScalableBorder(1))

ca =    rxd.Species([cyt, er], d=dCa, name='ca', charge=2, initial=c0)
And it looks as though the ca[er] initial concentration is then updated essentially between h.init() and h.run() - I haven't nailed down the exact timing though, because my code isn't quite running. After I define ca, I set up my fluxes: I have a simplified exchange between cytoplasm and ER, let's say a simple pump in each direction for now:

Code: Select all

Jer_out = vpump_fromER*(ca[er] - ca[cyt])**2  
Jer_in = np.divide(vpump_toER*ca[cyt]**2,ca[cyt]**2 + kp**2)
where the vpump and kp parameters are scalars. Then I use a multi-compartment reaction:

Code: Select all

dCa_er_dt = rxd.MultiCompartmentReaction(ca[er]<>ca[cyt],Jer_out,Jer_in, membrane = cyt_er_membrane) #,custom_dynamics=True)
I get an error at the h.run() line that traces back to rxd.generalizedReaction.py saying "ValueError: operands could not be broadcast together with shapes (276,) (552,)".
276 turns out to be my total number of segments, and 552 is obviously double that. However, I'm having trouble getting deeper than that. Where am I doubling one vector such that two vectors cannot be combined? It seems from the line in generalizedReaction.py (161) that it's an issue with the rates, which made me look at Jer_out and Jer_in, but it's a little opaque because I can only evaluate them before h.run(), at which point they're each of length 1.

In comparing with the three examples within cawave.py in Neymotin et. al, I wonder if it's because I have a bidirectional reaction AND have custom dynamics (you can see I tried "custom_dynamics = True" but commented it out - it looks like 4 is the max for arguments?). In cawave.py, I only see custom dynamics in one direction:

Code: Select all

serca = rxd.MultiCompartmentReaction(ca[cyt]>ca[er],gserca[cyt]*(1e3*ca[cyt])**2/(Kserca**2+(1e3*ca[cyt])**2),membrane=cyt_er_membrane,custom_dynamics=True)
OR bidirectional reactions with equivalent rates in both directions:

Code: Select all

ip3r = rxd.MultiCompartmentReaction(ca[er]<>ca[cyt], k, k, membrane=cyt_er_membrane)
leak = rxd.MultiCompartmentReaction(ca[er]<>ca[cyt], gleak[cyt], gleak[cyt], membrane=cyt_er_membrane)
Could that be related? If so, is there a way around this? Thanks!
bschneiders
Posts: 34
Joined: Thu Feb 02, 2017 11:30 am

Re: Different species parameters for different regions

Post by bschneiders »

A quick update: I tried splitting the Multi-Compartment Reaction into two separate ones so that I could use custom dynamics:

Code: Select all

ERin = rxd.MultiCompartmentReaction(ca[cyt]>ca[er],Jer_in,membrane=cyt_er_membrane,custom_dynamics=True)
ERout = rxd.MultiCompartmentReaction(ca[cyt]<ca[er],Jer_out,membrane=cyt_er_membrane,custom_dynamics=True)
This did not work either, so I think that rules out my last suspicion. I get the same error as before.
adamjhn
Posts: 54
Joined: Tue Apr 18, 2017 10:05 am

Re: Different species parameters for different regions

Post by adamjhn »

I was not able to reproduce the error using the current release NEURON 7.5 or development version NEURON 7.6 https://github.com/neuronsimulator/nrn. We may have already fixed this bug.

Could you try running the following code;

Code: Select all

#Additional code 
from neuron import h, rxd

h.load_file('stdrun.hoc')

sec = h.Section()
sec.nseg = 276

#Arbitrary parameter values 
er_cyt_ratio = 0.5 
dCa = 4.0
c0 = 1.0
vpump_toER = 10.0
vpump_fromER = 1.0
kp = 7.0
Kserca=12.0


#Code you posted on the forum
cyt = rxd.Region(h.allsec(), nrn_region='i', geometry=rxd.FractionalVolume(1-er_cyt_ratio, surface_fraction=1))
er = rxd.Region(h.allsec(), geometry=rxd.FractionalVolume(er_cyt_ratio))
cyt_er_membrane = rxd.Region(h.allsec(), geometry=rxd.ScalableBorder(1))

ca =    rxd.Species([cyt, er], d=dCa, name='ca', charge=2, initial=c0)

Jer_out = vpump_fromER*(ca[er] - ca[cyt])**2
Jer_in = vpump_toER*ca[cyt]**2/(ca[cyt]**2 + kp**2)
# I replaced np.divide with division, this did not change the outcome

dCa_er_dt = rxd.MultiCompartmentReaction(ca[er]<>ca[cyt],Jer_out,Jer_in, membrane = cyt_er_membrane) #,custom_dynamics=True)

h.init()
h.run()
If it gives the same error you could update your version of NEURON. Otherwise could you provide a complete example that gives rise to the ValueError.
bschneiders
Posts: 34
Joined: Thu Feb 02, 2017 11:30 am

Re: Different species parameters for different regions

Post by bschneiders »

Ok, I actually had to table that question (I have a workaround for a conference next week) but will definitely readdress it after to figure out if its my code or a more general neuron issue.

However, I have yet another related follow up question. Given that my spines are quite small, I have seen some sources say that the ER may only make it partway into the neck if at all.

Therefore, I'm trying to implement CICR in the dendrites but not the spines as a first pass. To do so, I need to add the Jer term into the dynamics for dendrite sections but not for spines. I have tried the following:

Code: Select all

dends_region = rxd.Region(all_dends, nrn_region='i')
spines_region = rxd.Region(all_spines, nrn_region='i')
r = rxd.Region(h.allsec(),nrn_region='i') 

ca =    rxd.Species(r, initial=c0,name='ca',charge=2, d = dCa)
ca_eff = rxd.Species(dends_region, initial=ceff0,name='caer', d = dCa) #effective ca_er - workaround for multi compartment for now

Jer = Jout - Jin
dcdt_d = rxd.Rate(ca, Jbuf + pump + Jer, regions = dends_region)
dcadt_s = rxd.Rate(ca, Jbuf + pump, regions = spines_region)
dcaer_dt = rxd.Rate(ca_eff, -Jer, regions = dends_region) 
This didn't work. It ran, and calcium changes, but ca_eff never changes from its initial value. Given that rxd.Rate objects add for effects on the same species, I also tried:

Code: Select all

dcdt = rxd.Rate(ca, Jbuf + pump)
dc_er_dt = rxd.Rate(ca, Jer, regions = dends_region)
dcaer_dt = rxd.Rate(ca_eff, -Jer, regions = dends_region) 
in the hopes that the Jer effect would just add to the dendrites with the second statement. However, I get the same result...ca_eff doesn't budge from its initial value. Am I implementing it incorrectly? Is there a better way to add a term to rate equations in one region but not another?
Post Reply