Using NEURON to model a system of endocrine cells

Extending NEURON to handle reaction-diffusion problems.

Moderators: hines, wwlytton, ramcdougal

Post Reply
Robert Claar
Posts: 23
Joined: Tue Jun 23, 2020 8:52 am

Using NEURON to model a system of endocrine cells

Post by Robert Claar »

Hello,

I am new to NEURON (and programming in general) and I am working with my lab partner to build a model of pancreatic islet cells using previous models not created in the NEURON framework. These previous models didn't have a spatial component and my lab partner and I want to utilize NEURON to add in this feature to the two existing models. Each cell type in the pancreatic islet secretes a hormone and that hormone affects the other cells. In the previous models, there was a single, well-mixed "sink" that held the concentrations of each of the hormones. In our model we want each cell to have it's own "sink" that is comprised of it's own internal hormone concentration, as well as the concentration of hormones secreted from other cells near it. If cell b secretes hormone i then we would like to have the concentration of i be greatest just outside the plasma membrane and decrease as you move radially away from the cell.

We have created mod file mechanisms for the secretion of each of the hormones as well as the ion channels in each of the cell types. This replicated the two models we are basing ours off of well. The next step we're trying to take is adding on the spatial relationships that I described above.

At first we tried to use NetCon because the weighting seemed like it would help us weight the effects of the hormones by distance, but the threshold argument was an issue because we needed a reading of the hormones coming from other cells at each time step, not just when a threshold was crossed. I was wondering if Reaction-Diffusion in NEURON could help us accomplish the individual "sinks" and spatial relationships I described above. If so, how would we set it up or where would we learn how to set it up? Any help is appreciated!
adamjhn
Posts: 54
Joined: Tue Apr 18, 2017 10:05 am

Re: Using NEURON to model a system of endocrine cells

Post by adamjhn »

Here is a simple example of two cells connected by diffusion through the extracellular space.

Code: Select all

from neuron import h, rxd
from neuron.rxd import v
from neuron.units import mV, s 
from matplotlib import pyplot

h.load_file('stdrun.hoc')

# create a couple of cells 
cell1 = h.Section('cell1')
cell1.pt3dclear()
cell1.pt3dadd(-20,0,0,10)
cell1.pt3dadd(-10,0,0,10)

cell2 = h.Section('cell2')
cell2.pt3dclear()
cell2.pt3dadd(10,0,0,10)
cell2.pt3dadd(20,0,0,10)

# Where?
# the intracellular spaces
cyt = rxd.Region(h.allsec(), name='cyt', nrn_region='i')

# plasma membrane 
mem = rxd.Region(h.allsec(), name='mem', geometry=rxd.membrane)

# the extracellular space
ecs = rxd.Extracellular(-35, -15, -15, 35, 15, 15, dx=10, volume_fraction=0.2, tortuosity=1.6)

# Who?
glucagon = rxd.Species([cyt, ecs], name='glucagon', d=1.0, initial=lambda nd: 1 if hasattr(nd, 'sec') and nd.segment in cell1 else 0)
gcyt = glucagon[cyt]
gecs = glucagon[ecs]

# interaction between intracellular and extracellular glucagon
R = 1e3
U = 1e5
rrate = R*gcyt*(-v) # release rate [molecules per square micron per ms]
urate = U*gecs      # uptake rate [molecules per square micron per ms]
glucagon_release = rxd.MultiCompartmentReaction(gcyt, gecs, rrate, urate,
                                                membrane=mem, 
                                                custom_dynamics=True)


# record the concentrations in the cells
t_vec = h.Vector()
t_vec.record(h._ref_t)
gl_cell1 = h.Vector().record(cell1(0.5)._ref_glucagoni)
gl_ecs = h.Vector().record(gecs.node_by_location(0,0,0)._ref_concentration)
gl_cell2 = h.Vector().record(cell2(0.5)._ref_glucagoni)


# run and plot the results
h.finitialize(-65 * mV)
h.continuerun(1 * s)

fig = pyplot.figure(dpi=200)
pyplot.plot(t_vec,gl_cell1, label='cell1')
pyplot.plot(t_vec,gl_ecs, label='ecs')
pyplot.plot(t_vec,gl_cell2, label='cell2')
pyplot.legend(frameon=False)
pyplot.xlabel('t (ms)')
pyplot.ylabel('glucagon (mM)')
pyplot.savefig("/tmp/300620_glucagon_example.png")
pyplot.show()
Image


In this simple example, if Gi is the intracellular glucagon and Go the extracellular glucagon (in the voxel that contains the segment) then release/update mechanism;
Gi’ = D∇²Gi - (R*Gi*(-v) - U*Go) * surface_area/ (intracellular_volume*NA)
Go’ = D∇²Go + (R*Gi*(-v) - U*Go) * surface_area/ (extracellular_volume*volume_fraction*NA)
I made release dependent on voltage, which I saw you mentioned in another post (although in this example the voltage doesn’t change).

