diff --git a/docs/source/Extensions/hidden_code_block.py b/docs/source/Extensions/hidden_code_block.py
new file mode 100644
index 0000000000000000000000000000000000000000..3782d47d582bdea3f50c9f0bd18af38c74dc667b
--- /dev/null
+++ b/docs/source/Extensions/hidden_code_block.py
@@ -0,0 +1,129 @@
+"""Simple, inelegant Sphinx extension which adds a directive for a
+highlighted code-block that may be toggled hidden and shown in HTML.
+This is possibly useful for teaching courses.
+
+The directive, like the standard code-block directive, takes
+a language argument and an optional linenos parameter. The
+hidden-code-block adds starthidden and label as optional
+parameters.
+
+Examples:
+
+.. hidden-code-block:: python
+ :starthidden: False
+
+ a = 10
+ b = a + 5
+
+.. hidden-code-block:: python
+ :label: --- SHOW/HIDE ---
+
+ x = 10
+ y = x + 5
+
+Thanks to http://www.javascriptkit.com/javatutors/dom3.shtml for
+inspiration on the javascript.
+
+Thanks to Milad 'animal' Fatenejad for suggesting this extension
+in the first place.
+
+Written by Anthony 'el Scopz' Scopatz, January 2012.
+
+Released under the WTFPL (http://sam.zoy.org/wtfpl/).
+"""
+
+from docutils import nodes
+from docutils.parsers.rst import directives
+from sphinx.directives.code import CodeBlock
+from sphinx.util.compat import make_admonition
+
+HCB_COUNTER = 0
+
+js_showhide = """\
+<script type="text/javascript">
+ function showhide(element){
+ if (!document.getElementById)
+ return
+
+ if (element.style.display == "block")
+ element.style.display = "none"
+ else
+ element.style.display = "block"
+ };
+</script>
+"""
+
+def nice_bool(arg):
+ tvalues = ('true', 't', 'yes', 'y')
+ fvalues = ('false', 'f', 'no', 'n')
+ arg = directives.choice(arg, tvalues + fvalues)
+ return arg in tvalues
+
+
+class hidden_code_block(nodes.General, nodes.FixedTextElement):
+ pass
+
+
+class HiddenCodeBlock(CodeBlock):
+ """Hidden code block is Hidden"""
+
+ option_spec = dict(starthidden=nice_bool,
+ label=str,
+ **CodeBlock.option_spec)
+
+ def run(self):
+ # Body of the method is more or less copied from CodeBlock
+ code = u'\n'.join(self.content)
+ hcb = hidden_code_block(code, code)
+ hcb['language'] = self.arguments[0]
+ hcb['linenos'] = 'linenos' in self.options
+ hcb['starthidden'] = self.options.get('starthidden', True)
+ hcb['label'] = self.options.get('label', '+ show/hide code')
+ hcb.line = self.lineno
+ return [hcb]
+
+
+def visit_hcb_html(self, node):
+ """Visit hidden code block"""
+ global HCB_COUNTER
+ HCB_COUNTER += 1
+
+ # We want to use the original highlighter so that we don't
+ # have to reimplement it. However it raises a SkipNode
+ # error at the end of the function call. Thus we intercept
+ # it and raise it again later.
+ try:
+ self.visit_literal_block(node)
+ except nodes.SkipNode:
+ pass
+
+ # The last element of the body should be the literal code
+ # block that was just made.
+ code_block = self.body[-1]
+
+ fill_header = {'divname': 'hiddencodeblock{0}'.format(HCB_COUNTER),
+ 'startdisplay': 'none' if node['starthidden'] else 'block',
+ 'label': node.get('label'),
+ }
+
+ divheader = ("""<a href="javascript:showhide(document.getElementById('{divname}'))">"""
+ """{label}</a><br />"""
+ '''<div id="{divname}" style="display: {startdisplay}">'''
+ ).format(**fill_header)
+
+ code_block = js_showhide + divheader + code_block + "</div>"
+
+ # reassign and exit
+ self.body[-1] = code_block
+ raise nodes.SkipNode
+
+
+def depart_hcb_html(self, node):
+ """Depart hidden code block"""
+ # Stub because of SkipNode in visit
+
+
+def setup(app):
+ app.add_directive('hidden-code-block', HiddenCodeBlock)
+ app.add_node(hidden_code_block, html=(visit_hcb_html, depart_hcb_html))
+
diff --git a/docs/source/conf.py b/docs/source/conf.py
index f3d787c8d8f668d52c5df341e587151d8e9eeb0e..4b6e833ef26bf1789124f6d221c566bc7235c8b4 100644
--- a/docs/source/conf.py
+++ b/docs/source/conf.py
@@ -11,6 +11,11 @@
# All configuration values have a default; values that are commented out
# serve to show the default.
+#Workaround to fix bug where extensions werent added
+from docutils.parsers.rst.directives.admonitions import BaseAdmonition
+from sphinx.util import compat
+compat.make_admonition = BaseAdmonition
+
import subprocess
import os
import sys
@@ -21,6 +26,7 @@ import mock
# 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.insert(0, os.path.abspath('./Extensions'))
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'))
@@ -41,6 +47,7 @@ extensions = ['sphinx.ext.autodoc',
'sphinx.ext.autosummary',
'sphinx.ext.todo',
'sphinx.ext.viewcode',
+ 'hidden_code_block'
]
todo_include_todos = True
@@ -110,6 +117,7 @@ todo_include_todos = True
# A list of ignored prefixes for module index sorting.
#modindex_common_prefix = []
+
# -- Options for HTML output ---------------------------------------------------
# The theme to use for HTML and HTML Help pages. See the documentation for
@@ -286,4 +294,5 @@ import subprocess, os
read_the_docs_build = os.environ.get('READTHEDOCS', None) == 'True'
if read_the_docs_build:
- subprocess.call('cd ../../doxygen; echo HELLO......................; doxygen Doxyfile', shell=True)
+ subprocess.call('cd doxygen; echo HELLO......................; doxygen Doxyfile', shell=True)
+
diff --git a/docs/source/docs/source/modules.rst b/docs/source/docs/source/modules.rst
deleted file mode 100644
index bf86f144c054261a95bfc75c1b58fc47b941f688..0000000000000000000000000000000000000000
--- a/docs/source/docs/source/modules.rst
+++ /dev/null
@@ -1,2 +0,0 @@
-projectdir
-==========
diff --git a/docs/source/doxygen/doxy.rst b/docs/source/doxygen/doxy.rst
index 4d62729abe5307f9c78297eb319bc899817d9c05..8d271e7f433c9a9013ac59d76972e0abc3bf83e9 100644
--- a/docs/source/doxygen/doxy.rst
+++ b/docs/source/doxygen/doxy.rst
@@ -2,3 +2,6 @@ Doxygen
-------
Here you can find all the references necessary for MOOSE.
+
+`Click here <../../../source/doxygen/cpp/html/index.html>`_
+
diff --git a/docs/source/images/FB.png b/docs/source/images/FB.png
new file mode 100644
index 0000000000000000000000000000000000000000..7e16810c5850f724ec01274487b71d95aedc2d88
Binary files /dev/null and b/docs/source/images/FB.png differ
diff --git a/docs/source/images/Kholodenko_tut.png b/docs/source/images/Kholodenko_tut.png
new file mode 100644
index 0000000000000000000000000000000000000000..590fada9e83bb999380aa4412e98fc90e24ed9d6
Binary files /dev/null and b/docs/source/images/Kholodenko_tut.png differ
diff --git a/docs/source/images/doseR.png b/docs/source/images/doseR.png
new file mode 100644
index 0000000000000000000000000000000000000000..480fb5ef1784b1e405fc83ef3652bc0015829772
Binary files /dev/null and b/docs/source/images/doseR.png differ
diff --git a/docs/source/images/findS.png b/docs/source/images/findS.png
new file mode 100644
index 0000000000000000000000000000000000000000..912216a323f4629272a8580109783ffe0e564c28
Binary files /dev/null and b/docs/source/images/findS.png differ
diff --git a/docs/source/images/mapkFB.png b/docs/source/images/mapkFB.png
new file mode 100644
index 0000000000000000000000000000000000000000..62a1d730e828df26443ac4ed16ae2cf77d14c359
Binary files /dev/null and b/docs/source/images/mapkFB.png differ
diff --git a/docs/source/images/mapkFB2.png b/docs/source/images/mapkFB2.png
new file mode 100644
index 0000000000000000000000000000000000000000..c8a3d9f51e867449497cb3635c5f6166ca16e86a
Binary files /dev/null and b/docs/source/images/mapkFB2.png differ
diff --git a/docs/source/images/propBis.gif b/docs/source/images/propBis.gif
new file mode 100644
index 0000000000000000000000000000000000000000..b500a387f0edba3d722c579da837c24d935fbb70
Binary files /dev/null and b/docs/source/images/propBis.gif differ
diff --git a/docs/source/images/propBis.png b/docs/source/images/propBis.png
new file mode 100644
index 0000000000000000000000000000000000000000..27ac875c33b02dc27db9b8a6ac505cbf98c936c5
Binary files /dev/null and b/docs/source/images/propBis.png differ
diff --git a/docs/source/images/relax.png b/docs/source/images/relax.png
new file mode 100644
index 0000000000000000000000000000000000000000..495ce2517af016bf25d7bf951cd9f4f008ff996f
Binary files /dev/null and b/docs/source/images/relax.png differ
diff --git a/docs/source/images/relaxOsc_tut.png b/docs/source/images/relaxOsc_tut.png
new file mode 100644
index 0000000000000000000000000000000000000000..154cd222592e52ff854f2160cac1a932382d745d
Binary files /dev/null and b/docs/source/images/relaxOsc_tut.png differ
diff --git a/docs/source/images/repressillatorOsc.png b/docs/source/images/repressillatorOsc.png
new file mode 100644
index 0000000000000000000000000000000000000000..6e9824fc920a31d5ed16c5b71adb3a682c257271
Binary files /dev/null and b/docs/source/images/repressillatorOsc.png differ
diff --git a/docs/source/images/repris.png b/docs/source/images/repris.png
new file mode 100644
index 0000000000000000000000000000000000000000..a1e080f63844596dfb0b7cef05d5a2520f2fd65e
Binary files /dev/null and b/docs/source/images/repris.png differ
diff --git a/docs/source/images/sV1.png b/docs/source/images/sV1.png
new file mode 100644
index 0000000000000000000000000000000000000000..4c66962f29987ea44658887b0edc40b8bba36765
Binary files /dev/null and b/docs/source/images/sV1.png differ
diff --git a/docs/source/images/sV2.png b/docs/source/images/sV2.png
new file mode 100644
index 0000000000000000000000000000000000000000..2ab09075496188a038b20d86dccf3f369cfaeac0
Binary files /dev/null and b/docs/source/images/sV2.png differ
diff --git a/docs/source/images/sV3.png b/docs/source/images/sV3.png
new file mode 100644
index 0000000000000000000000000000000000000000..4b0218aba4f6c66f0746e37f75bc1e92cc2afcb2
Binary files /dev/null and b/docs/source/images/sV3.png differ
diff --git a/docs/source/images/sV4.png b/docs/source/images/sV4.png
new file mode 100644
index 0000000000000000000000000000000000000000..eec6f0e22603b4c139f53a5b4cf374f373fadc0f
Binary files /dev/null and b/docs/source/images/sV4.png differ
diff --git a/docs/source/images/sV5.png b/docs/source/images/sV5.png
new file mode 100644
index 0000000000000000000000000000000000000000..dc699319207ab16065f7d15b738251677f6a3e03
Binary files /dev/null and b/docs/source/images/sV5.png differ
diff --git a/docs/source/images/sV6.png b/docs/source/images/sV6.png
new file mode 100644
index 0000000000000000000000000000000000000000..30f9c4ac37b10ef872837741984baf680de21b6a
Binary files /dev/null and b/docs/source/images/sV6.png differ
diff --git a/docs/source/images/sV7.png b/docs/source/images/sV7.png
new file mode 100644
index 0000000000000000000000000000000000000000..d1771eb68582e7f42642f0edca1f0716eed05028
Binary files /dev/null and b/docs/source/images/sV7.png differ
diff --git a/docs/source/images/simpleB.png b/docs/source/images/simpleB.png
new file mode 100644
index 0000000000000000000000000000000000000000..459b0b4ab912c7d2bd4f348eea4268e6557ec8a3
Binary files /dev/null and b/docs/source/images/simpleB.png differ
diff --git a/docs/source/images/strongB.png b/docs/source/images/strongB.png
new file mode 100644
index 0000000000000000000000000000000000000000..4116cdc9c3047f0ca0f752d0e6caf76821ebe8e9
Binary files /dev/null and b/docs/source/images/strongB.png differ
diff --git a/docs/source/images/strongBis.png b/docs/source/images/strongBis.png
new file mode 100644
index 0000000000000000000000000000000000000000..5968aa2c5506fd644531ba9e8dd9b975acbdaa06
Binary files /dev/null and b/docs/source/images/strongBis.png differ
diff --git a/docs/source/images/turing.png b/docs/source/images/turing.png
new file mode 100644
index 0000000000000000000000000000000000000000..c0e19eef74ba00da3c643b36daff60bd6d2eca0e
Binary files /dev/null and b/docs/source/images/turing.png differ
diff --git a/docs/source/images/turingPatternTut.png b/docs/source/images/turingPatternTut.png
new file mode 100644
index 0000000000000000000000000000000000000000..e45d84d15424f447af4259208a6a7aab50f5e9cb
Binary files /dev/null and b/docs/source/images/turingPatternTut.png differ
diff --git a/docs/source/index.rst b/docs/source/index.rst
index aa463b195ce811c66214430238e874299d41f6ea..62dfb486b3f34d5a187cbaaba46a51c3d1642ca0 100644
--- a/docs/source/index.rst
+++ b/docs/source/index.rst
@@ -28,13 +28,13 @@ MOOSE is a simulation environment, not just a numerical engine: It provides data
/install/index_install
/user/py/quickstart/index_qs
/user/py/cookbook/index_ckbk
+ /user/py/tutorials/index_tut
/user/py/graphics/index_graphics
/user/py/references/index_ref
/doxygen/doxy
.. toctree::
:maxdepth: 1
-
:titlesonly:
.. toctree::
diff --git a/docs/source/install/install.rst b/docs/source/install/install.rst
index 41a73c401862cb08207d44abe4d22c3b55c27170..6d83001c30d57a9cbe04f831d40f49bdcae2e714 100644
--- a/docs/source/install/install.rst
+++ b/docs/source/install/install.rst
@@ -25,11 +25,15 @@ to pick your distribution and follow instructions.
.. todo:: Packages for gentoo, Arch Linux
+Mac OSX
+^^^^^^^^
-MacOSX support is not complete yet. Python-scripting interface can be installed on MaxOSX using ``homebrew``
+
+MacOSX support for moose-gui is not complete yet. However, the python-scripting interface can be installed on OSX using ``homebrew``
::
- $ brew install homebrew/science/moose
+ $ brew tap BhallaLab/moose
+ $ brew install moose
Building MOOSE
@@ -37,11 +41,12 @@ Building MOOSE
In case your distribution is not listed on `our repository page
<https://software.opensuse.org/download.html?project=home:moose&package=moose>`_
-, or if you want to build the lastest development code, read on.
+, or if you want to build the latest development code, read on.
First, you need to get the source code. You can use ``git`` (clone the
repository) or download snapshot of github repo by clicking on `this link
-<https://github.com/BhallaLab/moose/archive/master.zip>`_.::
+<https://github.com/BhallaLab/moose/archive/master.zip>`__.
+::
$ git clone https://github.com/BhallaLab/moose
@@ -52,8 +57,8 @@ Or,
$ wget https://github.com/BhallaLab/moose/archive/master.zip
$ unzip master.zip
-If you don't want lasest snapshot of ``MOOSE``, you can download other released
-versions from `here <`https://github.com/BhallaLab/moose/releases>`_.
+If you don't want latest snapshot of ``MOOSE``, you can download other released
+versions from `here <https://github.com/BhallaLab/moose/releases>`__.
Install dependencies
^^^^^^^^^^^^^^^^^^^^
@@ -111,7 +116,7 @@ build moose
.. code-block:: bash
- $ cd /to/moose/source/code
+ $ cd /to/moose_directory/moose-core/
$ mkdir _build
$ cd _build
$ cmake ..
@@ -153,12 +158,15 @@ If you have installed the pre-built package, then you already have the GUI.
You can launch it by runnung `moosegui` command in terminal.
You can get the source of ``moose-gui`` from `here
-<https://github.com/BhallaLab/moose-gui>`_. You can download it either by
-clicking on `this link <https://github.com/BhallaLab/moose-gui/archive/master.zip>`_
+<https://github.com/BhallaLab/moose-gui>`__. You can download it either by
+clicking on `this link <https://github.com/BhallaLab/moose-gui/archive/master.zip>`__
or by using ``git`` ::
$ git clone https://github.com/BhallaLab/moose-gui
+
+Alternatively the moose-gui folder exists within the moose folder downloaded and built earlier in the installation process. It can be found under ``location_of_moose_folder/moose/moose-gui/``.
+
Below are packages which you may want to install to use MOOSE Graphical User Interface.
- Required:
@@ -196,7 +204,7 @@ Building moogli
---------------
``moogli`` is subproject of ``MOOSE`` for visualizing models. More details can
-be found `here <http://moose.ncbs.res.in/moogli>`_.
+be found `here <http://moose.ncbs.res.in/moogli>`__.
`Moogli` is part of `moose` package. Building moogli can be tricky because of
multiple depednecies it has.
@@ -205,7 +213,7 @@ multiple depednecies it has.
- OSG (3.2.x) For 3D rendering and simulation of neuronal models
- Qt4 (4.8.x) For C++ GUI of Moogli
-To get the latest source code of ``moogli``, click on `this link <https://github.com/BhallaLab/moogli/archive/master.zip>`_.
+To get the latest source code of ``moogli``, click on `this link <https://github.com/BhallaLab/moogli/archive/master.zip>`__.
Moogli depends on ``OpenSceneGraph`` (version 3.2.0 or higher) which may not
be easily available for your operating system.
diff --git a/docs/source/user/py/cookbook/Rd.rst b/docs/source/user/py/cookbook/Rd.rst
new file mode 100644
index 0000000000000000000000000000000000000000..625b8bc819cf5b695210c854050e7630d60c6ff0
--- /dev/null
+++ b/docs/source/user/py/cookbook/Rd.rst
@@ -0,0 +1,1186 @@
+**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/chem_GUI.rst b/docs/source/user/py/cookbook/chem_GUI.rst
index 390b56ee104ab200ada187c945987ffc3959d8a5..543f48c6b34abb7ecc4da7d71d4f7f822221f06b 100644
--- a/docs/source/user/py/cookbook/chem_GUI.rst
+++ b/docs/source/user/py/cookbook/chem_GUI.rst
@@ -355,6 +355,7 @@ variables can be input from pool object.
- **New**: **File -> New -> Model name**. This opens a empty widget for model building
- **Saving models**: **File -> Save Model -> select from dialog**.
- **Changing numerical methods**: **Preference->Chemical tab** item from Simulation Control. Currently supports:
+
1. Runge Kutta: This is the Runge-Kutta-Fehlberg implementation from the GNU Scientific Library (GSL). It is a fifth order variable timestep explicit method. Works well for most reaction systems except if they have very stiff reactions.
2. Gillespie: Optimized Gillespie stochastic systems algorithm, custom implementation. This uses variable timesteps internally. Note that it slows down with increasing numbers of molecules in each pool. It also slows down, but not so badly, if the number of reactions goes up.
3. Exponential Euler:This methods computes the solution of partial and ordinary differential equations.
diff --git a/docs/source/user/py/cookbook/elec_sim_eg.rst b/docs/source/user/py/cookbook/elec_sim_eg.rst
index 0506b94c093caaf775b3a3fc0e0399c931718772..224cf913e7c5429701ef78692b5a0fd31ab08862 100644
--- a/docs/source/user/py/cookbook/elec_sim_eg.rst
+++ b/docs/source/user/py/cookbook/elec_sim_eg.rst
@@ -53,8 +53,9 @@ Insert Spine heads
.. automodule:: insertSpinesWithoutRdesigneur
:members:
-Multi compartmental Neuron model
---------------------------------
+.. Multi compartmental Neuron model
+.. --------------------------------
-.. automodule:: testHsolve
+.. automodule: testHsolve
:members:
+ (testHsolve no longer exists in any of the appended paths)
diff --git a/docs/source/user/py/cookbook/multi_rdes.rst b/docs/source/user/py/cookbook/multi_rdes.rst
index dce3f35878901693b3f6a530a84c6b310e2ecd44..baf5719bd85b42dec18aac48e1e1a8d7fb06f091 100644
--- a/docs/source/user/py/cookbook/multi_rdes.rst
+++ b/docs/source/user/py/cookbook/multi_rdes.rst
@@ -1,6 +1,6 @@
-**********
-RDesigneur
-**********
+************************
+More Rdesigneur Examples
+************************
Building Chemical-Electrical Signalling Models
----------------------------------------------
diff --git a/docs/source/user/py/cookbook/multi_sim_eg.rst b/docs/source/user/py/cookbook/multi_sim_eg.rst
index 564bdb7fab6c5ae715e87a6e84405dc0fc36ca57..529ee4c3911a9729ba3db76dfe404330b6823cec 100644
--- a/docs/source/user/py/cookbook/multi_sim_eg.rst
+++ b/docs/source/user/py/cookbook/multi_sim_eg.rst
@@ -20,5 +20,6 @@ Modeling chemical reactions in neurons
.. automodule:: gssaRDspiny
:members:
-.. automodule:: rxdSpineSize
+.. automodule: rxdSpineSize
:members:
+ (rxdSpineSize does not exist in any appended folders -> gives error while building html)
diff --git a/docs/source/user/py/cookbook/multiscale.rst b/docs/source/user/py/cookbook/multiscale.rst
index a8a83165c1675bd3d697107ed76aa9be96a40915..d751d3e4bb807f031020ba03fd4d4dd99854e20f 100644
--- a/docs/source/user/py/cookbook/multiscale.rst
+++ b/docs/source/user/py/cookbook/multiscale.rst
@@ -6,4 +6,8 @@ MultiScale Modeling
:maxdepth: 2
multi_sim_eg
+ Rd
multi_rdes
+
+
+
diff --git a/docs/source/user/py/quickstart/classes_demos.rst b/docs/source/user/py/quickstart/classes_demos.rst
index d8e3bb3493861f7edbb7991735d0158fe0e53503..290c2eb70994a5b2a4857ff00500dfad83782f19 100644
--- a/docs/source/user/py/quickstart/classes_demos.rst
+++ b/docs/source/user/py/quickstart/classes_demos.rst
@@ -70,6 +70,7 @@ Function
.. automodule:: func
:members:
+ :noindex:
SymCompartment
^^^^^^^^^^^^^^
diff --git a/docs/source/user/py/quickstart/demos.rst b/docs/source/user/py/quickstart/demos.rst
index 390d05fb801b32c6a4742e49e65d2cabb9bcf021..7fe408f1f09e5ae1b75ebf740262bafe92fc3678 100644
--- a/docs/source/user/py/quickstart/demos.rst
+++ b/docs/source/user/py/quickstart/demos.rst
@@ -21,31 +21,13 @@ Start, Stop, and setting clocks
.. automodule:: stimtable
:members:
-The Hodgkin-Huxley demo
------------------------
-
-This is a self-contained graphical demo implemented by Subhasis Ray,
-closely based on the 'Squid' demo by Mark Nelson which ran in GENESIS.
-
-.. figure:: ../../../images/squid_demo.png
- :alt: Hodgkin-Huxley's squid giant axon experiment
-
-Simulation of Hodgkin-Huxley's experiment on squid giant axon
-showing action potentials generated by a step current injection.
-
-
-The demo has built-in documentation and may be run from the
-``moose-examples/squid`` subdirectory of MOOSE.
-
-.. _quickstart-python_from_moose:
-
Run Python from MOOSE
---------------------
.. automodule:: pyrun
:members:
-.. automodule:: pyrun1
- :members:
+.. automodule: pyrun1
+.. :members:
.. _quickstart-classes:
diff --git a/docs/source/user/py/quickstart/index_qs.rst b/docs/source/user/py/quickstart/index_qs.rst
new file mode 100644
index 0000000000000000000000000000000000000000..c75b441d2b44c1bb44d586ff115a765d74b323e7
--- /dev/null
+++ b/docs/source/user/py/quickstart/index_qs.rst
@@ -0,0 +1,15 @@
+.. MOOSE documentation master file, created by
+ sphinx-quickstart on Tue Jul 1 19:05:47 2014.
+ You can adapt this file completely to your liking, but it should at least
+ contain the root `toctree` directive.
+
+Quick Start
+===========
+
+.. toctree::
+ :maxdepth: 2
+
+ qs_GUI
+ moose_quickstart
+ demos
+ classes_demos
diff --git a/docs/source/user/py/quickstart/moose_quickstart.rst b/docs/source/user/py/quickstart/moose_quickstart.rst
index a446b679fe69ccdd4bfac9add3761a4d35140adf..55b30186fe2adcefe4cace233b7cc287a0cef017 100644
--- a/docs/source/user/py/quickstart/moose_quickstart.rst
+++ b/docs/source/user/py/quickstart/moose_quickstart.rst
@@ -2,13 +2,9 @@
Getting started with python scripting for MOOSE
***********************************************
-.. contents::
- :local:
- :depth: 1
-
.. figure:: ../../../images/Gallery_Moose_Multiscale.png
:alt: **multiple scales in moose**
- :scale: 50%
+ :scale: 100%
*Multiple scales can be modelled and simulated in MOOSE*
@@ -30,18 +26,58 @@ MOOSE is a simulation environment, not just a numerical engine: It provides data
Contents:
-.. toctree::
- :maxdepth: 2
- :numbered:
+.. contents::
+ :local:
+ :depth: 1
- moose_quickstart
+Coding basics and how to use this document
+==========================================
-Indices and tables
-==================
+This page acts as the first stepping stone for learning how moose works. The tutorials here are intended to be interactive, and are presented as python commands. Python commands are identifiable by the ``>>>`` before the command as opposed to ``$`` which identifies a command-line command.
-* :ref:`genindex`
-* :ref:`modindex`
-* :ref:`search`
+ >>> this_is_a_python_command
+
+You are encouraged to run a python shell while reading through this documentation and trying out each command for yourself. Python shells are environments within your terminal wherein everything you type is interpreted as a python command. They can be accessed by typing
+::
+
+ $ python
+
+in your command-line terminal.
+
+While individually typing lines of code in a python terminal is useful for practicing using moose and coding in general, note that once you close the python environment all the code you typed is gone and the moose models created are also lost. In order to 'save' models that you create, you would have to type your code in a text file with a ``.py`` extension. The easiest way to do this is to create a text file in command line, open it with a text editor (for example, gedit), and simply type your code in (make sure you indent correctly).
+::
+
+ $ touch code.py
+ $ gedit code.py
+
+Once you have written your code in the file, you can run it through your python environment.
+::
+
+ $ python code.py
+
+Note that apart from this section of the quickstart, most of the moose documentation is in the form of ``snippets``. These are basically ``.py`` files with code that demonstrates a certain functionality in moose. If you see a dialogue box like this one:
+
+.. automodule:: func
+ :members: main
+
+You can view the code by clicking the green source button on the left side of the box. Alternatively, the source code for all of the examples in the documentation can be found in ``moose/moose-examples/snippets``. Once you run each file in python, it is encouraged that you look through the code to understand how it works.
+
+In the quickstart, most of the snippets demonstrate the functionality of specific classes. However, snippets in later sections such as the cookbook show how to do specific things in moose such as creating networks, chemical models, and synaptic channels.
+
+
+Importing moose and accessing documentation
+===========================================
+
+
+In a python script you import modules to access the functionalities they provide. In order to use moose, you need to import it within a python environment or at the beginning of your python script.
+
+ >>> import moose
+
+This make the ``moose`` module available for use in Python. You can use Python's built-in ``help`` function to read the top-level documentation for the moose module.
+
+ >>> help(moose)
+
+This will give you an overview of the module. Press ``q`` to exit the pager and get back to the interpreter. You can also access the documentation for individual classes and functions this way.
>>> help(moose.connect)
@@ -69,6 +105,14 @@ Note that you need to put the class-name followed by dot followed by
field-name within quotes. Otherwise, ``moose.doc`` will receive the
field value as parameter and get confused.
+Alternatively, if you want to see a full list of classes, functions and their fields, you can browse through the following pages. This is especially helpful when going through snippets.
+
+* :ref:`genindex`
+* :ref:`modindex`
+
+..
+ :ref:`search`
+
.. _quickstart-creating:
Creating objects and traversing the object hierarchy
@@ -593,9 +637,9 @@ loop through the elements in an ``vec`` like a Python list ::
shows ::
- /model/comp[0] <type 'moose.melement'>
- /model/comp[1] <type 'moose.melement'>
- /model/comp[2] <type 'moose.melement'>
+ /model[0]/comp[0] <type 'moose.Compartment'>
+ /model[0]/comp[1] <type 'moose.Compartment'>
+ /model[0]/comp[2] <type 'moose.Compartment'>
Thus elements are instances of class ``melement``. All elements in an
``vec`` share the ``id_`` of the ``vec`` which can retrieved by
diff --git a/docs/source/user/py/tutorials/ChemicalBistables.rst b/docs/source/user/py/tutorials/ChemicalBistables.rst
new file mode 100644
index 0000000000000000000000000000000000000000..6ac360cb3b67e38a77ed7c33bd017e23615951c2
--- /dev/null
+++ b/docs/source/user/py/tutorials/ChemicalBistables.rst
@@ -0,0 +1,1393 @@
+******************
+Chemical Bistables
+******************
+
+A `bistable system <https://en.wikipedia.org/wiki/Bistability>`_ is a dynamic system that has two stable equilibrium states. The following examples can be used to teach and demonstrate different aspects of bistable systems or to learn how to model them using moose. Each example contains a short description, the model's code, and the output with default settings.
+
+Each example can be found as a python file within the main moose folder under
+::
+
+ (...)/moose/moose-examples/tutorials/ChemicalBistables
+
+In order to run the example, run the script
+::
+
+ python filename.py
+
+in command line, where ``filename.py`` is the name of the python file you would like to run. The filenames of each example are written in **bold** at the beginning of their respective sections, and the files themselves can be found in the aformentioned directory.
+
+In chemical bistable models that use solvers, there are optional arguments that allow you to specify which solver you would like to use.
+::
+
+ python filename.py [gsl | gssa | ee]
+
+Where ``gsl`` is Gnu Scientific Library's deterministic solver, ``gssa`` stands for Gillespie stochastic simulation algorithm, and ``ee`` is the exponential euler algorithm.
+
+All the following examples can be run with either of the three solvers, which in some cases produces a different outcome. However, simply running the file without the optional argument will by default use the ``gsl`` solver. These ``gsl`` outputs are the ones shown below.
+
+Simple Bistables
+================
+
+Filename: **simpleBis.py**
+
+
+This example shows the key property of a chemical bistable system: it
+has two stable states. Here we start out with the system settling rather
+quickly to the first stable state, where molecule A is high (blue) and
+the complementary molecule B (green) is low. At t = 100s, we deliver a
+perturbation, which is to move 90% of the A molecules into B. This
+triggers a state flip, which settles into a distinct stable state where
+there is more of B than of A. At t = 200s we reverse the flip by moving
+99% of B molecules back to A.
+
+If we run the simulation with the gssa option python simpleBis.py gssa
+
+we see exactly the same sequence of events, except now the switch is
+noisy. The calculations are now run with the Gillespie Stochastic
+Systems Algorithm (gssa) which incorporates probabilistic reaction
+events. The switch still switches but one can see that it might flip
+spontaneously once in a while.
+
+Things to do:
+
+1. Open a copy of the script file in an editor, and around
+line 124 and 129 you will see how the state flip is implemented while
+maintaining mass conservation. What happens if you flip over fewer
+molecules? What is the threshold for a successful flip? Why are these
+thresholds different for the different states?
+
+2. Try different volumes in line 31, and rerun using the gssa. Will you
+ see more or less noise if you increase the volume to 1e-20 m^3?
+
+**Code:**
+
+.. hidden-code-block:: python
+ :linenos:
+ :label: Show/Hide code
+
+ #########################################################################
+ ## This program is part of 'MOOSE', the
+ ## Messaging Object Oriented Simulation Environment.
+ ## Copyright (C) 2013 Upinder S. Bhalla. and NCBS
+ ## It is made available under the terms of the
+ ## GNU Lesser General Public License version 2.1
+ ## See the file COPYING.LIB for the full notice.
+ #########################################################################
+
+ # This example illustrates how to set up a kinetic solver and kinetic model
+ # using the scripting interface. Normally this would be done using the
+ # Shell::doLoadModel command, and normally would be coordinated by the
+ # SimManager as the base of the entire model.
+ # This example creates a bistable model having two enzymes and a reaction.
+ # One of the enzymes is autocatalytic.
+ # The model is set up to run using deterministic integration.
+ # If you pass in the argument 'gssa' it will run with the stochastic
+ # solver instead
+ # You may also find it interesting to change the volume.
+
+ import math
+ import pylab
+ import numpy
+ import moose
+ import sys
+
+ def makeModel():
+ # create container for model
+ model = moose.Neutral( 'model' )
+ compartment = moose.CubeMesh( '/model/compartment' )
+ compartment.volume = 1e-21 # m^3
+ # the mesh is created automatically by the compartment
+ mesh = moose.element( '/model/compartment/mesh' )
+
+ # create molecules and reactions
+ a = moose.Pool( '/model/compartment/a' )
+ b = moose.Pool( '/model/compartment/b' )
+ c = moose.Pool( '/model/compartment/c' )
+ enz1 = moose.Enz( '/model/compartment/b/enz1' )
+ enz2 = moose.Enz( '/model/compartment/c/enz2' )
+ cplx1 = moose.Pool( '/model/compartment/b/enz1/cplx' )
+ cplx2 = moose.Pool( '/model/compartment/c/enz2/cplx' )
+ reac = moose.Reac( '/model/compartment/reac' )
+
+ # connect them up for reactions
+ moose.connect( enz1, 'sub', a, 'reac' )
+ moose.connect( enz1, 'prd', b, 'reac' )
+ moose.connect( enz1, 'enz', b, 'reac' )
+ moose.connect( enz1, 'cplx', cplx1, 'reac' )
+
+ moose.connect( enz2, 'sub', b, 'reac' )
+ moose.connect( enz2, 'prd', a, 'reac' )
+ moose.connect( enz2, 'enz', c, 'reac' )
+ moose.connect( enz2, 'cplx', cplx2, 'reac' )
+
+ moose.connect( reac, 'sub', a, 'reac' )
+ moose.connect( reac, 'prd', b, 'reac' )
+
+ # connect them up to the compartment for volumes
+ #for x in ( a, b, c, cplx1, cplx2 ):
+ # moose.connect( x, 'mesh', mesh, 'mesh' )
+
+ # Assign parameters
+ a.concInit = 1
+ b.concInit = 0
+ c.concInit = 0.01
+ enz1.kcat = 0.4
+ enz1.Km = 4
+ enz2.kcat = 0.6
+ enz2.Km = 0.01
+ reac.Kf = 0.001
+ reac.Kb = 0.01
+
+ # Create the output tables
+ graphs = moose.Neutral( '/model/graphs' )
+ outputA = moose.Table ( '/model/graphs/concA' )
+ outputB = moose.Table ( '/model/graphs/concB' )
+
+ # connect up the tables
+ moose.connect( outputA, 'requestOut', a, 'getConc' );
+ moose.connect( outputB, 'requestOut', b, 'getConc' );
+
+ # Schedule the whole lot
+ moose.setClock( 4, 0.01 ) # for the computational objects
+ moose.setClock( 8, 1.0 ) # for the plots
+ # The wildcard uses # for single level, and ## for recursive.
+ moose.useClock( 4, '/model/compartment/##', 'process' )
+ moose.useClock( 8, '/model/graphs/#', 'process' )
+
+ def displayPlots():
+ for x in moose.wildcardFind( '/model/graphs/conc#' ):
+ t = numpy.arange( 0, x.vector.size, 1 ) #sec
+ pylab.plot( t, x.vector, label=x.name )
+ pylab.legend()
+ pylab.show()
+
+ def main():
+ solver = "gsl"
+ makeModel()
+ if ( len ( sys.argv ) == 2 ):
+ solver = sys.argv[1]
+ stoich = moose.Stoich( '/model/compartment/stoich' )
+ stoich.compartment = moose.element( '/model/compartment' )
+ if ( solver == 'gssa' ):
+ gsolve = moose.Gsolve( '/model/compartment/ksolve' )
+ stoich.ksolve = gsolve
+ else:
+ ksolve = moose.Ksolve( '/model/compartment/ksolve' )
+ stoich.ksolve = ksolve
+ stoich.path = "/model/compartment/##"
+ #solver.method = "rk5"
+ #mesh = moose.element( "/model/compartment/mesh" )
+ #moose.connect( mesh, "remesh", solver, "remesh" )
+ moose.setClock( 5, 1.0 ) # clock for the solver
+ moose.useClock( 5, '/model/compartment/ksolve', 'process' )
+
+ moose.reinit()
+ moose.start( 100.0 ) # Run the model for 100 seconds.
+
+ a = moose.element( '/model/compartment/a' )
+ b = moose.element( '/model/compartment/b' )
+
+ # move most molecules over to b
+ b.conc = b.conc + a.conc * 0.9
+ a.conc = a.conc * 0.1
+ moose.start( 100.0 ) # Run the model for 100 seconds.
+
+ # move most molecules back to a
+ a.conc = a.conc + b.conc * 0.99
+ b.conc = b.conc * 0.01
+ moose.start( 100.0 ) # Run the model for 100 seconds.
+
+ # Iterate through all plots, dump their contents to data.plot.
+ displayPlots()
+
+ quit()
+
+ # Run the 'main' if this script is executed standalone.
+ if __name__ == '__main__':
+ main()
+
+
+|
+
+**Output:**
+
+.. image:: ../../../images/simpleB.png
+
+
+Scale Volumes
+=============
+
+File name: **scaleVolumes.py**
+
+This script runs exactly the same model as in simpleBis.py, but it
+automatically scales the volumes from 1e-19 down to smaller values.
+
+Note how the simulation successively becomes noisier, until at very
+small volumes there are spontaneous state transitions.
+
+**Code:**
+
+.. hidden-code-block:: python
+ :linenos:
+ :label: Show/Hide code
+
+ #########################################################################
+ ## This program is part of 'MOOSE', the
+ ## Messaging Object Oriented Simulation Environment.
+ ## Copyright (C) 2013 Upinder S. Bhalla. and NCBS
+ ## It is made available under the terms of the
+ ## GNU Lesser General Public License version 2.1
+ ## See the file COPYING.LIB for the full notice.
+ #########################################################################
+
+ import math
+ import pylab
+ import numpy
+ import moose
+
+ def makeModel():
+ # create container for model
+ model = moose.Neutral( 'model' )
+ compartment = moose.CubeMesh( '/model/compartment' )
+ compartment.volume = 1e-20
+ # the mesh is created automatically by the compartment
+ mesh = moose.element( '/model/compartment/mesh' )
+
+ # create molecules and reactions
+ a = moose.Pool( '/model/compartment/a' )
+ b = moose.Pool( '/model/compartment/b' )
+ c = moose.Pool( '/model/compartment/c' )
+ enz1 = moose.Enz( '/model/compartment/b/enz1' )
+ enz2 = moose.Enz( '/model/compartment/c/enz2' )
+ cplx1 = moose.Pool( '/model/compartment/b/enz1/cplx' )
+ cplx2 = moose.Pool( '/model/compartment/c/enz2/cplx' )
+ reac = moose.Reac( '/model/compartment/reac' )
+
+ # connect them up for reactions
+ moose.connect( enz1, 'sub', a, 'reac' )
+ moose.connect( enz1, 'prd', b, 'reac' )
+ moose.connect( enz1, 'enz', b, 'reac' )
+ moose.connect( enz1, 'cplx', cplx1, 'reac' )
+
+ moose.connect( enz2, 'sub', b, 'reac' )
+ moose.connect( enz2, 'prd', a, 'reac' )
+ moose.connect( enz2, 'enz', c, 'reac' )
+ moose.connect( enz2, 'cplx', cplx2, 'reac' )
+
+ moose.connect( reac, 'sub', a, 'reac' )
+ moose.connect( reac, 'prd', b, 'reac' )
+
+ # connect them up to the compartment for volumes
+ #for x in ( a, b, c, cplx1, cplx2 ):
+ # moose.connect( x, 'mesh', mesh, 'mesh' )
+
+ # Assign parameters
+ a.concInit = 1
+ b.concInit = 0
+ c.concInit = 0.01
+ enz1.kcat = 0.4
+ enz1.Km = 4
+ enz2.kcat = 0.6
+ enz2.Km = 0.01
+ reac.Kf = 0.001
+ reac.Kb = 0.01
+
+ # Create the output tables
+ graphs = moose.Neutral( '/model/graphs' )
+ outputA = moose.Table ( '/model/graphs/concA' )
+ outputB = moose.Table ( '/model/graphs/concB' )
+
+ # connect up the tables
+ moose.connect( outputA, 'requestOut', a, 'getConc' );
+ moose.connect( outputB, 'requestOut', b, 'getConc' );
+
+ # Schedule the whole lot
+ moose.setClock( 4, 0.01 ) # for the computational objects
+ moose.setClock( 8, 1.0 ) # for the plots
+ # The wildcard uses # for single level, and ## for recursive.
+ moose.useClock( 4, '/model/compartment/##', 'process' )
+ moose.useClock( 8, '/model/graphs/#', 'process' )
+
+ def displayPlots():
+ for x in moose.wildcardFind( '/model/graphs/conc#' ):
+ t = numpy.arange( 0, x.vector.size, 1 ) #sec
+ pylab.plot( t, x.vector, label=x.name )
+
+ def main():
+
+ """
+ This example illustrates how to run a model at different volumes.
+ The key line is just to set the volume of the compartment::
+
+ compt.volume = vol
+
+ If everything
+ else is set up correctly, then this change propagates through to all
+ reactions molecules.
+
+ For a deterministic reaction one would not see any change in output
+ concentrations.
+ For a stochastic reaction illustrated here, one sees the level of
+ 'noise'
+ changing, even though the concentrations are similar up to a point.
+ This example creates a bistable model having two enzymes and a reaction.
+ One of the enzymes is autocatalytic.
+ This model is set up within the script rather than using an external
+ file.
+ The model is set up to run using the GSSA (Gillespie Stocahstic systems
+ algorithim) method in MOOSE.
+
+ To run the example, run the script
+
+ ``python scaleVolumes.py``
+
+ and close the plots every cycle to see the outcome of stochastic
+ calculations at ever smaller volumes, keeping concentrations the same.
+ """
+ makeModel()
+ moose.seed( 11111 )
+ gsolve = moose.Gsolve( '/model/compartment/gsolve' )
+ stoich = moose.Stoich( '/model/compartment/stoich' )
+ compt = moose.element( '/model/compartment' );
+ stoich.compartment = compt
+ stoich.ksolve = gsolve
+ stoich.path = "/model/compartment/##"
+ moose.setClock( 5, 1.0 ) # clock for the solver
+ moose.useClock( 5, '/model/compartment/gsolve', 'process' )
+ a = moose.element( '/model/compartment/a' )
+
+ for vol in ( 1e-19, 1e-20, 1e-21, 3e-22, 1e-22, 3e-23, 1e-23 ):
+ # Set the volume
+ compt.volume = vol
+ print('vol = {}, a.concInit = {}, a.nInit = {}'.format( vol, a.concInit, a.nInit))
+ print('Close graph to go to next plot\n')
+
+ moose.reinit()
+ moose.start( 100.0 ) # Run the model for 100 seconds.
+
+ a = moose.element( '/model/compartment/a' )
+ b = moose.element( '/model/compartment/b' )
+
+ # move most molecules over to b
+ b.conc = b.conc + a.conc * 0.9
+ a.conc = a.conc * 0.1
+ moose.start( 100.0 ) # Run the model for 100 seconds.
+
+ # move most molecules back to a
+ a.conc = a.conc + b.conc * 0.99
+ b.conc = b.conc * 0.01
+ moose.start( 100.0 ) # Run the model for 100 seconds.
+
+ # Iterate through all plots, dump their contents to data.plot.
+ displayPlots()
+ pylab.show()
+
+ quit()
+
+ # Run the 'main' if this script is executed standalone.
+ if __name__ == '__main__':
+ main()
+|
+**Output:**
+
+.. parsed-literal::
+
+ vol = 1e-19, a.concInit = 1.0, a.nInit = 60221.415
+
+
+
+
+.. image:: ../../../images/sV1.png
+
+
+.. parsed-literal::
+
+ vol = 1e-20, a.concInit = 1.0, a.nInit = 6022.1415
+
+
+
+
+.. image:: ../../../images/sV2.png
+
+
+.. parsed-literal::
+
+ vol = 1e-21, a.concInit = 1.0, a.nInit = 602.21415
+
+
+
+
+.. image:: ../../../images/sV3.png
+
+
+.. parsed-literal::
+
+ vol = 3e-22, a.concInit = 1.0, a.nInit = 180.664245
+
+
+
+
+.. image:: ../../../images/sV4.png
+
+
+.. parsed-literal::
+
+ vol = 1e-22, a.concInit = 1.0, a.nInit = 60.221415
+
+
+
+
+.. image:: ../../../images/sV5.png
+
+
+.. parsed-literal::
+
+ vol = 3e-23, a.concInit = 1.0, a.nInit = 18.0664245
+
+
+
+
+.. image:: ../../../images/sV6.png
+
+
+.. parsed-literal::
+
+ vol = 1e-23, a.concInit = 1.0, a.nInit = 6.0221415
+
+
+
+
+.. image:: ../../../images/sV7.png
+
+
+Strong Bistable System
+======================
+
+File name: **strongBis.py**
+
+This example illustrates a particularly strong, that is, parametrically
+robust bistable system. The model topology is symmetric between
+molecules **b** and **c**. We have both positive feedback of molecules
+**b** and **c** onto themselves, and also inhibition of **b** by **c**
+and vice versa.
+
+.. image:: ../../../images/strongBis.png
+
+Open the python file to see what is happening. The simulation starts at
+a symmetric point and the model settles at precisely the balance point
+where **a**, **b**, and **c** are at the same concentration. At t= 100
+we apply a small molecular 'tap' to push it over to a state where **c**
+is larger. This is stable. At t = 210 we apply a moderate push to show
+that it is now very stably in this state, and the system rebounds to its
+original levels. At t = 320 we apply a strong push to take it over to a
+state where **b** is larger. At t = 430 we give it a strong push to take
+it back to the **c** dominant state.
+
+**Code:**
+
+.. hidden-code-block:: python
+ :linenos:
+ :label: Show/Hide code
+
+ #########################################################################
+ ## This program is part of 'MOOSE', the
+ ## Messaging Object Oriented Simulation Environment.
+ ## Copyright (C) 2014 Upinder S. Bhalla. and NCBS
+ ## It is made available under the terms of the
+ ## GNU Lesser General Public License version 2.1
+ ## See the file COPYING.LIB for the full notice.
+ #########################################################################
+
+ import moose
+ import matplotlib.pyplot as plt
+ import matplotlib.image as mpimg
+ import pylab
+ import numpy
+ import sys
+
+ def main():
+
+ solver = "gsl" # Pick any of gsl, gssa, ee..
+ #solver = "gssa" # Pick any of gsl, gssa, ee..
+ #moose.seed( 1234 ) # Needed if stochastic.
+ mfile = '../../genesis/M1719.g'
+ runtime = 100.0
+ if ( len( sys.argv ) >= 2 ):
+ solver = sys.argv[1]
+ modelId = moose.loadModel( mfile, 'model', solver )
+ # Increase volume so that the stochastic solver gssa
+ # gives an interesting output
+ compt = moose.element( '/model/kinetics' )
+ compt.volume = 0.2e-19
+ r = moose.element( '/model/kinetics/equil' )
+
+ moose.reinit()
+ moose.start( runtime )
+ r.Kf *= 1.1 # small tap to break symmetry
+ moose.start( runtime/10 )
+ r.Kf = r.Kb
+ moose.start( runtime )
+
+ r.Kb *= 2.0 # Moderate push does not tip it back.
+ moose.start( runtime/10 )
+ r.Kb = r.Kf
+ moose.start( runtime )
+
+ r.Kb *= 5.0 # Strong push does tip it over
+ moose.start( runtime/10 )
+ r.Kb = r.Kf
+ moose.start( runtime )
+ r.Kf *= 5.0 # Strong push tips it back.
+ moose.start( runtime/10 )
+ r.Kf = r.Kb
+ moose.start( runtime )
+
+
+ # Display all plots.
+ img = mpimg.imread( 'strongBis.png' )
+ fig = plt.figure( figsize=(12, 10 ) )
+ png = fig.add_subplot( 211 )
+ imgplot = plt.imshow( img )
+ ax = fig.add_subplot( 212 )
+ x = moose.wildcardFind( '/model/#graphs/conc#/#' )
+ dt = moose.element( '/clock' ).tickDt[18]
+ t = numpy.arange( 0, x[0].vector.size, 1 ) * dt
+ ax.plot( t, x[0].vector, 'r-', label=x[0].name )
+ ax.plot( t, x[1].vector, 'g-', label=x[1].name )
+ ax.plot( t, x[2].vector, 'b-', label=x[2].name )
+ plt.ylabel( 'Conc (mM)' )
+ plt.xlabel( 'Time (seconds)' )
+ pylab.legend()
+ pylab.show()
+
+ # Run the 'main' if this script is executed standalone.
+ if __name__ == '__main__':
+ main()
+|
+**Output:**
+
+.. image:: ../../../images/strongB.png
+
+
+MAPK Feedback Model
+===================
+
+File name: **mapkFB.py**
+
+This example illustrates loading, and running a kinetic model for a much
+more complex bistable positive feedback system, defined in kkit format.
+This is based on Bhalla, Ram and Iyengar, Science 2002.
+
+The core of this model is a positive feedback loop comprising of the
+MAPK cascade, PLA2, and PKC. It receives PDGF and Ca2+ as inputs.
+
+.. image:: ../../../images/mapkFB.png
+
+This model is quite a large one and due to some stiffness in its
+equations, it takes about 30 seconds to execute. Note that this is still
+200 times faster than the events it models.
+
+The simulation illustrated here shows how the model starts out in a
+state of low activity. It is induced to 'turn on' when a a PDGF stimulus
+is given for 400 seconds, starting at t = 500s. After it has settled to
+the new 'on' state, the model is made to 'turn off' by setting the
+system calcium levels to zero. This inhibition starts at t = 2900 and
+goes on for 500 s.
+
+Note that this is a somewhat unphysiological manipulation! Following
+this the model settles back to the same 'off' state it was in
+originally.
+
+**Code:**
+
+.. hidden-code-block:: python
+ :linenos:
+ :label: Show/Hide code
+
+ #########################################################################
+ ## This program is part of 'MOOSE', the
+ ## Messaging Object Oriented Simulation Environment.
+ ## Copyright (C) 2014 Upinder S. Bhalla. and NCBS
+ ## It is made available under the terms of the
+ ## GNU Lesser General Public License version 2.1
+ ## See the file COPYING.LIB for the full notice.
+ #########################################################################
+
+ import moose
+ import matplotlib.pyplot as plt
+ import matplotlib.image as mpimg
+ import pylab
+ import numpy
+ import sys
+ import os
+
+ scriptDir = os.path.dirname( os.path.realpath( __file__ ) )
+
+ def main():
+ """
+ This example illustrates loading, and running a kinetic model
+ for a bistable positive feedback system, defined in kkit format.
+ This is based on Bhalla, Ram and Iyengar, Science 2002.
+
+ The core of this model is a positive feedback loop comprising of
+ the MAPK cascade, PLA2, and PKC. It receives PDGF and Ca2+ as
+ inputs.
+
+ This model is quite a large one and due to some stiffness in its
+ equations, it runs somewhat slowly.
+
+ The simulation illustrated here shows how the model starts out in
+ a state of low activity. It is induced to 'turn on' when a
+ a PDGF stimulus is given for 400 seconds.
+ After it has settled to the new 'on' state, model is made to
+ 'turn off'
+ by setting the system calcium levels to zero for a while. This
+ is a somewhat unphysiological manipulation!
+
+ """
+
+ solver = "gsl" # Pick any of gsl, gssa, ee..
+ #solver = "gssa" # Pick any of gsl, gssa, ee..
+ mfile = os.path.join( scriptDir, '..', '..', 'genesis' , 'acc35.g' )
+ runtime = 2000.0
+ if ( len( sys.argv ) == 2 ):
+ solver = sys.argv[1]
+ modelId = moose.loadModel( mfile, 'model', solver )
+ # Increase volume so that the stochastic solver gssa
+ # gives an interesting output
+ compt = moose.element( '/model/kinetics' )
+ compt.volume = 5e-19
+
+ moose.reinit()
+ moose.start( 500 )
+ moose.element( '/model/kinetics/PDGFR/PDGF' ).concInit = 0.0001
+ moose.start( 400 )
+ moose.element( '/model/kinetics/PDGFR/PDGF' ).concInit = 0.0
+ moose.start( 2000 )
+ moose.element( '/model/kinetics/Ca' ).concInit = 0.0
+ moose.start( 500 )
+ moose.element( '/model/kinetics/Ca' ).concInit = 0.00008
+ moose.start( 2000 )
+
+ # Display all plots.
+ img = mpimg.imread( 'mapkFB.png' )
+ fig = plt.figure( figsize=(12, 10 ) )
+ png = fig.add_subplot( 211 )
+ imgplot = plt.imshow( img )
+ ax = fig.add_subplot( 212 )
+ x = moose.wildcardFind( '/model/#graphs/conc#/#' )
+ t = numpy.arange( 0, x[0].vector.size, 1 ) * x[0].dt
+ ax.plot( t, x[0].vector, 'b-', label=x[0].name )
+ ax.plot( t, x[1].vector, 'c-', label=x[1].name )
+ ax.plot( t, x[2].vector, 'r-', label=x[2].name )
+ ax.plot( t, x[3].vector, 'm-', label=x[3].name )
+ plt.ylabel( 'Conc (mM)' )
+ plt.xlabel( 'Time (seconds)' )
+ pylab.legend()
+ pylab.show()
+
+ # Run the 'main' if this script is executed standalone.
+ if __name__ == '__main__':
+ main()
+|
+
+**Output:**
+
+.. image:: ../../../images/mapkFB2.png
+
+
+Propogation of a Bistable System
+================================
+
+File name: **propagationBis.py**
+
+All the above models have been well-mixed, that is point or non-spatial
+models. Bistables do interesting things when they are dispersed in
+space. This is illustrated in this example. Here we have a tapering
+cylinder, that is a pseudo 1-dimensional reaction-diffusion system.
+Every point in this cylinder has the bistable system from the strongBis
+example.
+
+.. image:: ../../../images/propBis.png
+
+The example has two stages. First it starts out with the model in the
+unstable transition point, and introduces a small symmetry-breaking
+perturbation at one end. This rapidly propagates through the entire
+length model, leaving molecule **b** at a higher value than **c**.
+
+At t = 100 we carry out a different manipulation. We flip the
+concentrations of molecules b and c for the left half of the model, and
+then just let it run. Now we have opposing bistable states on either
+half. In the middle, the two systems battle it out. Molecule **c** from
+the left side diffuses over to the right, and tries to inhibit **b**,
+and vice versa. However we have a small asymmetry due to the tapering of
+the cylinder. As there is a slightly larger volume on the left, the
+transition point gradually advances to the right, as molecule **b**
+yields to the slightly larger amounts of molecule **c**.
+
+**Code:**
+
+.. hidden-code-block:: python
+ :linenos:
+ :label: Show/Hide code
+
+ #########################################################################
+ ## This program is part of 'MOOSE', the
+ ## Messaging Object Oriented Simulation Environment.
+ ## Copyright (C) 2014 Upinder S. Bhalla. and NCBS
+ ## It is made available under the terms of the
+ ## GNU Lesser General Public License version 2.1
+ ## See the file COPYING.LIB for the full notice.
+ #########################################################################
+
+ """
+ This example illustrates propagation of state flips in a
+ linear 1-dimensional reaction-diffusion system. It uses a
+ bistable system loaded in from a kkit definition file, and
+ places this in a tapering cylinder for pseudo 1-dimentionsional
+ diffusion.
+
+ This example illustrates a number of features of reaction-diffusion
+ calculations.
+
+ First, it shows how to set up such systems. Key steps are to create
+ the compartment and define its voxelization, then create the Ksolve,
+ Dsolve, and Stoich. Then we assign stoich.compartment, ksolve and
+ dsolve in that order. Finally we assign the path of the Stoich.
+
+ For running the model, we start by introducing
+ a small symmetry-breaking increment of concInit
+ of the molecule **b** in the last compartment on the cylinder. The model
+ starts out with molecules at equal concentrations, so that the system would
+ settle to the unstable fixed point. This symmetry breaking leads
+ to the last compartment moving towards the state with an
+ increased concentration of **b**,
+ and this effect propagates to all other compartments.
+
+ Once the model has settled to the state where **b** is high throughout,
+ we simply exchange the concentrations of **b** with **c** in the left
+ half of the cylinder. This introduces a brief transient at the junction,
+ which soon settles to a smooth crossover.
+
+ Finally, as we run the simulation, the tapering geometry comes into play.
+ Since the left hand side has a larger diameter than the right, the
+ state on the left gradually wins over and the transition point slowly
+ moves to the right.
+
+ """
+
+ import math
+ import numpy
+ import matplotlib.pyplot as plt
+ import matplotlib.image as mpimg
+ import moose
+ import sys
+
+ def makeModel():
+ # create container for model
+ r0 = 1e-6 # m
+ r1 = 0.5e-6 # m. Note taper.
+ num = 200
+ diffLength = 1e-6 # m
+ comptLength = num * diffLength # m
+ diffConst = 20e-12 # m^2/sec
+ concA = 1 # millimolar
+ diffDt = 0.02 # for the diffusion
+ chemDt = 0.2 # for the reaction
+ mfile = '../../genesis/M1719.g'
+
+ model = moose.Neutral( 'model' )
+ compartment = moose.CylMesh( '/model/kinetics' )
+
+ # load in model
+ modelId = moose.loadModel( mfile, '/model', 'ee' )
+ a = moose.element( '/model/kinetics/a' )
+ b = moose.element( '/model/kinetics/b' )
+ c = moose.element( '/model/kinetics/c' )
+
+ ac = a.concInit
+ bc = b.concInit
+ cc = c.concInit
+
+ compartment.r0 = r0
+ compartment.r1 = r1
+ compartment.x0 = 0
+ compartment.x1 = comptLength
+ compartment.diffLength = diffLength
+ assert( compartment.numDiffCompts == num )
+
+ # Assign parameters
+ for x in moose.wildcardFind( '/model/kinetics/##[ISA=PoolBase]' ):
+ #print 'pools: ', x, x.name
+ x.diffConst = diffConst
+
+ # Make solvers
+ ksolve = moose.Ksolve( '/model/kinetics/ksolve' )
+ dsolve = moose.Dsolve( '/model/dsolve' )
+ # Set up clocks.
+ moose.setClock( 10, diffDt )
+ for i in range( 11, 17 ):
+ moose.setClock( i, chemDt )
+
+ stoich = moose.Stoich( '/model/kinetics/stoich' )
+ stoich.compartment = compartment
+ stoich.ksolve = ksolve
+ stoich.dsolve = dsolve
+ stoich.path = "/model/kinetics/##"
+ print(('dsolve.numPools, num = ', dsolve.numPools, num))
+ b.vec[num-1].concInit *= 1.01 # Break symmetry.
+
+ def main():
+ runtime = 100
+ displayInterval = 2
+ makeModel()
+ dsolve = moose.element( '/model/dsolve' )
+ moose.reinit()
+ #moose.start( runtime ) # Run the model for 10 seconds.
+
+ a = moose.element( '/model/kinetics/a' )
+ b = moose.element( '/model/kinetics/b' )
+ c = moose.element( '/model/kinetics/c' )
+
+ img = mpimg.imread( 'propBis.png' )
+ #imgplot = plt.imshow( img )
+ #plt.show()
+
+ plt.ion()
+ fig = plt.figure( figsize=(12,10) )
+ png = fig.add_subplot(211)
+ imgplot = plt.imshow( img )
+ ax = fig.add_subplot(212)
+ ax.set_ylim( 0, 0.001 )
+ plt.ylabel( 'Conc (mM)' )
+ plt.xlabel( 'Position along cylinder (microns)' )
+ pos = numpy.arange( 0, a.vec.conc.size, 1 )
+ line1, = ax.plot( pos, a.vec.conc, 'r-', label='a' )
+ line2, = ax.plot( pos, b.vec.conc, 'g-', label='b' )
+ line3, = ax.plot( pos, c.vec.conc, 'b-', label='c' )
+ timeLabel = plt.text(60, 0.0009, 'time = 0')
+ plt.legend()
+ fig.canvas.draw()
+
+ for t in range( displayInterval, runtime, displayInterval ):
+ moose.start( displayInterval )
+ line1.set_ydata( a.vec.conc )
+ line2.set_ydata( b.vec.conc )
+ line3.set_ydata( c.vec.conc )
+ timeLabel.set_text( "time = %d" % t )
+ fig.canvas.draw()
+
+ print('Swapping concs of b and c in half the cylinder')
+ for i in range( b.numData/2 ):
+ temp = b.vec[i].conc
+ b.vec[i].conc = c.vec[i].conc
+ c.vec[i].conc = temp
+
+ newruntime = 200
+ for t in range( displayInterval, newruntime, displayInterval ):
+ moose.start( displayInterval )
+ line1.set_ydata( a.vec.conc )
+ line2.set_ydata( b.vec.conc )
+ line3.set_ydata( c.vec.conc )
+ timeLabel.set_text( "time = %d" % (t + runtime) )
+ fig.canvas.draw()
+
+ print( "Hit 'enter' to exit" )
+ sys.stdin.read(1)
+
+
+
+ # Run the 'main' if this script is executed standalone.
+ if __name__ == '__main__':
+ main()
+
+|
+
+**Output:**
+
+.. image:: ../../../images/propBis.gif
+
+
+Steady-state Finder
+===================
+
+File name: **findSteadyState**
+
+This is an example of how to use an internal MOOSE solver to find steady
+states of a system very rapidly. The method starts from a random
+position in state space that obeys mass conservation. It then finds the
+nearest steady state and reports it. If it does this enough times it
+should find all the steady states.
+
+We illustrate this process for 50 attempts to find the steady states. It
+does find all of them. Each time it plots and prints the values, though
+the plotting is not necessary.
+
+The printout shows the concentrations of all molecules in the first 5
+columns. Then it prints the type of solution, and the numbers of
+negative and positive eigenvalues. In all cases the calculations are
+successful, though it takes different numbers of iterations to arrive at
+the steady state. In some models it would be necessary to put a cap on
+the number of iterations, if the system is not able to find a steady
+state.
+
+In this example we run the bistable model using the ODE solver right at
+the end, and manually enforce transitions to show where the target
+steady states are.
+
+For more information on the algorithm used, look in the comments within
+the main method of the code below.
+
+**Code:**
+
+.. hidden-code-block:: python
+ :linenos:
+ :label: Show/Hide code
+
+ #########################################################################
+ ## This program is part of 'MOOSE', the
+ ## Messaging Object Oriented Simulation Environment.
+ ## Copyright (C) 2013 Upinder S. Bhalla. and NCBS
+ ## It is made available under the terms of the
+ ## GNU Lesser General Public License version 2.1
+ ## See the file COPYING.LIB for the full notice.
+ #########################################################################
+
+ from __future__ import print_function
+
+ import math
+ import pylab
+ import numpy
+ import moose
+
+ def main():
+ """
+ This example sets up the kinetic solver and steady-state finder, on
+ a bistable model of a chemical system. The model is set up within the
+ script.
+ The algorithm calls the steady-state finder 50 times with different
+ (randomized) initial conditions, as follows:
+
+ * Set up the random initial condition that fits the conservation laws
+ * Run for 2 seconds. This should not be mathematically necessary, but
+ for obscure numerical reasons it makes it much more likely that the
+ steady state solver will succeed in finding a state.
+ * Find the fixed point
+ * Print out the fixed point vector and various diagnostics.
+ * Run for 10 seconds. This is completely unnecessary, and is done here
+ just so that the resultant graph will show what kind of state has
+ been found.
+
+ After it does all this, the program runs for 100 more seconds on the
+ last found fixed point (which turns out to be a saddle node), then
+ is hard-switched in the script to the first attractor basin from which
+ it runs for another 100 seconds till it settles there, and then
+ is hard-switched yet again to the second attractor and runs for 400
+ seconds.
+
+ Looking at the output you will see many features of note:
+
+ * the first attractor (stable point) and the saddle point (unstable
+ fixed point) are both found quite often. But the second
+ attractor is found just once.
+ It has a very small basin of attraction.
+ * The values found for each of the fixed points match well with the
+ values found by running the system to steady-state at the end.
+ * There are a large number of failures to find a fixed point. These are
+ found and reported in the diagnostics. They show up on the plot
+ as cases where the 10-second runs are not flat.
+
+ If you wanted to find fixed points in a production model, you would
+ not need to do the 10-second runs, and you would need to eliminate the
+ cases where the state-finder failed. Then you could identify the good
+ points and keep track of how many of each were found.
+
+ There is no way to guarantee that all fixed points have been found
+ using this algorithm! If there are points in an obscure corner of state
+ space (as for the singleton second attractor convergence in this
+ example) you may have to iterate very many times to find them.
+
+ You may wish to sample concentration space logarithmically rather than
+ linearly.
+ """
+ compartment = makeModel()
+ ksolve = moose.Ksolve( '/model/compartment/ksolve' )
+ stoich = moose.Stoich( '/model/compartment/stoich' )
+ stoich.compartment = compartment
+ stoich.ksolve = ksolve
+ stoich.path = "/model/compartment/##"
+ state = moose.SteadyState( '/model/compartment/state' )
+
+ moose.reinit()
+ state.stoich = stoich
+ state.showMatrices()
+ state.convergenceCriterion = 1e-6
+ moose.seed( 111 ) # Used when generating the samples in state space
+
+ for i in range( 0, 50 ):
+ getState( ksolve, state )
+
+ # Now display the states of the system at more length to compare.
+ moose.start( 100.0 ) # Run the model for 100 seconds.
+
+ a = moose.element( '/model/compartment/a' )
+ b = moose.element( '/model/compartment/b' )
+
+ # move most molecules over to b
+ b.conc = b.conc + a.conc * 0.9
+ a.conc = a.conc * 0.1
+ moose.start( 100.0 ) # Run the model for 100 seconds.
+
+ # move most molecules back to a
+ a.conc = a.conc + b.conc * 0.99
+ b.conc = b.conc * 0.01
+ moose.start( 400.0 ) # Run the model for 200 seconds.
+
+ # Iterate through all plots, dump their contents to data.plot.
+ displayPlots()
+
+ quit()
+
+ def makeModel():
+ """ This function creates a bistable reaction system using explicit
+ MOOSE calls rather than load from a file
+ """
+ # create container for model
+ model = moose.Neutral( 'model' )
+ compartment = moose.CubeMesh( '/model/compartment' )
+ compartment.volume = 1e-15
+ # the mesh is created automatically by the compartment
+ mesh = moose.element( '/model/compartment/mesh' )
+
+ # create molecules and reactions
+ a = moose.Pool( '/model/compartment/a' )
+ b = moose.Pool( '/model/compartment/b' )
+ c = moose.Pool( '/model/compartment/c' )
+ enz1 = moose.Enz( '/model/compartment/b/enz1' )
+ enz2 = moose.Enz( '/model/compartment/c/enz2' )
+ cplx1 = moose.Pool( '/model/compartment/b/enz1/cplx' )
+ cplx2 = moose.Pool( '/model/compartment/c/enz2/cplx' )
+ reac = moose.Reac( '/model/compartment/reac' )
+
+ # connect them up for reactions
+ moose.connect( enz1, 'sub', a, 'reac' )
+ moose.connect( enz1, 'prd', b, 'reac' )
+ moose.connect( enz1, 'enz', b, 'reac' )
+ moose.connect( enz1, 'cplx', cplx1, 'reac' )
+
+ moose.connect( enz2, 'sub', b, 'reac' )
+ moose.connect( enz2, 'prd', a, 'reac' )
+ moose.connect( enz2, 'enz', c, 'reac' )
+ moose.connect( enz2, 'cplx', cplx2, 'reac' )
+
+ moose.connect( reac, 'sub', a, 'reac' )
+ moose.connect( reac, 'prd', b, 'reac' )
+
+ # Assign parameters
+ a.concInit = 1
+ b.concInit = 0
+ c.concInit = 0.01
+ enz1.kcat = 0.4
+ enz1.Km = 4
+ enz2.kcat = 0.6
+ enz2.Km = 0.01
+ reac.Kf = 0.001
+ reac.Kb = 0.01
+
+ # Create the output tables
+ graphs = moose.Neutral( '/model/graphs' )
+ outputA = moose.Table2 ( '/model/graphs/concA' )
+ outputB = moose.Table2 ( '/model/graphs/concB' )
+ outputC = moose.Table2 ( '/model/graphs/concC' )
+ outputCplx1 = moose.Table2 ( '/model/graphs/concCplx1' )
+ outputCplx2 = moose.Table2 ( '/model/graphs/concCplx2' )
+
+ # connect up the tables
+ moose.connect( outputA, 'requestOut', a, 'getConc' );
+ moose.connect( outputB, 'requestOut', b, 'getConc' );
+ moose.connect( outputC, 'requestOut', c, 'getConc' );
+ moose.connect( outputCplx1, 'requestOut', cplx1, 'getConc' );
+ moose.connect( outputCplx2, 'requestOut', cplx2, 'getConc' );
+
+ return compartment
+
+ def displayPlots():
+ for x in moose.wildcardFind( '/model/graphs/conc#' ):
+ t = numpy.arange( 0, x.vector.size, 1 ) #sec
+ pylab.plot( t, x.vector, label=x.name )
+ pylab.legend()
+ pylab.show()
+
+ def getState( ksolve, state ):
+ """ This function finds a steady state starting from a random
+ initial condition that is consistent with the stoichiometry rules
+ and the original model concentrations.
+ """
+ scale = 1.0 / ( 1e-15 * 6.022e23 )
+ state.randomInit() # Randomize init conditions, subject to stoichiometry
+ moose.start( 2.0 ) # Run the model for 2 seconds.
+ state.settle() # This function finds the steady states.
+ for x in ksolve.nVec[0]:
+ print( "{:.2f}".format( x * scale ), end=' ')
+
+ print( "Type={} NegEig={} PosEig={} status={} {} Iter={:2d}".format( state.stateType, state.nNegEigenvalues, state.nPosEigenvalues, state.solutionStatus, state.status, state.nIter))
+ moose.start( 10.0 ) # Run model for 10 seconds, just for display
+
+
+ # Run the 'main' if this script is executed standalone.
+ if __name__ == '__main__':
+ main()
+
+|
+
+**Output:**
+
+.. parsed-literal::
+
+ 0.92 0.05 0.00 0.01 0.01 Type=2 NegEig=2 PosEig=1 status=0 success Iter=16
+ 0.92 0.05 0.00 0.01 0.01 Type=2 NegEig=2 PosEig=1 status=0 success Iter=29
+ 0.18 0.75 0.00 0.03 0.01 Type=0 NegEig=3 PosEig=0 status=0 success Iter=10
+ 0.18 0.75 0.00 0.03 0.01 Type=0 NegEig=3 PosEig=0 status=0 success Iter=26
+ 0.18 0.75 0.00 0.03 0.01 Type=0 NegEig=3 PosEig=0 status=0 success Iter=27
+ 0.92 0.05 0.00 0.01 0.01 Type=2 NegEig=2 PosEig=1 status=0 success Iter=30
+ 0.18 0.75 0.00 0.03 0.01 Type=0 NegEig=3 PosEig=0 status=0 success Iter=12
+ 0.92 0.05 0.00 0.01 0.01 Type=2 NegEig=2 PosEig=1 status=0 success Iter=29
+ 0.18 0.75 0.00 0.03 0.01 Type=0 NegEig=3 PosEig=0 status=0 success Iter=12
+ 0.92 0.05 0.00 0.01 0.01 Type=2 NegEig=2 PosEig=1 status=0 success Iter=41
+ 0.18 0.75 0.00 0.03 0.01 Type=0 NegEig=3 PosEig=0 status=0 success Iter=29
+ 0.18 0.75 0.00 0.03 0.01 Type=0 NegEig=3 PosEig=0 status=0 success Iter=18
+ 0.18 0.75 0.00 0.03 0.01 Type=0 NegEig=3 PosEig=0 status=0 success Iter=27
+ 0.18 0.75 0.00 0.03 0.01 Type=0 NegEig=3 PosEig=0 status=0 success Iter=14
+ 0.18 0.75 0.00 0.03 0.01 Type=0 NegEig=3 PosEig=0 status=0 success Iter=12
+ 0.18 0.75 0.00 0.03 0.01 Type=0 NegEig=3 PosEig=0 status=0 success Iter=19
+ 0.18 0.75 0.00 0.03 0.01 Type=0 NegEig=3 PosEig=0 status=0 success Iter= 6
+ 0.18 0.75 0.00 0.03 0.01 Type=0 NegEig=3 PosEig=0 status=0 success Iter=14
+ 0.18 0.75 0.00 0.03 0.01 Type=0 NegEig=3 PosEig=0 status=0 success Iter=23
+ 0.18 0.75 0.00 0.03 0.01 Type=0 NegEig=3 PosEig=0 status=0 success Iter=25
+ 0.18 0.75 0.00 0.03 0.01 Type=0 NegEig=3 PosEig=0 status=0 success Iter=16
+ 0.18 0.75 0.00 0.03 0.01 Type=0 NegEig=3 PosEig=0 status=0 success Iter= 5
+ 0.92 0.05 0.00 0.01 0.01 Type=2 NegEig=2 PosEig=1 status=0 success Iter=43
+ 0.92 0.05 0.00 0.01 0.01 Type=2 NegEig=2 PosEig=1 status=0 success Iter= 9
+ 0.92 0.05 0.00 0.01 0.01 Type=2 NegEig=2 PosEig=1 status=0 success Iter=43
+ 0.92 0.05 0.00 0.01 0.01 Type=2 NegEig=2 PosEig=1 status=0 success Iter=29
+ 0.18 0.75 0.00 0.03 0.01 Type=0 NegEig=3 PosEig=0 status=0 success Iter=27
+ 0.18 0.75 0.00 0.03 0.01 Type=0 NegEig=3 PosEig=0 status=0 success Iter= 9
+ 0.18 0.75 0.00 0.03 0.01 Type=0 NegEig=3 PosEig=0 status=0 success Iter=12
+ 0.92 0.05 0.00 0.01 0.01 Type=2 NegEig=2 PosEig=1 status=0 success Iter=24
+ 0.92 0.05 0.00 0.01 0.01 Type=2 NegEig=2 PosEig=1 status=0 success Iter=26
+ 0.18 0.75 0.00 0.03 0.01 Type=0 NegEig=3 PosEig=0 status=0 success Iter=14
+ 0.92 0.05 0.00 0.01 0.01 Type=2 NegEig=2 PosEig=1 status=0 success Iter=14
+ 0.18 0.75 0.00 0.03 0.01 Type=0 NegEig=3 PosEig=0 status=0 success Iter=10
+ 0.18 0.75 0.00 0.03 0.01 Type=0 NegEig=3 PosEig=0 status=0 success Iter=13
+ 0.18 0.75 0.00 0.03 0.01 Type=0 NegEig=3 PosEig=0 status=0 success Iter=26
+ 0.18 0.75 0.00 0.03 0.01 Type=0 NegEig=3 PosEig=0 status=0 success Iter=21
+ 0.18 0.75 0.00 0.03 0.01 Type=0 NegEig=3 PosEig=0 status=0 success Iter=26
+ 0.18 0.75 0.00 0.03 0.01 Type=0 NegEig=3 PosEig=0 status=0 success Iter=24
+ 0.18 0.75 0.00 0.03 0.01 Type=0 NegEig=3 PosEig=0 status=0 success Iter=24
+ 0.92 0.05 0.00 0.01 0.01 Type=2 NegEig=2 PosEig=1 status=0 success Iter=18
+ 0.18 0.75 0.00 0.03 0.01 Type=0 NegEig=3 PosEig=0 status=0 success Iter=26
+ 0.18 0.75 0.00 0.03 0.01 Type=5 NegEig=4 PosEig=0 status=0 success Iter=13
+ 0.18 0.75 0.00 0.03 0.01 Type=0 NegEig=3 PosEig=0 status=0 success Iter=23
+ 0.92 0.05 0.00 0.01 0.01 Type=2 NegEig=2 PosEig=1 status=0 success Iter=24
+ 0.18 0.75 0.00 0.03 0.01 Type=0 NegEig=3 PosEig=0 status=0 success Iter= 8
+ 0.18 0.75 0.00 0.03 0.01 Type=0 NegEig=3 PosEig=0 status=0 success Iter=18
+ 0.18 0.75 0.00 0.03 0.01 Type=0 NegEig=3 PosEig=1 status=0 success Iter=21
+ 0.99 0.00 0.01 0.00 0.00 Type=0 NegEig=3 PosEig=0 status=0 success Iter=15
+ 0.92 0.05 0.00 0.01 0.01 Type=2 NegEig=2 PosEig=1 status=0 success Iter=29
+
+.. image:: ../../../images/findS.png
+
+Dose Response (Under construction)
+==================================
+
+File name: **doseResponse.py**
+
+This example generates a doseResponse plot for a bistable system,
+against a control parameter (dose) that takes the system in and out
+again from the bistable regime. Like the previous example, it uses the
+steady-state solver to find the stable points for each value of the
+control parameter. Unfortunately it doesn't work right now. Seems like
+the kcat scaling isn't being registered.
+
+**Code:**
+
+.. hidden-code-block:: python
+ :linenos:
+ :label: Show/Hide code
+
+ ## Makes and plots the dose response curve for bistable models
+ ## Author: Sahil Moza
+ ## June 26, 2014
+
+ import moose
+ import pylab
+ import numpy as np
+ from matplotlib import pyplot as plt
+
+ def setupSteadyState(simdt,plotDt):
+
+ ksolve = moose.Ksolve( '/model/kinetics/ksolve' )
+ stoich = moose.Stoich( '/model/kinetics/stoich' )
+ stoich.compartment = moose.element('/model/kinetics')
+
+ stoich.ksolve = ksolve
+ #ksolve.stoich = stoich
+ stoich.path = "/model/kinetics/##"
+ state = moose.SteadyState( '/model/kinetics/state' )
+
+ #### Set clocks here
+ #moose.useClock(4, "/model/kinetics/##[]", "process")
+ #moose.setClock(4, float(simdt))
+ #moose.setClock(5, float(simdt))
+ #moose.useClock(5, '/model/kinetics/ksolve', 'process' )
+ #moose.useClock(8, '/model/graphs/#', 'process' )
+ #moose.setClock(8, float(plotDt))
+
+ moose.reinit()
+
+ state.stoich = stoich
+ state.showMatrices()
+ state.convergenceCriterion = 1e-8
+
+ return ksolve, state
+
+ def parseModelName(fileName):
+ pos1=fileName.rfind('/')
+ pos2=fileName.rfind('.')
+ directory=fileName[:pos1]
+ prefix=fileName[pos1+1:pos2]
+ suffix=fileName[pos2+1:len(fileName)]
+ return directory, prefix, suffix
+
+ # Solve for the steady state
+ def getState( ksolve, state, vol):
+ scale = 1.0 / ( vol * 6.022e23 )
+ moose.reinit
+ state.randomInit() # Removing random initial condition to systematically make Dose reponse curves.
+ moose.start( 2.0 ) # Run the model for 2 seconds.
+ state.settle()
+
+ vector = []
+ a = moose.element( '/model/kinetics/a' ).conc
+ for x in ksolve.nVec[0]:
+ vector.append( x * scale)
+ moose.start( 10.0 ) # Run model for 10 seconds, just for display
+ failedSteadyState = any([np.isnan(x) for x in vector])
+
+ if not (failedSteadyState):
+ return state.stateType, state.solutionStatus, a, vector
+
+
+ def main():
+ # Setup parameters for simulation and plotting
+ simdt= 1e-2
+ plotDt= 1
+
+ # Factors to change in the dose concentration in log scale
+ factorExponent = 10 ## Base: ten raised to some power.
+ factorBegin = -20
+ factorEnd = 21
+ factorStepsize = 1
+ factorScale = 10.0 ## To scale up or down the factors
+
+ # Load Model and set up the steady state solver.
+ # model = sys.argv[1] # To load model from a file.
+ model = './19085.cspace'
+ modelPath, modelName, modelType = parseModelName(model)
+ outputDir = modelPath
+
+ modelId = moose.loadModel(model, 'model', 'ee')
+ dosePath = '/model/kinetics/b/DabX' # The dose entity
+
+ ksolve, state = setupSteadyState( simdt, plotDt)
+ vol = moose.element( '/model/kinetics' ).volume
+ iterInit = 100
+ solutionVector = []
+ factorArr = []
+
+ enz = moose.element(dosePath)
+ init = enz.kcat # Dose parameter
+
+ # Change Dose here to .
+ for factor in range(factorBegin, factorEnd, factorStepsize ):
+ scale = factorExponent ** (factor/factorScale)
+ enz.kcat = init * scale
+ print( "scale={:.3f}\tkcat={:.3f}".format( scale, enz.kcat) )
+ for num in range(iterInit):
+ stateType, solStatus, a, vector = getState( ksolve, state, vol)
+ if solStatus == 0:
+ #solutionVector.append(vector[0]/sum(vector))
+ solutionVector.append(a)
+ factorArr.append(scale)
+
+ joint = np.array([factorArr, solutionVector])
+ joint = joint[:,joint[1,:].argsort()]
+
+ # Plot dose response.
+ fig0 = plt.figure()
+ pylab.semilogx(joint[0,:],joint[1,:],marker="o",label = 'concA')
+ pylab.xlabel('Dose')
+ pylab.ylabel('Response')
+ pylab.suptitle('Dose-Reponse Curve for a bistable system')
+
+ pylab.legend(loc=3)
+ #plt.savefig(outputDir + "/" + modelName +"_doseResponse" + ".png")
+ plt.show()
+ #plt.close(fig0)
+ quit()
+
+
+
+ if __name__ == '__main__':
+ main()
+|
+**Output:**
+
+.. parsed-literal::
+
+ scale=0.010 kcat=0.004
+ scale=0.013 kcat=0.005
+ scale=0.016 kcat=0.006
+ scale=0.020 kcat=0.007
+ scale=0.025 kcat=0.009
+ scale=0.032 kcat=0.011
+ scale=0.040 kcat=0.014
+ scale=0.050 kcat=0.018
+ scale=0.063 kcat=0.023
+ scale=0.079 kcat=0.029
+ scale=0.100 kcat=0.036
+ scale=0.126 kcat=0.045
+ scale=0.158 kcat=0.057
+ scale=0.200 kcat=0.072
+ scale=0.251 kcat=0.091
+ scale=0.316 kcat=0.114
+ scale=0.398 kcat=0.144
+ scale=0.501 kcat=0.181
+ scale=0.631 kcat=0.228
+ scale=0.794 kcat=0.287
+ scale=1.000 kcat=0.361
+ scale=1.259 kcat=0.454
+ scale=1.585 kcat=0.572
+ scale=1.995 kcat=0.720
+ scale=2.512 kcat=0.907
+ scale=3.162 kcat=1.142
+ scale=3.981 kcat=1.437
+ scale=5.012 kcat=1.809
+ scale=6.310 kcat=2.278
+ scale=7.943 kcat=2.868
+ scale=10.000 kcat=3.610
+ scale=12.589 kcat=4.545
+ scale=15.849 kcat=5.722
+ scale=19.953 kcat=7.203
+ scale=25.119 kcat=9.068
+ scale=31.623 kcat=11.416
+ scale=39.811 kcat=14.372
+ scale=50.119 kcat=18.093
+ scale=63.096 kcat=22.778
+ scale=79.433 kcat=28.676
+ scale=100.000 kcat=36.101
+
+
+.. image:: ../../../images/doseR.png
+
+
diff --git a/docs/source/user/py/tutorials/ChemicalOscillators.rst b/docs/source/user/py/tutorials/ChemicalOscillators.rst
new file mode 100644
index 0000000000000000000000000000000000000000..c356ab6aabcd389356d929c5db7c9d2e1af0f290
--- /dev/null
+++ b/docs/source/user/py/tutorials/ChemicalOscillators.rst
@@ -0,0 +1,558 @@
+********************
+Chemical Oscillators
+********************
+
+`Chemical Oscillators <https://en.wikipedia.org/wiki/Chemical_clock>`_, also known as chemical clocks, are chemical systems in which the concentrations of one or more reactants undergoes periodic changes.
+
+These Oscillatory reactions can be modelled using moose. The examples below demonstrate different types of chemical oscillators, as well as how they can be simulated using moose. Each example has a short description, the code used in the simulation, and the default (gsl solver) output of the code.
+
+Each example can be found as a python file within the main moose folder under
+::
+
+ (...)/moose/moose-examples/tutorials/ChemicalOscillators
+
+In order to run the example, run the script
+::
+
+ python filename.py
+
+in command line, where ``filename.py`` is the name of the python file you would like to run. The filenames of each example are written in **bold** at the beginning of their respective sections, and the files themselves can be found in the aformentioned directory.
+
+In chemical models that use solvers, there are optional arguments that allow you to specify which solver you would like to use.
+::
+
+ python filename.py [gsl | gssa | ee]
+
+Where ``gsl`` is Gnu Scientific Library's deterministic solver, ``gssa`` stands for Gillespie stochastic simulation algorithm, and ``ee`` is the exponential euler algorithm.
+
+All the following examples can be run with either of the three solvers, which in some cases produces a different outcome. However, simply running the file without the optional argument will by default use the ``gsl`` solver. These ``gsl`` outputs are the ones shown below.
+
+|
+|
+
+Slow Feedback Oscillator
+========================
+
+File name: **slowFbOsc.py**
+
+This example illustrates loading, and running a kinetic model for a
+delayed -ve feedback oscillator, defined in kkit format. The model is
+one by Boris N. Kholodenko from Eur J Biochem. (2000) 267(6):1583-8
+
+.. image:: ../../../images/Kholodenko_tut.png
+
+This model has a high-gain MAPK stage, whose effects are visible whem
+one looks at the traces from successive stages in the plots. The
+upstream pools have small early peaks, and the downstream pools have
+large delayed ones. The negative feedback step is mediated by a simple
+binding reaction of the end-product of oscillation with an upstream
+activator.
+
+We use the gsl solver here. The model already defines some plots and
+sets the runtime to 4000 seconds. The model does not really play nicely
+with the GSSA solver, since it involves some really tiny amounts of the
+MAPKKK.
+
+Things to do with the model:
+
+::
+
+ - Look at model once it is loaded in::
+
+ moose.le( '/model' )
+ moose.showfields( '/model/kinetics/MAPK/MAPK' )
+
+ - Behold the amplification properties of the cascade. Could do this by blocking the feedback step and giving a small pulse input.
+ - Suggest which parameters you would alter to change the period of the oscillator:
+ - Concs of various molecules, for example::
+
+ ras_MAPKKKK = moose.element( '/model/kinetics/MAPK/Ras_dash_MKKKK' )
+ moose.showfields( ras_MAPKKKK )
+ ras_MAPKKKK.concInit = 1e-5
+ - Feedback reaction rates
+ - Rates of all the enzymes::
+
+ for i in moose.wildcardFind( '/##[ISA=EnzBase]'):
+ i.kcat *= 10.0
+
+**Code:**
+
+.. hidden-code-block:: python
+ :linenos:
+ :label: Show/Hide code
+
+ #########################################################################
+ ## This program is part of 'MOOSE', the
+ ## Messaging Object Oriented Simulation Environment.
+ ## Copyright (C) 2014 Upinder S. Bhalla. and NCBS
+ ## It is made available under the terms of the
+ ## GNU Lesser General Public License version 2.1
+ ## See the file COPYING.LIB for the full notice.
+ #########################################################################
+
+ import moose
+ import matplotlib.pyplot as plt
+ import matplotlib.image as mpimg
+ import pylab
+ import numpy
+ import sys
+
+ def main():
+
+ solver = "gsl"
+ mfile = '../../genesis/Kholodenko.g'
+ runtime = 5000.0
+ if ( len( sys.argv ) >= 2 ):
+ solver = sys.argv[1]
+ modelId = moose.loadModel( mfile, 'model', solver )
+ dt = moose.element( '/clock' ).tickDt[18]
+ moose.reinit()
+ moose.start( runtime )
+
+ # Display all plots.
+ img = mpimg.imread( 'Kholodenko_tut.png' )
+ fig = plt.figure( figsize=( 12, 10 ) )
+ png = fig.add_subplot( 211 )
+ imgplot = plt.imshow( img )
+ ax = fig.add_subplot( 212 )
+ x = moose.wildcardFind( '/model/#graphs/conc#/#' )
+ t = numpy.arange( 0, x[0].vector.size, 1 ) * dt
+ ax.plot( t, x[0].vector * 100, 'b-', label='Ras-MKKK * 100' )
+ ax.plot( t, x[1].vector, 'y-', label='MKKK-P' )
+ ax.plot( t, x[2].vector, 'm-', label='MKK-PP' )
+ ax.plot( t, x[3].vector, 'r-', label='MAPK-PP' )
+ plt.ylabel( 'Conc (mM)' )
+ plt.xlabel( 'Time (seconds)' )
+ pylab.legend()
+ pylab.show()
+
+ # Run the 'main' if this script is executed standalone.
+ if __name__ == '__main__':
+ main()
+
+|
+**Output:**
+
+
+.. image:: ../../../images/FB.png
+
+|
+|
+
+
+Turing Pattern Oscillator in One Dimension
+==========================================
+
+File name: **TuringOneDim.py**
+
+This example illustrates how to set up a oscillatory Turing pattern in
+1-D using reaction diffusion calculations. Reaction system 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.
+
+in sum, **a** has a positive feedback onto itself and also forms **b**.
+**b** has a negative feedback onto **a**. Finally, the diffusion
+constant for **a** is 1/10 that of **b**.
+
+.. image:: ../../../images/turingPatternTut.png
+
+This chemical system is present in a 1-dimensional (cylindrical)
+compartment. The entire reaction-diffusion system is set up within the
+script.
+
+**Code:**
+
+.. hidden-code-block:: python
+ :linenos:
+ :label: Show/Hide code
+
+ #########################################################################
+ ## This program is part of 'MOOSE', the
+ ## Messaging Object Oriented Simulation Environment.
+ ## Copyright (C) 2014 Upinder S. Bhalla. and NCBS
+ ## It is made available under the terms of the
+ ## GNU Lesser General Public License version 2.1
+ ## See the file COPYING.LIB for the full notice.
+ #########################################################################
+
+ import math
+ import numpy
+ import matplotlib.pyplot as plt
+ import matplotlib.image as mpimg
+ import moose
+
+ def makeModel():
+
+ # create container for model
+ r0 = 1e-6 # m
+ r1 = 1e-6 # m
+ num = 100
+ diffLength = 1e-6 # m
+ len = num * diffLength # m
+ diffConst = 5e-12 # m^2/sec
+ motorRate = 1e-6 # m/sec
+ concA = 1 # millimolar
+ dt4 = 0.02 # for the diffusion
+ dt5 = 0.2 # for the reaction
+
+ model = moose.Neutral( 'model' )
+ compartment = moose.CylMesh( '/model/compartment' )
+ compartment.r0 = r0
+ compartment.r1 = r1
+ compartment.x0 = 0
+ compartment.x1 = len
+ compartment.diffLength = diffLength
+
+ assert( compartment.numDiffCompts == num )
+
+ # create molecules and reactions
+ a = moose.Pool( '/model/compartment/a' )
+ b = moose.Pool( '/model/compartment/b' )
+ s = moose.Pool( '/model/compartment/s' )
+ e1 = moose.MMenz( '/model/compartment/e1' )
+ e2 = moose.MMenz( '/model/compartment/e2' )
+ e3 = moose.MMenz( '/model/compartment/e3' )
+ r1 = moose.Reac( '/model/compartment/r1' )
+ moose.connect( e1, 'sub', s, 'reac' )
+ moose.connect( e1, 'prd', a, 'reac' )
+ moose.connect( a, 'nOut', e1, 'enzDest' )
+ e1.Km = 1
+ e1.kcat = 1
+
+ moose.connect( e2, 'sub', s, 'reac' )
+ moose.connect( e2, 'prd', b, 'reac' )
+ moose.connect( a, 'nOut', e2, 'enzDest' )
+ e2.Km = 1
+ e2.kcat = 0.5
+
+ moose.connect( e3, 'sub', a, 'reac' )
+ moose.connect( e3, 'prd', s, 'reac' )
+ moose.connect( b, 'nOut', e3, 'enzDest' )
+ e3.Km = 0.1
+ e3.kcat = 1
+
+ moose.connect( r1, 'sub', b, 'reac' )
+ moose.connect( r1, 'prd', s, 'reac' )
+ r1.Kf = 0.3 # 1/sec
+ r1.Kb = 0 # 1/sec
+
+ # Assign parameters
+ a.diffConst = diffConst/10
+ b.diffConst = diffConst
+ s.diffConst = 0
+
+ # Make solvers
+ ksolve = moose.Ksolve( '/model/compartment/ksolve' )
+ dsolve = moose.Dsolve( '/model/dsolve' )
+ # Set up clocks. The dsolver to know before assigning stoich
+ moose.setClock( 4, dt4 )
+ moose.setClock( 5, dt5 )
+ moose.useClock( 4, '/model/dsolve', 'process' )
+ # Ksolve must be scheduled after dsolve.
+ moose.useClock( 5, '/model/compartment/ksolve', 'process' )
+
+ stoich = moose.Stoich( '/model/compartment/stoich' )
+ stoich.compartment = compartment
+ stoich.ksolve = ksolve
+ stoich.dsolve = dsolve
+ stoich.path = "/model/compartment/##"
+ assert( dsolve.numPools == 3 )
+ a.vec.concInit = [0.1]*num
+ a.vec[0].concInit *= 1.2 # slight perturbation at one end.
+ b.vec.concInit = [0.1]*num
+ s.vec.concInit = [1]*num
+
+ def displayPlots():
+ a = moose.element( '/model/compartment/a' )
+ b = moose.element( '/model/compartment/b' )
+ pos = numpy.arange( 0, a.vec.conc.size, 1 )
+ pylab.plot( pos, a.vec.conc, label='a' )
+ pylab.plot( pos, b.vec.conc, label='b' )
+ pylab.legend()
+ pylab.show()
+
+ def main():
+ runtime = 400
+ displayInterval = 2
+ makeModel()
+ dsolve = moose.element( '/model/dsolve' )
+ moose.reinit()
+ #moose.start( runtime ) # Run the model for 10 seconds.
+
+ a = moose.element( '/model/compartment/a' )
+ b = moose.element( '/model/compartment/b' )
+ s = moose.element( '/model/compartment/s' )
+
+ img = mpimg.imread( 'turingPatternTut.png' )
+ #imgplot = plt.imshow( img )
+ #plt.show()
+
+ plt.ion()
+ fig = plt.figure( figsize=(12,10) )
+ png = fig.add_subplot(211)
+ imgplot = plt.imshow( img )
+ ax = fig.add_subplot(212)
+ ax.set_ylim( 0, 0.5 )
+ plt.ylabel( 'Conc (mM)' )
+ plt.xlabel( 'Position along cylinder (microns)' )
+ pos = numpy.arange( 0, a.vec.conc.size, 1 )
+ line1, = ax.plot( pos, a.vec.conc, label='a' )
+ line2, = ax.plot( pos, b.vec.conc, label='b' )
+ timeLabel = plt.text(60, 0.4, 'time = 0')
+ plt.legend()
+ fig.canvas.draw()
+
+ for t in range( displayInterval, runtime, displayInterval ):
+ moose.start( displayInterval )
+ line1.set_ydata( a.vec.conc )
+ line2.set_ydata( b.vec.conc )
+ timeLabel.set_text( "time = %d" % t )
+ fig.canvas.draw()
+
+ print( "Hit 'enter' to exit" )
+ raw_input( )
+
+
+
+ # Run the 'main' if this script is executed standalone.
+ if __name__ == '__main__':
+ main()
+
+|
+
+**Output:**
+
+.. image:: ../../../images/turing.png
+
+|
+|
+
+
+Relaxation Oscillator
+=====================
+
+File name: **relaxationOsc.py**
+
+This example illustrates a **Relaxation Oscillator**. This is an
+oscillator built around a switching reaction, which tends to flip into
+one or other state and stay there. The relaxation bit comes in because
+once it is in state 1, a slow (relaxation) process begins which
+eventually flips it into state 2, and vice versa.
+
+.. image:: ../../../images/relaxOsc_tut.png
+
+The model is based on Bhalla, Biophys J. 2011. It is defined in kkit
+format. It uses the deterministic gsl solver by default. You can specify
+the stochastic Gillespie solver on the command line
+
+::
+
+ ``python relaxationOsc.py gssa``
+
+Things to do with the model:
+
+::
+
+ * Figure out what determines its frequency. You could change
+ the initial concentrations of various model entities::
+
+ ma = moose.element( '/model/kinetics/A/M' )
+ ma.concInit *= 1.5
+
+ Alternatively, you could scale the rates of molecular traffic
+ between the compartments::
+
+ exo = moose.element( '/model/kinetics/exo' )
+ endo = moose.element( '/model/kinetics/endo' )
+ exo.Kf *= 1.0
+ endo.Kf *= 1.0
+
+ * Play with stochasticity. The standard thing here is to scale the
+ volume up and down::
+
+ compt.volume = 1e-18
+ compt.volume = 1e-20
+ compt.volume = 1e-21
+
+**Code:**
+
+.. hidden-code-block:: python
+ :linenos:
+ :label: Show/Hide code
+
+ #########################################################################
+ ## This program is part of 'MOOSE', the
+ ## Messaging Object Oriented Simulation Environment.
+ ## Copyright (C) 2014 Upinder S. Bhalla. and NCBS
+ ## It is made available under the terms of the
+ ## GNU Lesser General Public License version 2.1
+ ## See the file COPYING.LIB for the full notice.
+ #########################################################################
+
+ import moose
+ import matplotlib.pyplot as plt
+ import matplotlib.image as mpimg
+ import pylab
+ import numpy
+ import sys
+
+ def main():
+
+ solver = "gsl" # Pick any of gsl, gssa, ee..
+ #solver = "gssa" # Pick any of gsl, gssa, ee..
+ mfile = '../../genesis/OSC_Cspace.g'
+ runtime = 4000.0
+ if ( len( sys.argv ) >= 2 ):
+ solver = sys.argv[1]
+ modelId = moose.loadModel( mfile, 'model', solver )
+ # Increase volume so that the stochastic solver gssa
+ # gives an interesting output
+ compt = moose.element( '/model/kinetics' )
+ compt.volume = 1e-19
+ dt = moose.element( '/clock' ).tickDt[18] # 18 is the plot clock.
+
+ moose.reinit()
+ moose.start( runtime )
+
+ # Display all plots.
+ img = mpimg.imread( 'relaxOsc_tut.png' )
+ fig = plt.figure( figsize=(12, 10 ) )
+ png = fig.add_subplot( 211 )
+ imgplot = plt.imshow( img )
+ ax = fig.add_subplot( 212 )
+ x = moose.wildcardFind( '/model/#graphs/conc#/#' )
+ t = numpy.arange( 0, x[0].vector.size, 1 ) * dt
+ ax.plot( t, x[0].vector, 'b-', label=x[0].name )
+ ax.plot( t, x[1].vector, 'c-', label=x[1].name )
+ ax.plot( t, x[2].vector, 'r-', label=x[2].name )
+ ax.plot( t, x[3].vector, 'm-', label=x[3].name )
+ plt.ylabel( 'Conc (mM)' )
+ plt.xlabel( 'Time (seconds)' )
+ pylab.legend()
+ pylab.show()
+
+ # Run the 'main' if this script is executed standalone.
+ if __name__ == '__main__':
+ main()
+
+|
+
+**Output:**
+
+.. image:: ../../../images/relax.png
+
+
+|
+|
+
+Repressilator
+=============
+
+File name: **repressilator.py**
+
+This example illustrates the classic **Repressilator** model, based on
+Elowitz and Liebler, Nature 2000. The model has the basic architecture
+
+.. image:: ../../../images/repressillatorOsc.png
+
+where **TetR**, **Lac**, and **Lcl** are genes whose products repress
+eachother. The circle symbol indicates inhibition. The model uses the
+Gillespie (stochastic) method by default but you can run it using a
+deterministic method by saying ``python repressillator.py gsl``
+
+Good things to do with this model include:
+
+::
+
+ * Ask what it would take to change period of repressillator:
+
+ * Change inhibitor rates::
+
+ inhib = moose.element( '/model/kinetics/TetR_gene/inhib_reac' )
+ moose.showfields( inhib )
+ inhib.Kf *= 0.1
+
+ * Change degradation rates::
+
+ degrade = moose.element( '/model/kinetics/TetR_gene/TetR_degradation' )
+ degrade.Kf *= 10.0
+ * Run in stochastic mode:
+
+ * Change volumes, figure out how many molecules are present::
+
+ lac = moose.element( '/model/kinetics/lac_gene/lac' )
+ print lac.n``
+
+ * Find when it becomes hopelessly unreliable with small volumes.
+
+**Code:**
+
+.. hidden-code-block:: python
+ :linenos:
+ :label: Show/Hide code
+
+ #########################################################################
+ ## This program is part of 'MOOSE', the
+ ## Messaging Object Oriented Simulation Environment.
+ ## Copyright (C) 2014 Upinder S. Bhalla. and NCBS
+ ## It is made available under the terms of the
+ ## GNU Lesser General Public License version 2.1
+ ## See the file COPYING.LIB for the full notice.
+ #########################################################################
+
+ import moose
+ import matplotlib.pyplot as plt
+ import matplotlib.image as mpimg
+ import pylab
+ import numpy
+ import sys
+
+ def main():
+
+ #solver = "gsl" # Pick any of gsl, gssa, ee..
+ solver = "gssa" # Pick any of gsl, gssa, ee..
+ mfile = '../../genesis/Repressillator.g'
+ runtime = 6000.0
+ if ( len( sys.argv ) >= 2 ):
+ solver = sys.argv[1]
+ modelId = moose.loadModel( mfile, 'model', solver )
+ # Increase volume so that the stochastic solver gssa
+ # gives an interesting output
+ compt = moose.element( '/model/kinetics' )
+ compt.volume = 1e-19
+ dt = moose.element( '/clock' ).tickDt[18]
+
+ moose.reinit()
+ moose.start( runtime )
+
+ # Display all plots.
+ img = mpimg.imread( 'repressillatorOsc.png' )
+ fig = plt.figure( figsize=(12, 10 ) )
+ png = fig.add_subplot( 211 )
+ imgplot = plt.imshow( img )
+ ax = fig.add_subplot( 212 )
+ x = moose.wildcardFind( '/model/#graphs/conc#/#' )
+ plt.ylabel( 'Conc (mM)' )
+ plt.xlabel( 'Time (seconds)' )
+ for x in moose.wildcardFind( '/model/#graphs/conc#/#' ):
+ t = numpy.arange( 0, x.vector.size, 1 ) * dt
+ pylab.plot( t, x.vector, label=x.name )
+ pylab.legend()
+ pylab.show()
+
+ # Run the 'main' if this script is executed standalone.
+ if __name__ == '__main__':
+ main()
+
+|
+
+**Output:**
+
+.. image:: ../../../images/repris.png
+
+
diff --git a/docs/source/user/py/tutorials/Squid.rst b/docs/source/user/py/tutorials/Squid.rst
new file mode 100644
index 0000000000000000000000000000000000000000..60657e5c4b6f2a486a7a9822e9616c29907a0eb6
--- /dev/null
+++ b/docs/source/user/py/tutorials/Squid.rst
@@ -0,0 +1,35 @@
+****************
+Squid giant axon
+****************
+
+This tutorial is an interactive graphical simulation of a squid giant axon,
+closely based on the 'Squid' demo by Mark Nelson which ran in GENESIS.
+
+The `squid giant axon <https://en.wikipedia.org/wiksi/Squid_giant_axon>`_ is a very large axon that plays a role in the water jet propulsion systems of squid.
+
+Alan Hodgkin, Andrew Huxley, and John Eccles won the nobel prize in physiology or medicine for their pioneering work on the squid axon. Hodgin and Huxley were the first to qualitatively describe action potentials within neurons. The large diameter of the squid giant axon (0.5 mm to 1 mm) allowed them to affix electrodes and voltage clamps to precisely measure the action potential as it travelled through the axon. They later went on to mathematically describe this in an equation that paved the road for mathematical and computational biology's development.
+
+This tutorial models the Hodglin-Huxley equation within a neat graphical interface that allows one to change different parameters and see how it affects the resulting action potential.
+
+.. figure:: ../../../images/squid_demo.png
+ :alt: Hodgkin-Huxley's squid giant axon experiment
+
+ The GUI of the simulation.
+
+The tutorial can be run from within the ``.../moose/moose-examples/squid`` directory by running
+::
+
+ python squid_demo.py
+
+in command line from within the directory.
+
+Once the model loads, you can access the inbuilt documentation by clicking on ``Help running`` in the top right of the window as shown below
+
+.. figure:: ../../../images/Help.png
+ :alt: Help running
+
+ The "Help running" tab is highlighted in red
+
+The page that pops up will take you through using the GUI, changing the parameters and understanding the model.
+
+For more details on the biophysics behind the model, you can click on ``Help biophysics`` tab to the immediate right of the ``Help running`` tab (note that for smaller default window sizes, the tab might not be visible and can be accessed by clicking the ``>>`` on the top right corner of the GUI).
diff --git a/docs/source/user/py/tutorials/index_tut.rst b/docs/source/user/py/tutorials/index_tut.rst
new file mode 100644
index 0000000000000000000000000000000000000000..9bcc8a87d49db6125ee1a9f06ad4312e47be4c26
--- /dev/null
+++ b/docs/source/user/py/tutorials/index_tut.rst
@@ -0,0 +1,12 @@
+Teaching Tutorials
+==================
+
+These tutorials use the moose simulation environment to illustrate various scientific concepts, phenomena, and research. As such the materials in this section are useful for both teaching and learning about these topics.
+
+.. toctree::
+ :maxdepth: 2
+
+ ChemicalBistables
+ ChemicalOscillators
+ Squid
+
diff --git a/moose-examples/snippets/compartment_net.py b/moose-examples/snippets/compartment_net.py
index 55b95f546fb4c434bf58cf02e6a4b5e9116a5b8e..b77b25db8de95ab18c4fa24227d2e1a056c5ae22 100644
--- a/moose-examples/snippets/compartment_net.py
+++ b/moose-examples/snippets/compartment_net.py
@@ -177,5 +177,10 @@ def main():
plt.show()
+
+if __name__ == '__main__':
+ main()
+
+
#
# compartment_net.py ends here
diff --git a/moose-examples/snippets/wildcard.py b/moose-examples/snippets/wildcard.py
index ee2aee82c34485fc6be99dcd76cce4603b164cc4..ba2712e45a798ceee0d394e7104bc0cab513c955 100644
--- a/moose-examples/snippets/wildcard.py
+++ b/moose-examples/snippets/wildcard.py
@@ -91,13 +91,15 @@ def wildcard_test():
for element in moose.wildcardFind(wildcard):
print(('\t', element.path))
- # `#` can be used only once and matches all subsequent characters in name
+
+ # `?` can be used any number of times but substitutes a single character
wildcard = '/alfa/bravo/charl?e'
print(('\nElements Matching:', wildcard))
for element in moose.wildcardFind(wildcard):
print(('\t', element.path))
- # `?` can be used any number of times but substitutes a single character
+
+ # `#` can be used only once and matches all subsequent characters in name
wildcard = '/alfa/bravo/fox#'
print(('\nElements Matching:', wildcard))
for element in moose.wildcardFind(wildcard):