.. A cookbook for MOOSE .. Lists all the snippets in Demos/snippets directory MOOSE Cookbook ============== The MOOSE Cookbook contains recipes showing you how to do specific tasks in MOOSE. Loading and running models -------------------------- This section of the documentation explains how to load and run predefined models in MOOSE. Hello, MOOSE: Load, run and display existing models ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. automodule:: helloMoose :members: .. _squid: 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 ``Demos/squid`` subdirectory of MOOSE. If you want to read a simpler implementation of the same (without the code for setting up GUI), check out :ref:`hhmodel`. Start, Stop, and setting clocks ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. automodule:: startstop :members: Run Python from MOOSE ^^^^^^^^^^^^^^^^^^^^^ .. automodule:: pyrun :members: Accessing and tweaking parameters ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. automodule:: tweakingParameters :members: .. figure:: images/tweakingParameters.png :alt: Three oscillation patterns after tweaking model parameters. Storing simulation output ^^^^^^^^^^^^^^^^^^^^^^^^^ Here we'll show how to store and dump from a table and also using HDF5. .. automodule:: tabledemo :members: .. automodule:: hdfdemo :members: Computing arbitrary functions on the fly ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Sometimes you want to calculate arbitrary function of the state variables of one or more elements and feed the result into another element during a simulation. The Function class is useful for this. .. figure:: images/function.png :alt: Outputs of Function object calculating z = c0 * exp(c1 * x) * cos(y) :scale: 50% .. automodule:: function :members: Solving arbitrary differential equations ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ The Function classes are also useful for setting up arbitrary differential equations. MOOSE supports solution of such equation systems, with the caveat that it thinks they should be chemical systems and therefore interprets the variables as concentrations. MOOSE does not permit the concentrations to go negative. .. automodule:: diffEqSolution :members: This limitation on positive solutions can be overcome with offsets to the variables. .. automodule:: funcRateHarmonicOsc :members: It is also possible to set up other forms of differential equations, where instead of directly controlling the rate of change of a pool, the equations take the form of modifying the rate of a reaction involving a pool. .. automodule:: funcReacLotkaVolterra :members: The solutions of such systems are much more accurate and faster using the chemical kinetic solver than with the basic exponential Euler method. However, it is important to point out that the use of arbitrary differential equations to represent chemical systems is discouraged in MOOSE. This is for the simple reason that they are abstractions, frequently with serious flaws, of the underlying chemistry. At the current time the stochastic solver cannot handle systems formulated with functions rather than chemical objects. Chemical Signaling models ------------------------- This section of the documentation explains how to do operations specific to the chemical signaling. Running with different numerical methods ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. automodule:: switchKineticSolvers :members: Changing volumes ^^^^^^^^^^^^^^^^ .. automodule:: scaleVolumes :members: Feeding tabulated input to a model ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. automodule:: analogStimTable :members: Finding steady states ^^^^^^^^^^^^^^^^^^^^^ .. automodule:: findChemSteadyState :members: Making a dose-response curve ^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. figure:: images/chemDoseResponse.png :alt: Dose-response curve example for a bistable system. .. automodule:: chemDoseResponse :members: Building a chemical model from parts ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Disclaimer: Avoid doing this for all but the very simplest models. This is error-prone, tedious, and non-portable. For preference use one of the standard model formats like SBML, which MOOSE and many other tools can read and write. Nevertheless, it is useful to see how these models are set up. There are several tutorials and snippets that build the entire chemical model system using the basic MOOSE calls. The sequence of steps is typically: #. Create container (chemical compartment) for model. This is typically a CubeMesh, a CylMesh, and if you really know what you are doing, a NeuroMesh. #. Create the reaction components: pools of molecules **moose.Pool**; reactions **moose.Reac**; and enzymes **moose.Enz**. Note that when creating an enzyme, one must also create a molecule beneath it to serve as the enzyme-substrate complex. Other less-used components include Michaelis-Menten enzymes **moose.MMenz**, input tables, pulse generators and so on. These are illustrated in other examples. #. Assign parameters for the components. * Compartments have a **volume**, and each subtype will have quite elaborate options for partitioning the compartment into voxels. * **Pool** s have one key parameter, the initial concentration **concInit**. * **Reac** tions have two parameters: **Kf** and **Kb**. * **Enz** ymes have two primary parameters **kcat** and **Km**. That is enough for **MMenz** ymes. Regular **Enz** ymes have an additional parameter **k2** which by default is set to 4.0 times **kcat**, but you may also wish to explicitly assign it if you know its value. #. Connect up the reaction system using moose messaging. #. Create and connect up input and output tables as needed. #. Create and connect up the solvers as needed. This has to be done in a specific order. Examples are linked below, but briefly the order is: a. Make the compartment and reaction system. b. Make the Ksolve or Gsolve. c. Make the Stoich. d. Assign **stoich.compartment** to the compartment e. Assign **stoich.ksolve** to either the Ksolve or Gsolve. f. Assign **stoich.path** to finally fill in the reaction system. There is an additional step if diffusion is also present, see `Reaction-diffusion in a cylinder`_. Some examples of doing this are in: * `Making a dose-response curve`_ , which defines a small bistable system including three **Pool** s, two **Enz** ymes and a **Reac** tion. * `Feeding tabulated input to a model`_, which shows how to connect up a **StimulusTable** object to a simple 2-molecule reaction. * `Reaction-diffusion in a cylinder`_, which defines a simple binding reaction and embeds it in a 1-dimensional diffusive volume of a cylinder. The recommended way to build a chemical model, of course, is to load it in from a file format specific to such models. MOOSE understands **SBML**, **kkit.g** (a legacy GENESIS format), and **cspace** (a very compact format used in a large study of bistables from Ramakrishnan and Bhalla, PLoS Comp. Biol 2008). One key concept is that in MOOSE the components, messaging, and access to model components is identical regardless of whether the model was built from parts, or loaded in from a file. All that the file loaders do is to use the file to automate the steps above. Thus the model components and their fields are completely accessible from the script even if the model has been loaded from a file. See `Accessing and tweaking parameters`_ for an example of this. Oscillation models ^^^^^^^^^^^^^^^^^^ There are several chemical oscillators defined in the ``Demos/tutorials/ChemkcalOscillators`` directory. These include: 1. Slow Feedback Oscillator based on a model by Boris Kholdenko .. automodule:: slowFbOsc :members: 2. Repressilator, based on Elowitz and Liebler, Nature 2000. .. automodule:: repressillator :members: 3. Relaxation oscillator. .. automodule:: relaxationOsc :members: Bistability models ^^^^^^^^^^^^^^^^^^ There are several bistable models defined in the ``Demos/tutorials/ChemkcalBistables`` directory. These include: 1. MAPK feedback loop model. .. automodule:: mapkFB :members: 2. Simple minimal bistable model, run stochastically at different volumes to illustrate the effects of chemical noise. .. automodule:: scaleVolumes :members: 3. Strongly bistable model. .. automodule:: strongBis :members: Reaction-diffusion models ------------------------- The MOOSE design for reaction-diffusion is to specify one or more cellular 'compartments', and embed reaction systems in each of them. A 'compartment', in the context of reaction-diffusion, is used in the cellular sense of a biochemically defined, volume restricted subpart of a cell. Many but not all compartments are bounded by a cell membrane, but biochemically the membrane itself may form a compartment. Note that this interpretation differs from that of a 'compartment' in detailed electrical models of neurons. A reaction system can be loaded in from any of the supported MOOSE formats, or built within Python from MOOSE parts. The computations for such models are done by a set of objects: Stoich, Ksolve and Dsolve. Respectively, these handle the model reactions and stoichiometry matrix, the reaction computations for each voxel, and the diffusion between voxels. The 'Compartment' specifies how the model should be spatially discretized. Reaction-diffusion + transport in a tapering cylinder ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. automodule:: cylinderDiffusion :members: A Turing model ^^^^^^^^^^^^^^ .. automodule:: TuringOneDim :members: A spatial bistable model ^^^^^^^^^^^^^^^^^^^^^^^^ .. automodule:: propagationBis Reaction-diffusion in neurons ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Reaction-diffusion systems can easily be embedded into neuronal geometries. MOOSE does so by treating each neuron as a pseudo 1-dimensional object. This means that diffusion only happens along the axis of dendritic segments, not radially from inside to outside a dendrite, nor tangentially around the dendrite circumference. Here we illustrate two cases. The simple case treats the entire neuron as a single, chemically equivalent reaction-diffusion system in a binary branching neuronal tree. The more complex example shows how to set up three chemically distinct kinds of subdivisions within the neuron: the dendritic tree, the dendritic spine heads, and the postsynaptic densities. In both examples we embed a simple Turing-like spatial oscillator in every compartment of the model neurons, so as to see nice oscillations and animations. The first example has a particularly striking pseudo-3D rendition of the neuron and the molecular spatial oscillations within it. .. figure:: images/reacDiffBranchingNeuron.png :alt: Pseudo-3-D rendition of branching neuron and the concs in it. .. automodule:: reacDiffBranchingNeuron :members: .. automodule:: reacDiffSpinyNeuron :members: Transport in branching dendritic tree ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. automodule:: transportBranchingNeuron :members: Cross-compartment reaction systems ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Frequently reaction systems span cellular (chemical) compartments. For example, a membrane-bound molecule may be phosphorylated by a cytosolic kinase, using soluble ATP as one of the substrates. Here the membrane and the cytsol are different chemical compartments. MOOSE supports such reactions. The following snippets illustrate cross-compartment chemistry. Note that the interpretation of the rates of enzymes and reactions does depend on which compartment they reside in. .. automodule:: crossComptSimpleReac :members: .. automodule:: crossComptOscillator :members: .. automodule:: crossComptNeuroMesh :members: Single neuron models -------------------- Neurons are modelled as equivalent electrical circuits. The morphology of a neuron can be broken into isopotential compartments connected by axial resistances R\ :sub:`a`\ denoting the cytoplasmic resistance. In each compartment, the neuronal membrane is represented as a capacitance C\ :sub:`m`\ with a shunt leak resistance R\ :sub:`m`\ . Electrochemical gradient (due to ion pumps) across the leaky membrane causes a voltage drive E\ :sub:`m`\ , that hyperpolarizes the inside of the cell membrane compared to the outside. Each voltage dependent ion channel, present on the membrane, is modelled as a voltage dependent conductance G\ :sub:`k`\ with gating kinetics, in series with an electrochemical voltage drive (battery) E\ :sub:`k`\ , across the membrane capacitance C\ :sub:`m`\ , as in the figure below. -------------- .. figure:: ../../images/neuroncompartment.png :align: center :alt: **Equivalent circuit of neuronal compartments** **Equivalent circuit of neuronal compartments** -------------- Neurons fire action potentials / spikes (sharp rise and fall of membrane potential V\ :sub:`m`\ ) due to voltage dependent channels. These result in opening of excitatory / inhibitory synaptic channels (conductances with batteries, similar to voltage gated channels) on other connected neurons in the network. MOOSE can handle large networks of detailed neurons, each with complicated channel dynamics. Further, MOOSE can integrate chemical signalling with electrical activity. Presently, creating and simulating these requires PyMOOSE scripting, but these will be incorporated into the GUI in the future. To understand channel kinetics and neuronal action potentials, run the Squid Axon demo installed along with MOOSEGUI and consult its help/tutorial. Read more about compartmental modelling in the first few chapters of the `Book of Genesis `_. Models can be defined in `NeuroML `_, an XML format which is mostly supported across simulators. Channels, neuronal morphology (compartments), and networks can be specified using various levels of NeuroML, namely ChannelML, MorphML and NetworkML. Importing of cell models in the `GENESIS `_ .p format is supported for backwards compatibitility. Modeling details ^^^^^^^^^^^^^^^^^ Some salient properties of neuronal building blocks in MOOSE are described below. Variables that are updated at every simulation time step are are listed **dynamical**. Rest are parameters. - **Compartment** When you select a compartment, you can view and edit its properties in the right pane. V\ :sub:`m`\ and I\ :sub:`m`\ are plot-able. - V\ :sub:`m`\ membrane potential (across C\ :sub:`m`\ ) in Volts. It is a dynamical variable. - C\ :sub:`m`\ membrane capacitance in Farads. - E\ :sub:`m`\ membrane leak potential in Volts due to the electrochemical gradient setup by ion pumps. - I\ :sub:`m`\ current in Amperes across the membrane via leak resistance R\ :sub:`m`\ . - inject current in Amperes injected externally into the compartment. - initVm initial V\ :sub:`m`\ in Volts. - R\ :sub:`m`\ membrane leak resistance in Ohms due to leaky channels. - diameter diameter of the compartment in metres. - length length of the compartment in metres. - **HHChannel** Hodgkin-Huxley channel with voltage dependent dynamical gates. - Gbar peak channel conductance in Siemens. - E\ :sub:`k`\ reversal potential of the channel, due to electrochemical gradient of the ion(s) it allows. - G\ :sub:`k`\ conductance of the channel in Siemens. G\ :sub:`k`\ (t) = Gbar × X(t)\ :sup:`Xpower`\ × Y(t)\ :sup:`Ypower`\ × Z(t)\ :sup:`Zpower`\ - I\ :sub:`k`\ current through the channel into the neuron in Amperes. I\ :sub:`k`\ (t) = G\ :sub:`k`\ (t) × (E\ :sub:`k`\ -V\ :sub:`m`\ (t)) - X, Y, Z gating variables (range 0.0 to 1.0) that may turn on or off as voltage increases with different time constants. dX(t)/dt = X\ :sub:`inf`\ /τ - X(t)/τ Here, X\ :sub:`inf`\ and τ are typically sigmoidal/linear/linear-sigmoidal functions of membrane potential V\ :sub:`m`\ , which are described in a ChannelML file and presently not editable from MOOSEGUI. Thus, a gate may open (X\ :sub:`inf`\ (V\ :sub:`m`\ ) → 1) or close (X\ :sub:`inf`\ (V\ :sub:`m`\ ) → 0) on increasing V\ :sub:`m`\ , with time constant τ(V\ :sub:`m`\ ). - Xpower, Ypower, Zpower powers to which gates are raised in the G\ :sub:`k`\ (t) formula above. - **HHChannel2D** The Hodgkin-Huxley channel2D can have the usual voltage dependent dynamical gates, and also gates that depend on voltage and an ionic concentration, as for say Ca-dependent K conductance. It has the properties of HHChannel above, and a few more, similar to in the `GENESIS tab2Dchannel reference `_. - **CaConc** This is a pool of Ca ions in each compartment, in a shell volume under the cell membrane. The dynamical Ca concentration increases when Ca channels open, and decays back to resting with a specified time constant τ. Its concentration controls Ca-dependent K channels, etc. - Ca Ca concentration in the pool in units mM ( i.e., mol/m\ :sup:`3`\ ). d[Ca\ :sup:`2+`\ ]/dt = B × I\ :sub:`Ca`\ - [Ca\ :sup:`2+`\ ]/τ - CaBasal/Ca_base Base Ca concentration to which the Ca decays - tau time constant with which the Ca concentration decays to the base Ca level. - B constant in the [Ca\ :sup:`2+`\ ] equation above. - thick thickness of the Ca shell within the cell membrane which is used to calculate B (see Chapter 19 of `Book of GENESIS `_.) Neuronal simulations in MOOSEGUI ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Neuronal models in various formats can be loaded and simulated in the **MOOSE Graphical User Interface**. The GUI displays the neurons in 3D, and allows visual selection and editing of neuronal properties. Plotting and visualization of activity proceeds concurrently with the simulation. Support for creating and editing channels, morphology and networks is planned for the future. MOOSEGUI uses SI units throughout. Demos ^^^^^ - **Cerebellar granule cell** **File -> Load ->** ~/moose/Demos/neuroml/GranuleCell/GranuleCell.net.xml This is a single compartment Cerebellar granule cell with a variety of channels `Maex, R. and De Schutter, E., 1997 `_ (exported from http://www.neuroconstruct.org/). Click on its soma, and **See children** for its list of channels. Vary the Gbar of these channels to obtain regular firing, adapting and bursty behaviour (may need to increase tau of the Ca pool). - **Pyloric rhythm generator in the stomatogastric ganglion of lobster** **File -> Load ->** ~/moose/Demos/neuroml/pyloric/Generated.net.xml - **Purkinje cell** **File -> Load ->** ~/moose/Demos/neuroml/PurkinjeCell/Purkinje.net.xml This is a purely passive cell, but with extensive morphology [De Schutter, E. and Bower, J. M., 1994] (exported from http://www.neuroconstruct.org/). The channel specifications are in an obsolete ChannelML format which MOOSE does not support. - **Olfactory bulb subnetwork** **File -> Load ->** ~/moose/Demos/neuroml/OlfactoryBulb/numgloms2_seed100.0_decimated.xml This is a pruned and decimated version of a detailed network model of the Olfactory bulb [Gilra A. and Bhalla U., in preparation] without channels and synaptic connections. We hope to post the ChannelML specifications of the channels and synapses soon. - **All channels cell** **File -> Load ->** ~/moose/Demos/neuroml/allChannelsCell/allChannelsCell.net.xml This is the Cerebellar granule cell as above, but with loads of channels from various cell types (exported from http://www.neuroconstruct.org/). Play around with the channel properties to see what they do. You can also edit the ChannelML files in ~/moose/Demos/neuroml/allChannelsCell/cells_channels/ to experiment further. - **NeuroML python scripts** In directory ~/moose/Demos/neuroml/GranuleCell, you can run python FvsI_Granule98.py which plots firing rate vs injected current for the granule cell. Consult this python script to see how to read in a NeuroML model and to set up simulations. There are ample snippets in ~/moose/Demos/snippets too. Loading, modifying, saving ^^^^^^^^^^^^^^^^^^^^^^^^^^ Explicit vs. implict methods ^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Integrate-and-fire models ^^^^^^^^^^^^^^^^^^^^^^^^^ .. automodule:: IntegrateFireZoo :members: The HH model ^^^^^^^^^^^^ .. _hhmodel: This is a standalone script for simulating the Hodgkin-Huxley squid axon experiment with a step current injection. The graphical version of the same is :ref:`squid`. .. automodule:: ionchannel :members: Analyzing spike trains ^^^^^^^^^^^^^^^^^^^^^^ Network models -------------- Connecting two cells via a synapse ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Below is the connectivity diagram for setting up a synaptic connection from one neuron to another. The PulseGen object is there for stimulating the presynaptic cell as part of experimental setup. The cells are defined as single-compartments with Hodgkin-Huxley type Na+ and K+ channels (see :ref:`hhmodel`) .. figure:: images/twoCells.png :scale: 50% :alt: Two cells connected via synapse .. automodule:: twocells :members: Plastic synapse: STDP ^^^^^^^^^^^^^^^^^^^^^ .. automodule:: STDP :members: Network with Ca-based plasticity ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. automodule:: ExcInhNet_HigginsGraupnerBrunel2014 :members: Providing random input to a cell ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. automodule:: randomspike :members: .. figure:: images/randomSpike.png :scale: 50% :alt: Random spike input to a cell Recurrent integrate-and-fire network ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. automodule:: ExcInhNet_Ostojic2014_Brunel2000 :members: Recurrent integrate-and-fire network with plasticity ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ A feed-forward network with random input ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Using compartmental models in networks ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ **Pyloric rhythm generator in the stomatogastric ganglion of lobster** .. automodule:: STG_net :members: Multiscale models ----------------- Single-compartment multiscale model ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. automodule:: multiscaleOneCompt :members: Multi-compartment multiscale model ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Graphics -------- Displaying time-series plots ^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Animation of values along an axis ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Using MOOGLI widgets to display a neuron ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^