From 84aaf322f9d9bd6e05c002f87fe381ccad0c876b Mon Sep 17 00:00:00 2001
From: Dhruva Gowda Storz <dhruvastorz@gmail.com>
Date: Tue, 19 Dec 2017 15:54:13 +0530
Subject: [PATCH] Recovered rdesigneur tutorial. Modified file by adding
extended contents, commented out incomplete sections. Modified rdesigneur
snippets page title to 'more rdesigneur examples' instead of just
'Rdesigneur'
---
docs/source/user/py/cookbook/Rd.rst | 1188 +++++++++++++++++++
docs/source/user/py/cookbook/multi_rdes.rst | 6 +-
docs/source/user/py/cookbook/multiscale.rst | 4 +
3 files changed, 1195 insertions(+), 3 deletions(-)
create mode 100644 docs/source/user/py/cookbook/Rd.rst
diff --git a/docs/source/user/py/cookbook/Rd.rst b/docs/source/user/py/cookbook/Rd.rst
new file mode 100644
index 00000000..04836369
--- /dev/null
+++ b/docs/source/user/py/cookbook/Rd.rst
@@ -0,0 +1,1188 @@
+**Rdesigneur: Building multiscale models**
+==========================================
+
+.. Upi Bhalla
+
+.. Aug 26 2016.
+
+.. --------------
+
+Contents
+--------
+
+.. contents::
+ :depth: 3
+
+Introduction
+------------
+
+**Rdesigneur** (Reaction Diffusion and Electrical SIGnaling in NEURons)
+is an interface to the multiscale modeling capabilities in MOOSE. It is
+designed to build models incorporating biochemical signaling pathways in
+dendrites and spines, coupled to electrical events in neurons.
+Rdesigneur assembles models from predefined parts: it delegates the
+details to specialized model definition formats. Rdesigneur combines one
+or more of the following cell parts to build models:
+
+- Neuronal morphology
+- Dendritic spines
+- Ion channels
+- Reaction systems
+
+Rdesigneur's main role is to specify how these are put together,
+including assigning parameters to do so. Rdesigneur also helps with
+setting up the simulation input and output.
+
+Quick Start
+-----------
+
+Here we provide a few use cases, building up from a minimal model to a
+reasonably complete multiscale model spanning chemical and electrical
+signaling.
+
+Bare Rdesigneur: single passive compartment
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+If we don't provide any arguments at all to the Rdesigneur, it makes a
+model with a single passive electrical compartment in the MOOSE path
+``/model/elec/soma``. Here is how to do this:
+
+::
+
+ import moose
+ import rdesigneur as rd
+ rdes = rd.rdesigneur()
+ rdes.buildModel()
+
+To confirm that it has made a compartment with some default values we
+can add a line:
+
+::
+
+ moose.showfields( rdes.soma )
+
+This should produce the output:
+
+::
+
+ [ /model[0]/elec[0]/soma[0] ]
+ diameter = 0.0005
+ fieldIndex = 0
+ Ra = 7639437.26841
+ y0 = 0.0
+ Rm = 424413.177334
+ index = 0
+ numData = 1
+ inject = 0.0
+ initVm = -0.065
+ Em = -0.0544
+ y = 0.0
+ numField = 1
+ path = /model[0]/elec[0]/soma[0]
+ dt = 0.0
+ tick = -2
+ z0 = 0.0
+ name = soma
+ Cm = 7.85398163398e-09
+ x0 = 0.0
+ Vm = -0.06
+ className = ZombieCompartment
+ idValue = 465
+ length = 0.0005
+ Im = 1.3194689277e-08
+ x = 0.0005
+ z = 0.0
+
+Simulate and display current pulse to soma
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+A more useful script would run and display the model. Rdesigneur can
+help with the stimulus and the plotting. This simulation has the same
+passive compartment, and current is injected as the simulation runs.
+This script displays the membrane potential of the soma as it charges
+and discharges.
+
+::
+
+ import moose
+ import rdesigneur as rd
+ rdes = rd.rdesigneur(
+ stimList = [['soma', '1', '.', 'inject', '(t>0.1 && t<0.2) * 2e-8']],
+ plotList = [['soma', '1', '.', 'Vm', 'Soma membrane potential']],
+ )
+ rdes.buildModel()
+ moose.reinit()
+ moose.start( 0.3 )
+ rdes.display()
+
+The *stimList* defines a stimulus. Each entry has five arguments:
+
+::
+
+ `[region_in_cell, region_expression, moose_object, parameter, expression_string]`
+
+- ``region_in_cell`` specifies the objects to stimulate. Here it is
+ just the soma.
+- ``region_expression`` specifies a geometry based calculation to
+ decide whether to apply the stimulus. The value must be >0 for the
+ stimulus to be present. Here it is just 1. ``moose_object`` specifies
+ the simulation object to operate upon during the stimulus. Here the
+ ``.`` means that it is the soma itself. In other models it might be a
+ channel on the soma, or a synapse, and so on.
+- ``parameter`` specifies the simulation parameter on the moose object
+ that the stimulus will modify. Here it is the injection current to
+ the soma compartment.
+- ``expression_string`` calculates the value of the parameter,
+ typically as a function of time. Here we use the function
+ ``(t>0.1 && t<0.2) * 2e-8`` which evaluates as 2e-8 between the times
+ of 0.1 and 0.2 seconds.
+
+To summarise this, the *stimList* here means *inject a current of 20nA
+to the soma between the times of 0.1 and 0.2 s*.
+
+The *plotList* defines what to plot. It has a similar set of arguments:
+
+::
+
+ `[region_in_cell, region_expression, moose_object, parameter, title_of_plot]`
+
+These mean the same thing as for the stimList except for the title of
+the plot.
+
+The *rdes.display()* function causes the plots to be displayed.
+
+.. figure:: ../../../../images/rdes2_passive_squid.png
+ :alt: Plot for current input to passive compartment
+
+ Plot for current input to passive compartment
+
+When we run this we see an initial depolarization as the soma settles
+from its initial -65 mV to a resting Em = -54.4 mV. These are the
+original HH values, see the example above. At t = 0.1 seconds there is
+another depolarization due to the current injection, and at t = 0.2
+seconds this goes back to the resting potential.
+
+HH Squid model in a single compartment
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Here we put the Hodgkin-Huxley squid model channels into a passive
+compartment. The HH channels are predefined as prototype channels for
+Rdesigneur,
+
+::
+
+ import moose
+ import pylab
+ import rdesigneur as rd
+ rdes = rd.rdesigneur(
+ chanProto = [['make_HH_Na()', 'Na'], ['make_HH_K()', 'K']],
+ chanDistrib = [
+ ['Na', 'soma', 'Gbar', '1200' ],
+ ['K', 'soma', 'Gbar', '360' ]],
+ stimList = [['soma', '1', '.', 'inject', '(t>0.1 && t<0.2) * 1e-8' ]],
+ plotList = [['soma', '1', '.', 'Vm', 'Membrane potential']]
+ )
+
+ rdes.buildModel()
+ moose.reinit()
+ moose.start( 0.3 )
+ rdes.display()
+
+Here we introduce two new model specification lines:
+
+- **chanProto**: This specifies which ion channels will be used in the
+ model. Each entry here has two fields: the source of the channel
+ definition, and (optionally) the name of the channel. In this example
+ we specify two channels, an Na and a K channel using the original
+ Hodgkin-Huxley parameters. As the source of the channel definition we
+ use the name of the Python function that builds the channel. The
+ *make\_HH\_Na()* and *make\_HH\_K()* functions are predefined but we
+ can also specify our own functions for making prototypes. We could
+ also have specified the channel prototype using the name of a channel
+ definition file in ChannelML (a subset of NeuroML) format.
+- **chanDistrib**: This specifies *where* the channels should be placed
+ over the geometry of the cell. Each entry in the chanDistrib list
+ specifies the distribution of parameters for one channel using four
+ entries:
+
+ ``[object_name, region_in_cell, parameter, expression_string]``
+
+ In this case the job is almost trivial, since we just have a single
+ compartment named *soma*. So the line
+
+ ``['Na', 'soma', 'Gbar', '1200' ]``
+
+ means *Put the Na channel in the soma, and set its maximal
+ conductance density (Gbar) to 1200 Siemens/m^2*.
+
+As before we apply a somatic current pulse. Since we now have HH
+channels in the model, this generates action potentials.
+
+.. figure:: ../../../../images/rdes3_squid.png
+ :alt: Plot for HH squid simulation
+
+ Plot for HH squid simulation
+
+HH Squid model in an axon
+~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Here we put the Hodgkin-Huxley squid model into a long compartment that
+is subdivided into many segments, so that we can watch action potentials
+propagate. Most of this example is boilerplate code to build a spiral
+axon. There is a short *rdesigneur* segment that takes the spiral axon
+prototype and populates it with channels, and sets up the display. Later
+examples will show you how to read morphology files to specify the
+neuronal geometry.
+
+::
+
+ import numpy as np
+ import moose
+ import pylab
+ import rdesigneur as rd
+
+ numAxonSegments = 200
+ comptLen = 10e-6
+ comptDia = 1e-6
+ RM = 1.0
+ RA = 10.0
+ CM = 0.01
+
+ def makeAxonProto():
+ axon = moose.Neuron( '/library/axon' )
+ prev = rd.buildCompt( axon, 'soma', RM = RM, RA = RA, CM = CM, dia = 10e-6, x=0, dx=comptLen)
+ theta = 0
+ x = comptLen
+ y = 0.0
+
+ for i in range( numAxonSegments ):
+ dx = comptLen * np.cos( theta )
+ dy = comptLen * np.sin( theta )
+ r = np.sqrt( x * x + y * y )
+ theta += comptLen / r
+ compt = rd.buildCompt( axon, 'axon' + str(i), RM = RM, RA = RA, CM = CM, x = x, y = y, dx = dx, dy = dy, dia = comptDia )
+ moose.connect( prev, 'axial', compt, 'raxial' )
+ prev = compt
+ x += dx
+ y += dy
+
+ return axon
+
+ moose.Neutral( '/library' )
+ makeAxonProto()
+
+ rdes = rd.rdesigneur(
+ chanProto = [['make_HH_Na()', 'Na'], ['make_HH_K()', 'K']],
+ cellProto = [['elec','axon']],
+ chanDistrib = [
+ ['Na', '#', 'Gbar', '1200' ],
+ ['K', '#', 'Gbar', '360' ]],
+ stimList = [['soma', '1', '.', 'inject', '(t>0.01 && t<0.2) * 2e-11' ]],
+ plotList = [['soma', '1', '.', 'Vm', 'Membrane potential']],
+ moogList = [['#', '1', '.', 'Vm', 'Vm (mV)']]
+ )
+
+ rdes.buildModel()
+ moose.reinit()
+
+ rdes.displayMoogli( 0.00005, 0.05, 0.0 )
+
+.. figure:: ../../../../images/rdes3.1_axon.png
+ :alt: Axon with propagating action potential
+
+ Axon with propagating action potential
+
+Note how we explicitly create the prototype axon on '/library', and then
+specify it using the *cellProto* line in the rdesigneur. The moogList
+specifies the 3-D display. See below for how to set up and use these
+displays.
+
+HH Squid model in a myelinated axon
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+This is a curious cross-species chimera model, where we embed the HH
+equations into a myelinated example model. As for the regular axon
+above, most of the example is boilerplate setup code. Note how we
+restrict the HH channels to the nodes of Ranvier using a conditional
+test for the diameter of the axon segment.
+
+::
+
+ import numpy as np
+ import moose
+ import pylab
+ import rdesigneur as rd
+
+ numAxonSegments = 405
+ nodeSpacing = 100
+ comptLen = 10e-6
+ comptDia = 2e-6 # 2x usual
+ RM = 100.0 # 10x usual
+ RA = 5.0
+ CM = 0.001 # 0.1x usual
+
+ nodeDia = 1e-6
+ nodeRM = 1.0
+ nodeCM = 0.01
+
+ def makeAxonProto():
+ axon = moose.Neuron( '/library/axon' )
+ x = 0.0
+ y = 0.0
+ prev = rd.buildCompt( axon, 'soma', RM = RM, RA = RA, CM = CM, dia = 10e-6, x=0, dx=comptLen)
+ theta = 0
+ x = comptLen
+
+ for i in range( numAxonSegments ):
+ r = comptLen
+ dx = comptLen * np.cos( theta )
+ dy = comptLen * np.sin( theta )
+ r = np.sqrt( x * x + y * y )
+ theta += comptLen / r
+ if i % nodeSpacing == 0:
+ compt = rd.buildCompt( axon, 'axon' + str(i), RM = nodeRM, RA = RA, CM = nodeCM, x = x, y = y, dx = dx, dy = dy, dia = nodeDia )
+ else:
+ compt = rd.buildCompt( axon, 'axon' + str(i), RM = RM, RA = RA, CM = CM, x = x, y = y, dx = dx, dy = dy, dia = comptDia )
+ moose.connect( prev, 'axial', compt, 'raxial' )
+ prev = compt
+ x += dx
+ y += dy
+
+ return axon
+
+ moose.Neutral( '/library' )
+ makeAxonProto()
+
+ rdes = rd.rdesigneur(
+ chanProto = [['make_HH_Na()', 'Na'], ['make_HH_K()', 'K']],
+ cellProto = [['elec','axon']],
+ chanDistrib = [
+ ['Na', '#', 'Gbar', '12000 * (dia < 1.5e-6)' ],
+ ['K', '#', 'Gbar', '3600 * (dia < 1.5e-6)' ]],
+ stimList = [['soma', '1', '.', 'inject', '(t>0.01 && t<0.2) * 1e-10' ]],
+ plotList = [['soma,axon100,axon200,axon300,axon400', '1', '.', 'Vm', 'Membrane potential']],
+ moogList = [['#', '1', '.', 'Vm', 'Vm (mV)']]
+ )
+
+ rdes.buildModel()
+
+ for i in moose.wildcardFind( "/model/elec/#/Na" ):
+ print i.parent.name, i.Gbar
+
+ moose.reinit()
+
+ rdes.displayMoogli( 0.00005, 0.05, 0.0 )
+
+When you run the example, keep an eye out for a few things:
+
+- **saltatory conduction:** This is the way the action potential jumps
+ from one node of Ranvier to the next. Between the nodes it is just
+ passive propagation.
+- **Failure to propagate:** Observe that the second and fourth action
+ potentials fails to trigger propagation along the axon. Here we have
+ specially tuned the model properties so that this happens. With a
+ larger RA of 10.0, the model will be more reliable.
+- **Speed:** Compare the propagation speed with the previous,
+ unmyelinated axon. Note that the current model is larger!
+
+.. figure:: ../../../../images/rdes3.2_myelinated_axon.png
+ :alt: Myelinated axon with propagating action potential
+
+ Myelinated axon with propagating action potential
+
+Reaction system in a single compartment
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Here we use the compartment as a place in which to embed a chemical
+model. The chemical oscillator model is predefined in the rdesigneur
+prototypes.
+
+::
+
+ import moose
+ import pylab
+ import rdesigneur as rd
+ rdes = rd.rdesigneur(
+ turnOffElec = True,
+ diffusionLength = 1e-3, # Default diffusion length is 2 microns
+ chemProto = [['makeChemOscillator()', 'osc']],
+ chemDistrib = [['osc', 'soma', 'install', '1' ]],
+ plotList = [['soma', '1', 'dend/a', 'conc', 'a Conc'],
+ ['soma', '1', 'dend/b', 'conc', 'b Conc']]
+ )
+ rdes.buildModel()
+ b = moose.element( '/model/chem/dend/b' )
+ b.concInit *= 5
+ moose.reinit()
+ moose.start( 200 )
+
+ rdes.display()
+
+In this special case we set the turnOffElec flag to True, so that
+Rdesigneur only sets up chemical and not electrical calculations. This
+makes the calculations much faster, since we disable electrical
+calculations and delink chemical calculations from them.
+
+We also have a line which sets the ``diffusionLength`` to 1 mm, so that
+it is bigger than the 0.5 mm squid axon segment in the default
+compartment. If you don't do this the system will subdivide the
+compartment into the default 2 micron voxels for the purposes of putting
+in a reaction-diffusion system. We discuss this case below.
+
+Note how the *plotList* is done here. To remind you, each entry has five
+arguments
+
+::
+
+ [region_in_cell, region_expression, moose_object, parameter, title_of_plot]
+
+The change from the earlier usage is that the ``moose_object`` now
+refers to a chemical entity, in this example the molecule *dend/a*. The
+simulator builds a default chemical compartment named *dend* to hold the
+reactions defined in the *chemProto*. What we do in this plot is to
+select molecule *a* sitting in *dend*, and plot its concentration. Then
+we do this again for molecule *b*.
+
+After the model is built, we add a couple of lines to change the initial
+concentration of the molecular pool *b*. Note its full path within
+MOOSE: */model/chem/dend/b*. It is scaled up 5x to give rise to slowly
+decaying oscillations.
+
+.. figure:: ../../../../images/rdes4_osc.png
+ :alt: Plot for single-compartment reaction simulation
+
+ Plot for single-compartment reaction simulation
+
+Reaction-diffusion system
+~~~~~~~~~~~~~~~~~~~~~~~~~
+
+In order to see what a reaction-diffusion system looks like, delete the
+``diffusionLength`` expression in the previous example and add a couple
+of lines to set up 3-D graphics for the reaction-diffusion product:
+
+::
+
+ import moose
+ import pylab
+ import rdesigneur as rd
+ rdes = rd.rdesigneur(
+ turnOffElec = True,
+ chemProto = [['makeChemOscillator()', 'osc']],
+ chemDistrib = [['osc', 'soma', 'install', '1' ]],
+ plotList = [['soma', '1', 'dend/a', 'conc', 'Concentration of a'],
+ ['soma', '1', 'dend/b', 'conc', 'Concentration of b']],
+ moogList = [['soma', '1', 'dend/a', 'conc', 'a Conc', 0, 360 ]]
+ )
+
+ rdes.buildModel()
+ bv = moose.vec( '/model/chem/dend/b' )
+ bv[0].concInit *= 2
+ bv[-1].concInit *= 2
+ moose.reinit()
+
+ rdes.displayMoogli( 1, 400, 0.001 )
+
+This is the line we deleted.
+
+::
+
+ diffusionLength = 1e-3,
+
+With this change we permit *rdesigneur* to use the default diffusion
+length of 2 microns. The 500-micron axon segment is now subdivided into
+250 voxels, each of which has a reaction system and diffusing molecules.
+To make it more picturesque, we have added a line after the plotList, to
+display the outcome in 3-D:
+
+::
+
+ moogList = [['soma', '1', 'dend/a', 'conc', 'a Conc', 0, 360 ]]
+
+This line says: take the model compartments defined by ``soma`` as the
+region to display, do so throughout the the geometry (the ``1``
+signifies this), and over this range find the chemical entity defined by
+``dend/a``. For each ``a`` molecule, find the ``conc`` and dsiplay it.
+There are two optional arguments, ``0`` and ``360``, which specify the
+low and high value of the displayed variable.
+
+In order to initially break the symmetry of the system, we change the
+initial concentration of molecule b at each end of the cylinder:
+
+::
+
+ bv[0].concInit *= 2
+ bv[-1].concInit *= 2
+
+If we didn't do this the entire system would go through a few cycles of
+decaying oscillation and then reach a boring, spatially uniform, steady
+state. Try putting an initial symmetry break elsewhere to see what
+happens.
+
+To display the concenctration changes in the 3-D soma as the simulation
+runs, we use the line
+
+::
+
+ `rdes.displayMoogli( 1, 400, 0.001 )`
+
+The arguments mean: *displayMoogli( frametime, runtime, rotation )*
+Here,
+
+::
+
+ frametime = time by which simulation advances between display updates
+ runtime = Total simulated time
+ rotation = angle by which display rotates in each frame, in radians.
+
+When we run this, we first get a 3-D display with the oscillating
+reaction-diffusion system making its way inward from the two ends. After
+the simulation ends the plots for all compartments for the whole run
+come up.
+
+.. figure:: ../../../../images/rdes5_reacdiff.png
+ :alt: Display for oscillatory reaction-diffusion simulation
+
+ Display for oscillatory reaction-diffusion simulation
+
+Primer on using the 3-D MOOGLI display
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Here is a short primer on the 3-D display controls.
+
+- *Roll, pitch, and yaw*: Use the letters *r*, *p*, and *y*. To rotate
+ backwards, use capitals.
+- *Zoom out and in*: Use the *,* and *.* keys, or their upper-case
+ equivalents, *<* and *>*. Easier to remember if you think in terms of
+ the upper-case.
+- *Left/right/up/down*: Arrow keys.
+- *Quit*: control-q or control-w.
+- You can also use the mouse or trackpad to control most of the above.
+- By default rdesigneur gives Moogli a small rotation each frame. It is
+ the *rotation* argument in the line:
+
+ ``displayMoogli( frametime, runtime, rotation )``
+
+These controls operate over and above this rotation, but the rotation
+continues. If you set the rotation to zero you can, with a suitable
+flick of the mouse, get the image to rotate in any direction you choose
+as long as the window is updating.
+
+Make a toy multiscale model with electrical and chemical signaling.
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Now we put together chemical and electrical models. In this toy model we
+have an HH-squid type single compartment electrical model, cohabiting
+with a chemical oscillator. The chemical oscillator regulates K+ channel
+amounts, and the average membrane potential regulates the amounts of a
+reactant in the chemical oscillator. This is a recipe for some strange
+firing patterns.
+
+::
+
+ import moose
+ import pylab
+ import rdesigneur as rd
+ rdes = rd.rdesigneur(
+ # We want just one compartment so we set diffusion length to be
+ # bigger than the 0.5 mm HH axon compartment default.
+ diffusionLength = 1e-3,
+ chanProto = [['make_HH_Na()', 'Na'], ['make_HH_K()', 'K']],
+ chanDistrib = [
+ ['Na', 'soma', 'Gbar', '1200' ],
+ ['K', 'soma', 'Gbar', '360' ]],
+ chemProto = [['makeChemOscillator()', 'osc']],
+ chemDistrib = [['osc', 'soma', 'install', '1' ]],
+ # These adaptor parameters give interesting-looking but
+ # not particularly physiological behaviour.
+ adaptorList = [
+ [ 'dend/a', 'conc', 'Na', 'modulation', 1, -5.0 ],
+ [ 'dend/b', 'conc', 'K', 'modulation', 1, -0.2],
+ [ 'dend/b', 'conc', '.', 'inject', -1.0e-7, 4e-7 ],
+ [ '.', 'Vm', 'dend/s', 'conc', 2.5, 20.0 ]
+ ],
+ plotList = [['soma', '1', 'dend/a', 'conc', 'a Conc'],
+ ['soma', '1', 'dend/b', 'conc', 'b Conc'],
+ ['soma', '1', 'dend/s', 'conc', 's Conc'],
+ ['soma', '1', 'Na', 'Gk', 'Na Gk'],
+ ['soma', '1', '.', 'Vm', 'Membrane potential']
+ ]
+ )
+
+ rdes.buildModel()
+ moose.reinit()
+ moose.start( 250 ) # Takes a few seconds to run this.
+
+ rdes.display()
+
+We've already modeled the HH squid model and the oscillator
+individually, and you should recognize the parts of those models above.
+The new section that makes this work the *adaptorList* which specifies
+how the electrical and chemical parts talk to each other. This entirely
+fictional set of interactions goes like this:
+
+::
+
+ [ 'dend/a', 'conc', 'Na', 'modulation', 1, -5.0 ]
+
+- *dend/a*: The originating variable comes from the 'a' pool on the
+ 'dend' compartment.
+
+ *conc*: This is the originating variable name on the 'a' pool.
+
+ *Na*: This is the target variable
+
+ *modulation*: scale the Gbar of Na up and down. Use 'modulation'
+ rather than direct assignment of Gbar since Gbar is different for
+ each differently-sized compartment.
+
+ *1*: This is the initial offset
+
+ *-5.0*: This is the scaling from the input to the parameter updated
+ in the simulation.
+
+A similar set of adaptor entries couple the molecule *dend/b* to the K
+channel, *dend/b* again to the current injection into the soma, and the
+membrane potential to the concentration of *dend/s*.
+
+.. figure:: ../../../../images/rdes6_multiscale.png
+ :alt: Plot for toy multiscale model
+
+ Plot for toy multiscale model
+
+Morphology: Load .swc morphology file and view it
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Here we build a passive model using a morphology file in the .swc file
+format (as used by NeuroMorpho.org). The morphology file is predefined
+for Rdesigneur and resides in the directory ``./cells``. We apply a
+somatic current pulse, and view the somatic membrane potential in a
+plot, as before. To make things interesting we display the morphology in
+3-D upon which we represent the membrane potential as colors.
+
+::
+
+ import moose
+ import rdesigneur as rd
+ rdes = rd.rdesigneur(
+ cellProto = [['./cells/h10.CNG.swc', 'elec']],
+ stimList = [['soma', '1', '.', 'inject', 't * 25e-9' ]],
+ plotList = [['#', '1', '.', 'Vm', 'Membrane potential'],
+ ['#', '1', 'Ca_conc', 'Ca', 'Ca conc (uM)']],
+ moogList = [['#', '1', '.', 'Vm', 'Soma potential']]
+ )
+
+ rdes.buildModel()
+
+ moose.reinit()
+ rdes.displayMoogli( 0.0002, 0.1 )
+
+Here the new concept is the cellProto line, which loads in the specified
+cell model:
+
+::
+
+ `[ filename, cellname ]`
+
+The system recognizes the filename extension and builds a model from the
+swc file. It uses the cellname **elec** in this example.
+
+We use a similar line as in the reaction-diffusion example, to build up
+a Moogli display of the cell model:
+
+::
+
+ `moogList = [['#', '1', '.', 'Vm', 'Soma potential']]`
+
+Here we have:
+
+::
+
+ *#*: the path to use for selecting the compartments to display.
+ This wildcard means use all compartments.
+ *1*: The expression to use for the compartments. Again, `1` means use
+ all of them.
+ *.*: Which object in the compartment to display. Here we are using the
+ compartment itself, so it is just a dot.
+ *Vm*: Field to display
+ *Soma potential*: Title for display.
+
+.. figure:: ../../../../images/rdes7_passive.png
+ :alt: 3-D display for passive neuron
+
+ 3-D display for passive neuron
+
+Build an active neuron model by putting channels into a morphology file
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+We load in a morphology file and distribute voltage-gated ion channels
+over the neuron. Here the voltage-gated channels are obtained from a
+number of channelML files, located in the ``./channels`` subdirectory.
+Since we have a spatially extended neuron, we need to specify the
+spatial distribution of channel densities too.
+
+::
+
+ import moose
+ import rdesigneur as rd
+ rdes = rd.rdesigneur(
+ chanProto = [
+ ['./chans/hd.xml'],
+ ['./chans/kap.xml'],
+ ['./chans/kad.xml'],
+ ['./chans/kdr.xml'],
+ ['./chans/na3.xml'],
+ ['./chans/nax.xml'],
+ ['./chans/CaConc.xml'],
+ ['./chans/Ca.xml']
+ ],
+ cellProto = [['./cells/h10.CNG.swc', 'elec']],
+ chanDistrib = [ \
+ ["hd", "#dend#,#apical#", "Gbar", "50e-2*(1+(p*3e4))" ],
+ ["kdr", "#", "Gbar", "p < 50e-6 ? 500 : 100" ],
+ ["na3", "#soma#,#dend#,#apical#", "Gbar", "850" ],
+ ["nax", "#soma#,#axon#", "Gbar", "1250" ],
+ ["kap", "#axon#,#soma#", "Gbar", "300" ],
+ ["kap", "#dend#,#apical#", "Gbar",
+ "300*(H(100-p*1e6)) * (1+(p*1e4))" ],
+ ["Ca_conc", "#", "tau", "0.0133" ],
+ ["kad", "#soma#,#dend#,#apical#", "Gbar", "50" ],
+ ["Ca", "#", "Gbar", "50" ]
+ ],
+ stimList = [['soma', '1', '.', 'inject', '(t>0.02) * 1e-9' ]],
+ plotList = [['#', '1', '.', 'Vm', 'Membrane potential'],
+ ['#', '1', 'Ca_conc', 'Ca', 'Ca conc (uM)']],
+ moogList = [['#', '1', 'Ca_conc', 'Ca', 'Calcium conc (uM)', 0, 120],
+ ['#', '1', '.', 'Vm', 'Soma potential']]
+ )
+
+ rdes.buildModel()
+
+ moose.reinit()
+ rdes.displayMoogli( 0.0002, 0.052 )
+
+Here we make more extensive use of two concepts which we've already seen
+from the single compartment squid model:
+
+1. *chanProto*: This defines numerous channels, each of which is of the
+ form:
+
+ ``[ filename ]``
+
+ or
+
+ ``[ filename, channelname ]``
+
+If the *channelname* is not specified the system uses the last part of
+the channel name, before the filetype suffix.
+
+2. *chanDistrib*: This defines the spatial distribution of each channel
+ type. Each line is of a form that should be familiar now:
+
+ ``[channelname, region_in_cell, parameter, expression_string]``
+
+- The *channelname* is the name of the prototype from *chanproto*. This
+ is usually an ion channel, but in the example above you can also see
+ a calcium concentration pool defined.
+- The *region\_in\_cell* is typically defined using wildcards, so that
+ it generalizes to any cell morphology. For example, the plain
+ wildcard ``#`` means to consider all cell compartments. The wildcard
+ ``#dend#`` means to consider all compartments with the string
+ ``dend`` somewhere in the name. Wildcards can be comma-separated, so
+ ``#soma#,#dend#`` means consider all compartments with either soma or
+ dend in their name. The naming in MOOSE is defined by the model file.
+ Importantly, in **.swc** files MOOSE generates names that respect the
+ classification of compartments into axon, soma, dendrite, and apical
+ dendrite compartments respectively. SWC files generate compartment
+ names such as:
+
+ ::
+
+ soma_<number>
+ dend_<number>
+ apical_<number>
+ axon_<number>
+
+where the number is automatically assigned by the reader. In order to
+select all dendritic compartments, for example, one would use *"#dend#"*
+where the *"#"* acts as a wildcard to accept any string. - The
+*parameter* is usually Gbar, the channel conductance density in *S/m^2*.
+If *Gbar* is zero or less, then the system economizes by not
+incorporating this channel mechanism in this part of the cell.
+Similarly, for calcium pools, if the *tau* is below zero then the
+calcium pool object is simply not inserted into this part of the cell. -
+The *expression\_string* defines the value of the parameter, such as
+Gbar. This is typically a function of position in the cell. The
+expression evaluator knows about several parameters of cell geometry.
+All units are in metres:
+
+- *x*, *y* and *z* coordinates.
+- *g*, the geometrical distance from the soma
+- *p*, the path length from the soma, measured along the dendrites.
+- *dia*, the diameter of the dendrite.
+- *L*, The electrotonic length from the soma (no units).
+
+Along with these geometrical arguments, we make liberal use of the
+Heaviside function H(x) to set up the channel distributions. The
+expression evaluator also knows about pretty much all common algebraic,
+trignometric, and logarithmic functions, should you wish to use these.
+
+Also note the two Moogli displays. The first is the calcium
+concentration. The second is the membrane potential in each compartment.
+Easy!
+
+.. figure:: ../../../../images/rdes8_active.png
+ :alt: 3-D display for active neuron
+
+ 3-D display for active neuron
+
+Build a spiny neuron from a morphology file and put active channels in it.
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+This model is one step elaborated from the previous one, in that we now
+also have dendritic spines. MOOSE lets one decorate a bare neuronal
+morphology file with dendritic spines, specifying various geometric
+parameters of their location. As before, we use an swc file for the
+morphology, and the same ion channels and distribution.
+
+::
+
+ import moose
+ import pylab
+ import rdesigneur as rd
+ rdes = rd.rdesigneur(
+ chanProto = [
+ ['./chans/hd.xml'],
+ ['./chans/kap.xml'],
+ ['./chans/kad.xml'],
+ ['./chans/kdr.xml'],
+ ['./chans/na3.xml'],
+ ['./chans/nax.xml'],
+ ['./chans/CaConc.xml'],
+ ['./chans/Ca.xml']
+ ],
+ cellProto = [['./cells/h10.CNG.swc', 'elec']],
+ spineProto = [['makeActiveSpine()', 'spine']],
+ chanDistrib = [
+ ["hd", "#dend#,#apical#", "Gbar", "50e-2*(1+(p*3e4))" ],
+ ["kdr", "#", "Gbar", "p < 50e-6 ? 500 : 100" ],
+ ["na3", "#soma#,#dend#,#apical#", "Gbar", "850" ],
+ ["nax", "#soma#,#axon#", "Gbar", "1250" ],
+ ["kap", "#axon#,#soma#", "Gbar", "300" ],
+ ["kap", "#dend#,#apical#", "Gbar",
+ "300*(H(100-p*1e6)) * (1+(p*1e4))" ],
+ ["Ca_conc", "#", "tau", "0.0133" ],
+ ["kad", "#soma#,#dend#,#apical#", "Gbar", "50" ],
+ ["Ca", "#", "Gbar", "50" ]
+ ],
+ spineDistrib = [['spine', '#dend#,#apical#', '20e-6', '1e-6']],
+ stimList = [['soma', '1', '.', 'inject', '(t>0.02) * 1e-9' ]],
+ plotList = [['#', '1', '.', 'Vm', 'Membrane potential'],
+ ['#', '1', 'Ca_conc', 'Ca', 'Ca conc (uM)']],
+ moogList = [['#', '1', 'Ca_conc', 'Ca', 'Calcium conc (uM)', 0, 120],
+ ['#', '1', '.', 'Vm', 'Soma potential']]
+ )
+
+ rdes.buildModel()
+
+ moose.reinit()
+ rdes.displayMoogli( 0.0002, 0.023 )
+
+Spines are set up in a familiar way: we first define one (or more)
+prototype spines, and then distribute these around the cell. Here is the
+prototype string:
+
+::
+
+ [spine_proto, spinename]
+
+*spineProto*: This is typically a function. One can define one's own,
+but there are several predefined ones in rdesigneur. All these define a
+spine with the following parameters:
+
+- head diameter 0.5 microns
+- head length 0.5 microns
+- shaft length 1 micron
+- shaft diameter of 0.2 microns
+- RM = 1.0 ohm-metre square
+- RA = 1.0 ohm-meter
+- CM = 0.01 Farads per square metre.
+
+Here are the predefined spine prototypes:
+
+- *makePassiveSpine()*: This just makes a passive spine with the
+ default parameters
+- *makeExcSpine()*: This makes a spine with NMDA and glu receptors, and
+ also a calcium pool. The NMDA channel feeds the Ca pool.
+- *makeActiveSpine()*: This adds a Ca channel to the exc\_spine. and
+ also a calcium pool.
+
+The spine distributions are specified in a familiar way for the first
+few arguments, and then there are multiple (optional) spine-specific
+parameters:
+
+*[spinename, region\_in\_cell, spacing, spacing\_distrib, size,
+size\_distrib, angle, angle\_distrib ]*
+
+Only the first two arguments are mandatory.
+
+- *spinename*: The prototype name
+- *region\_in\_cell*: Usual wildcard specification of names of
+ compartments in which to put the spines.
+- *spacing*: Math expression to define spacing between spines. In the
+ current implementation this evaluates to
+ ``1/probability_of_spine_per_unit_length``. Defaults to 10 microns.
+ Thus, there is a 10% probability of a spine insertion in every
+ micron. This evaluation method has the drawback that it is possible
+ to space spines rather too close to each other. If spacing is zero or
+ less, no spines are inserted.
+- *spacing\_distrib*: Math expression for distribution of spacing. In
+ the current implementation, this specifies the interval at which the
+ system samples from the spacing probability above. Defaults to 1
+ micron.
+- *size*: Linear scale factor for size of spine. All dimensions are
+ scaled by this factor. The default spine head here is 0.5 microns in
+ diameter and length. If the scale factor were to be 2, the volume
+ would be 8 times as large. Defaults to 1.0.
+- *size\_distrib*: Range for size of spine. A random number R is
+ computed in the range 0 to 1, and the final size used is
+ ``size + (R - 0.5) * size_distrib``. Defaults to 0.5
+- *angle*: This specifies the initial angle at which the spine sticks
+ out of the dendrite. If all angles were zero, they would all point
+ away from the soma. Defaults to 0 radians.
+- *angle\_distrib*: Specifies a random number to add to the initial
+ angle. Defaults to 2 PI radians, so the spines come out in any
+ direction.
+
+One may well ask why we are not using a Python dictionary to handle all
+these parameters. Short answer is: terseness. Longer answer is that the
+rdesigneur format is itself meant to be an intermediate form for an
+eventual high-level, possibly XML-based multiscale modeling format.
+
+.. figure:: ../../../../images/rdes9_spiny_active.png
+ :alt: 3-D display for spiny active neuron
+
+ 3-D display for spiny active neuron
+
+Put a spine on a cylindrical compartment, give it synaptic input, and watch Ca diffuse.
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Calcium enters spines during strong synaptic input, how far does it
+spread along the dendrites? This model is simple conceptually but
+illustrates several things:
+
+- Setting up random (Poisson) synaptic input to one or more spines
+- Setting up a reaction-diffusion system (Ca) coupled to an electrical
+ system (Ca influx through channels). This uses a separate chemical
+ definition script. One can replace this with more elaborate chemical
+ schemes simply by changing the name of the script.
+- User definitions of prototypes (the soma in this example).
+- Assigning random number seeds.
+- The difference between electrical and chemical length scales. For
+ numerical reasons, the discretization of reaction-diffusion systems
+ into voxels normally happens on a smaller length scale (microns) than
+ for electrical systems (tens to hundreds of microns). In this example
+ there is just one electrical compartment, but 50 chemical
+ subdivisions.
+
+Most of the script is setting up the input and the prototypes.
+Rdesigneur is compact, as usual. First, the headers and parameter list:
+
+::
+
+ import moose
+ import numpy as np
+ import rdesigneur as rd
+
+ params = {
+ 'diffusionLength':1.0e-6, # Diffusion characteristic length, used as voxel length too.
+ 'dendDiameter': 1e-6, # Diameter of section of dendrite in model
+ 'dendLength': 50e-6, # Length of section of dendrite in model
+ 'spineSpacing': 30e-6, # mean spacing between spines.
+ 'diffConstCa':20.0e-12, # Diffusion constant of Ca, m^2/sec
+ 'spineFreq': 1, # Frequencey of input to spines
+ 'gluWeight': 100, # Weight for glutamate receptor
+ 'nmdaWeight': 100, # weight for NMDA receptor
+ 'chemModel':'spineCa_diffn.g', # Chem model definition.
+ 'RA': 1.0, # Axial resistivity of compartment, ohms.metre
+ 'RM': 1.0, # membrane resistivity of compartment, ohms.metre^2
+ 'CM': 0.01, # Specific capacitance of membrane, Farads/metre^2
+ 'runtime': 3, # Simulation run time, sec.
+ }
+
+Then, we define the prototypes for the soma compartment and the CaConc
+object that handles conversion of calcium current into calcium
+concentration:
+
+::
+
+ def makePassiveSoma( name, length, diameter ):
+ elecid = moose.Neuron( '/library/' + name )
+ dend = moose.Compartment( elecid.path + '/soma' )
+ dend.diameter = diameter
+ dend.length = length
+ dend.Ra = params['RA'] * length * 4.0 / (diameter * diameter * np.pi)
+ dend.Rm = params['RM'] / (length * diameter * np.pi)
+ dend.Cm = params['CM'] * length * diameter * np.pi
+ dend.x = length
+ return elecid
+
+ def makeCaConc( name ):
+ conc = moose.CaConc( '/library/' + name )
+ conc.tau = 0.0133333
+ conc.B = 17.402e12 # Conversion from Amps to milliMolar for soma
+ conc.Ca_base = 0.0
+
+Then, we define the stimulus including the poisson spike generator:
+
+::
+
+ def attachStimulus():
+ numSpine = len( moose.wildcardFind( '/model/elec/head#' ) )
+ spikeInput = moose.RandSpike( '/model/elec/spineInput', numSpine )
+ spikeVec = spikeInput.vec
+ spikeVec.rate = params['spineFreq']
+
+ j = 0
+ for i in moose.wildcardFind( '/model/elec/head#' ):
+ sh = moose.element( i.path + '/glu/sh' )
+ sh.numSynapses = 1
+ sh.synapse[0].weight = params['gluWeight']
+ moose.connect( spikeVec[j], 'spikeOut', sh.synapse[0], 'addSpike')
+ sh = moose.element( i.path + '/NMDA/sh' )
+ sh.numSynapses = 1
+ sh.synapse[0].weight = params['nmdaWeight']
+ moose.connect( spikeVec[j], 'spikeOut', sh.synapse[0], 'addSpike')
+ j += 1
+
+Having specified the framework for the model, here is the actual
+rdesigneur setup:
+
+::
+
+ moose.seed( 123 ) # One seed for the layout
+ library = moose.Neutral( '/library' )
+ makePassiveSoma( 'cell', params['dendLength'], params['dendDiameter'] )
+ makeCaConc( 'Ca_conc' )
+
+ rdes = rd.rdesigneur(
+ chemPlotDt = 0.02,
+ diffusionLength = params['diffusionLength'],
+ spineProto = [['makeExcSpine()', 'spine']],
+ spineDistrib = [['spine', '#', str( params['spineSpacing'] ),'1e-7']],
+ chanDistrib = [["Ca_conc", "#", "tau", "0.0133", "thick","0.1e-6"]],
+ cellProto = [['cell', 'elec']],
+ chemProto = [['../chem/' + params['chemModel'], 'chem']],
+ chemDistrib = [['chem', '#soma#', 'install', '1' ]],
+ plotList = [
+ ['soma', '1', '.', 'Vm', 'soma Vm'],
+ ['soma', '1', 'dend/DEND/Ca', 'conc', '[dend Ca]'],
+ ['#head#', '1', 'spine/Ca', 'conc', '[Spine Ca]'],
+ ['#head#', '1', 'psd/Ca', 'conc', '[PSD Ca]'],
+ ],
+ moogList = [['#', '1', 'dend/DEND/Ca', 'conc', 'dend Ca', 0, 0.5]],
+ adaptorList = [
+ [ 'Ca_conc', 'Ca', 'psd/Ca_input', 'concInit', 2e-6, 0.1 ],
+ [ 'Ca_conc', 'Ca','dend/DEND/Ca_input','concInit',2e-6,0.001],
+ ]
+ )
+ for ca in moose.wildcardFind( '/library/##/Ca' ):
+ ca.diffConst = params['diffConstCa']
+ rdes.buildModel()
+ attachStimulus()
+ moose.reinit()
+ moose.seed( 3 ) # Another seed because the reinit does a reseed.
+ rdes.displayMoogli( 0.01, params['runtime'], 0.0 )
+
+You will additionally need to copy over the chemical models for the
+calcium. These reside in the moose-examples/genesis directory. The
+simple model ``spineCa_diffn.g`` has calcium in spine PSD, spine head,
+and dendrite pools, with reactions for controlling input from the
+*Ca\_conc* object.
+
+With this done you can run the script and watch calcium spreading from
+the location of the spine. There are three pulses of calcium, the first
+being quite weak.
+
+.. figure:: ../../../../images/rdes10_CaSpread.png
+ :alt: Calcium influx from spine in middle of cylindrical compartment,
+ spreading into the dendrite and then axially in both directions.
+
+ Calcium influx from spine in middle of cylindrical compartment,
+ spreading into the dendrite and then axially in both directions.
+
+Once the simulation completes you'll also see a number of plots, so that
+you can figure out what the calcium influx was doing. Here, we have a
+single spine so there is just one trace for it. We have 50 chemical
+voxels along the dendrite, so there are 50 traces for the chemical
+time-course.
+
+.. figure:: ../../../../images/rdes10_CaTimecourse.png
+ :alt: Time-course of Calcium buildup in spine and dendrite.
+
+ Time-course of Calcium buildup in spine and dendrite.
+
+Note the explicit assignment of random seeds using
+``moose.seed( 123 )``. By default, MOOSE generates a reasonably random
+seed using system information. Here, however, we want to be sure that
+the simulation always gives the same result. So we explicitly set the
+seed to a known number. Note also that we set the seed at two places:
+First, before setup, so that we ensure that the spines are going to come
+in the same number and place each time. Second, after ``moose.reinit()``
+to make sure that the same pseudo-random sequence of synaptic input and
+chemical stochastic calculations happens each run.
+
+Note also that this is run using stochastic methods for the chemical
+calculations. This is the default. One can alter this using the
+following line in rdesigneur:
+
+::
+
+ useGssa = False
+
+As an exercise for the user, we also have a plug-in replaceable model
+for Ca influx, this one having a lot of calmodulin to act as a buffer.
+Replace the original line in the params dictionary as follows:
+
+::
+
+ 'chemModel':'spineCa_diffn.g'
+
+with
+
+::
+
+ 'chemModel':'spineCa_CaM_diffn.g'
+
+The resultant model has reduced free Ca++ buildup.
+
+Some other interesting things to try are:
+
+- Increase the diffusion constant for calcium. You might expect that
+ this would lead to faster flow of Ca from the spine to the dendrite,
+ and hence a higher peak in the dendrite. Try it and see.
+- Change the diameter of the dendrite.
+
+.. Build a spiny neuron from a morphology file and put a reaction-diffusion system in it.
+.. ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+.. Rdesigneur is specially designed to take reaction systems with a
+ dendrite, a spine head, and a spine PSD compartment, and embed these
+ systems into neuronal morphologies. This example shows how this is done.
+
+ The dendritic molecules diffuse along the dendrite in the region
+ specified by the *chemDistrib* keyword. In this case they are placed on
+ all apical and basal dendrites, but only at distances over 500 microns
+ from the soma. The spine head and PSD reaction systems are inserted only
+ into spines within this same *chemDistrib* zone. Diffusion coupling
+ between dendrite, and each spine head and PSD is also set up. It takes a
+ predefined chemical model file for Rdesigneur, which resides in the
+ ``./chem`` subdirectory. As in an earlier example, we turn off the
+ electrical calculations here as they are not needed. Here we plot out
+ the number of receptors on every single spine as a function of time.
+
+ (Documentation still to come here)
+
+.. Make a full multiscale model with complex spiny morphology and electrical and chemical signaling.
+.. ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+.. (Documentation still to come here)
diff --git a/docs/source/user/py/cookbook/multi_rdes.rst b/docs/source/user/py/cookbook/multi_rdes.rst
index dce3f358..baf5719b 100644
--- a/docs/source/user/py/cookbook/multi_rdes.rst
+++ b/docs/source/user/py/cookbook/multi_rdes.rst
@@ -1,6 +1,6 @@
-**********
-RDesigneur
-**********
+************************
+More Rdesigneur Examples
+************************
Building Chemical-Electrical Signalling Models
----------------------------------------------
diff --git a/docs/source/user/py/cookbook/multiscale.rst b/docs/source/user/py/cookbook/multiscale.rst
index a8a83165..d751d3e4 100644
--- a/docs/source/user/py/cookbook/multiscale.rst
+++ b/docs/source/user/py/cookbook/multiscale.rst
@@ -6,4 +6,8 @@ MultiScale Modeling
:maxdepth: 2
multi_sim_eg
+ Rd
multi_rdes
+
+
+
--
GitLab