diff --git a/docs/images/ex2.1_vclamp_a.png b/docs/images/ex2.1_vclamp_a.png
new file mode 100644
index 0000000000000000000000000000000000000000..bfb68c12d83563c3bd8b6f0ea281b713479fe94e
Binary files /dev/null and b/docs/images/ex2.1_vclamp_a.png differ
diff --git a/docs/images/ex2.1_vclamp_b.png b/docs/images/ex2.1_vclamp_b.png
new file mode 100644
index 0000000000000000000000000000000000000000..234dd59f6e3d74e8c35343af9e7c36df6a7ff421
Binary files /dev/null and b/docs/images/ex2.1_vclamp_b.png differ
diff --git a/docs/images/ex3.1_squid_vclamp.png b/docs/images/ex3.1_squid_vclamp.png
new file mode 100644
index 0000000000000000000000000000000000000000..0748b411e3b149792eac238e6860b0c13253f4c6
Binary files /dev/null and b/docs/images/ex3.1_squid_vclamp.png differ
diff --git a/docs/images/ex3.2_axon_propagating_AP.png b/docs/images/ex3.2_axon_propagating_AP.png
new file mode 100644
index 0000000000000000000000000000000000000000..649c380249e1bf7cc9586086ee0fe44d1ec251bc
Binary files /dev/null and b/docs/images/ex3.2_axon_propagating_AP.png differ
diff --git a/docs/images/ex3.3_AP_collision.png b/docs/images/ex3.3_AP_collision.png
new file mode 100644
index 0000000000000000000000000000000000000000..2c7b6ef81c20fbd756c51cc6b3de00a267c37a58
Binary files /dev/null and b/docs/images/ex3.3_AP_collision.png differ
diff --git a/docs/images/ex4.0_scaledSoma.png b/docs/images/ex4.0_scaledSoma.png
new file mode 100644
index 0000000000000000000000000000000000000000..1d5f861af2ce95d3195dc54564bdf4090366ee38
Binary files /dev/null and b/docs/images/ex4.0_scaledSoma.png differ
diff --git a/docs/images/ex4.1_ballAndStick.png b/docs/images/ex4.1_ballAndStick.png
new file mode 100644
index 0000000000000000000000000000000000000000..dcd4fb64bd6e8dac13322b1aa738b3792ef381ed
Binary files /dev/null and b/docs/images/ex4.1_ballAndStick.png differ
diff --git a/docs/images/ex4.2_sine_stim.png b/docs/images/ex4.2_sine_stim.png
new file mode 100644
index 0000000000000000000000000000000000000000..f592f8f3b70fb9fc1e211a5309624ace913c11d4
Binary files /dev/null and b/docs/images/ex4.2_sine_stim.png differ
diff --git a/docs/images/ex4.2_spiking.png b/docs/images/ex4.2_spiking.png
new file mode 100644
index 0000000000000000000000000000000000000000..33b3ab4aec5b97e1823c974250d33b0db3708060
Binary files /dev/null and b/docs/images/ex4.2_spiking.png differ
diff --git a/docs/images/ex5.0_random_syn_input.png b/docs/images/ex5.0_random_syn_input.png
new file mode 100644
index 0000000000000000000000000000000000000000..e7ef4c2683fd6764c60fbaae18a63be35e8b8520
Binary files /dev/null and b/docs/images/ex5.0_random_syn_input.png differ
diff --git a/docs/images/ex5.1_periodic_syn_input.png b/docs/images/ex5.1_periodic_syn_input.png
new file mode 100644
index 0000000000000000000000000000000000000000..e9503c3395615bf15e0306a2e0240f102c72c103
Binary files /dev/null and b/docs/images/ex5.1_periodic_syn_input.png differ
diff --git a/docs/images/ex7.0_spatial_chem_osc.png b/docs/images/ex7.0_spatial_chem_osc.png
new file mode 100644
index 0000000000000000000000000000000000000000..66927b396a5b500f602d81b86527e6cf4899b45b
Binary files /dev/null and b/docs/images/ex7.0_spatial_chem_osc.png differ
diff --git a/docs/images/ex7.1_diffusive_gradient.png b/docs/images/ex7.1_diffusive_gradient.png
new file mode 100644
index 0000000000000000000000000000000000000000..3ebb5c1dc63e78ca8017fa266fe5ca1516a6cc9c
Binary files /dev/null and b/docs/images/ex7.1_diffusive_gradient.png differ
diff --git a/docs/images/ex8.0_multiscale_Ca.png b/docs/images/ex8.0_multiscale_Ca.png
new file mode 100644
index 0000000000000000000000000000000000000000..eacdebe336fed7174a2ba1204e66a56cdd0e8a86
Binary files /dev/null and b/docs/images/ex8.0_multiscale_Ca.png differ
diff --git a/docs/images/ex8.0_multiscale_KA_conc.png b/docs/images/ex8.0_multiscale_KA_conc.png
new file mode 100644
index 0000000000000000000000000000000000000000..d396c781a9fe59be7c293c649d758d2a5cb2bbcc
Binary files /dev/null and b/docs/images/ex8.0_multiscale_KA_conc.png differ
diff --git a/docs/images/ex8.0_multiscale_cell_spiking.png b/docs/images/ex8.0_multiscale_cell_spiking.png
new file mode 100644
index 0000000000000000000000000000000000000000..27cb4cf92fea542c5b82cce696dd41c304407575
Binary files /dev/null and b/docs/images/ex8.0_multiscale_cell_spiking.png differ
diff --git a/docs/images/ex8.0_multiscale_currInj.png b/docs/images/ex8.0_multiscale_currInj.png
new file mode 100644
index 0000000000000000000000000000000000000000..6701d6425a2923a0d9dbb74e0aa7ce8a41e3e1f8
Binary files /dev/null and b/docs/images/ex8.0_multiscale_currInj.png differ
diff --git a/docs/images/ex8.2_Ca_dend.png b/docs/images/ex8.2_Ca_dend.png
new file mode 100644
index 0000000000000000000000000000000000000000..079abb90ccd16410a7fb66ec87e3f43e3710e18c
Binary files /dev/null and b/docs/images/ex8.2_Ca_dend.png differ
diff --git a/docs/images/ex8.2_Ca_spine.png b/docs/images/ex8.2_Ca_spine.png
new file mode 100644
index 0000000000000000000000000000000000000000..fca5539669053deb7a91af5c04ce826bc5a122eb
Binary files /dev/null and b/docs/images/ex8.2_Ca_spine.png differ
diff --git a/docs/images/ex8.2_Vm.png b/docs/images/ex8.2_Vm.png
new file mode 100644
index 0000000000000000000000000000000000000000..97abd1b9f74a0de685f8943ebd7a71666cd6e044
Binary files /dev/null and b/docs/images/ex8.2_Vm.png differ
diff --git a/docs/images/ex8.2_active_CaMKII.png b/docs/images/ex8.2_active_CaMKII.png
new file mode 100644
index 0000000000000000000000000000000000000000..4a22f8766fbf25e3ee9e19a1c259d6a4356966ab
Binary files /dev/null and b/docs/images/ex8.2_active_CaMKII.png differ
diff --git a/docs/images/ex8.3_CaMKII_spine.png b/docs/images/ex8.3_CaMKII_spine.png
new file mode 100644
index 0000000000000000000000000000000000000000..ba132db6374787848d89f720dca84d4678af98d5
Binary files /dev/null and b/docs/images/ex8.3_CaMKII_spine.png differ
diff --git a/docs/images/ex8.3_Vm.png b/docs/images/ex8.3_Vm.png
new file mode 100644
index 0000000000000000000000000000000000000000..cead72861c9f7a67c2e8aed44ae3bf8222f343c1
Binary files /dev/null and b/docs/images/ex8.3_Vm.png differ
diff --git a/docs/images/ex8.3_chan_p.png b/docs/images/ex8.3_chan_p.png
new file mode 100644
index 0000000000000000000000000000000000000000..092a508d683496207eed9e6c6c0368c425c5667a
Binary files /dev/null and b/docs/images/ex8.3_chan_p.png differ
diff --git a/docs/images/ex8.3_gluR.png b/docs/images/ex8.3_gluR.png
new file mode 100644
index 0000000000000000000000000000000000000000..f2a0da86489987a18799a80ff0ea45404f20cb1a
Binary files /dev/null and b/docs/images/ex8.3_gluR.png differ
diff --git a/docs/images/ex9.0_passive_cell_morpho.png b/docs/images/ex9.0_passive_cell_morpho.png
new file mode 100644
index 0000000000000000000000000000000000000000..0d00e176ffc1d9591b0ee3de560aaa982bb46735
Binary files /dev/null and b/docs/images/ex9.0_passive_cell_morpho.png differ
diff --git a/docs/source/conf.py b/docs/source/conf.py
index 4bd3770f99af1e814c6f891139a6ba8542dbb47a..43a10e7d8b472fa41b4fe9a57d1c353495244ef2 100644
--- a/docs/source/conf.py
+++ b/docs/source/conf.py
@@ -27,13 +27,13 @@ conf_dir_ = os.path.dirname( __file__ )
# If extensions (or modules to document with autodoc) are in another directory,
# add these directories to sys.path here. If the directory is relative to the
# documentation root, use os.path.abspath to make it absolute, like shown here.
-sys.path.insert(0, os.path.abspath('../python'))
-sys.path.append(os.path.abspath('../../moose-examples/snippets'))
-sys.path.append(os.path.abspath('../../moose-examples/tutorials/ChemicalOscillators'))
-sys.path.append(os.path.abspath('../../moose-examples/tutorials/ChemicalBistables'))
-sys.path.append(os.path.abspath('../../moose-examples/tutorials/ExcInhNet'))
-sys.path.append(os.path.abspath('../../moose-examples/neuroml/lobster_pyloric'))
-sys.path.append(os.path.abspath('../../moose-examples/tutorials/ExcInhNetCaPlasticity'))
+sys.path.insert(0, os.path.abspath('/home/harsha/MOOSE/FullMoose/moose-core/python'))
+sys.path.append(os.path.abspath('/home/harsha/MOOSE/moose-examples-july/snippets'))
+sys.path.append(os.path.abspath('/home/harsha/MOOSE/moose-examples-july/tutorials/ChemicalOscillators'))
+sys.path.append(os.path.abspath('/home/harsha/MOOSE/moose-examples-july/tutorials/ChemicalBistables'))
+sys.path.append(os.path.abspath('/home/harsha/MOOSE/moose-examples-july/tutorials/ExcInhNet'))
+sys.path.append(os.path.abspath('/home/harsha/MOOSE/moose-examples-july/neuroml/lobster_pyloric'))
+sys.path.append(os.path.abspath('/home/harsha/MOOSE/moose-examples-july/tutorials/ExcInhNetCaPlasticity'))
sys.path.append(os.path.join(conf_dir_, 'Extensions') )
# -- General configuration -----------------------------------------------------
@@ -301,6 +301,6 @@ exclude_patterns = ['/docs/source/user/py/references/*.rst']
import subprocess, os
read_the_docs_build = os.environ.get('READTHEDOCS', None) == 'True'
-if not read_the_docs_build:
- subprocess.call('cd doxygen; echo HELLO......................; doxygen Doxyfile', shell=True)
+#if not read_the_docs_build:
+# subprocess.call('cd doxygen; echo HELLO......................; doxygen Doxyfile', shell=True)
diff --git a/docs/source/user/py/cookbook/Rd.rst b/docs/source/user/py/cookbook/Rd.rst
deleted file mode 100644
index 625b8bc819cf5b695210c854050e7630d60c6ff0..0000000000000000000000000000000000000000
--- a/docs/source/user/py/cookbook/Rd.rst
+++ /dev/null
@@ -1,1186 +0,0 @@
-**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/multiscale.rst b/docs/source/user/py/cookbook/multiscale.rst
index d751d3e4bb807f031020ba03fd4d4dd99854e20f..799ee8045213ca6c3f1d28300b74bb4ec888fd46 100644
--- a/docs/source/user/py/cookbook/multiscale.rst
+++ b/docs/source/user/py/cookbook/multiscale.rst
@@ -6,8 +6,6 @@ MultiScale Modeling
:maxdepth: 2
multi_sim_eg
- Rd
- multi_rdes
diff --git a/docs/source/user/py/rdesigneur/index_rd.rst b/docs/source/user/py/rdesigneur/index_rd.rst
new file mode 100644
index 0000000000000000000000000000000000000000..a729cb8b0a2f2a9d357f566bd8821085078f7571
--- /dev/null
+++ b/docs/source/user/py/rdesigneur/index_rd.rst
@@ -0,0 +1,13 @@
+.. MOOSE documentation master file, created by
+ sphinx-quickstart on Tue Jul 31 19:05:47 2018.
+ You can adapt this file completely to your liking, but it should at least
+ contain the root `toctree` directive.
+
+Rdesignuer
+===========
+
+.. toctree::
+ :maxdepth: 2
+
+ rdes
+ multi_rdes
diff --git a/docs/source/user/py/cookbook/multi_rdes.rst b/docs/source/user/py/rdesigneur/multi_rdes.rst
similarity index 98%
rename from docs/source/user/py/cookbook/multi_rdes.rst
rename to docs/source/user/py/rdesigneur/multi_rdes.rst
index 2015743a75e5e827cd0b40b1a46f40c9a4024f9b..d240f1b1cbb2412ab0d9331eff0ccb4f25d136c7 100644
--- a/docs/source/user/py/cookbook/multi_rdes.rst
+++ b/docs/source/user/py/rdesigneur/multi_rdes.rst
@@ -1,5 +1,5 @@
************************
-More Rdesigneur Examples
+Rdesigneur Examples
************************
.. hidden-code-block:: reStructuredText
diff --git a/docs/source/user/py/rdesigneur/rdes.rst b/docs/source/user/py/rdesigneur/rdes.rst
new file mode 100644
index 0000000000000000000000000000000000000000..7d77f8fb9f06e9b28c1119b58bd155319e54f9b8
--- /dev/null
+++ b/docs/source/user/py/rdesigneur/rdes.rst
@@ -0,0 +1,1962 @@
+**Rdesigneur: Building multiscale models**
+==========================================
+
+.. Upi Bhalla
+
+.. Aug 26 2016. Updated July 2018
+
+.. --------------
+
+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
+- Adaptors that couple between these for multiscale models
+
+It also folds in simulation input and output
+
+- Time-series stimuli for molecular concentration change and reaction rates
+- Current and voltage clamp
+- Synaptic input.
+- Time-series plots
+- File dumps
+- 3-D neuronal graphics
+
+Rdesigneur's main role is to specify how these are put together,
+including assigning parameters for the model. Using Rdesigneur one can compactly
+and quickly put together quite complex multiscale models.
+
+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. The files for these examples are also available in
+``moose-examples/tutorials/Rdesigneur``, and the file names are mentioned
+as we go along.
+
+Bare Rdesigneur: single passive compartment
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+*ex1_minimalModel.py*
+
+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
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+*ex2.0_currentPulse.py*
+
+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.
+
+Simulate and display voltage clamp stimulus to soma
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+*ex2.1_vclamp.py*
+
+This model introduces the voltage clamp stimulus on a passive compartment.
+As before, we add a few lines to define the stimulus and plot.
+This script displays both the membrane potential, and the holding current
+of the voltage clamp circuit as
+it charges and discharges the passive compartment model.
+
+::
+
+ import moose
+ import rdesigneur as rd
+ rdes = rd.rdesigneur(
+ stimList = [['soma', '1', '.', 'vclamp', '-0.065 + (t>0.1 && t<0.2) * 0.02' ]],
+ plotList = [
+ ['soma', '1', '.', 'Vm', 'Soma membrane potential'],
+ ['soma', '1', 'vclamp', 'current', 'Soma holding current'],
+ ]
+ )
+ rdes.buildModel()
+ moose.reinit()
+ moose.start( 0.3 )
+ rdes.display()
+
+Here the *stimList* line tells the system to deliver a voltage clamp (vclamp)
+on the soma, starting at -65 mV and jumping up by 20 mV between 0.1 and 0.2
+seconds. The *plotList* now includes two entries, and will generate two plots.
+The first is for plotting the soma membrane potential, just to be sure that
+the voltage clamp is doing its job.
+
+.. figure:: ../../../../images/ex2.1_vclamp_a.png
+ :alt: Plot for membrane potential in voltage clamp
+
+ Plot for membrane potential in voltage clamp
+
+The second graph plots the holding current. Note the capacitive transients.
+
+.. figure:: ../../../../images/ex2.1_vclamp_b.png
+ :alt: Plot for holding current for voltage clamp
+
+ Plot for holding current for voltage clamp
+
+HH Squid model in a single compartment
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+*ex3.0_squid_currentPulse.py*
+
+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
+
+There are several interesting things to do with the model by varying stimulus
+parameters:
+
+ - Change injection current.
+ - Put in a protocol to get rebound action potential.
+ - Put in a current ramp, and run it for a different duration
+ - Put in a frequency chirp, and see how the squid model is tuned
+ to a certain frequency range.
+ - Modify channel or passive parameters. See if it still fires.
+ - Try the frequency chirp on the cell with parameters changed. Does
+ the tuning change?
+
+
+HH Squid model with voltage clamp
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+*ex3.1_squid_vclamp.py*
+
+This is the same squid model, but now we add a voltage clamp to the squid
+and monitor the holding current. This stimulus line is identical to ex2.1.
+
+::
+
+ 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', '.', 'vclamp', '-0.065 + (t>0.1 && t<0.2) * 0.02' ]],
+ plotList = [
+ ['soma', '1', '.', 'Vm', 'Membrane potential'],
+ ['soma', '1', 'vclamp', 'current', 'Soma holding current']
+ ]
+ )
+ rdes.buildModel()
+ moose.reinit()
+ moose.start( 0.3 )
+ rdes.display()
+
+Here we see the classic HH current response, a downward brief deflection due to
+the Na channel, and a slower upward sustained current due to the K delayed
+rectifier.
+
+.. figure:: ../../../../images/ex3.1_squid_vclamp.png
+ :alt: Plot for HH squid voltage clamp pulse.
+
+ Plot for HH squid voltage clamp pulse.
+
+Here are some suggestions for further exploration:
+
+ - Monitor individual channel currents through additional plots.
+ - Convert this into a voltage clamp series. Easiest way to do this is
+ to complete the rdes.BuildModel, then delete the Function object
+ on the */model/elec/soma/vclamp*. Now you can simply set the 'command'
+ field of the vclamp in a for loop, going from -ve to +ve voltages.
+ Remember, SI units. You may wish to capture the plot vectors each
+ cycle. The plot vectors are accessed by something like
+
+ ``moose.element( '/model/graphs/plot1' ).vector``
+
+
+HH Squid model in an axon
+~~~~~~~~~~~~~~~~~~~~~~~~~
+
+*ex3.2_squid_axon_propgn.py*
+
+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/ex3.2_axon_propagating_AP.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.
+
+Action potential collision in HH Squid axon model
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+*ex3.3_AP_collision.py*
+
+This is identical to the previous example, except that now we deliver current
+injection at at two points, the soma and a point along the axon. The modified
+stimulus line is:
+
+::
+
+ ...
+ stimList = [['soma', '1', '.', 'inject', '(t>0.01 && t<0.2) * 2e-11' ],
+ ['axon100', '1', '.', 'inject', '(t>0.01 && t<0.2) * 3e-11' ]],
+ ...
+
+Watch how the AP is triggered bidirectionally from the stimulus point on the
+100th segment of the axon, and observe what happens when two action potentials
+bump into each other.
+
+.. figure:: ../../../../images/ex3.3_AP_collision.png
+ :alt: Colliding action potentials
+
+ Colliding action potentials
+
+
+
+HH Squid model in a myelinated axon
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+*ex3.4_myelinated_axon.py*
+
+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
+
+Alternate (non-squid) way to define soma
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+*ex4.0_scaledSoma.py*
+
+The default HH-squid axon is not a very convincing soma. Rdesigneur offers a
+somewhat more general way to define the soma in the cell prototype line.
+
+::
+
+ import moose
+ import pylab
+ import rdesigneur as rd
+ rdes = rd.rdesigneur(
+ # cellProto syntax: ['somaProto', 'name', dia, length]
+ cellProto = [['somaProto', 'soma', 20e-6, 200e-6]],
+ chanProto = [['make_HH_Na()', 'Na'], ['make_HH_K()', 'K']],
+ chanDistrib = [
+ ['Na', 'soma', 'Gbar', '1200' ],
+ ['K', 'soma', 'Gbar', '360' ]],
+ stimList = [['soma', '1', '.', 'inject', '(t>0.01 && t<0.05) * 1e-9' ]],
+ plotList = [['soma', '1', '.', 'Vm', 'Membrane potential']],
+ moogList = [['#', '1', '.', 'Vm', 'Vm (mV)']]
+ )
+
+ rdes.buildModel()
+ soma = moose.element( '/model/elec/soma' )
+ print( 'Soma dia = {}, length = {}'.format( soma.diameter, soma.length ) )
+ moose.reinit()
+
+ rdes.displayMoogli( 0.0005, 0.06, 0.0 )
+
+Here the crucial line is the *cellProto* line. There are four arguments here:
+
+ ``['somaProto', 'name', dia, length]``
+
+ - The first argument tells the system to use a prototype soma, that is
+ a single cylindrical compartment.
+ - The second argument is the name to give the cell.
+ - The third argument is the diameter. Note that this is a double,
+ not a string.
+ - The fourth argument is the length of the cylinder that makes up the
+ soma. This too is a double, not a string.
+ The cylinder is oriented along the x axis, with one end at (0,0,0)
+ and the other end at (length, 0, 0).
+
+This is what the soma looks like:
+
+.. figure:: ../../../../images/ex4.0_scaledSoma.png
+ :alt: Image of soma.
+
+ Image of soma.
+
+It a somewhat elongated soma, being a cylinder 10 times as long as it is wide.
+
+Ball-and-stick model of a neuron
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+*ex4.1_ballAndStick.py*
+
+A somewhat more electrically reasonable model of a neuron has a soma and a
+single dendrite, which can itself be subdivided into segments so that it
+can exhibit voltage gradients, have channel and receptor distributions,
+and so on. This is accomplished in *rdesigneur* using a variant of the
+cellProto syntax.
+
+::
+
+ import moose
+ import pylab
+ import rdesigneur as rd
+ rdes = rd.rdesigneur(
+ # cellProto syntax: ['ballAndStick', 'name', somaDia, somaLength, dendDia, dendLength, numDendSegments ]
+ # The numerical arguments are all optional
+ cellProto = [['ballAndStick', 'soma', 20e-6, 20e-6, 4e-6, 500e-6, 10]],
+ chanProto = [['make_HH_Na()', 'Na'], ['make_HH_K()', 'K']],
+ chanDistrib = [
+ ['Na', 'soma', 'Gbar', '1200' ],
+ ['K', 'soma', 'Gbar', '360' ],
+ ['Na', 'dend#', 'Gbar', '400' ],
+ ['K', 'dend#', 'Gbar', '120' ]
+ ],
+ stimList = [['soma', '1', '.', 'inject', '(t>0.01 && t<0.05) * 1e-9' ]],
+ plotList = [['soma', '1', '.', 'Vm', 'Membrane potential']],
+ moogList = [['#', '1', '.', 'Vm', 'Vm (mV)']]
+ )
+ rdes.buildModel()
+ soma = moose.element( '/model/elec/soma' )
+ moose.reinit()
+ rdes.displayMoogli( 0.0005, 0.06, 0.0 )
+
+As before, the *cellProto* line plays a key role. Here, because we have a long
+dendrite, we have a few more numerical arguments. All of the numerical
+arguments are optional.
+
+ ``['ballAndStick', 'name', somaDia, somaLength, dendDia, dendLength, numDendSegments ]``
+
+ - The first argument specifies a ballAndStick model: soma + dendrite.
+ The length of the dendrite is along the x axis. The soma is a single
+ segment, the dendrite can be more than one.
+ - The second argument is the name to give the cell.
+ - Arg 3 is the soma diameter, as a double.
+ - Arg 4 is the length of the soma, as a double.
+ - Arg 5 is the diameter of the dendrite, as a double.
+ - Arg 6 is the length of the dendrite, as a double.
+ - Arg 7 is the number of segments into which the dendrite should be
+ divided. This is a positive integer greater than 0.
+
+This is what the ball-and-stick cell looks like:
+
+.. figure:: ../../../../images/ex4.1_ballAndStick.png
+ :alt: Image of ball and stick cell.
+
+ Image of ball and stick cell.
+
+In this version of the 3-D display, the soma is displayed as a bit blocky
+rather than round.
+Note that we have populated the dendrite with Na and K channels and it has
+10 segments, so it supports action potential propagation. The snapshot
+illustrates this.
+
+Here are some things to try:
+
+ - Change the length of the dendrite
+ - Change the number of segments. Explore what it does to accuracy. How
+ will you know that you have an accurate model?
+
+Benchmarking simulation speed
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+*ex4.2_ballAndStickSpeed.py*
+
+The ball-and-stick model gives us an opportunity to check out your system
+and how computation scales with model size. While we're at it we'll deliver
+a sine-wave stimulus just to see how it can be done. The test model is
+very similar to the previous one, ex4.1:
+
+::
+
+ import moose
+ import pylab
+ import rdesigneur as rd
+ import time
+ rdes = rd.rdesigneur(
+ cellProto = [['ballAndStick', 'soma', 20e-6, 20e-6, 4e-6, 500e-6, 10]],
+ chanProto = [['make_HH_Na()', 'Na'], ['make_HH_K()', 'K']],
+ chanDistrib = [
+ ['Na', 'soma', 'Gbar', '1200' ],
+ ['K', 'soma', 'Gbar', '360' ],
+ ['Na', 'dend#', 'Gbar', '400' ],
+ ['K', 'dend#', 'Gbar', '120' ]
+ ],
+ stimList = [['soma', '1', '.', 'inject', '(1+cos(t/10))*(t>31.4 && t<94) * 0
+ .2e-9' ]],
+ plotList = [
+ ['soma', '1', '.', 'Vm', 'Membrane potential'],
+ ['soma', '1', '.', 'inject', 'Stimulus current']
+ ],
+ )
+ rdes.buildModel()
+ runtime = 100
+ moose.reinit()
+ t0= time.time()
+ moose.start( runtime )
+ print "Real time to run {} simulated seconds = {} seconds".format( runtime, time
+ .time() - t0 )
+
+ rdes.display()
+
+While the real point of this simulation is to check speed, it does illustrate
+how to deliver a stimulus shaped like a sine wave:
+
+.. figure:: ../../../../images/ex4.2_sine_stim.png
+ :alt: Sine-wave shaped stimulus.
+
+ Sine-wave shaped stimulus.
+
+We can see that the cell has a peculiar response to this. Not surprising, as
+the cell uses HH channels which are not good at rate coding.
+
+.. figure:: ../../../../images/ex4.2_spiking.png
+ :alt: Spiking response to sine-wave shaped stimulus.
+
+ Spiking response to sine-wave shaped stimulus.
+
+As a reference point, on a fast 2018 laptop this benchmark runs in 5.4 seconds.
+Some more things to try for benchmarking:
+
+ - How slow does it get if you turn on the 3-D moogli display?
+ - Is it costlier to run 2 compartments for 1000 seconds, or
+ 200 compartments for 10 seconds?
+
+Synaptic stimulus: random (Possion)
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+*ex5.0_random_syn_input.py*
+
+In this example we introduce synaptic inputs: both the receptor channels
+and a means for stimulating the channels. We do this in a passive model.
+
+::
+
+ import moose
+ import rdesigneur as rd
+ rdes = rd.rdesigneur(
+ cellProto = [['somaProto', 'soma', 20e-6, 200e-6]],
+ chanProto = [['make_glu()', 'glu']],
+ chanDistrib = [['glu', 'soma', 'Gbar', '1' ]],
+ stimList = [['soma', '0.5', 'glu', 'randsyn', '50' ]],
+ # Deliver stimulus to glu synapse on soma, at mean 50 Hz Poisson.
+ plotList = [['soma', '1', '.', 'Vm', 'Soma membrane potential']]
+ )
+ rdes.buildModel()
+ moose.reinit()
+ moose.start( 0.3 )
+ rdes.display()
+
+Most of the rdesigneur setup uses familiar syntax.
+
+Novelty 1: we use the default built-in glutamate receptor model, in chanProto.
+We just put it in the soma at a max conductance of 1 Siemen/sq metre.
+
+Novelty 2: We specify a new kind of stimulus in the stimList:
+
+ ``['soma', '0.5', 'glu', 'randsyn', '50' ]``
+
+Most of this is similar to previous stimLists.
+
+ - arg0: 'soma': the named compartments in the cell to populate with
+ the *glu* receptor
+ - arg1: '0.5': Tell the system to use a uniform synaptic weight of 0.5.
+ This argument could be a more complicated expression incorporating
+ spatial arguments. Here it is just uniform.
+ - arg2: 'glu': Which receptor to stimulate
+ - arg3: 'randsyn': Apply random (Poisson) synaptic input.
+ - arg4: '50': Mean firing rate of the Poisson input. Note that this last
+ argument could be a function of time and hence is quite versatile.
+
+As the model has no voltage-gated channels, we do not see spiking.
+
+.. figure:: ../../../../images/ex5.0_random_syn_input.png
+ :alt: Random synaptic input with a Poisson distribution.
+
+ Random synaptic input with a Poisson distribution.
+
+Things to try: Vary the rate and the weight of the synaptic input.
+
+Synaptic stimulus: periodic
+~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+*ex5.1_periodic_syn_input.py*
+
+This is almost identical to 5.0, except that the input is now perfectly
+periodic. The one change is of an argument in the stimList to say
+``periodicsyn`` rather than ``randsyn``.
+
+::
+
+ import moose
+ import rdesigneur as rd
+ rdes = rd.rdesigneur(
+ cellProto = [['somaProto', 'soma', 20e-6, 200e-6]],
+ chanProto = [['make_glu()', 'glu']],
+ chanDistrib = [['glu', 'soma', 'Gbar', '1' ]],
+
+ # Deliver stimulus to glu synapse on soma, periodically at 50 Hz.
+ stimList = [['soma', '0.5', 'glu', 'periodicsyn', '50' ]],
+ plotList = [['soma', '1', '.', 'Vm', 'Soma membrane potential']]
+ )
+ rdes.buildModel()
+ moose.reinit()
+ moose.start( 0.3 )
+ rdes.display()
+
+As designed, we get periodically firing synaptic input.
+
+.. figure:: ../../../../images/ex5.1_periodic_syn_input.png
+ :alt: Periodic synaptic input
+
+ Periodic synaptic input
+
+
+Reaction system in a single compartment
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+*ex6_chem_osc.py*
+
+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. Its general form is:
+
+::
+
+ s ---a---> a // s goes to a, catalyzed by a.
+ s ---a---> b // s goes to b, catalyzed by a.
+ a ---b---> s // a goes to s, catalyzed by b.
+ b -------> s // b is degraded irreversibly to s
+
+Here is the script:
+
+::
+
+ 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
+~~~~~~~~~~~~~~~~~~~~~~~~~
+
+*ex7.0_spatial_chem_osc.py*
+
+In order to see what a reaction-diffusion system looks like, we assign the
+``diffusionLength`` expression in the previous example to a much shorter
+length, 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,
+ #This subdivides the length of the soma into 2 micron voxels
+ diffusionLength = 2e-6,
+ 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, rotation = 0, azim = np.pi/2, elev = 0.0 )
+
+This is the new value for diffusion length.
+
+::
+
+ diffusionLength = 2e-3,
+
+With this change we tell *rdesigneur* to use the diffusion length of 2 microns.
+This happens to be the default too. 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, rotation = 0, azim = np.pi/2, elev = 0.0 )
+
+The arguments mean: *displayMoogli( frametime, runtime, rotation, azimuth, elevation )*
+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.
+ azimuth = Azimuth angle of view point, in radians
+ elevation = elevation angle of view point, 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
+
+For those who would rather use the much simpler matplotlib 3-D display option,
+this is what the same simulation looks like:
+
+.. figure:: ../../../../images/ex7.0_spatial_chem_osc.png
+ :alt: Display for oscillatory reac-diff simulation using matplotlib
+
+ Display for oscillatory reac-diff simulation using matplotlib
+
+Primer on using the 3-D MOOGLI display
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+There are two variants of the MOOGLI display. The first, named Moogli,
+uses OpenGL and OpenSceneGraph. It is fast to display, slow to load, and
+difficult to compile. It produces much better looking 3-D graphics.
+The second is a fallback interface using mplot3d, which is a library of
+Matplotlib and so should be generally available. It is slower to display,
+faster to load, but needs no special compilation. It uses stick graphics
+and though it conveys much the same information, isn't as nice to look at
+as the original Moogli. Its controls are more or less the same but less
+smooth than the original Moogli.
+
+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.
+
+Diffusion of a single molecule
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+*ex7.1_diffusive_gradient.py*
+
+This is simply a test model to confirm that simple diffusion happens as
+expected. While the model is just that of a single pool, we spend a few lines
+taking snapshots of the spatial profile of this pool.
+
+::
+
+ import moose
+ import pylab
+ import re
+ import rdesigneur as rd
+ import matplotlib.pyplot as plt
+ import numpy as np
+
+ moose.Neutral( '/library' )
+ moose.Neutral( '/library/diffn' )
+ moose.CubeMesh( '/library/diffn/dend' )
+ A = moose.Pool( '/library/diffn/dend/A' )
+ A.diffConst = 1e-10
+
+ rdes = rd.rdesigneur(
+ turnOffElec = True,
+ diffusionLength = 1e-6,
+ chemProto = [['diffn', 'diffn']],
+ chemDistrib = [['diffn', 'soma', 'install', '1' ]],
+ moogList = [
+ ['soma', '1', 'dend/A', 'conc', 'A Conc', 0, 360 ]
+ ]
+ )
+ rdes.buildModel()
+
+ rdes.displayMoogli( 1, 2, rotation = 0, azim = -np.pi/2, elev = 0.0, block = False )
+ av = moose.vec( '/model/chem/dend/A' )
+ for i in range(10):
+ av[i].concInit = 1
+ moose.reinit()
+ plist = []
+ for i in range( 20 ):
+ plist.append( av.conc[:200] )
+ moose.start( 2 )
+ fig = plt.figure( figsize = ( 10, 12 ) )
+ plist = np.array( plist ).T
+ plt.plot( range( 0, 200 ), plist )
+ plt.xlabel( "position ( microns )" )
+ plt.ylabel( "concentration ( mM )" )
+ plt.show( block = True )
+
+
+Here are the snapshots, overlaid in a single plot:
+
+.. figure:: ../../../../images/ex7.1_diffusive_gradient.png
+ :alt: Display of how a molecule A spreads through the inter
+
+ Display for simple time-series of spread of a diffusing molecule
+ using matplotlib
+
+
+Multiscale models: single compartment
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+*ex8.0_multiscale_KA_phosph.py*
+
+The next few examples are for the multiscale modeling that is the main purpose
+of rdesigneur and MOOSE as a whole. These are 'toy' examples in that the
+chemical and electrical signaling is simplified, but they exhibit dynamics
+that are of real interest.
+
+The first example is of a bistable system where the feedback loop comprises of
+
+`calcium influx -> chemical activity -> channel modulation -> electrical activity -> calcium influx.`
+
+Calcium enters through voltage gated calcium channels, leads to enzyme
+activation and phosphorylation of a KA channel, which depolarizes the cell,
+so it spikes more, so more calcium enters.
+
+::
+
+ import moose
+ import pylab
+ import rdesigneur as rd
+ rdes = rd.rdesigneur(
+ elecDt = 50e-6,
+ chemDt = 0.002,
+ chemPlotDt = 0.002,
+ # cellProto syntax: ['somaProto', 'name', dia, length]
+ cellProto = [['somaProto', 'soma', 12e-6, 12e-6]],
+ chemProto = [['./chem/chanPhosphByCaMKII.g', 'chem']],
+ chanProto = [
+ ['make_Na()', 'Na'],
+ ['make_K_DR()', 'K_DR'],
+ ['make_K_A()', 'K_A' ],
+ ['make_Ca()', 'Ca' ],
+ ['make_Ca_conc()', 'Ca_conc' ]
+ ],
+ # Some changes to the default passive properties of the cell.
+ passiveDistrib = [['.', 'soma', 'CM', '0.03', 'Em', '-0.06']],
+ chemDistrib = [['chem', 'soma', 'install', '1' ]],
+ chanDistrib = [
+ ['Na', 'soma', 'Gbar', '300' ],
+ ['K_DR', 'soma', 'Gbar', '250' ],
+ ['K_A', 'soma', 'Gbar', '200' ],
+ ['Ca_conc', 'soma', 'tau', '0.0333' ],
+ ['Ca', 'soma', 'Gbar', '40' ]
+ ],
+ adaptorList = [
+ [ 'dend/chan', 'conc', 'K_A', 'modulation', 0.0, 70 ],
+ [ 'Ca_conc', 'Ca', 'dend/Ca', 'conc', 0.00008, 2 ]
+ ],
+ # Give a + pulse from 5 to 7s, and a - pulse from 20 to 21.
+ stimList = [['soma', '1', '.', 'inject', '((t>5 && t<7) - (t>20 && t<21)) * 1.0e-12' ]],
+ plotList = [
+ ['soma', '1', '.', 'Vm', 'Membrane potential'],
+ ['soma', '1', '.', 'inject', 'current inj'],
+ ['soma', '1', 'K_A', 'Ik', 'K_A current'],
+ ['soma', '1', 'dend/chan', 'conc', 'Unphosph K_A conc'],
+ ['soma', '1', 'dend/Ca', 'conc', 'Chem Ca'],
+ ],
+ )
+
+ rdes.buildModel()
+ moose.reinit()
+ moose.start( 30 )
+
+ rdes.display()
+
+There is only one fundamentally new element in this script:
+
+**adaptor List:** `[source, sourceField, dest, destField, offset, scale]`
+The adaptor list maps between molecular, electrical or even structural
+quantities in the simulation. At present it is linear mapping, in due course
+it may evolve to an arbitrary function.
+
+The two adaptorLists in the above script do the following:
+
+ ``[ 'dend/chan', 'conc', 'K_A', 'modulation', 0.0, 70 ]``:
+
+Use the concentration of the 'chan' molecule in the 'dend' compartment,
+to modulate the conductance of the 'K_A' channel such that the basal
+conductance is zero and 1 millimolar of 'chan' results in a conductance that is
+70 times greater than the baseline conductance of the channel, *Gbar*.
+
+It is advisable to use the field *'modulation'* on channels undergoing scaling,
+rather than to directly assign the conductance *'Gbar'*. This is because
+*Gbar* is an absolute conductance, and therefore it is scaled to the area of
+the electrical segment. This makes it difficult to keep track of. *Modulation*
+is a simple multiplier term onto *Gbar*, and is therefore easier to work with.
+
+ ``[ 'Ca_conc', 'Ca', 'dend/Ca', 'conc', 0.00008, 2 ]``:
+
+Use the concentration of *Ca* as computed in the electrical model, to assign
+the concentration of molecule *Ca* on the dendrite compartment. There is a
+basal level of 80 nanomolar, and every unit of electrical *Ca* maps to 2
+millimolar of chemical *Ca*.
+
+The arguments in the adaptorList are:
+
+ * **Source and Dest**: Strings. These can be either a molecular or an
+ electrical object. To identify a molecular object, it should be
+ prefixed with the name of the chemical compartment, which is one
+ of *dend, spine, psd*. Thus *dend/chan* specifies a molecule
+ named *'chan'* sitting in the *'dend'* compartment.
+
+ To identify an electrical object, just pass in its path,
+ such as '.' or *'Ca_conc'*.
+
+ Note that the adaptors do **not** need to know anything about the
+ location. It is assumed that the adaptors do their job wherever
+ the specified source and dest coexist. There is a subtlety here
+ due to the different length and time scales. The rule of thumb
+ is that the adaptor averages whichever one is subdivided more finely.
+
+ - Example 1: Molecules are typically spatially partitioned into
+ short voxels (micron-scale) compared to typical 100-micron
+ electrical
+ segments. So an adaptor going from molecules to, say, channel
+ conductance, would average all the molecular voxels that fit
+ in the electrical segment.
+ - Example 2: Electrical activity is typically much faster than
+ chemical.
+ So an adaptor going from an electrical entity (Ca computed from
+ channel opening) to molecules (Chemical Ca concentration) would
+ average all the time-steps between updates to the molecule.
+
+ * **Fields**: Strings. These are simply the field names on the
+ objects coupled by the adaptors.
+
+ * **offset and scale**: Doubles. At present the adaptor is just a
+ straight-line conversion, obeying ``y = mx + c``. The computed
+ output is *y*, averaged input is *x*, offset is *c* and scale is *m*.
+
+There is a handy new line to specify cellular passive properties:
+
+**passiveDistrib:** `['.', path, field, value, field, value, ... ]`,
+
+ * '.': This is just a placeholder.
+ * path: String. Specifies the object whose parameters are to be changed.
+ * field: String. Name of the field on the object.
+ * value: String, that is the value has to be enclosed in quotes. The
+ value to be assigned to the object.
+
+With these in place, the model behavior is rather neat. It starts out silent,
+then we apply 2 seconds of +ve current injection.
+
+.. figure:: ../../../../images/ex8.0_multiscale_currInj.png
+ :alt: Current injection stimuli for multiscale model.
+
+ Current injection stimuli for multiscale model.
+
+The cell fires briskly, and keeps firing even when the current injection
+drops to zero.
+
+.. figure:: ../../../../images/ex8.0_multiscale_cell_spiking.png
+ :alt: Firing responses of cell with multiscale signaling.
+
+ Firing responses of cell with multiscale signaling.
+
+The firing of the neuron leads to Ca influx.
+
+.. figure:: ../../../../images/ex8.0_multiscale_Ca.png
+ :alt: Calcium buildup in cell due to firing.
+
+ Calcium buildup in cell due to firing.
+
+The chemical reactions downstream of Ca lead to phosphorylation of the K_A
+channel. Only the unphosphorylated K_A channel is active, so the net effect
+is to reduce K_A conductance while the Ca influx persists.
+
+.. figure:: ../../../../images/ex8.0_multiscale_KA_conc.png
+ :alt: Removal of KA channel due to phosphorylation.
+
+ Removal of KA channel due to phosphorylation.
+
+
+Since the phosphorylated form has low conductance, the cell becomes more
+excitable and keeps firing even when the current injection is stopped. It takes
+a later, -ve current injection to turn the firing off again.
+
+Suggestions for things to do with the model:
+
+ - Vary the adaptor settings, which couple electrical to chemical
+ signaling and vice versa.
+ - Play with the channel densities
+ - Open the chem model in moosegui and vary its parameters too.
+
+Multiscale model spanning PSD, spine head and dendrite
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+*ex8.2_multiscale_glurR_phosph_3compt.py*
+
+This is another multiscale model on similar lines to 8.0. It is structurally
+and computationally more complicated, because the action is distributed between
+spines and dendrites, but formally it does the same thing: it turns on and
+stays on after a strong stimulus, due to phosphorylation of a (receptor)
+channel leading to greater excitability.
+
+`calcium influx -> chemical activity -> channel modulation -> electrical activity -> calcium influx.`
+
+The model is bistable as long as synaptic input keeps coming along at a basal
+rate, in this case 1 Hz.
+
+Here we have two new lines, to do with addition of spines. These are discussed
+in detail in a later example. For now it is enough to know that the
+**spineProto** line defines one of the prototype spines to be used to put into
+the model, and the **spineDistrib** line tells the system where to put them,
+and how widely to space them.
+
+::
+
+ import moose
+ import rdesigneur as rd
+ rdes = rd.rdesigneur(
+ elecDt = 50e-6,
+ chemDt = 0.002,
+ diffDt = 0.002,
+ chemPlotDt = 0.02,
+ useGssa = False,
+ # cellProto syntax: ['ballAndStick', 'name', somaDia, somaLength, dendDia, d
+ endLength, numDendSegments ]
+ cellProto = [['ballAndStick', 'soma', 12e-6, 12e-6, 4e-6, 100e-6, 2 ]],
+ chemProto = [['./chem/chanPhosph3compt.g', 'chem']],
+ spineProto = [['makeActiveSpine()', 'spine']],
+ chanProto = [
+ ['make_Na()', 'Na'],
+ ['make_K_DR()', 'K_DR'],
+ ['make_K_A()', 'K_A' ],
+ ['make_Ca()', 'Ca' ],
+ ['make_Ca_conc()', 'Ca_conc' ]
+ ],
+ passiveDistrib = [['.', 'soma', 'CM', '0.01', 'Em', '-0.06']],
+ spineDistrib = [['spine', '#dend#', '50e-6', '1e-6']],
+ chemDistrib = [['chem', '#', 'install', '1' ]],
+ chanDistrib = [
+ ['Na', 'soma', 'Gbar', '300' ],
+ ['K_DR', 'soma', 'Gbar', '250' ],
+ ['K_A', 'soma', 'Gbar', '200' ],
+ ['Ca_conc', 'soma', 'tau', '0.0333' ],
+ ['Ca', 'soma', 'Gbar', '40' ]
+ ],
+ adaptorList = [
+ [ 'psd/chan_p', 'n', 'glu', 'modulation', 0.1, 1.0 ],
+ [ 'Ca_conc', 'Ca', 'spine/Ca', 'conc', 0.00008, 8 ]
+ ],
+ # Syn input basline 1 Hz, and 40Hz burst for 1 sec at t=20. Syn weight
+ # is 0.5, specified in 2nd argument as a special case stimLists.
+ stimList = [['head#', '0.5','glu', 'periodicsyn', '1 + 40*(t>10 && t<11)']],
+ plotList = [
+ ['soma', '1', '.', 'Vm', 'Membrane potential'],
+ ['#', '1', 'spine/Ca', 'conc', 'Ca in Spine'],
+ ['#', '1', 'dend/DEND/Ca', 'conc', 'Ca in Dend'],
+ ['#', '1', 'spine/Ca_CaM', 'conc', 'Ca_CaM'],
+ ['head#', '1', 'psd/chan_p', 'conc', 'Phosph gluR'],
+ ['head#', '1', 'psd/Ca_CaM_CaMKII', 'conc', 'Active CaMKII'],
+ ]
+ )
+ moose.seed(123)
+ rdes.buildModel()
+ moose.reinit()
+ moose.start( 25 )
+ rdes.display()
+
+
+This is how it works:
+
+This is a ball-and-stick model with a couple of spines sitting on the dendrite.
+The spines get synaptic input onto NMDARs and gluRs. There is a baseline
+input rate of 1 Hz thoughout, and there is a burst at 40 Hz for 1 second at
+t = 10s.
+
+.. figure:: ../../../../images/ex8.2_Vm.png
+ :alt: Membrane potential responses of cell with synaptic input and multiscale signaling
+
+ Membrane potential responses of cell with synaptic input and multiscale signaling
+
+
+At baseline, we just have small EPSPs and little Ca influx. A burst of
+strong synaptic input causes Ca entry into the spine via NMDAR.
+
+.. figure:: ../../../../images/ex8.2_Ca_spine.png
+ :alt: Calcium influx into spine.
+
+ Calcium influx into spine.
+
+Ca diffuses from the spine into the dendrite and spreads. In the graph below
+we see how Calcium goes into the 50-odd voxels of the dendrite.
+
+.. figure:: ../../../../images/ex8.2_Ca_dend.png
+ :alt: Calcium influx and diffusion in dendrite.
+
+ Calcium influx and diffusion in dendrite.
+
+
+The Ca influx into the spine
+triggers activation of CaMKII and its translocation to the PSD, where
+it phosphorylates and increases the conductance of gluR. We have two spines
+with slightly different geometry, so the CaMKII activity differs slightly.
+
+.. figure:: ../../../../images/ex8.2_active_CaMKII.png
+ :alt: Activation of CaMKII and translocation to PSD
+
+ Activation of CaMKII and translocation to PSD
+
+
+Now that gluR has a greater weight, the baseline synaptic input keeps
+Ca trickling in enough to keep the CaMKII active.
+
+Here are the reactions:
+
+::
+
+ Ca+CaM <===> Ca_CaM; Ca_CaM + CaMKII <===> Ca_CaM_CaMKII (all in
+ spine head, except that the Ca_CaM_CaMKII translocates to the PSD)
+
+ chan ------Ca_CaM_CaMKII-----> chan_p; chan_p ------> chan (all in PSD)
+
+Suggestions:
+
+ - Add GABAR using make_GABA(), put it on soma or dendrite. Stimulate it
+ after 20 s to see if you can turn off the sustained activation
+ - Replace the 'periodicsyn' in stimList with 'randsyn'. This gives
+ Poisson activity at the specified mean frequency. Does the switch
+ remain reliable?
+ - What are the limits of various parameters for this switching? You
+ could try basal synaptic rate, burst rate, the various scaling factors
+ for the adaptors, the densities of various channels, synaptic weight,
+ and so on.
+ - In real life an individual synaptic EPSP is tiny, under a millivolt.
+ How many synapses would you need to achieve this kind of switching?
+ You can play with # of synapses by altering the spacing between
+ spines as the third argument of spineDistrib.
+
+Multiscale model in which spine geometry changes due to signaling
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+*ex8.3_spine_vol_change.py*
+
+This model is very similar to 8.2. The main design difference is that
+*adaptor*, instead of just modulating the gluR conductance, scales the
+entire spine cross-section area, with all sorts of electrical and chemical
+ramifications. There are a lot of plots, to illustrate some of these outcomes.
+
+::
+
+ import moose
+ import rdesigneur as rd
+ rdes = rd.rdesigneur(
+ elecDt = 50e-6,
+ chemDt = 0.002,
+ diffDt = 0.002,
+ chemPlotDt = 0.02,
+ useGssa = False,
+ stealCellFromLibrary = True, # Simply move library model to use for sim
+ cellProto = [['ballAndStick', 'soma', 12e-6, 12e-6, 4e-6, 100e-6, 2 ]],
+ chemProto = [['./chem/chanPhosph3compt.g', 'chem']],
+ spineProto = [['makeActiveSpine()', 'spine']],
+ chanProto = [
+ ['make_Na()', 'Na'],
+ ['make_K_DR()', 'K_DR'],
+ ['make_K_A()', 'K_A' ],
+ ['make_Ca()', 'Ca' ],
+ ['make_Ca_conc()', 'Ca_conc' ]
+ ],
+ passiveDistrib = [['.', 'soma', 'CM', '0.01', 'Em', '-0.06']],
+ spineDistrib = [['spine', '#dend#', '50e-6', '1e-6']],
+ chemDistrib = [['chem', '#', 'install', '1' ]],
+ chanDistrib = [
+ ['Na', 'soma', 'Gbar', '300' ],
+ ['K_DR', 'soma', 'Gbar', '250' ],
+ ['K_A', 'soma', 'Gbar', '200' ],
+ ['Ca_conc', 'soma', 'tau', '0.0333' ],
+ ['Ca', 'soma', 'Gbar', '40' ]
+ ],
+ adaptorList = [
+ # This scales the psdArea of the spine by # of chan_p. Note that
+ # the cross-section area of the spine head is identical to psdArea.
+ [ 'psd/chan_p', 'n', 'spine', 'psdArea', 0.1e-12, 0.01e-12 ],
+ [ 'Ca_conc', 'Ca', 'spine/Ca', 'conc', 0.00008, 8 ]
+ ],
+ # Syn input basline 1 Hz, and 40Hz burst for 1 sec at t=20. Syn wt=10
+ stimList = [['head#', '10','glu', 'periodicsyn', '1 + 40*(t>10 && t<11)']],
+ plotList = [
+ ['soma', '1', '.', 'Vm', 'Membrane potential'],
+ ['#', '1', 'spine/Ca', 'conc', 'Ca in Spine'],
+ ['#', '1', 'dend/DEND/Ca', 'conc', 'Ca in Dend'],
+ ['head#', '1', 'psd/chan_p', 'n', 'Amount of Phospho-chan'],
+ ['head#', '1', 'spine/CaMKII', 'conc', 'Conc of CaMKII in spine'],
+ ['head#', '1', '.', 'Cm', 'Capacitance of spine head'],
+ ['head#', '1', '.', 'Rm', 'Membrane res of spine head'],
+ ['head#', '1', '.', 'Ra', 'Axial res of spine head'],
+ ['head#', '1', 'glu', 'Gbar', 'Conductance of gluR'],
+ ['head#', '1', 'NMDA', 'Gbar', 'Conductance of NMDAR'],
+ ]
+ )
+ moose.seed(123)
+ rdes.buildModel()
+ moose.reinit()
+ moose.start( 25 )
+ rdes.display()
+
+
+The key *adaptor* line is as follows:
+
+``[ 'psd/chan_p', 'n', 'spine', 'psdArea', 0.1e-12, 0.01e-12 ]``
+
+Here, we use the phosphorylated *chan_p* molecule in the PSD as a proxy for
+processes that control spine size. We operate on a special object called
+*spine* which manages many aspects of spines in the model (see below). Here
+we control the *psdArea*, which defines the cross-section area of the spine
+head and by extension of the PSD too. We keep a minimum spine area of 0.1 um^2,
+and a scaling factor of 0.01um^2 per phosphorylated molecule.
+
+The reaction system is identical to the one in *ex8.2*:
+
+::
+
+ Ca+CaM <===> Ca_CaM; Ca_CaM + CaMKII <===> Ca_CaM_CaMKII (all in
+ spine head, except that the Ca_CaM_CaMKII translocates to the PSD)
+
+ chan ------Ca_CaM_CaMKII-----> chan_p; chan_p ------> chan (all in PSD)
+
+Rather than list all the 10 plots, here are a few to show what is going on.
+
+First, just the spiking activity of the cell. Here the burst of activity is
+followed by a few seconds of enhanced synaptic weight, followed by subthreshold
+EPSPs:
+
+.. figure:: ../../../../images/ex8.3_Vm.png
+ :alt: Membrane potential and spiking.
+
+ Membrane potential and spiking.
+
+Then, we fast-forward to the amount of *chan_p* which is the molecule that
+controls spine size scaling:
+
+.. figure:: ../../../../images/ex8.3_chan_p.png
+ :alt: Molecule that controles spine size
+
+ Molecule that controles spine size
+
+This causes some obvious outcomes. One of them is to increase the synaptic
+conductance of the glutamate receptor. The system assumes that the conductance
+of all channels in the PSD scales linearly with the psdArea.
+
+.. figure:: ../../../../images/ex8.3_gluR.png
+ :alt: Conductance of glutamate receptor
+
+ Conductance of glutamate receptor
+
+Here is one of several non-intuitive outcomes. Because the spine volume has
+increased, the concentration of molecules in the spine is diluted out. So
+the concentration of active CaMKII actually falls when the spine gets bigger.
+In a more detailed model, this would be a race between the increase in spine
+size and the time taken for diffusion and further reactions to replenish
+CaMKII. In the current model we don't have a diffusive coupling of CaMKII to
+the dendrite, so this replenishment doesn't happen.
+
+.. figure:: ../../../../images/ex8.3_CaMKII_spine.png
+ :alt: Concentration of CaMKII in the spine
+
+ Concentration of CaMKII in the spine
+
+In the simulation we display several other electrical and chemical properties
+that change with spine size. The diffusion properties also change since the
+cross-section areas are altered. This is harder to visualize but has large
+effects on coupling to the dendrite,
+especially if the *shaftDiameter* is the parameter scaled by the signaling.
+
+
+Suggestions:
+
+ - The Spine class (instance: spine) manages several possible scaling
+ targets on the spine geometry: shaftLength, shaftDiameter,
+ headLength, headDiameter, psdArea, headVolume, totalLength. Try them
+ out. Think about mechanisms by which molecular concentrations might
+ affect each.
+ - When volume changes, we assume that the molecular numbers stay
+ fixed, so concentration changes. Except for buffered molecules, where
+ we assume concentration remains fixed. Use this to design a bistable
+ simply relying on molecules and spine geometry terms.
+ - Even more interesting, use it to design an oscillator. You could look
+ at Bhalla, BiophysJ 2011 for some ideas.
+
+
+
+Morphology: Load .swc morphology file and view it
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+*ex9.0_load_neuronal_morphology_file.py*
+
+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 sys
+ import moose
+ import rdesigneur as rd
+
+ if len( sys.argv ) > 1:
+ fname = sys.argv[1]
+ else:
+ fname = './cells/h10.CNG.swc'
+ rdes = rd.rdesigneur(
+ cellProto = [[fname, '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.001, 0.1, rotation = 0.02 )
+
+
+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/ex9.0_passive_cell_morpho.png
+ :alt: 3-D display for passive neuron
+
+ 3-D display for passive neuron
+
+Suggestions:
+
+ - The tutorial directory already has a number of pre-loaded files from
+ NeuroMorpho. Pass them in to ex9.0 on the command line:
+
+ `python ex9.0_load_neuronal_morphology_file.py <morpho.swc>`
+ - Grab other morphology files from NeuroMorpho.org, try them out.
+
+Build an active neuron model by putting channels into a morphology file
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+*ex9.1_chans_in_neuronal_morph.py*
+
+Here we load in a morphology file and distribute voltage-gated ion channels
+over the neuron. 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 ]``
+
+ or
+
+ ``[ channelFunction(), 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
+ternary expressions like *p < 50e-6 ? 500 : 100* or multiplying a channel
+density with a logical function or 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
+
+Suggestions:
+
+ - Try another morphology file.
+ - Try different channel distributions by editing the chanDistrib lines.
+ - There are numerous predefined channels available within Rdesigneur.
+ These can be defined using the following chanProto options:
+
+ ::
+
+ ['make_HH_Na()', 'HH_Na']
+ ['make_HH_K_DR()', 'HH_K']
+ ['make_Na()', 'Na']
+ ['make_K_DR()', 'K_DR']
+ ['make_K_A()', 'K_A']
+ ['make_K_AHP()', 'K_AHP']
+ ['make_K_C()', 'K_C']
+ ['make_Ca()', 'Ca']
+ ['make_Ca_conc()', 'Ca_conc']
+ ['make_glu()', 'glu']
+ ['make_GABA()', 'GABA']
+
+ Then the chanDistrib can refer to these channels instead.
+ - Deliver stimuli on the dendrites rather than the soma.
+
+
+Build a spiny neuron from a morphology file and put active channels in it.
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+*ex9.2_spines_in_neuronal_morpho.py*
+
+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.
+
+.. figure:: ../../../../images/rdes9_spiny_active.png
+ :alt: 3-D display for spiny active neuron
+
+ 3-D display for spiny active neuron
+
+Suggestions:
+
+ - Try different spine settings. Warning: if you put in too many spines
+ it will take much longer to load and run!
+ - Try different spine geometry layouts.
+ - See if you can deliver the current injection to the spine. Hint: the
+ name of the spine compartments is 'head#' where # is the index of the
+ spine.
+
+
diff --git a/docs/source/user/py/rdesigneur/rdesigneur/index_rd.rst b/docs/source/user/py/rdesigneur/rdesigneur/index_rd.rst
new file mode 100644
index 0000000000000000000000000000000000000000..a729cb8b0a2f2a9d357f566bd8821085078f7571
--- /dev/null
+++ b/docs/source/user/py/rdesigneur/rdesigneur/index_rd.rst
@@ -0,0 +1,13 @@
+.. MOOSE documentation master file, created by
+ sphinx-quickstart on Tue Jul 31 19:05:47 2018.
+ You can adapt this file completely to your liking, but it should at least
+ contain the root `toctree` directive.
+
+Rdesignuer
+===========
+
+.. toctree::
+ :maxdepth: 2
+
+ rdes
+ multi_rdes
diff --git a/docs/source/user/py/rdesigneur/rdesigneur/multi_rdes.rst b/docs/source/user/py/rdesigneur/rdesigneur/multi_rdes.rst
new file mode 100644
index 0000000000000000000000000000000000000000..d240f1b1cbb2412ab0d9331eff0ccb4f25d136c7
--- /dev/null
+++ b/docs/source/user/py/rdesigneur/rdesigneur/multi_rdes.rst
@@ -0,0 +1,52 @@
+************************
+Rdesigneur Examples
+************************
+
+.. hidden-code-block:: reStructuredText
+ :label: How to run these examples
+
+ Each of the following examples can be run by clicking on the green source button
+ on the right side of each example, and running from within a ``.py`` python file
+ on a computer where moose is installed.
+
+ Alternatively, all the files mentioned on this page can be found in the main
+ moose directory. They can be found under
+
+ (...)/moose/moose-examples/snippets
+
+ They can be run by typing
+
+ $ python filename.py
+
+ in your command line, where filename.py is the python file you want to run.
+
+ All of the following examples show one or more methods within each python file.
+ For example, in the ``cubeMeshSigNeur`` section, there are two blue tabs
+ describing the ``cubeMeshSigNeur.createSquid()`` and ``cubeMeshSigNeur.main()``
+ methods.
+
+ The filename is the bit that comes before the ``.`` in the blue boxes, with
+ ``.py`` added at the end of it. In this case, the file name would be
+ ``cubeMeshSigNeur.py``.
+|
+
+Building Chemical-Electrical Signalling Models
+----------------------------------------------
+
+Building a compartment
+^^^^^^^^^^^^^^^^^^^^^^
+
+.. automodule:: testRdesigneur
+ :members:
+
+Inserting Spines and viewing
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+.. automodule:: insertSpines
+ :members:
+
+Proceeding with Spines
+^^^^^^^^^^^^^^^^^^^^^^
+
+.. automodule:: testWigglySpines
+ :members:
diff --git a/docs/source/user/py/rdesigneur/rdesigneur/rdes.rst b/docs/source/user/py/rdesigneur/rdesigneur/rdes.rst
new file mode 100644
index 0000000000000000000000000000000000000000..161ff2f97e15771e76ca1d6f81315e8eb3ca4743
--- /dev/null
+++ b/docs/source/user/py/rdesigneur/rdesigneur/rdes.rst
@@ -0,0 +1,1964 @@
+==========================================
+Rdesigneur: Building multiscale models
+==========================================
+
+| Author: Upi Bhalla
+| Date: Aug 26 2016,
+| Last-Updated: July 31 2018
+| By: Upi Bhalla
+
+------------------------------
+
+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
+- Adaptors that couple between these for multiscale models
+
+It also folds in simulation input and output
+
+- Time-series stimuli for molecular concentration change and reaction rates
+- Current and voltage clamp
+- Synaptic input.
+- Time-series plots
+- File dumps
+- 3-D neuronal graphics
+
+Rdesigneur's main role is to specify how these are put together,
+including assigning parameters for the model. Using Rdesigneur one can compactly
+and quickly put together quite complex multiscale models.
+
+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. The files for these examples are also available in
+``moose-examples/tutorials/Rdesigneur``, and the file names are mentioned
+as we go along.
+
+Bare Rdesigneur: single passive compartment
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+*ex1_minimalModel.py*
+
+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
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+*ex2.0_currentPulse.py*
+
+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.
+
+Simulate and display voltage clamp stimulus to soma
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+*ex2.1_vclamp.py*
+
+This model introduces the voltage clamp stimulus on a passive compartment.
+As before, we add a few lines to define the stimulus and plot.
+This script displays both the membrane potential, and the holding current
+of the voltage clamp circuit as
+it charges and discharges the passive compartment model.
+
+::
+
+ import moose
+ import rdesigneur as rd
+ rdes = rd.rdesigneur(
+ stimList = [['soma', '1', '.', 'vclamp', '-0.065 + (t>0.1 && t<0.2) * 0.02' ]],
+ plotList = [
+ ['soma', '1', '.', 'Vm', 'Soma membrane potential'],
+ ['soma', '1', 'vclamp', 'current', 'Soma holding current'],
+ ]
+ )
+ rdes.buildModel()
+ moose.reinit()
+ moose.start( 0.3 )
+ rdes.display()
+
+Here the *stimList* line tells the system to deliver a voltage clamp (vclamp)
+on the soma, starting at -65 mV and jumping up by 20 mV between 0.1 and 0.2
+seconds. The *plotList* now includes two entries, and will generate two plots.
+The first is for plotting the soma membrane potential, just to be sure that
+the voltage clamp is doing its job.
+
+.. figure:: ../../../../images/ex2.1_vclamp_a.png
+ :alt: Plot for membrane potential in voltage clamp
+
+ Plot for membrane potential in voltage clamp
+
+The second graph plots the holding current. Note the capacitive transients.
+
+.. figure:: ../../../../images/ex2.1_vclamp_b.png
+ :alt: Plot for holding current for voltage clamp
+
+ Plot for holding current for voltage clamp
+
+HH Squid model in a single compartment
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+*ex3.0_squid_currentPulse.py*
+
+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
+
+There are several interesting things to do with the model by varying stimulus
+parameters:
+
+ - Change injection current.
+ - Put in a protocol to get rebound action potential.
+ - Put in a current ramp, and run it for a different duration
+ - Put in a frequency chirp, and see how the squid model is tuned
+ to a certain frequency range.
+ - Modify channel or passive parameters. See if it still fires.
+ - Try the frequency chirp on the cell with parameters changed. Does
+ the tuning change?
+
+
+HH Squid model with voltage clamp
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+*ex3.1_squid_vclamp.py*
+
+This is the same squid model, but now we add a voltage clamp to the squid
+and monitor the holding current. This stimulus line is identical to ex2.1.
+
+::
+
+ 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', '.', 'vclamp', '-0.065 + (t>0.1 && t<0.2) * 0.02' ]],
+ plotList = [
+ ['soma', '1', '.', 'Vm', 'Membrane potential'],
+ ['soma', '1', 'vclamp', 'current', 'Soma holding current']
+ ]
+ )
+ rdes.buildModel()
+ moose.reinit()
+ moose.start( 0.3 )
+ rdes.display()
+
+Here we see the classic HH current response, a downward brief deflection due to
+the Na channel, and a slower upward sustained current due to the K delayed
+rectifier.
+
+.. figure:: ../../../../images/ex3.1_squid_vclamp.png
+ :alt: Plot for HH squid voltage clamp pulse.
+
+ Plot for HH squid voltage clamp pulse.
+
+Here are some suggestions for further exploration:
+
+ - Monitor individual channel currents through additional plots.
+ - Convert this into a voltage clamp series. Easiest way to do this is
+ to complete the rdes.BuildModel, then delete the Function object
+ on the */model/elec/soma/vclamp*. Now you can simply set the 'command'
+ field of the vclamp in a for loop, going from -ve to +ve voltages.
+ Remember, SI units. You may wish to capture the plot vectors each
+ cycle. The plot vectors are accessed by something like
+
+ ``moose.element( '/model/graphs/plot1' ).vector``
+
+
+HH Squid model in an axon
+~~~~~~~~~~~~~~~~~~~~~~~~~
+
+*ex3.2_squid_axon_propgn.py*
+
+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/ex3.2_axon_propagating_AP.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.
+
+Action potential collision in HH Squid axon model
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+*ex3.3_AP_collision.py*
+
+This is identical to the previous example, except that now we deliver current
+injection at at two points, the soma and a point along the axon. The modified
+stimulus line is:
+
+::
+
+ ...
+ stimList = [['soma', '1', '.', 'inject', '(t>0.01 && t<0.2) * 2e-11' ],
+ ['axon100', '1', '.', 'inject', '(t>0.01 && t<0.2) * 3e-11' ]],
+ ...
+
+Watch how the AP is triggered bidirectionally from the stimulus point on the
+100th segment of the axon, and observe what happens when two action potentials
+bump into each other.
+
+.. figure:: ../../../../images/ex3.3_AP_collision.png
+ :alt: Colliding action potentials
+
+ Colliding action potentials
+
+
+
+HH Squid model in a myelinated axon
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+*ex3.4_myelinated_axon.py*
+
+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
+
+Alternate (non-squid) way to define soma
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+*ex4.0_scaledSoma.py*
+
+The default HH-squid axon is not a very convincing soma. Rdesigneur offers a
+somewhat more general way to define the soma in the cell prototype line.
+
+::
+
+ import moose
+ import pylab
+ import rdesigneur as rd
+ rdes = rd.rdesigneur(
+ # cellProto syntax: ['somaProto', 'name', dia, length]
+ cellProto = [['somaProto', 'soma', 20e-6, 200e-6]],
+ chanProto = [['make_HH_Na()', 'Na'], ['make_HH_K()', 'K']],
+ chanDistrib = [
+ ['Na', 'soma', 'Gbar', '1200' ],
+ ['K', 'soma', 'Gbar', '360' ]],
+ stimList = [['soma', '1', '.', 'inject', '(t>0.01 && t<0.05) * 1e-9' ]],
+ plotList = [['soma', '1', '.', 'Vm', 'Membrane potential']],
+ moogList = [['#', '1', '.', 'Vm', 'Vm (mV)']]
+ )
+
+ rdes.buildModel()
+ soma = moose.element( '/model/elec/soma' )
+ print( 'Soma dia = {}, length = {}'.format( soma.diameter, soma.length ) )
+ moose.reinit()
+
+ rdes.displayMoogli( 0.0005, 0.06, 0.0 )
+
+Here the crucial line is the *cellProto* line. There are four arguments here:
+
+ ``['somaProto', 'name', dia, length]``
+
+ - The first argument tells the system to use a prototype soma, that is
+ a single cylindrical compartment.
+ - The second argument is the name to give the cell.
+ - The third argument is the diameter. Note that this is a double,
+ not a string.
+ - The fourth argument is the length of the cylinder that makes up the
+ soma. This too is a double, not a string.
+ The cylinder is oriented along the x axis, with one end at (0,0,0)
+ and the other end at (length, 0, 0).
+
+This is what the soma looks like:
+
+.. figure:: ../../../../images/ex4.0_scaledSoma.png
+ :alt: Image of soma.
+
+ Image of soma.
+
+It a somewhat elongated soma, being a cylinder 10 times as long as it is wide.
+
+Ball-and-stick model of a neuron
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+*ex4.1_ballAndStick.py*
+
+A somewhat more electrically reasonable model of a neuron has a soma and a
+single dendrite, which can itself be subdivided into segments so that it
+can exhibit voltage gradients, have channel and receptor distributions,
+and so on. This is accomplished in *rdesigneur* using a variant of the
+cellProto syntax.
+
+::
+
+ import moose
+ import pylab
+ import rdesigneur as rd
+ rdes = rd.rdesigneur(
+ # cellProto syntax: ['ballAndStick', 'name', somaDia, somaLength, dendDia, dendLength, numDendSegments ]
+ # The numerical arguments are all optional
+ cellProto = [['ballAndStick', 'soma', 20e-6, 20e-6, 4e-6, 500e-6, 10]],
+ chanProto = [['make_HH_Na()', 'Na'], ['make_HH_K()', 'K']],
+ chanDistrib = [
+ ['Na', 'soma', 'Gbar', '1200' ],
+ ['K', 'soma', 'Gbar', '360' ],
+ ['Na', 'dend#', 'Gbar', '400' ],
+ ['K', 'dend#', 'Gbar', '120' ]
+ ],
+ stimList = [['soma', '1', '.', 'inject', '(t>0.01 && t<0.05) * 1e-9' ]],
+ plotList = [['soma', '1', '.', 'Vm', 'Membrane potential']],
+ moogList = [['#', '1', '.', 'Vm', 'Vm (mV)']]
+ )
+ rdes.buildModel()
+ soma = moose.element( '/model/elec/soma' )
+ moose.reinit()
+ rdes.displayMoogli( 0.0005, 0.06, 0.0 )
+
+As before, the *cellProto* line plays a key role. Here, because we have a long
+dendrite, we have a few more numerical arguments. All of the numerical
+arguments are optional.
+
+ ``['ballAndStick', 'name', somaDia, somaLength, dendDia, dendLength, numDendSegments ]``
+
+ - The first argument specifies a ballAndStick model: soma + dendrite.
+ The length of the dendrite is along the x axis. The soma is a single
+ segment, the dendrite can be more than one.
+ - The second argument is the name to give the cell.
+ - Arg 3 is the soma diameter, as a double.
+ - Arg 4 is the length of the soma, as a double.
+ - Arg 5 is the diameter of the dendrite, as a double.
+ - Arg 6 is the length of the dendrite, as a double.
+ - Arg 7 is the number of segments into which the dendrite should be
+ divided. This is a positive integer greater than 0.
+
+This is what the ball-and-stick cell looks like:
+
+.. figure:: ../../../../images/ex4.1_ballAndStick.png
+ :alt: Image of ball and stick cell.
+
+ Image of ball and stick cell.
+
+In this version of the 3-D display, the soma is displayed as a bit blocky
+rather than round.
+Note that we have populated the dendrite with Na and K channels and it has
+10 segments, so it supports action potential propagation. The snapshot
+illustrates this.
+
+Here are some things to try:
+
+ - Change the length of the dendrite
+ - Change the number of segments. Explore what it does to accuracy. How
+ will you know that you have an accurate model?
+
+Benchmarking simulation speed
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+*ex4.2_ballAndStickSpeed.py*
+
+The ball-and-stick model gives us an opportunity to check out your system
+and how computation scales with model size. While we're at it we'll deliver
+a sine-wave stimulus just to see how it can be done. The test model is
+very similar to the previous one, ex4.1:
+
+::
+
+ import moose
+ import pylab
+ import rdesigneur as rd
+ import time
+ rdes = rd.rdesigneur(
+ cellProto = [['ballAndStick', 'soma', 20e-6, 20e-6, 4e-6, 500e-6, 10]],
+ chanProto = [['make_HH_Na()', 'Na'], ['make_HH_K()', 'K']],
+ chanDistrib = [
+ ['Na', 'soma', 'Gbar', '1200' ],
+ ['K', 'soma', 'Gbar', '360' ],
+ ['Na', 'dend#', 'Gbar', '400' ],
+ ['K', 'dend#', 'Gbar', '120' ]
+ ],
+ stimList = [['soma', '1', '.', 'inject', '(1+cos(t/10))*(t>31.4 && t<94) * 0
+ .2e-9' ]],
+ plotList = [
+ ['soma', '1', '.', 'Vm', 'Membrane potential'],
+ ['soma', '1', '.', 'inject', 'Stimulus current']
+ ],
+ )
+ rdes.buildModel()
+ runtime = 100
+ moose.reinit()
+ t0= time.time()
+ moose.start( runtime )
+ print "Real time to run {} simulated seconds = {} seconds".format( runtime, time
+ .time() - t0 )
+
+ rdes.display()
+
+While the real point of this simulation is to check speed, it does illustrate
+how to deliver a stimulus shaped like a sine wave:
+
+.. figure:: ../../../../images/ex4.2_sine_stim.png
+ :alt: Sine-wave shaped stimulus.
+
+ Sine-wave shaped stimulus.
+
+We can see that the cell has a peculiar response to this. Not surprising, as
+the cell uses HH channels which are not good at rate coding.
+
+.. figure:: ../../../../images/ex4.2_spiking.png
+ :alt: Spiking response to sine-wave shaped stimulus.
+
+ Spiking response to sine-wave shaped stimulus.
+
+As a reference point, on a fast 2018 laptop this benchmark runs in 5.4 seconds.
+Some more things to try for benchmarking:
+
+ - How slow does it get if you turn on the 3-D moogli display?
+ - Is it costlier to run 2 compartments for 1000 seconds, or
+ 200 compartments for 10 seconds?
+
+Synaptic stimulus: random (Possion)
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+*ex5.0_random_syn_input.py*
+
+In this example we introduce synaptic inputs: both the receptor channels
+and a means for stimulating the channels. We do this in a passive model.
+
+::
+
+ import moose
+ import rdesigneur as rd
+ rdes = rd.rdesigneur(
+ cellProto = [['somaProto', 'soma', 20e-6, 200e-6]],
+ chanProto = [['make_glu()', 'glu']],
+ chanDistrib = [['glu', 'soma', 'Gbar', '1' ]],
+ stimList = [['soma', '0.5', 'glu', 'randsyn', '50' ]],
+ # Deliver stimulus to glu synapse on soma, at mean 50 Hz Poisson.
+ plotList = [['soma', '1', '.', 'Vm', 'Soma membrane potential']]
+ )
+ rdes.buildModel()
+ moose.reinit()
+ moose.start( 0.3 )
+ rdes.display()
+
+Most of the rdesigneur setup uses familiar syntax.
+
+Novelty 1: we use the default built-in glutamate receptor model, in chanProto.
+We just put it in the soma at a max conductance of 1 Siemen/sq metre.
+
+Novelty 2: We specify a new kind of stimulus in the stimList:
+
+ ``['soma', '0.5', 'glu', 'randsyn', '50' ]``
+
+Most of this is similar to previous stimLists.
+
+ - arg0: 'soma': the named compartments in the cell to populate with
+ the *glu* receptor
+ - arg1: '0.5': Tell the system to use a uniform synaptic weight of 0.5.
+ This argument could be a more complicated expression incorporating
+ spatial arguments. Here it is just uniform.
+ - arg2: 'glu': Which receptor to stimulate
+ - arg3: 'randsyn': Apply random (Poisson) synaptic input.
+ - arg4: '50': Mean firing rate of the Poisson input. Note that this last
+ argument could be a function of time and hence is quite versatile.
+
+As the model has no voltage-gated channels, we do not see spiking.
+
+.. figure:: ../../../../images/ex5.0_random_syn_input.png
+ :alt: Random synaptic input with a Poisson distribution.
+
+ Random synaptic input with a Poisson distribution.
+
+Things to try: Vary the rate and the weight of the synaptic input.
+
+Synaptic stimulus: periodic
+~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+*ex5.1_periodic_syn_input.py*
+
+This is almost identical to 5.0, except that the input is now perfectly
+periodic. The one change is of an argument in the stimList to say
+``periodicsyn`` rather than ``randsyn``.
+
+::
+
+ import moose
+ import rdesigneur as rd
+ rdes = rd.rdesigneur(
+ cellProto = [['somaProto', 'soma', 20e-6, 200e-6]],
+ chanProto = [['make_glu()', 'glu']],
+ chanDistrib = [['glu', 'soma', 'Gbar', '1' ]],
+
+ # Deliver stimulus to glu synapse on soma, periodically at 50 Hz.
+ stimList = [['soma', '0.5', 'glu', 'periodicsyn', '50' ]],
+ plotList = [['soma', '1', '.', 'Vm', 'Soma membrane potential']]
+ )
+ rdes.buildModel()
+ moose.reinit()
+ moose.start( 0.3 )
+ rdes.display()
+
+As designed, we get periodically firing synaptic input.
+
+.. figure:: ../../../../images/ex5.1_periodic_syn_input.png
+ :alt: Periodic synaptic input
+
+ Periodic synaptic input
+
+
+Reaction system in a single compartment
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+*ex6_chem_osc.py*
+
+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. Its general form is:
+
+::
+
+ s ---a---> a // s goes to a, catalyzed by a.
+ s ---a---> b // s goes to b, catalyzed by a.
+ a ---b---> s // a goes to s, catalyzed by b.
+ b -------> s // b is degraded irreversibly to s
+
+Here is the script:
+
+::
+
+ 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
+~~~~~~~~~~~~~~~~~~~~~~~~~
+
+*ex7.0_spatial_chem_osc.py*
+
+In order to see what a reaction-diffusion system looks like, we assign the
+``diffusionLength`` expression in the previous example to a much shorter
+length, 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,
+ #This subdivides the length of the soma into 2 micron voxels
+ diffusionLength = 2e-6,
+ 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, rotation = 0, azim = np.pi/2, elev = 0.0 )
+
+This is the new value for diffusion length.
+
+::
+
+ diffusionLength = 2e-3,
+
+With this change we tell *rdesigneur* to use the diffusion length of 2 microns.
+This happens to be the default too. 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, rotation = 0, azim = np.pi/2, elev = 0.0 )
+
+The arguments mean: *displayMoogli( frametime, runtime, rotation, azimuth, elevation )*
+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.
+ azimuth = Azimuth angle of view point, in radians
+ elevation = elevation angle of view point, 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
+
+For those who would rather use the much simpler matplotlib 3-D display option,
+this is what the same simulation looks like:
+
+.. figure:: ../../../../images/ex7.0_spatial_chem_osc.png
+ :alt: Display for oscillatory reac-diff simulation using matplotlib
+
+ Display for oscillatory reac-diff simulation using matplotlib
+
+Primer on using the 3-D MOOGLI display
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+There are two variants of the MOOGLI display. The first, named Moogli,
+uses OpenGL and OpenSceneGraph. It is fast to display, slow to load, and
+difficult to compile. It produces much better looking 3-D graphics.
+The second is a fallback interface using mplot3d, which is a library of
+Matplotlib and so should be generally available. It is slower to display,
+faster to load, but needs no special compilation. It uses stick graphics
+and though it conveys much the same information, isn't as nice to look at
+as the original Moogli. Its controls are more or less the same but less
+smooth than the original Moogli.
+
+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.
+
+Diffusion of a single molecule
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+*ex7.1_diffusive_gradient.py*
+
+This is simply a test model to confirm that simple diffusion happens as
+expected. While the model is just that of a single pool, we spend a few lines
+taking snapshots of the spatial profile of this pool.
+
+::
+
+ import moose
+ import pylab
+ import re
+ import rdesigneur as rd
+ import matplotlib.pyplot as plt
+ import numpy as np
+
+ moose.Neutral( '/library' )
+ moose.Neutral( '/library/diffn' )
+ moose.CubeMesh( '/library/diffn/dend' )
+ A = moose.Pool( '/library/diffn/dend/A' )
+ A.diffConst = 1e-10
+
+ rdes = rd.rdesigneur(
+ turnOffElec = True,
+ diffusionLength = 1e-6,
+ chemProto = [['diffn', 'diffn']],
+ chemDistrib = [['diffn', 'soma', 'install', '1' ]],
+ moogList = [
+ ['soma', '1', 'dend/A', 'conc', 'A Conc', 0, 360 ]
+ ]
+ )
+ rdes.buildModel()
+
+ rdes.displayMoogli( 1, 2, rotation = 0, azim = -np.pi/2, elev = 0.0, block = False )
+ av = moose.vec( '/model/chem/dend/A' )
+ for i in range(10):
+ av[i].concInit = 1
+ moose.reinit()
+ plist = []
+ for i in range( 20 ):
+ plist.append( av.conc[:200] )
+ moose.start( 2 )
+ fig = plt.figure( figsize = ( 10, 12 ) )
+ plist = np.array( plist ).T
+ plt.plot( range( 0, 200 ), plist )
+ plt.xlabel( "position ( microns )" )
+ plt.ylabel( "concentration ( mM )" )
+ plt.show( block = True )
+
+
+Here are the snapshots, overlaid in a single plot:
+
+.. figure:: ../../../../images/ex7.1_diffusive_gradient.png
+ :alt: Display of how a molecule A spreads through the inter
+
+ Display for simple time-series of spread of a diffusing molecule
+ using matplotlib
+
+
+Multiscale models: single compartment
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+*ex8.0_multiscale_KA_phosph.py*
+
+The next few examples are for the multiscale modeling that is the main purpose
+of rdesigneur and MOOSE as a whole. These are 'toy' examples in that the
+chemical and electrical signaling is simplified, but they exhibit dynamics
+that are of real interest.
+
+The first example is of a bistable system where the feedback loop comprises of
+
+`calcium influx -> chemical activity -> channel modulation -> electrical activity -> calcium influx.`
+
+Calcium enters through voltage gated calcium channels, leads to enzyme
+activation and phosphorylation of a KA channel, which depolarizes the cell,
+so it spikes more, so more calcium enters.
+
+::
+
+ import moose
+ import pylab
+ import rdesigneur as rd
+ rdes = rd.rdesigneur(
+ elecDt = 50e-6,
+ chemDt = 0.002,
+ chemPlotDt = 0.002,
+ # cellProto syntax: ['somaProto', 'name', dia, length]
+ cellProto = [['somaProto', 'soma', 12e-6, 12e-6]],
+ chemProto = [['./chem/chanPhosphByCaMKII.g', 'chem']],
+ chanProto = [
+ ['make_Na()', 'Na'],
+ ['make_K_DR()', 'K_DR'],
+ ['make_K_A()', 'K_A' ],
+ ['make_Ca()', 'Ca' ],
+ ['make_Ca_conc()', 'Ca_conc' ]
+ ],
+ # Some changes to the default passive properties of the cell.
+ passiveDistrib = [['.', 'soma', 'CM', '0.03', 'Em', '-0.06']],
+ chemDistrib = [['chem', 'soma', 'install', '1' ]],
+ chanDistrib = [
+ ['Na', 'soma', 'Gbar', '300' ],
+ ['K_DR', 'soma', 'Gbar', '250' ],
+ ['K_A', 'soma', 'Gbar', '200' ],
+ ['Ca_conc', 'soma', 'tau', '0.0333' ],
+ ['Ca', 'soma', 'Gbar', '40' ]
+ ],
+ adaptorList = [
+ [ 'dend/chan', 'conc', 'K_A', 'modulation', 0.0, 70 ],
+ [ 'Ca_conc', 'Ca', 'dend/Ca', 'conc', 0.00008, 2 ]
+ ],
+ # Give a + pulse from 5 to 7s, and a - pulse from 20 to 21.
+ stimList = [['soma', '1', '.', 'inject', '((t>5 && t<7) - (t>20 && t<21)) * 1.0e-12' ]],
+ plotList = [
+ ['soma', '1', '.', 'Vm', 'Membrane potential'],
+ ['soma', '1', '.', 'inject', 'current inj'],
+ ['soma', '1', 'K_A', 'Ik', 'K_A current'],
+ ['soma', '1', 'dend/chan', 'conc', 'Unphosph K_A conc'],
+ ['soma', '1', 'dend/Ca', 'conc', 'Chem Ca'],
+ ],
+ )
+
+ rdes.buildModel()
+ moose.reinit()
+ moose.start( 30 )
+
+ rdes.display()
+
+There is only one fundamentally new element in this script:
+
+**adaptor List:** `[source, sourceField, dest, destField, offset, scale]`
+The adaptor list maps between molecular, electrical or even structural
+quantities in the simulation. At present it is linear mapping, in due course
+it may evolve to an arbitrary function.
+
+The two adaptorLists in the above script do the following:
+
+ ``[ 'dend/chan', 'conc', 'K_A', 'modulation', 0.0, 70 ]``:
+
+Use the concentration of the 'chan' molecule in the 'dend' compartment,
+to modulate the conductance of the 'K_A' channel such that the basal
+conductance is zero and 1 millimolar of 'chan' results in a conductance that is
+70 times greater than the baseline conductance of the channel, *Gbar*.
+
+It is advisable to use the field *'modulation'* on channels undergoing scaling,
+rather than to directly assign the conductance *'Gbar'*. This is because
+*Gbar* is an absolute conductance, and therefore it is scaled to the area of
+the electrical segment. This makes it difficult to keep track of. *Modulation*
+is a simple multiplier term onto *Gbar*, and is therefore easier to work with.
+
+ ``[ 'Ca_conc', 'Ca', 'dend/Ca', 'conc', 0.00008, 2 ]``:
+
+Use the concentration of *Ca* as computed in the electrical model, to assign
+the concentration of molecule *Ca* on the dendrite compartment. There is a
+basal level of 80 nanomolar, and every unit of electrical *Ca* maps to 2
+millimolar of chemical *Ca*.
+
+The arguments in the adaptorList are:
+
+ * **Source and Dest**: Strings. These can be either a molecular or an
+ electrical object. To identify a molecular object, it should be
+ prefixed with the name of the chemical compartment, which is one
+ of *dend, spine, psd*. Thus *dend/chan* specifies a molecule
+ named *'chan'* sitting in the *'dend'* compartment.
+
+ To identify an electrical object, just pass in its path,
+ such as '.' or *'Ca_conc'*.
+
+ Note that the adaptors do **not** need to know anything about the
+ location. It is assumed that the adaptors do their job wherever
+ the specified source and dest coexist. There is a subtlety here
+ due to the different length and time scales. The rule of thumb
+ is that the adaptor averages whichever one is subdivided more finely.
+
+ - Example 1: Molecules are typically spatially partitioned into
+ short voxels (micron-scale) compared to typical 100-micron
+ electrical
+ segments. So an adaptor going from molecules to, say, channel
+ conductance, would average all the molecular voxels that fit
+ in the electrical segment.
+ - Example 2: Electrical activity is typically much faster than
+ chemical.
+ So an adaptor going from an electrical entity (Ca computed from
+ channel opening) to molecules (Chemical Ca concentration) would
+ average all the time-steps between updates to the molecule.
+
+ * **Fields**: Strings. These are simply the field names on the
+ objects coupled by the adaptors.
+
+ * **offset and scale**: Doubles. At present the adaptor is just a
+ straight-line conversion, obeying ``y = mx + c``. The computed
+ output is *y*, averaged input is *x*, offset is *c* and scale is *m*.
+
+There is a handy new line to specify cellular passive properties:
+
+**passiveDistrib:** `['.', path, field, value, field, value, ... ]`,
+
+ * '.': This is just a placeholder.
+ * path: String. Specifies the object whose parameters are to be changed.
+ * field: String. Name of the field on the object.
+ * value: String, that is the value has to be enclosed in quotes. The
+ value to be assigned to the object.
+
+With these in place, the model behavior is rather neat. It starts out silent,
+then we apply 2 seconds of +ve current injection.
+
+.. figure:: ../../../../images/ex8.0_multiscale_currInj.png
+ :alt: Current injection stimuli for multiscale model.
+
+ Current injection stimuli for multiscale model.
+
+The cell fires briskly, and keeps firing even when the current injection
+drops to zero.
+
+.. figure:: ../../../../images/ex8.0_multiscale_cell_spiking.png
+ :alt: Firing responses of cell with multiscale signaling.
+
+ Firing responses of cell with multiscale signaling.
+
+The firing of the neuron leads to Ca influx.
+
+.. figure:: ../../../../images/ex8.0_multiscale_Ca.png
+ :alt: Calcium buildup in cell due to firing.
+
+ Calcium buildup in cell due to firing.
+
+The chemical reactions downstream of Ca lead to phosphorylation of the K_A
+channel. Only the unphosphorylated K_A channel is active, so the net effect
+is to reduce K_A conductance while the Ca influx persists.
+
+.. figure:: ../../../../images/ex8.0_multiscale_KA_conc.png
+ :alt: Removal of KA channel due to phosphorylation.
+
+ Removal of KA channel due to phosphorylation.
+
+
+Since the phosphorylated form has low conductance, the cell becomes more
+excitable and keeps firing even when the current injection is stopped. It takes
+a later, -ve current injection to turn the firing off again.
+
+Suggestions for things to do with the model:
+
+ - Vary the adaptor settings, which couple electrical to chemical
+ signaling and vice versa.
+ - Play with the channel densities
+ - Open the chem model in moosegui and vary its parameters too.
+
+Multiscale model spanning PSD, spine head and dendrite
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+*ex8.2_multiscale_glurR_phosph_3compt.py*
+
+This is another multiscale model on similar lines to 8.0. It is structurally
+and computationally more complicated, because the action is distributed between
+spines and dendrites, but formally it does the same thing: it turns on and
+stays on after a strong stimulus, due to phosphorylation of a (receptor)
+channel leading to greater excitability.
+
+`calcium influx -> chemical activity -> channel modulation -> electrical activity -> calcium influx.`
+
+The model is bistable as long as synaptic input keeps coming along at a basal
+rate, in this case 1 Hz.
+
+Here we have two new lines, to do with addition of spines. These are discussed
+in detail in a later example. For now it is enough to know that the
+**spineProto** line defines one of the prototype spines to be used to put into
+the model, and the **spineDistrib** line tells the system where to put them,
+and how widely to space them.
+
+::
+
+ import moose
+ import rdesigneur as rd
+ rdes = rd.rdesigneur(
+ elecDt = 50e-6,
+ chemDt = 0.002,
+ diffDt = 0.002,
+ chemPlotDt = 0.02,
+ useGssa = False,
+ # cellProto syntax: ['ballAndStick', 'name', somaDia, somaLength, dendDia, d
+ endLength, numDendSegments ]
+ cellProto = [['ballAndStick', 'soma', 12e-6, 12e-6, 4e-6, 100e-6, 2 ]],
+ chemProto = [['./chem/chanPhosph3compt.g', 'chem']],
+ spineProto = [['makeActiveSpine()', 'spine']],
+ chanProto = [
+ ['make_Na()', 'Na'],
+ ['make_K_DR()', 'K_DR'],
+ ['make_K_A()', 'K_A' ],
+ ['make_Ca()', 'Ca' ],
+ ['make_Ca_conc()', 'Ca_conc' ]
+ ],
+ passiveDistrib = [['.', 'soma', 'CM', '0.01', 'Em', '-0.06']],
+ spineDistrib = [['spine', '#dend#', '50e-6', '1e-6']],
+ chemDistrib = [['chem', '#', 'install', '1' ]],
+ chanDistrib = [
+ ['Na', 'soma', 'Gbar', '300' ],
+ ['K_DR', 'soma', 'Gbar', '250' ],
+ ['K_A', 'soma', 'Gbar', '200' ],
+ ['Ca_conc', 'soma', 'tau', '0.0333' ],
+ ['Ca', 'soma', 'Gbar', '40' ]
+ ],
+ adaptorList = [
+ [ 'psd/chan_p', 'n', 'glu', 'modulation', 0.1, 1.0 ],
+ [ 'Ca_conc', 'Ca', 'spine/Ca', 'conc', 0.00008, 8 ]
+ ],
+ # Syn input basline 1 Hz, and 40Hz burst for 1 sec at t=20. Syn weight
+ # is 0.5, specified in 2nd argument as a special case stimLists.
+ stimList = [['head#', '0.5','glu', 'periodicsyn', '1 + 40*(t>10 && t<11)']],
+ plotList = [
+ ['soma', '1', '.', 'Vm', 'Membrane potential'],
+ ['#', '1', 'spine/Ca', 'conc', 'Ca in Spine'],
+ ['#', '1', 'dend/DEND/Ca', 'conc', 'Ca in Dend'],
+ ['#', '1', 'spine/Ca_CaM', 'conc', 'Ca_CaM'],
+ ['head#', '1', 'psd/chan_p', 'conc', 'Phosph gluR'],
+ ['head#', '1', 'psd/Ca_CaM_CaMKII', 'conc', 'Active CaMKII'],
+ ]
+ )
+ moose.seed(123)
+ rdes.buildModel()
+ moose.reinit()
+ moose.start( 25 )
+ rdes.display()
+
+
+This is how it works:
+
+This is a ball-and-stick model with a couple of spines sitting on the dendrite.
+The spines get synaptic input onto NMDARs and gluRs. There is a baseline
+input rate of 1 Hz thoughout, and there is a burst at 40 Hz for 1 second at
+t = 10s.
+
+.. figure:: ../../../../images/ex8.2_Vm.png
+ :alt: Membrane potential responses of cell with synaptic input and multiscale signaling
+
+ Membrane potential responses of cell with synaptic input and multiscale signaling
+
+
+At baseline, we just have small EPSPs and little Ca influx. A burst of
+strong synaptic input causes Ca entry into the spine via NMDAR.
+
+.. figure:: ../../../../images/ex8.2_Ca_spine.png
+ :alt: Calcium influx into spine.
+
+ Calcium influx into spine.
+
+Ca diffuses from the spine into the dendrite and spreads. In the graph below
+we see how Calcium goes into the 50-odd voxels of the dendrite.
+
+.. figure:: ../../../../images/ex8.2_Ca_dend.png
+ :alt: Calcium influx and diffusion in dendrite.
+
+ Calcium influx and diffusion in dendrite.
+
+
+The Ca influx into the spine
+triggers activation of CaMKII and its translocation to the PSD, where
+it phosphorylates and increases the conductance of gluR. We have two spines
+with slightly different geometry, so the CaMKII activity differs slightly.
+
+.. figure:: ../../../../images/ex8.2_active_CaMKII.png
+ :alt: Activation of CaMKII and translocation to PSD
+
+ Activation of CaMKII and translocation to PSD
+
+
+Now that gluR has a greater weight, the baseline synaptic input keeps
+Ca trickling in enough to keep the CaMKII active.
+
+Here are the reactions:
+
+::
+
+ Ca+CaM <===> Ca_CaM; Ca_CaM + CaMKII <===> Ca_CaM_CaMKII (all in
+ spine head, except that the Ca_CaM_CaMKII translocates to the PSD)
+
+ chan ------Ca_CaM_CaMKII-----> chan_p; chan_p ------> chan (all in PSD)
+
+Suggestions:
+
+ - Add GABAR using make_GABA(), put it on soma or dendrite. Stimulate it
+ after 20 s to see if you can turn off the sustained activation
+ - Replace the 'periodicsyn' in stimList with 'randsyn'. This gives
+ Poisson activity at the specified mean frequency. Does the switch
+ remain reliable?
+ - What are the limits of various parameters for this switching? You
+ could try basal synaptic rate, burst rate, the various scaling factors
+ for the adaptors, the densities of various channels, synaptic weight,
+ and so on.
+ - In real life an individual synaptic EPSP is tiny, under a millivolt.
+ How many synapses would you need to achieve this kind of switching?
+ You can play with # of synapses by altering the spacing between
+ spines as the third argument of spineDistrib.
+
+Multiscale model in which spine geometry changes due to signaling
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+*ex8.3_spine_vol_change.py*
+
+This model is very similar to 8.2. The main design difference is that
+*adaptor*, instead of just modulating the gluR conductance, scales the
+entire spine cross-section area, with all sorts of electrical and chemical
+ramifications. There are a lot of plots, to illustrate some of these outcomes.
+
+::
+
+ import moose
+ import rdesigneur as rd
+ rdes = rd.rdesigneur(
+ elecDt = 50e-6,
+ chemDt = 0.002,
+ diffDt = 0.002,
+ chemPlotDt = 0.02,
+ useGssa = False,
+ stealCellFromLibrary = True, # Simply move library model to use for sim
+ cellProto = [['ballAndStick', 'soma', 12e-6, 12e-6, 4e-6, 100e-6, 2 ]],
+ chemProto = [['./chem/chanPhosph3compt.g', 'chem']],
+ spineProto = [['makeActiveSpine()', 'spine']],
+ chanProto = [
+ ['make_Na()', 'Na'],
+ ['make_K_DR()', 'K_DR'],
+ ['make_K_A()', 'K_A' ],
+ ['make_Ca()', 'Ca' ],
+ ['make_Ca_conc()', 'Ca_conc' ]
+ ],
+ passiveDistrib = [['.', 'soma', 'CM', '0.01', 'Em', '-0.06']],
+ spineDistrib = [['spine', '#dend#', '50e-6', '1e-6']],
+ chemDistrib = [['chem', '#', 'install', '1' ]],
+ chanDistrib = [
+ ['Na', 'soma', 'Gbar', '300' ],
+ ['K_DR', 'soma', 'Gbar', '250' ],
+ ['K_A', 'soma', 'Gbar', '200' ],
+ ['Ca_conc', 'soma', 'tau', '0.0333' ],
+ ['Ca', 'soma', 'Gbar', '40' ]
+ ],
+ adaptorList = [
+ # This scales the psdArea of the spine by # of chan_p. Note that
+ # the cross-section area of the spine head is identical to psdArea.
+ [ 'psd/chan_p', 'n', 'spine', 'psdArea', 0.1e-12, 0.01e-12 ],
+ [ 'Ca_conc', 'Ca', 'spine/Ca', 'conc', 0.00008, 8 ]
+ ],
+ # Syn input basline 1 Hz, and 40Hz burst for 1 sec at t=20. Syn wt=10
+ stimList = [['head#', '10','glu', 'periodicsyn', '1 + 40*(t>10 && t<11)']],
+ plotList = [
+ ['soma', '1', '.', 'Vm', 'Membrane potential'],
+ ['#', '1', 'spine/Ca', 'conc', 'Ca in Spine'],
+ ['#', '1', 'dend/DEND/Ca', 'conc', 'Ca in Dend'],
+ ['head#', '1', 'psd/chan_p', 'n', 'Amount of Phospho-chan'],
+ ['head#', '1', 'spine/CaMKII', 'conc', 'Conc of CaMKII in spine'],
+ ['head#', '1', '.', 'Cm', 'Capacitance of spine head'],
+ ['head#', '1', '.', 'Rm', 'Membrane res of spine head'],
+ ['head#', '1', '.', 'Ra', 'Axial res of spine head'],
+ ['head#', '1', 'glu', 'Gbar', 'Conductance of gluR'],
+ ['head#', '1', 'NMDA', 'Gbar', 'Conductance of NMDAR'],
+ ]
+ )
+ moose.seed(123)
+ rdes.buildModel()
+ moose.reinit()
+ moose.start( 25 )
+ rdes.display()
+
+
+The key *adaptor* line is as follows:
+
+``[ 'psd/chan_p', 'n', 'spine', 'psdArea', 0.1e-12, 0.01e-12 ]``
+
+Here, we use the phosphorylated *chan_p* molecule in the PSD as a proxy for
+processes that control spine size. We operate on a special object called
+*spine* which manages many aspects of spines in the model (see below). Here
+we control the *psdArea*, which defines the cross-section area of the spine
+head and by extension of the PSD too. We keep a minimum spine area of 0.1 um^2,
+and a scaling factor of 0.01um^2 per phosphorylated molecule.
+
+The reaction system is identical to the one in *ex8.2*:
+
+::
+
+ Ca+CaM <===> Ca_CaM; Ca_CaM + CaMKII <===> Ca_CaM_CaMKII (all in
+ spine head, except that the Ca_CaM_CaMKII translocates to the PSD)
+
+ chan ------Ca_CaM_CaMKII-----> chan_p; chan_p ------> chan (all in PSD)
+
+Rather than list all the 10 plots, here are a few to show what is going on.
+
+First, just the spiking activity of the cell. Here the burst of activity is
+followed by a few seconds of enhanced synaptic weight, followed by subthreshold
+EPSPs:
+
+.. figure:: ../../../../images/ex8.3_Vm.png
+ :alt: Membrane potential and spiking.
+
+ Membrane potential and spiking.
+
+Then, we fast-forward to the amount of *chan_p* which is the molecule that
+controls spine size scaling:
+
+.. figure:: ../../../../images/ex8.3_chan_p.png
+ :alt: Molecule that controles spine size
+
+ Molecule that controles spine size
+
+This causes some obvious outcomes. One of them is to increase the synaptic
+conductance of the glutamate receptor. The system assumes that the conductance
+of all channels in the PSD scales linearly with the psdArea.
+
+.. figure:: ../../../../images/ex8.3_gluR.png
+ :alt: Conductance of glutamate receptor
+
+ Conductance of glutamate receptor
+
+Here is one of several non-intuitive outcomes. Because the spine volume has
+increased, the concentration of molecules in the spine is diluted out. So
+the concentration of active CaMKII actually falls when the spine gets bigger.
+In a more detailed model, this would be a race between the increase in spine
+size and the time taken for diffusion and further reactions to replenish
+CaMKII. In the current model we don't have a diffusive coupling of CaMKII to
+the dendrite, so this replenishment doesn't happen.
+
+.. figure:: ../../../../images/ex8.3_CaMKII_spine.png
+ :alt: Concentration of CaMKII in the spine
+
+ Concentration of CaMKII in the spine
+
+In the simulation we display several other electrical and chemical properties
+that change with spine size. The diffusion properties also change since the
+cross-section areas are altered. This is harder to visualize but has large
+effects on coupling to the dendrite,
+especially if the *shaftDiameter* is the parameter scaled by the signaling.
+
+
+Suggestions:
+
+ - The Spine class (instance: spine) manages several possible scaling
+ targets on the spine geometry: shaftLength, shaftDiameter,
+ headLength, headDiameter, psdArea, headVolume, totalLength. Try them
+ out. Think about mechanisms by which molecular concentrations might
+ affect each.
+ - When volume changes, we assume that the molecular numbers stay
+ fixed, so concentration changes. Except for buffered molecules, where
+ we assume concentration remains fixed. Use this to design a bistable
+ simply relying on molecules and spine geometry terms.
+ - Even more interesting, use it to design an oscillator. You could look
+ at Bhalla, BiophysJ 2011 for some ideas.
+
+
+
+Morphology: Load .swc morphology file and view it
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+*ex9.0_load_neuronal_morphology_file.py*
+
+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 sys
+ import moose
+ import rdesigneur as rd
+
+ if len( sys.argv ) > 1:
+ fname = sys.argv[1]
+ else:
+ fname = './cells/h10.CNG.swc'
+ rdes = rd.rdesigneur(
+ cellProto = [[fname, '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.001, 0.1, rotation = 0.02 )
+
+
+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/ex9.0_passive_cell_morpho.png
+ :alt: 3-D display for passive neuron
+
+ 3-D display for passive neuron
+
+Suggestions:
+
+ - The tutorial directory already has a number of pre-loaded files from
+ NeuroMorpho. Pass them in to ex9.0 on the command line:
+
+ `python ex9.0_load_neuronal_morphology_file.py <morpho.swc>`
+ - Grab other morphology files from NeuroMorpho.org, try them out.
+
+Build an active neuron model by putting channels into a morphology file
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+*ex9.1_chans_in_neuronal_morph.py*
+
+Here we load in a morphology file and distribute voltage-gated ion channels
+over the neuron. 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 ]``
+
+ or
+
+ ``[ channelFunction(), 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
+ternary expressions like *p < 50e-6 ? 500 : 100* or multiplying a channel
+density with a logical function or 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
+
+Suggestions:
+
+ - Try another morphology file.
+ - Try different channel distributions by editing the chanDistrib lines.
+ - There are numerous predefined channels available within Rdesigneur.
+ These can be defined using the following chanProto options:
+
+ ::
+
+ ['make_HH_Na()', 'HH_Na']
+ ['make_HH_K_DR()', 'HH_K']
+ ['make_Na()', 'Na']
+ ['make_K_DR()', 'K_DR']
+ ['make_K_A()', 'K_A']
+ ['make_K_AHP()', 'K_AHP']
+ ['make_K_C()', 'K_C']
+ ['make_Ca()', 'Ca']
+ ['make_Ca_conc()', 'Ca_conc']
+ ['make_glu()', 'glu']
+ ['make_GABA()', 'GABA']
+
+ Then the chanDistrib can refer to these channels instead.
+ - Deliver stimuli on the dendrites rather than the soma.
+
+
+Build a spiny neuron from a morphology file and put active channels in it.
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+*ex9.2_spines_in_neuronal_morpho.py*
+
+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.
+
+.. figure:: ../../../../images/rdes9_spiny_active.png
+ :alt: 3-D display for spiny active neuron
+
+ 3-D display for spiny active neuron
+
+Suggestions:
+
+ - Try different spine settings. Warning: if you put in too many spines
+ it will take much longer to load and run!
+ - Try different spine geometry layouts.
+ - See if you can deliver the current injection to the spine. Hint: the
+ name of the spine compartments is 'head#' where # is the index of the
+ spine.
+
+