It’s similar to an example I posted here, which makes use of a mod file rather than a multicompartment reaction, to make rxd work with the mod file you just have to use the same name and charge for the rxd.Specie and USEION in the mod file.


There are also more rxd examples here https://neuron.yale.edu/neuron/docs/reaction-diffusion
There is an upcoming rxd tutorial on July 18; https://ocns.memberclicks.net/cns-2020-tutorials#T2
Robert Claar
Posts: 23
Joined: Tue Jun 23, 2020 8:52 am

Re: Using NEURON to model a system of endocrine cells

Post by Robert Claar »

Thank you so much, this helps a lot. I will try to create my own script inserting the mechanisms I have created.
Robert Claar
Posts: 23
Joined: Tue Jun 23, 2020 8:52 am

Re: Using NEURON to model a system of endocrine cells

Post by Robert Claar »

Thank you in advance for reading through this monster of a post/question.

I am still working on the same project I posted about above a year ago and am now stuck on something different. We are again using NEURON to recreate and add to a model created using XPPAUT. We have replicated the results from the model we are building upon in NEURON by first starting with a simple setup of putting all equations for all of our cells into one mod file, recording all of our variables, etc. to essentially test NEURON as a differential equation solver and make sure the results matched the ones from the XPPAUT model. This was done successfully.

Our next step will be to leverage NEURON's rxd capabilities to have our hormones (insulin, glucagon, somatostatin) be "created" by rxd statements rather than by statements in the mod file. Using insulin as an example, insulin concentration is governed by this equation in our mod file I'=JIS/vc-fb*I. Where JIS is a range variable that represents the rate at which insulin is secreted into the hormone sink in the XPPAUT model, vc is a constant representing the volume of that sink, and fb represents the rate of insulin clearance from the sink.

To start small and simple, rxd.Rate() seemed like the best place to start since it's assumptions match the assumptions of the XPPAUT model (from rxd tutorial "Reactions conserve mass. This property is especially important for stochastic simulations. Often, however, one might wish to model a source or a sink (e.g. protein degradation), or describe the dynamics of a state variable that is not a concentration. In this case, use the non-conservative rxd.Rate.").

Code: Select all

I_conc = rxd.Rate(Ins, JIS/vc - fb*Ins)
We are attempting to simulate this reaction in extracellular space, initializing insulin to match the initial concentration in the XPPAUT model. Since the JIS variable is a range variable in our mod file, we first tried to access this variable directly and use it in the Rate equation as shown below

Code: Select all

my_cell = h.Section(name = "my_cell")
# one is the name of our one big mod file containing all the models equations
my_cell.insert("one")

# Make extracellular space. Think about rubix cube to visualize its shape.
# 3 rows of 9 10x10x10 cubes
ecs = rxd.Extracellular(0,0,0,30,30,30, dx = (10, 10, 10))

# Leave d = 0 so no diffusion meaning no insulin should flow from one voxel to another
# The concentration curves of all voxels should be the same
# I do not know why the developers of the model use so many decimal places
Ins = rxd.Species(ecs, name = "Ins", d = 0, initial = 48.04298494148047)

# Constants for sink volume and insulin clearance from sink
vc = 1e-13
fb = 2000

# Access JIS value directly to use in Rate reaction
# Set it to initial value in mod file
my_cell(0.5).one.JIS = 9.608596975100945e-09
JIS = my_cell(0.5).one.JIS
Ins_conc = rxd.Rate(Ins, JIS/vc-fb*Ins, regions = ecs)
This did not match the results we got when running a simulation using only our one big mod file and reading the state variable I which is odd because insulin concentration should be governed by the same differential equation in both setups. I'=JIS/vc-fb*I in the mod file and then rxd.Rate(Ins, JIS/vc-fb*Ins, regions = ecs) should generate the same equation based on the documentation for rxd.Rate. We tried changing the integration method in the mod file along with both variable and fixed timestep methods to see if that would affect our results but they were the same for all.

We then tried to use JIS as a NEURON Vector in the rxd.Rate equation by keeping everything the same as above except

Code: Select all

JIS = h.Vector()
JIS.record(my_cell(0.5).one._ref_JIS)

Ins_conc = rxd.Rate(Ins, JIS/vc-fb*Ins, regions = ecs)
And when doing this we get a compilation error that has to do with some of the rxddll files when running an rxd simulation.

"FileNotFoundError: Could not find module ... rxddll3ecdb187-038c-11ec-857b-c858c0244ed2.so' (or one of its dependencies)"

The file rxddll3ecdb187-038c-11ec-857b-c858c0244ed2.c is in my directory, but I couldn't find much info on these rxddll files in the documentation to do any troubleshooting myself. So I have two questions:

1.) Is the use of a recording vector in this way not permitted? If not, is there a better way to access range variables at each time step to use in rxd simulations?

2.) Given that insulin concentration in our sink is governed by the same differential equation in both the mod file and the rxd setup, do you know what could be causing the difference?
Post Reply