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!