Skip to content
Snippets Groups Projects
Commit fe2d8940 authored by Dilawar Singh's avatar Dilawar Singh
Browse files

Deleted temp files.

parent 88e9ca95
Branches
Tags
No related merge requests found
Showing
with 0 additions and 2433 deletions
# -*- coding: utf-8 -*-
#
# MOOSE documentation build configuration file, created by
# sphinx-quickstart on Thu Feb 9 14:56:13 2017.
#
# This file is execfile()d with the current directory set to its
# containing dir.
#
# Note that not all possible configuration values are present in this
# autogenerated file.
#
# All configuration values have a default; values that are commented out
# serve to show the default.
# If extensions (or modules to document with autodoc) are in another directory,
# add these directories to sys.path here. If the directory is relative to the
# documentation root, use os.path.abspath to make it absolute, like shown here.
#
import os
import sys
import sphinx_rtd_theme
# sys.path.insert(0, os.path.abspath('.'))
sys.path.insert(0, os.path.abspath('../python'))
sys.path.append(os.path.abspath('../../moose-examples/snippets'))
sys.path.append(os.path.abspath('../../moose-examples/tutorials/ChemicalOscillators'))
sys.path.append(os.path.abspath('../../moose-examples/tutorials/ChemicalBistables'))
sys.path.append(os.path.abspath('../../moose-examples/tutorials/ExcInhNet'))
sys.path.append(os.path.abspath('../../moose-examples/neuroml/lobster_pyloric'))
sys.path.append(os.path.abspath('../../moose-examples/tutorials/ExcInhNetCaPlasticity'))
# -- General configuration ------------------------------------------------
# If your documentation needs a minimal Sphinx version, state it here.
#
# needs_sphinx = '1.0'
# Add any Sphinx extension module names here, as strings. They can be
# extensions coming with Sphinx (named 'sphinx.ext.*') or your custom
# ones.
extensions = ['sphinx.ext.autodoc']
# Add any paths that contain templates here, relative to this directory.
templates_path = ['_templates']
# The suffix(es) of source filenames.
# You can specify multiple suffix as a list of string:
#
# source_suffix = ['.rst', '.md']
source_suffix = '.rst'
# The master toctree document.
master_doc = 'index'
# General information about the project.
project = u'MOOSE'
copyright = u'2017, Upinder S. Bhall, Subhasis Ray, Aditya Gilra, Dilawar Singh, Harsha Rani, Aviral Goel, Malav Shah'
author = u'Upinder S. Bhall, Subhasis Ray, Aditya Gilra, Dilawar Singh, Harsha Rani, Aviral Goel, Malav Shah'
# The version info for the project you're documenting, acts as replacement for
# |version| and |release|, also used in various other places throughout the
# built documents.
#
# The short X.Y version.
version = u'3.1'
# The full version, including alpha/beta/rc tags.
release = u'ChamCham (3.1.1)'
# The language for content autogenerated by Sphinx. Refer to documentation
# for a list of supported languages.
#
# This is also used if you do content translation via gettext catalogs.
# Usually you set "language" from the command line for these cases.
language = None
# List of patterns, relative to source directory, that match files and
# directories to ignore when looking for source files.
# This patterns also effect to html_static_path and html_extra_path
exclude_patterns = []
# The name of the Pygments (syntax highlighting) style to use.
pygments_style = 'sphinx'
# If true, `todo` and `todoList` produce output, else they produce nothing.
todo_include_todos = False
# -- Options for HTML output ----------------------------------------------
# The theme to use for HTML and HTML Help pages. See the documentation for
# a list of builtin themes.
#
html_theme = "sphinx_rtd_theme"
html_theme_path = [sphinx_rtd_theme.get_html_theme_path()]
# Theme options are theme-specific and customize the look and feel of a theme
# further. For a list of options available for each theme, see the
# documentation.
#
# html_theme_options = {}
# Add any paths that contain custom static files (such as style sheets) here,
# relative to this directory. They are copied after the builtin static files,
# so a file named "default.css" will overwrite the builtin "default.css".
html_static_path = ['nstatic']
# -- Options for HTMLHelp output ------------------------------------------
html_logo = 'images/moose_logo.png'
# Output file base name for HTML help builder.
htmlhelp_basename = 'MOOSEdoc'
# -- Options for LaTeX output ---------------------------------------------
latex_elements = {
# The paper size ('letterpaper' or 'a4paper').
#
# 'papersize': 'letterpaper',
# The font size ('10pt', '11pt' or '12pt').
#
# 'pointsize': '10pt',
# Additional stuff for the LaTeX preamble.
#
# 'preamble': '',
# Latex figure (float) alignment
#
# 'figure_align': 'htbp',
}
# Grouping the document tree into LaTeX files. List of tuples
# (source start file, target name, title,
# author, documentclass [howto, manual, or own class]).
latex_documents = [
(master_doc, 'MOOSE.tex', u'MOOSE Documentation',
u'Upinder S. Bhall, Subhasis Ray, Aditya Gilra, Dilawar Singh, Harsha Rani, Aviral Goel, Malav Shah', 'manual'),
]
# -- Options for manual page output ---------------------------------------
# One entry per manual page. List of tuples
# (source start file, name, description, authors, manual section).
man_pages = [
(master_doc, 'moose', u'MOOSE Documentation',
[author], 1)
]
# -- Options for Texinfo output -------------------------------------------
# Grouping the document tree into Texinfo files. List of tuples
# (source start file, target name, title, author,
# dir menu entry, description, category)
texinfo_documents = [
(master_doc, 'MOOSE', u'MOOSE Documentation',
author, 'MOOSE', 'One line description of project.',
'Miscellaneous'),
]
****************
Chemical Aspects
****************
.. contents::
:local:
:depth: 1
.. _ckbk-chem-load-run:
Load - Run - Save models
========================
Load a Kinetic Model
--------------------
.. automodule:: loadKineticModel
:members:
Load an SBML Model
-------------------
.. automodule:: loadSbmlmodel
:members:
Load a CSpace Model
-------------------
.. automodule:: loadCspaceModel
:members:
Save a model into SBML format
-----------------------------
.. automodule:: convert_Genesis2Sbml
:members:
Save a model
------------
.. automodule:: savemodel
:members:
.. _ckbk-chem-snippets:
Simple Examples
===============================
Set-up a kinetic solver and model
---------------------------------
with Scripting
^^^^^^^^^^^^^^
.. automodule:: scriptGssaSolver
:members:
With something else
^^^^^^^^^^^^^^^^^^^
.. automodule:: changeFuncExpression
:members:
Building a chemical Model from Parts
------------------------------------
.. automodule:: scriptKineticSolver
:members:
Cross-Compartment Reaction Systems
----------------------------------
.. automodule:: crossComptSimpleReac
:members:
.. automodule:: crossComptOscillator
:members:
.. automodule:: crossComptNeuroMesh
:members:
.. automodule:: crossComptSimpleReacGSSA
:members:
Tweaking Parameters
-------------------
.. automodule:: tweakingParameters
:members:
Models' Demonstration
---------------------
Oscillation Model
^^^^^^^^^^^^^^^^^
.. automodule:: slowFbOsc
:members:
Bistability Models
^^^^^^^^^^^^^^^^^^
MAPK feedback loop model
~~~~~~~~~~~~~~~~~~~~~~~~
.. automodule:: mapkFB
:members:
Simple minimal bistable model
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
.. automodule:: scaleVolumes
:members:
Strongly bistable Model
~~~~~~~~~~~~~~~~~~~~~~~
.. automodule:: strongBis
:members:
Model of bidirectional synaptic plasticity
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[showing bistable chemical switch]
.. automodule:: bidirectionalPlasticity
:members:
Reaction Diffusion Models
-------------------------
[Reaction-diffusion + transport in a tapering cylinder]
.. automodule:: cylinderDiffusion
:members:
.. automodule:: cylinderMotor
:members:
.. automodule:: gssaCylinderDiffusion
:members:
Neuronal Diffusion Reaction
^^^^^^^^^^^^^^^^^^^^^^^^^^^
.. automodule:: rxdFuncDiffusion
:members:
.. automodule:: rxdReacDiffusion
:members:
.. automodule:: rxdFuncDiffusionStoch
:members:
A Turing Model
--------------
.. automodule:: TuringOneDim
:members:
A Spatial Bistable Model
------------------------
.. automodule::
:members:
Reaction Diffusion in Neurons
-----------------------------
.. automodule:: reacDiffConcGradient
:members:
.. automodule:: reacDiffBranchingNeuron
:members:
.. automodule:: reacDiffSpinyNeuron
:members:
.. automodule:: diffSpinyNeuron
:members:
Manipulating Chemical Models
----------------------------
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
^^^^^^^^^^^^^^^^^^^^^^^^^^^^
.. automodule:: chemDoseResponse
:members:
Transport in branching dendritic tree
-------------------------------------
.. automodule:: transportBranchingNeuron
:members:
.. _ckbk-chem-gui:
Interface for chemical kinetic models in MOOSEGUI
=================================================
Upinder Bhalla, Harsha Rani
Nov 8 2016.
--------------
- `Introduction <#introduction>`__
- `**TODO** What are chemical kinetic
models? <#todo-what-are-chemical-kinetic-models>`__
- `Levels of model <#levels-of-model>`__
- `Numerical methods <#numerical-methods>`__
- `Using Kinetikit 12 <#using-kinetikit-12>`__
- `Overview <#overview>`__
- `Model layout and icons <#model-layout-and-icons>`__
- `Compartment <#compartment>`__
- `Pool <#pool>`__
- `Buffered pools <#buffered-pools>`__
- `Reaction <#reaction>`__
- `Mass-action enzymes <#mass-action-enzymes>`__
- `Michaelis-Menten Enzymes <#michaelis-menten-enzymes>`__
- `Summation <#summation>`__
- `Model operations <#model-operations>`__
- `Model Building <#model-building>`__
`Introduction <#TOC>`__
-----------------------
Kinetikit 12 is a graphical interface for doing chemical kinetic modeling in MOOSE. It is derived in part from Kinetikit, which was the
graphical interface used in GENESIS for similar models. Kinetikit, also known as kkit, was at version 11 with GENESIS. Here we start with
Kinetikit 12.
`**TODO** What are chemical kinetic models? <#TOC>`__
-----------------------------------------------------
Much of neuronal computation occurs through chemical signaling. For
example, many forms of synaptic plasticity begin with calcium influx
into the synapse, followed by calcium binding to calmodulin, and then
calmodulin activation of numerous enzymes. These events can be
represented in chemical terms:
4 Ca2+ + CaM <===> Ca4.CaM
Such chemical equations can be modeled through standard Ordinary
Differential Equations, if we ignore space::
d[Ca]/dt = −4Kf ∗ [Ca]4 ∗ [CaM] + 4Kb ∗ [Ca4.CaM] d[CaM]/dt = −Kf ∗ [Ca]4 ∗ [CaM] + Kb ∗ [Ca4.CaM] d[Ca4.CaM]/dt = Kf ∗ [Ca]4 ∗ [CaM] − Kb ∗ [Ca4.CaM]
MOOSE models these chemical systems. This help document describes how to do such modelling using the graphical interface, Kinetikit 12.
`Levels of model <#TOC>`__
--------------------------
Chemical kinetic models can be simple well-stirred (or point) models, or
they could have multiple interacting compartments, or they could include
space explicitly using reaction-diffusion. In addition such models could
be solved either deterministically, or using a stochastic formulation.
At present Kinetikit handles compartmental models but does not compute
diffusion within the compartments, though MOOSE itself can do this at
the script level. Kkit12 will do deterministic as well as stochastic
chemical calculations.
`Numerical methods <#TOC>`__
----------------------------
- **Deterministic**: Adaptive timestep 5th order Runge-Kutta-Fehlberg from the GSL (GNU Scientific Library).
- **Stochastic**: Optimized Gillespie Stochastic Systems Algorithm, custom implementation.
`Using Kinetikit 12 <#TOC>`__
-----------------------------
`Overview <#TOC>`__
^^^^^^^^^^^^^^^^^^^
- Load models using **'File -> Load model'**. A reaction schematic for the chemical system appears in the **'Editor view'** tab.
- From **'Editor view'** tab
- View parameters by clicking on icons, and looking at entries in **'Properties'** table to the right.
- Edit parameters by changing their values in the **'Properties'** table.
- From **'Run View'**
- Pools can be plotted by clicking on their icons and dragging the icons onto the plot Window. Presently only concentration v/s time is plottable.
- Select simulation, diffusin dt's along updateInterval for plot and Gui with numerical method using options under **'Preferences'** button in simulation control.
- Run model using **'Run'** button.
- Save plots image using the icons at the top of the **'Plot Window'** or right click on plot to Export to csv.
Most of these operations are detailed in other sections, and are shared
with other aspects of the MOOSE simulation interface. Here we focus on
the Kinetikit-specific items.
`Model layout and icons <#TOC>`__
---------------------------------
When you are in the **'Editor View'** tab you will see a collection of
icons, arrows, and grey boxes surrounding these. This is a schematic of
the reaction scheme being modeled. You can view and change parameters,
and change the layout of the model.
.. figure:: ../../../images/Moose1.png
:alt:
Resizing the model layout and icons:
- **Zoom**:Â Â Comma and period keys. Alternatively, the mouse scroll wheel or vertical scroll line on the track pad will cause the display to zoom in and out.
- **Pan**:Â Â The arrow keys move the display left, right, up, and down.
- **Entire Model View**:Â Â Pressing the **'a'** key will fit the entire model into the entire field of view.
- **Resize Icons**:Â Â Angle bracket keys, that is, **'<'** and **'>'** or **'+'** and **'-'**. This resizes the icons while leaving their positions on the screen layout more or less the same.
- **Original Model View**:Â Â Pressing the **'A'** key (capital 'A') will revert to the original model view including the original icon scaling.
`Compartment <#TOC>`__
^^^^^^^^^^^^^^^^^^^^^^
The *compartment* in moose is usually a contiguous domain in which a
certain set of chemical reactions and molecular species occur. The
definition is very closely related to that of a cell-biological
compartment. Examples include the extracellular space, the cell
membrane, the cytosol, and the nucleus. Compartments can be nested, but
of course you cannot put a bigger compartment into a smaller one.
- **Icon**: Grey boundary around a set of reactions.
- **Moving Compartments**: Click and drag on the boundary.
- **Resizing Compartment boundary**: Happens automatically when contents are repositioned, so that the boundary just contains contents.
- **Compartment editable parameters**:
- **'name'**: The name of the compartment.
- **'size'**: This is the volume, surface area or length of the compartment, depending on its type.
- **Compartment fixed parameters**:
- **'numDimensions'**: This specifies whether the compartment is a volume, a 2-D surface, or if it is just being represented as a length.
`Pool <#TOC>`__
^^^^^^^^^^^^^^^
This is the set of molecules of a given species within a compartment.
Different chemical states of the same molecule are in different pools.
- **Icon**: |image0| Colored rectangle with pool name in it.
- **Moving pools**: Click and drag.
- **Pool editable parameters**:
- **name**: Name of the pool
- **n**: Number of molecules in the pool
- **nInit**: Initial number of molecules in the pool. 'n' gets set
to this value when the 'reinit' operation is done.
- **conc**: Concentration of the molecules in the pool.
``conc = n * unit_scale_factor / (N<sub>A</sub> * vol)``
- **concInit**: Initial concentration of the molecules in the pool.
'conc' is set to this value when the 'reinit' operation is done.
``concInit = nInit * unit_scale_factor / (N<sub>A</sub> * vol)``
- **Pool fixed parameters**
- **size**: Derived from the compartment that holds the pool.
Specifies volume, surface area or length of the holding
compartment.
`Buffered pools <#TOC>`__
^^^^^^^^^^^^^^^^^^^^^^^^^
Some pools are set to a fixed 'n', that is number of molecules, and
therefore a fixed concentration, throughout a simulation. These are
buffered pools.
- **Icon**: |image1| Colored rectangle with pool name in it.
- **Moving Buffered pools**: Click and drag.
- **Buffered Pool editable parameters**
- **name**: Name of the pool
- **nInit**: Fixed number of molecules in the pool. 'n' gets set to
this value throughout the run.
- **concInit**: Fixed concentration of the molecules in the pool.
'conc' is set to this value throughout the run.
``concInit = nInit * unit_scale_factor / (N<sub>A</sub> * vol)``
- **Pool fixed parameters**:
- **n**: Number of molecules in the pool. Derived from 'nInit'.
- **conc**: Concentration of molecules in the pool. Derived from
'concInit'.
- **size**: Derived from the compartment that holds the pool.
Specifies volume, surface area or length of the holding
compar'tment.
`Reaction <#TOC>`__
^^^^^^^^^^^^^^^^^^^
These are conversion reactions between sets of pools. They are
reversible, but you can set either of the rates to zero to get
irreversibility. In the illustration below, **'D'** and **'A'** are
substrates, and **'B'** is the product of the reaction. This is
indicated by the direction of the green arrow.
.. figure:: ../../../images/KkitReaction.png
:alt:
- **Icon**: |image2| Reversible reaction arrow.
- **Moving Reactions**: Click and drag.
- **Reaction editable parameters**:
- **Name** : Name of reaction
- **K\ :sub:`f`** : 'Forward rate' of reaction, in
'concentration/time' units. This is the normal way to express and
manipulate the reaction rate.
- **k\ :sub:`f`** : Forward rate of reaction, in 'number/time'
units. This is used internally for computations, but is
volume-dependent and should not be used to manipulate the reaction
rate unless you really know what you are doing.
- **K\ :sub:`b`** : Backward rate' of reaction, in
'concentration/time' units. This is the normal way to express and
manipulate the reaction rate.
- **k\ :sub:`b`** : Backward rate of reaction, in 'number/time'
units. This is used internally for computations, but is
volume-dependent and should not be used to manipulate the reaction
rate unless you really know what you are doing.
- **Reaction fixed parameters**:
- **numSubstrates**: Number of substrates molecules.
- **numProducts**: Number of product molecules.
`Mass-action enzymes <#TOC>`__
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
These are enzymes that model the chemical equation's
E + S <===> E.S -> E + P
Note that the second reaction is irreversible. Note also that
mass-action enzymes include a pool to represent the **'E.S'**
(enzyme-substrate) complex. In the example below, the enzyme pool is
named **'MassActionEnz'**, the substrate is **'C'**, and the product is
**'E'**. The direction of the enzyme reaction is indicated by the red
arrows.
.. figure:: ../../../images/MassActionEnzReac.png
:alt:
- **Icon**: |image3| Colored ellipse atop a small square. The ellipse represents the enzyme. The small square represents **'E.S'**, the enzyme-substrate complex. The ellipse icon has the same color as the enzyme pool **'E'**. It is connected to the enzyme pool **'E'** with a straight line of the same color.
The ellipse icon sits on a continuous, typically curved arrow in red,
from the substrate to the product.
A given enzyme pool can have any number of enzyme activities, since
the same enzyme might catalyze many reactions.
- **Moving Enzymes**: Click and drag on the ellipse.
- **Enzyme editable parameters**
- **name** : Name of enzyme.
- **K\ :sub:`m`** : Michaelis-Menten value for enzyme, in
'concentration' units.
- **k\ :sub:`cat`** : Production rate of enzyme, in '1/time' units.
Equal to k\ :sub:`3`, the rate of the second, irreversible
reaction.
- **k1** : Forward rate of the **E+S** reaction, in number and
'1/time' units. This is what is used in the internal calculations.
- **k2** : Backward rate of the **E+S** reaction, in '1/time' units.
Used in internal calculations.
- **k3** : Forward rate of the **E.S -> E + P** reaction, in
'1/time' units. Equivalent to k\ :sub:`cat`. Used in internal
calculations.
- **ratio** : This is equal to k\ :sub:`2`/k:sub:`3`. Needed to
define the internal rates in terms of K\ :sub:`m` and
k\ :sub:`cat`. I usually use a value of 4.
- **Enzyme-substrate-complex editable parameters**: These are identica'l to those of any other pool.
- **name**: Name of the **E.S** complex. Defaults to \*\*\_cplx\*\*.
- **n**: Number of molecules in the pool
- **nInit**: Initial number of molecules in the complex. 'n' gets set to this value when the 'reinit' operation is done.
- **conc**: Concentration of the molecules in the pool.
``conc = n * unit_scale_factor / (N<sub>A</sub> * vol)``
- **concInit**: Initial concentration of the molecules in the pool.
'conc' is set to this value when the 'reinit' operation is done.
``concI'nit = nInit * unit_scale_factor / (N<sub>A</sub> * vol)``
- **Enzyme-substrate-complex fixed parameters**:
- **size**: Derived from the compartment that holds the pool.
Specifies volume, surface area or length of the holding
compartment. Note that the Enzyme-substrate-complex is assumed to
be in the same compartment as the enzyme molecule.
`Michaelis-Menten Enzymes <#TOC>`__
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
These are enzymes that obey the Michaelis-Menten equation
``V = V<sub>max</sub> * [S] / ( K<sub>m</sub> + [S] ) = k<sub>cat</sub> * [Etot] * [S] / ( K<sub>m</sub> + [S] )``
where
- V\ :sub:`max` is the maximum rate of the enzyme
- [Etot] is the total amount of the enzyme
- K\ :sub:`m` is the Michaelis-Menten constant
- S is the substrate.
Nominally these enzymes model the same chemical equation as the mass-action enzyme':
``E + S <===> E.S -> E + P``
but they make the assumption that the **E.S** is in a quasi-steady-state
with **E** and **S**, and they also ignore sequestration of the enzyme
into the complex. So there is no representation of the **E.S** complex.
In the example below, the enzyme pool is named **MM\_Enz**, the
substrate is **E**, and the product is **P**. The direction of the
enzyme reaction is indicated by the red arrows.
.. figure:: ../../../images/MM_EnzReac.png
:alt:
- **Icon**: |image4| Colored ellipse. The ellipse represents the enzyme
The ellipse icon has the same color as the enzyme **'MM\_Enz'**. It
is connected to the enzyme pool **'MM\_Enz'** with a straight line of
the same color. The ellipse icon sits on a continuous, typically
curved arrow in red, from the substrate to the product. A given
enzyme pool can have any number of enzyme activities, since the same
enzyme might catalyze many reactions.
- **Moving Enzymes**: Click and drag.
- **Enzyme editable parameters**:
- **name**: Name of enzyme.
- K\ :sub:`m`: Michaelis-Menten value for enzyme, in 'concentration'
units.
- k\ :sub:`cat`: Production rate of enzyme, in '1/time' units. Equal to k\ :sub:`3`, the rate of the second, irreversible reaction.
`Summation <#TOC>`__
^^^^^^^^^^^^^^^^^^^^
Summation object can be used to add specified variable values. The
variables can be input from pool object.
- **Icon**: This is **Σ** in the example image below. The input pools
**'A'** and **'B'** connect to the **Σ** with blue arrows. The
function ouput's to BuffPool |image5|
`Model operations <#TOC>`__
----------------------------
- **Loading models**: **'File -> Load Model -> select from dialog'**.
This operation makes the previously loaded model disable and loads newly selected models in **'Model View''**
- **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:
- 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.
- 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.
- Exponential Euler:This methods computes the solution of partial
and ordinary differential equations.
`Model building <#TOC>`__
-------------------------
- The **Edit Widget** includes various menu options and model icons on
the top. Use the mouse buttton to click and drag icons from toolbar
to Edit Widget, two things will happen, **icon** will appear in the
editor widget and a **object editor** will pop up with lots of
parameters with respect to moose object.
.. figure:: ../../../images/chemical_CS.png
:alt:
**Rules**:
::
* Compartment has to be created firstly(At present only single compartment model is allowed)
* Enzyme should be dropped on a pool as parent
* function should be dropped on buffPool for output
**Note**:
::
* Drag in pool's and reaction on to the editor widget, now one can set up a reaction.
* Click on mooseObject one can find a little arrow on the top right corner of the object, drag from this little arrow to any object for connection. <br>E.g pool to reaction and reaction to pool. Specific connection type gets specific colored arrow. E.g. Green color arrow for specifying connection between reactant and product for reaction.
* Clicking on the object one can rearrange object for clean layout.
* Second order reaction can also be done by repeating the connection over again
* Each connection can be deleted and using rubberband selection each moose object can be deleted
- From **run widget**, pools are draggable to plot window for plotting.
(Currently **conc** is plotted as default field) Plots are
color-coded as per in model.
.. figure:: ../../../images/Chemical_run.png
:alt:
- Model can be run by clicking **start** button. One can stop button in
mid-stream and start up again without affectiong the calculations.
The reset button clears the simulation.
.. |image0| image:: ../../../images/Pool.png
.. |image1| image:: ../../../images/BufPool.png
.. |image2| image:: ../../../images/KkitReacIcon.png
.. |image3| image:: ../../../images/MassActionEnzIcon.png
.. |image4| image:: ../../../images/MM_EnzIcon.png
.. |image5| image:: ../../../images/func.png
.. _ckbk-chem-rde:
RDesigneur
==========
.. _ckbk-chem-tut:
Tutorials
=========
Deterministic Simulation
------------------------
Stochastic Simulation
---------------------
Finding Steady State (CSpace)
-----------------------------
.. automodule:: cspaceSteadyState
:members:
Building Simple Reaction Model
------------------------------
Define a kinetic model using the scripting
------------------------------------------
.. automodule:: scriptKineticModel
:members:
.. _ckbk-chem-adv:
Advanced Examples
=================
# analogStimTable.py ---
#
# Filename: analogStimTable.py
# Description:
# Author: Upi Bhalla
#
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public License as
# published by the Free Software Foundation; either version 3, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; see the file COPYING. If not, write to
# the Free Software Foundation, Inc., 51 Franklin Street, Fifth
# Floor, Boston, MA 02110-1301, USA.
#
import numpy as np
import pylab
import moose
from moose import utils
def analogStimTable():
simtime = 150
simdt = 0.1
model = moose.Neutral('/model')
data = moose.Neutral('/data')
# This is the stimulus generator
stimtable = moose.StimulusTable('/model/stim')
a = moose.BufPool( '/model/a' )
b = moose.Pool( '/model/b' )
reac = moose.Reac( '/model/reac' )
reac.Kf = 0.1
reac.Kb = 0.1
moose.connect( stimtable, 'output', a, 'setConcInit' )
moose.connect( reac, 'sub', a, 'reac' )
moose.connect( reac, 'prd', b, 'reac' )
aPlot = moose.Table('/data/aPlot')
moose.connect(aPlot, 'requestOut', a, 'getConc')
bPlot = moose.Table('/data/bPlot')
moose.connect(bPlot, 'requestOut', b, 'getConc')
moose.setClock( stimtable.tick, simdt )
moose.setClock( a.tick, simdt )
moose.setClock( aPlot.tick, simdt )
####################################################
# Here we set up the stimulus table. It is half a sine-wave.
stim = [ np.sin(0.01 * float(i) ) for i in range( 314 )]
stimtable.vector = stim
stimtable.stepSize = 0 # This forces use of current time as x value
# The table will interpolate its contents over the time start to stop:
# At values less than startTime, it emits the first value in table
stimtable.startTime = 5
# At values more than stopTime, it emits the last value in table
stimtable.stopTime = 60
stimtable.doLoop = 1 # Enable repeat playbacks.
stimtable.loopTime = 10 + simtime / 2.0 # Repeat playback over this time
moose.reinit()
moose.start(simtime)
t = [ x * simdt for x in range( len( aPlot.vector ) )]
pylab.plot( t, aPlot.vector, label='Stimulus waveform' )
pylab.plot( t, bPlot.vector, label='Reaction product b' )
pylab.legend()
pylab.show()
def main():
"""
Example of using a StimulusTable as an analog signal source in
a reaction system. It could be used similarly to give other
analog inputs to a model, such as a current or voltage clamp signal.
This demo creates a StimulusTable and assigns it half a sine wave.
Then we assign the start time and period over which to emit the wave.
The output of the StimTable is sent to a pool **a**, which participates
in a trivial reaction::
table ----> a <===> b
The output of **a** and **b** are recorded in a regular table
for plotting.
"""
analogStimTable()
if __name__ == '__main__':
main
#
# stimtable.py ends here
import math
import pylab
import numpy
import matplotlib.pyplot as plt
import sys
import moose
print(( '[INFO] Using moose from %s' % moose.__file__ ))
def makeCompt( name, parent, dx, dy, dia ):
RM = 1.0
RA = 1.0
CM = 0.01
EM = -0.065
pax = 0
pay = 0
if ( parent.className == "Compartment" ):
pax = parent.x
pay = parent.y
compt = moose.Compartment( name )
compt.x0 = pax
compt.y0 = pay
compt.z0 = 0
compt.x = pax + dx
compt.y = pay + dy
compt.z = 0
compt.diameter = dia
clen = numpy.sqrt( dx * dx + dy * dy )
compt.length = clen
compt.Rm = RM / (numpy.pi * dia * clen)
compt.Ra = RA * 4.0 * numpy.pi * clen / ( dia * dia )
compt.Cm = CM * numpy.pi * dia * clen
if ( parent.className == "Compartment" ):
moose.connect( parent, 'raxial', compt, 'axial' )
return compt
def makeNeuron( numSeg ):
segmentLength = 1e-6
segmentDia = 1e-6
shaftLength = 1e-6
shaftDia = 0.2e-6
headLength = 0.5e-6
headDia = 0.5e-6
cell = moose.Neutral( '/model/cell' )
model = moose.element( '/model' )
prev = makeCompt( '/model/cell/soma',
model, 0.0, segmentLength, segmentDia )
dend = prev
for i in range( 0, numSeg ):
name = '/model/cell/dend' + str( i )
dend = makeCompt( name, dend, 0.0, segmentLength, segmentDia )
name = '/model/cell/shaft' + str( i )
shaft = makeCompt( name, dend, shaftLength, 0.0, shaftDia )
name = '/model/cell/head' + str( i )
head = makeCompt( name, shaft, headLength, 0.0, headDia )
return cell
def makeModel():
numSeg = 5
diffConst = 0.0
# create container for model
model = moose.Neutral( 'model' )
compt0 = moose.NeuroMesh( '/model/compt0' )
compt0.separateSpines = 1
compt0.geometryPolicy = 'cylinder'
compt1 = moose.SpineMesh( '/model/compt1' )
moose.connect( compt0, 'spineListOut', compt1, 'spineList', 'OneToOne' )
compt2 = moose.PsdMesh( '/model/compt2' )
moose.connect( compt0, 'psdListOut', compt2, 'psdList', 'OneToOne' )
# create molecules and reactions
a = moose.Pool( '/model/compt0/a' )
b = moose.Pool( '/model/compt1/b' )
c = moose.Pool( '/model/compt2/c' )
reac0 = moose.Reac( '/model/compt0/reac0' )
reac1 = moose.Reac( '/model/compt1/reac1' )
# connect them up for reactions
moose.connect( reac0, 'sub', a, 'reac' )
moose.connect( reac0, 'prd', b, 'reac' )
moose.connect( reac1, 'sub', b, 'reac' )
moose.connect( reac1, 'prd', c, 'reac' )
# Assign parameters
a.diffConst = diffConst
b.diffConst = diffConst
c.diffConst = diffConst
a.concInit = 1
b.concInit = 12.1
c.concInit = 1
reac0.Kf = 1
reac0.Kb = 1
reac1.Kf = 1
reac1.Kb = 1
# Create a 'neuron' with a dozen spiny compartments.
elec = makeNeuron( numSeg )
# assign geometry to mesh
compt0.diffLength = 10e-6
#compt0.cell = elec
compt0.subTreePath = elec.path + "/##"
# Build the solvers. No need for diffusion in this version.
ksolve0 = moose.Ksolve( '/model/compt0/ksolve' )
ksolve1 = moose.Ksolve( '/model/compt1/ksolve' )
ksolve2 = moose.Ksolve( '/model/compt2/ksolve' )
stoich0 = moose.Stoich( '/model/compt0/stoich' )
stoich1 = moose.Stoich( '/model/compt1/stoich' )
stoich2 = moose.Stoich( '/model/compt2/stoich' )
# Configure solvers
stoich0.compartment = compt0
stoich1.compartment = compt1
stoich2.compartment = compt2
stoich0.ksolve = ksolve0
stoich1.ksolve = ksolve1
stoich2.ksolve = ksolve2
stoich0.path = '/model/compt0/#'
stoich1.path = '/model/compt1/#'
stoich2.path = '/model/compt2/#'
assert( stoich0.numVarPools == 1 )
assert( stoich0.numProxyPools == 1 )
assert( stoich0.numRates == 1 )
assert( stoich1.numVarPools == 1 )
assert( stoich1.numProxyPools == 1 )
assert( stoich1.numRates == 1 )
assert( stoich2.numVarPools == 1 )
assert( stoich2.numProxyPools == 0 )
assert( stoich2.numRates == 0 )
stoich0.buildXreacs( stoich1 )
stoich1.buildXreacs( stoich2 )
stoich0.filterXreacs()
stoich1.filterXreacs()
stoich2.filterXreacs()
print((a.vec.volume, b.vec.volume, c.vec.volume))
a.vec.concInit = list(range( numSeg + 1, 0, -1))
b.vec.concInit = [5.0 * ( 1 + x ) for x in range( numSeg )]
c.vec.concInit = list(range( 1, numSeg + 1))
print((a.vec.concInit, b.vec.concInit, c.vec.concInit))
# 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' )
# connect up the tables
a1 = moose.element( '/model/compt0/a[' + str( numSeg )+ ']')
b1 = moose.element( '/model/compt1/b[' +str(numSeg - 1)+']')
c1 = moose.element( '/model/compt2/c[' +str(numSeg - 1)+']')
moose.connect( outputA, 'requestOut', a1, 'getConc' );
moose.connect( outputB, 'requestOut', b1, 'getConc' );
moose.connect( outputC, 'requestOut', c1, 'getConc' );
def main():
"""
This example illustrates how to define a kinetic model embedded in
a NeuroMesh, and undergoing cross-compartment reactions. It is
completely self-contained and does not use any external model definition
files. Normally one uses standard model formats like
SBML or kkit to concisely define kinetic and neuronal models.
This example creates a simple reaction::
a <==> b <==> c
in which
**a, b**, and **c** are in the dendrite, spine head, and PSD
respectively.
The model is set up to run using the Ksolve for integration. Although
a diffusion solver is set up, the diff consts here are set to zero.
The display has two parts:
Above is a line plot of concentration against compartment#.
Below is a time-series plot that appears after # the simulation has
ended. The plot is for the last (rightmost) compartment.
Concs of **a**, **b**, **c** are plotted for both graphs.
"""
simdt = 0.01
plotdt = 0.01
makeModel()
# Schedule the whole lot
for i in range( 10, 19):
moose.setClock( i, simdt ) # for the compute objects
moose.reinit()
display()
quit()
def display():
dt = 0.01
runtime = 1
a = moose.element( '/model/compt0/a' )
b = moose.element( '/model/compt1/b' )
c = moose.element( '/model/compt2/c' )
plt.ion()
fig = plt.figure( figsize=(12,10))
timeseries = fig.add_subplot( 212 )
spatial = fig.add_subplot( 211)
spatial.set_ylim(0, 15)
pos = numpy.arange( 0, a.vec.conc.size, 1 )
line1, = spatial.plot( pos, a.vec.conc, 'b-', label='a' )
line2, = spatial.plot( pos[1:], b.vec.conc, 'g-', label='b' )
line3, = spatial.plot( pos[1:], c.vec.conc, 'r-', label='c' )
timeLabel = plt.text( 3, 12, 'time = 0' )
plt.legend()
fig.canvas.draw()
for t in numpy.arange( dt, runtime, dt ):
line1.set_ydata( a.vec.conc )
line2.set_ydata( b.vec.conc )
line3.set_ydata( c.vec.conc )
#print b.vec.volume
#print a.vec.n, b.vec.n, c.vec.n
timeLabel.set_text( "time = %f" % t )
fig.canvas.draw()
#raw_input()
moose.start( dt )
timeseries.set_ylim( 0, 20 )
for x in moose.wildcardFind( '/model/graphs/conc#' ):
t = numpy.arange( 0, x.vector.size *dt , dt ) # sec
line4, = timeseries.plot( t, x.vector, label=x.name )
plt.legend()
fig.canvas.draw()
outfile = '%s.png' % sys.argv[0]
# print( "Hit 'enter' to exit" )
# raw_input()
plt.savefig( outfile )
print(('[INFO] Results are saved to %s' % outfile ))
# Run the 'main' if this script is executed standalone.
if __name__ == '__main__':
main()
#########################################################################
## 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 displayPlots():
for x in moose.wildcardFind( '/model/graphs/#' ):
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 ):
state.randomInit()
moose.start( 0.1 ) # Run the model for 2 seconds.
state.settle()
'''
scale = 1.0 / ( 1e-15 * 6.022e23 )
for x in ksolve.nVec[0]:
print x * scale,
# print ksolve.nVec[0]
print state.nIter, state.status, state.stateType, state.nNegEigenvalues, state.nPosEigenvalues, state.solutionStatus
'''
moose.start( 20.0 ) # Run model for 10 seconds, just for display
def main():
"""
This example sets up the kinetic solver and steady-state finder, on
a bistable model.
It looks for the fixed points 100 times, 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 100
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. 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!
You may wish to sample concentration space logarithmically rather than
linearly.
"""
# The wildcard uses # for single level, and ## for recursive.
#compartment = makeModel()
moose.loadModel( '../genesis/M1719.cspace', '/model', 'ee' )
compartment = moose.element( 'model/kinetics' )
compartment.name = 'compartment'
ksolve = moose.Ksolve( '/model/compartment/ksolve' )
stoich = moose.Stoich( '/model/compartment/stoich' )
stoich.compartment = compartment
stoich.ksolve = ksolve
#ksolve.stoich = stoich
stoich.path = "/model/compartment/##"
state = moose.SteadyState( '/model/compartment/state' )
moose.reinit()
state.stoich = stoich
#state.showMatrices()
state.convergenceCriterion = 1e-7
moose.le( '/model/graphs' )
a = moose.element( '/model/compartment/a' )
b = moose.element( '/model/compartment/b' )
c = moose.element( '/model/compartment/c' )
for i in range( 0, 100 ):
getState( ksolve, state )
moose.start( 100.0 ) # Run the model for 100 seconds.
b = moose.element( '/model/compartment/b' )
c = moose.element( '/model/compartment/c' )
# move most molecules over to b
b.conc = b.conc + c.conc * 0.95
c.conc = c.conc * 0.05
moose.start( 100.0 ) # Run the model for 100 seconds.
# move most molecules back to a
c.conc = c.conc + b.conc * 0.95
b.conc = b.conc * 0.05
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()
#########################################################################
## This program is part of 'MOOSE', the
## Messaging Object Oriented Simulation Environment.
## Copyright (C) 2015 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 matplotlib.pyplot as plt
import moose
diffConst = 1e-11
chemdt = 0.001 # Tested various dts, this is reasonable.
diffdt = 0.001
plotdt = 0.01
animationdt = 0.01
runtime = 1
useGssa = False
def makeModel():
model = moose.Neutral( '/model' )
# Make neuronal model. It has no channels, just for geometry
cell = moose.loadModel( './spinyNeuron.p', '/model/cell', 'Neutral' )
# We don't want the cell to do any calculations. Disable everything.
for i in moose.wildcardFind( '/model/cell/##' ):
i.tick = -1
# create container for model
model = moose.element( '/model' )
chem = moose.Neutral( '/model/chem' )
# The naming of the compartments is dicated by the places that the
# chem model expects to be loaded.
compt0 = moose.NeuroMesh( '/model/chem/compt0' )
compt0.separateSpines = 1
compt0.geometryPolicy = 'cylinder'
compt1 = moose.SpineMesh( '/model/chem/compt1' )
moose.connect( compt0, 'spineListOut', compt1, 'spineList', 'OneToOne' )
compt2 = moose.PsdMesh( '/model/chem/compt2' )
moose.connect( compt0, 'psdListOut', compt2, 'psdList', 'OneToOne' )
#reacSystem = moose.loadModel( 'simpleOsc.g', '/model/chem', 'ee' )
makeChemModel( compt0, True ) # Populate all 3 compts with the chem system.
makeChemModel( compt1, False )
makeChemModel( compt2, True )
compt0.diffLength = 2e-6 # This will be over 100 compartments.
# This is the magic command that configures the diffusion compartments.
compt0.subTreePath = cell.path + "/#"
moose.showfields( compt0 )
# Build the solvers. No need for diffusion in this version.
ksolve0 = moose.Ksolve( '/model/chem/compt0/ksolve' )
if useGssa:
ksolve1 = moose.Gsolve( '/model/chem/compt1/ksolve' )
ksolve2 = moose.Gsolve( '/model/chem/compt2/ksolve' )
else:
ksolve1 = moose.Ksolve( '/model/chem/compt1/ksolve' )
ksolve2 = moose.Ksolve( '/model/chem/compt2/ksolve' )
dsolve0 = moose.Dsolve( '/model/chem/compt0/dsolve' )
dsolve1 = moose.Dsolve( '/model/chem/compt1/dsolve' )
dsolve2 = moose.Dsolve( '/model/chem/compt2/dsolve' )
stoich0 = moose.Stoich( '/model/chem/compt0/stoich' )
stoich1 = moose.Stoich( '/model/chem/compt1/stoich' )
stoich2 = moose.Stoich( '/model/chem/compt2/stoich' )
# Configure solvers
stoich0.compartment = compt0
stoich1.compartment = compt1
stoich2.compartment = compt2
stoich0.ksolve = ksolve0
stoich1.ksolve = ksolve1
stoich2.ksolve = ksolve2
stoich0.dsolve = dsolve0
stoich1.dsolve = dsolve1
stoich2.dsolve = dsolve2
stoich0.path = '/model/chem/compt0/#'
stoich1.path = '/model/chem/compt1/#'
stoich2.path = '/model/chem/compt2/#'
assert( stoich0.numVarPools == 1 )
assert( stoich0.numProxyPools == 0 )
assert( stoich0.numRates == 1 )
assert( stoich1.numVarPools == 1 )
assert( stoich1.numProxyPools == 0 )
if useGssa:
assert( stoich1.numRates == 2 )
assert( stoich2.numRates == 2 )
else:
assert( stoich1.numRates == 1 )
assert( stoich2.numRates == 1 )
assert( stoich2.numVarPools == 1 )
assert( stoich2.numProxyPools == 0 )
dsolve0.buildNeuroMeshJunctions( dsolve1, dsolve2 )
stoich0.buildXreacs( stoich1 )
stoich1.buildXreacs( stoich2 )
stoich0.filterXreacs()
stoich1.filterXreacs()
stoich2.filterXreacs()
Ca_input_dend = moose.vec( '/model/chem/compt0/Ca_input' )
print((len( Ca_input_dend )))
for i in range( 60 ):
Ca_input_dend[ 3 + i * 3 ].conc = 2.0
Ca_input_PSD = moose.vec( '/model/chem/compt2/Ca_input' )
print((len( Ca_input_PSD )))
for i in range( 5 ):
Ca_input_PSD[ 2 + i * 2].conc = 1.0
# Create the output tables
num = compt0.numDiffCompts - 1
graphs = moose.Neutral( '/model/graphs' )
makeTab( 'Ca_soma', '/model/chem/compt0/Ca[0]' )
makeTab( 'Ca_d1', '/model/chem/compt0/Ca[1]' )
makeTab( 'Ca_d2', '/model/chem/compt0/Ca[2]' )
makeTab( 'Ca_d3', '/model/chem/compt0/Ca[3]' )
makeTab( 'Ca_s3', '/model/chem/compt1/Ca[3]' )
makeTab( 'Ca_s4', '/model/chem/compt1/Ca[4]' )
makeTab( 'Ca_s5', '/model/chem/compt1/Ca[5]' )
makeTab( 'Ca_p3', '/model/chem/compt2/Ca[3]' )
makeTab( 'Ca_p4', '/model/chem/compt2/Ca[4]' )
makeTab( 'Ca_p5', '/model/chem/compt2/Ca[5]' )
def makeTab( plotname, molpath ):
tab = moose.Table2( '/model/graphs/' + plotname ) # Make output table
# connect up the tables
moose.connect( tab, 'requestOut', moose.element( molpath ), 'getConc' );
def makeDisplay():
plt.ion()
fig = plt.figure( figsize=(10,12) )
dend = fig.add_subplot( 411 )
plt.ylabel( 'Conc (mM)' )
plt.xlabel( 'Dend voxel #' )
plt.legend()
timeLabel = plt.text(200, 0.5, 'time = 0')
spine = fig.add_subplot( 412 )
plt.ylabel( 'Conc (mM)' )
plt.xlabel( 'Spine voxel #' )
plt.legend()
psd = fig.add_subplot( 413 )
plt.ylabel( 'Conc (mM)' )
plt.xlabel( 'PSD voxel #' )
plt.legend()
timeSeries = fig.add_subplot( 414 )
timeSeries.set_ylim( 0, 2 )
plt.ylabel( 'Conc (mM)' )
plt.xlabel( 'time (seconds)' )
plt.legend()
Ca = moose.vec( '/model/chem/compt0/Ca' )
Ca_input = moose.vec( '/model/chem/compt0/Ca_input' )
line1, = dend.plot( list(range( len( Ca ))), Ca.conc, label='Ca' )
line2, = dend.plot( list(range( len( Ca_input ))), Ca_input.conc, label='Ca_input' )
dend.set_ylim( 0, 2 )
Ca = moose.vec( '/model/chem/compt1/Ca' )
line3, = spine.plot( list(range( len( Ca ))), Ca.conc, label='Ca' )
spine.set_ylim( 0, 1 )
Ca = moose.vec( '/model/chem/compt2/Ca' )
Ca_input = moose.vec( '/model/chem/compt2/Ca_input' )
line4, = psd.plot( list(range( len( Ca ))), Ca.conc, label='Ca' )
line5, = psd.plot( list(range( len( Ca_input ))), Ca_input.conc, label='Ca_input' )
psd.set_ylim( 0, 1 )
fig.canvas.draw()
return ( timeSeries, dend, spine, psd, fig, line1, line2, line3, line4, line5, timeLabel )
def updateDisplay( plotlist ):
Ca = moose.vec( '/model/chem/compt0/Ca' )
Ca_input = moose.vec( '/model/chem/compt0/Ca_input' )
plotlist[5].set_ydata( Ca.conc )
plotlist[6].set_ydata( Ca_input.conc )
Ca = moose.vec( '/model/chem/compt1/Ca' )
plotlist[7].set_ydata( Ca.conc )
Ca = moose.vec( '/model/chem/compt2/Ca' )
Ca_input = moose.vec( '/model/chem/compt2/Ca_input' )
plotlist[8].set_ydata( Ca.conc )
plotlist[9].set_ydata( Ca_input.conc )
plotlist[4].canvas.draw()
def finalizeDisplay( plotlist, cPlotDt ):
for x in moose.wildcardFind( '/model/graphs/#[ISA=Table2]' ):
pos = numpy.arange( 0, x.vector.size, 1 ) * cPlotDt
line1, = plotlist[0].plot( pos, x.vector, label=x.name )
plotlist[4].canvas.draw()
print( "Hit '0' to exit" )
eval(str(input()))
def makeChemModel( compt, doInput ):
"""
This function setus up a simple chemical system in which Ca input
comes to the dend and to selected PSDs. There is diffusion between
PSD and spine head, and between dend and spine head.
:: Ca_input ------> Ca // in dend and spine head only.
"""
# create molecules and reactions
Ca = moose.Pool( compt.path + '/Ca' )
Ca.concInit = 0.08*1e-3
Ca.diffConst = diffConst
if doInput:
Ca_input = moose.BufPool( compt.path + '/Ca_input' )
Ca_input.concInit = 0.08*1e-3
Ca_input.diffConst = diffConst
rInput = moose.Reac( compt.path + '/rInput' )
moose.connect( rInput, 'sub', Ca, 'reac' )
moose.connect( rInput, 'prd', Ca_input, 'reac' )
rInput.Kf = 100 # 1/sec
rInput.Kb = 100 # 1/sec
else:
Ca_sink = moose.BufPool( compt.path + '/Ca_sink' )
Ca_sink.concInit = 0.08*1e-3
rSink = moose.Reac( compt.path + '/rSink' )
moose.connect( rSink, 'sub', Ca, 'reac' )
moose.connect( rSink, 'prd', Ca_sink, 'reac' )
rSink.Kf = 10 # 1/sec
rSink.Kb = 10 # 1/sec
def main():
"""
This example illustrates and tests diffusion embedded in
the branching pseudo-1-dimensional geometry of a neuron.
An input pattern of Ca stimulus is applied in a periodic manner both
on the dendrite and on the PSDs of the 13 spines. The Ca levels in
each of the dend, the spine head, and the spine PSD are monitored.
Since the same molecule name is used for Ca in the three compartments,
these are automagially connected up for diffusion. The simulation
shows the outcome of this diffusion.
This example uses an external electrical model file with basal
dendrite and three branches on
the apical dendrite. One of those branches has the 13 spines.
The model is set up to run using the Ksolve for integration and the
Dsolve for handling diffusion.
The timesteps here are not the defaults. It turns out that the
chem reactions and diffusion in this example are sufficiently fast
that the chemDt has to be smaller than default. Note that this example
uses rates quite close to those used in production models.
The display has four parts:
a. animated line plot of concentration against main compartment#.
b. animated line plot of concentration against spine compartment#.
c. animated line plot of concentration against psd compartment#.
d. time-series plot that appears after the simulation has
ended.
"""
makeModel()
plotlist = makeDisplay()
# Schedule the whole lot - autoscheduling already does this.
for i in range( 11, 17 ):
moose.setClock( i, chemdt ) # for the chem objects
moose.setClock( 10, diffdt ) # for the diffusion
moose.setClock( 18, plotdt ) # for the output tables.
'''
'''
moose.reinit()
for i in numpy.arange( 0, runtime, animationdt ):
moose.start( animationdt )
plotlist[10].set_text( "time = %d" % i )
updateDisplay( plotlist )
finalizeDisplay( plotlist, plotdt )
# Run the 'main' if this script is executed standalone.
if __name__ == '__main__':
main()
'''
FieldElements.
>>> import moose
on node 0, numNodes = 1, numCores = 8
Info: Time to define moose classes:0
Info: Time to initialize module:0.05
>>> a = moose.IntFire('a')
Created 123 path=a numData=1 isGlobal=0 baseType=IntFire
>>> b = moose.Synapse('a/b')
Created 125 path=a/b numData=1 isGlobal=0 baseType=Synapse
>>> b.numSynapse = 10
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
AttributeError: 'moose.Synapse' object has no attribute 'numSynapse'
>>> a.numSynapse = 10
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
AttributeError: 'moose.IntFire' object has no attribute 'numSynapse'
>>> dir(a)
['__class__', '__delattr__', '__doc__', '__eq__', '__format__', '__ge__', '__getattribute__', '__gt__', '__hash__', '__init__', '__le__', '__lt__', '__ne__', '__new__', '__reduce__', '__reduce_ex__', '__repr__', '__setattr__', '__sizeof__', '__str__', '__subclasshook__', 'connect', 'ematrix', 'getDataIndex', 'getField', 'getFieldIndex', 'getFieldNames', 'getFieldType', 'getId', 'getLookupField', 'getNeighbors', 'neighbors', 'parentMsg', 'process', 'reinit', 'setDestField', 'setField', 'setLookupField', 'synapse']
>>> a.synapse
/a.synapse
>>> a.synapse[0]
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
IndexError: moose.ElementField.getItem: index out of bounds.
>>> a.synapse.num = 10
>>> a.synapse[0]
<moose.Synapse: id=124, dataId=0, path=/a/synapse[0]>
>>> a.synapse[1]
<moose.Synapse: id=124, dataId=1, path=/a/synapse[0]>
>>>
----------------------
Subha, Tue Jan 7 11:58:10 IST 2014
Documentation/discussion on accessing synapses (ElementField):
A `synapse` is an ElementField of an IntFire element. The following
example is taken from Upi and modified to reflect the current API::
network = moose.IntFire('/network', size)
network.vec.buffferTime = [delayMax * 2] * size # `vec` allows vectorized access to fields on all `element`s
IntFire elements have an ElementField called `synapse`. You can access
the first `synapse` element on the first IntFire element with the
following::
synapse = moose.element('/network/synapse')
This is equivalent to::
synapse = moose.element('/network[0]/synapse[0]')
`synapse` is just another element but by default its length is 0
unless you set the number of elements in it either explicitly::
synapse.num = 10
or in an indirect way::
mid = moose.connect( network, 'spike', synapse, 'addSpike', 'Sparse' )
mid.setRandomConnectivity( connectionProbability, 5489 )
network.vec.Vm = [(Vmax*random.random()) for r in range(size)]
network.vec.thresh = thresh
network.vec.refractoryPeriod = refractoryPeriod
numSynVec = network.numSynapses
numTotSym = sum( numSynVec )
netvec = network.vec
for i in range( size ):
synvec = netvec[i].synapse.vec
synvec.weight = [ (random.random() * weightMax) for r in range( synvec.len )]
synvec.delay = [ (delayMin + random.random() * delayMax) for r in range( synvec.len )]
moose.useClock( '/network', 'process', 0 )
moose.setClock( 0, timestep )
moose.setClock( 9, timestep )
moose.reinit()
network.vec.Vm = [(Vmax*random.random()) for r in range(size)]
moose.start(runtime)
print network.vec[100].Vm, network.vec[900].Vm
'''
import os
import moose
# Create an IntFire vec containing 10 elements, a refers to alpha[0]
a = moose.IntFire('alpha', 10)
print(('a=', a))
for i in range( 10 ):
syn = moose.SimpleSynHandler( 'alpha[' + str(i) + ']/sh' )
moose.connect( syn, 'activationOut', a.vec[i], 'activation' )
syn = moose.SimpleSynHandler( 'alpha/sh', 10 )
moose.connect( syn, 'activationOut', a, 'activation', 'OneToOne' )
###############################
# FieldElement identity
###############################
x = syn.synapse # x is an ElementField alpha[0].synapse
print('x=',x)
print('x.num=', x.num) # Initially there are no synapses, so this will be 0
syn.synapse.num = 3 # We set number of field elements to 3
print('x.num=', x.num) # x refers to a.synapse, so this should be 3
b = moose.element('alpha[0]/sh/synapse[1]') # We access x[1]
print('b=',b)
print('x[1]=', x[1])
print('b==x[1]?', b == x[1])
###############################
# Check fieldIndex and dataId
###############################
print('syn.synapse[0]=', syn.synapse[0])
print('syn.synapse[1]=', syn.synapse[1]) # The fieldIndex should change, not dataId
#########################
# setVec call example
#########################
print('alpha[0].synapse.delay=', x.delay)
x.delay = [1.0, 2.0, 3.0] # This sets `delay` field in all elements via setVec call
print('alpha[0].synapse.delay=', x.delay)
x.delay = [1.141592] * len(x) # This is a Pythonic way of replicating the values in a list - ensures same length
print('alpha[0].synapse.delay=', x.delay)
#####################################################
# Play a little more with ObjId, FieldElement, Id
#####################################################
print('Length of alpha[1]/synapse=', len(moose.element('/alpha[1]/sh').synapse))
c = moose.element('alpha[1]/sh/synapse[2]') # This should throw an error - alpha[1] does not have 3 synapses.
print('b=', b, 'numData=', b.numData)
print('c=', c, 'numData=', c.numData)
try:
print('len(c)=', len(c))
except TypeError as e:
print(e)
d = moose.element('/alpha[1]/sh')
try:
print(d.synapse[1])
except IndexError as e:
print(e)
else:
print('Expected an IndexError. Length of synapse=', d.synapse.len)
# The fieldIndex should change, not dataId
x = moose.element(a.vec, 0, 1)
y = moose.element(a.vec, 1, 2)
print(x, y)
#########################################################################
## 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
import sys
def makeModel():
if len( sys.argv ) == 1:
useGsolve = True
else:
useGsolve = ( sys.argv[1] == 'True' )
# create container for model
model = moose.Neutral( 'model' )
compartment = moose.CubeMesh( '/model/compartment' )
compartment.volume = 1e-22
# the mesh is created automatically by the compartment
moose.le( '/model/compartment' )
mesh = moose.element( '/model/compartment/mesh' )
# create molecules and reactions
a = moose.Pool( '/model/compartment/a' )
b = moose.Pool( '/model/compartment/b' )
# create functions of time
f1 = moose.Function( '/model/compartment/f1' )
f2 = moose.Function( '/model/compartment/f2' )
# connect them up for reactions
moose.connect( f1, 'valueOut', a, 'setConc' )
moose.connect( f2, 'valueOut', b, 'increment' )
# Assign parameters
a.concInit = 0
b.concInit = 1
#f1.numVars = 1
#f2.numVars = 1
f1.expr = '1 + sin(t)'
f2.expr = '10 * cos(t)'
# Create the output tables
graphs = moose.Neutral( '/model/graphs' )
outputA = moose.Table2 ( '/model/graphs/nA' )
outputB = moose.Table2 ( '/model/graphs/nB' )
# connect up the tables
moose.connect( outputA, 'requestOut', a, 'getN' );
moose.connect( outputB, 'requestOut', b, 'getN' );
# Set up the solvers
if useGsolve:
gsolve = moose.Gsolve( '/model/compartment/gsolve' )
gsolve.useClockedUpdate = True
else:
gsolve = moose.Ksolve( '/model/compartment/gsolve' )
stoich = moose.Stoich( '/model/compartment/stoich' )
stoich.compartment = compartment
stoich.ksolve = gsolve
stoich.path = '/model/compartment/##'
'''
'''
# We need a finer timestep than the default 0.1 seconds,
# in order to get numerical accuracy.
for i in range (10, 19 ):
moose.setClock( i, 0.1 ) # for computational objects
def main():
"""
This example describes the special (and discouraged) use case where
functions provide input to a reaction system. Here we have two functions of
time which control the pool # and pool rate of change, respectively::
number of molecules of a = 1 + sin(t)
rate of change of number of molecules of b = 10 * cos(t)
In the stochastic case one must set a special flag *useClockedUpdate*
in order to achieve clock-triggered updates from the functions. This is
needed because the functions do not have reaction events to trigger them,
and even if there were reaction events they might not be frequent enough to
track the periodic updates. The use of this flag slows down the calculations,
so try to use a table to control a pool instead.
To run in stochastic mode::
''python funcInputToPools''
To run in deterministic mode::
''python funcInputToPools false''
"""
makeModel()
moose.seed()
moose.reinit()
moose.start( 50.0 ) # Run the model for 100 seconds.
a = moose.element( '/model/compartment/a' )
b = moose.element( '/model/compartment/b' )
# Iterate through all plots, dump their contents to data.plot.
for x in moose.wildcardFind( '/model/graphs/n#' ):
#x.xplot( 'scriptKineticModel.plot', x.name )
t = numpy.arange( 0, x.vector.size, 1 ) * x.dt # sec
pylab.plot( t, x.vector, label=x.name )
pylab.legend()
pylab.show()
quit()
# Run the 'main' if this script is executed standalone.
if __name__ == '__main__':
main()
# function.py ---
#
# Filename: function.py
# Description:
# Author: Subhasis Ray
# Maintainer:
# Created: Tue Sep 9 17:59:50 2014 (+0530)
# Version:
# Last-Updated: Sun Dec 20 00:02:50 2015 (-0500)
# By: subha
# Update #: 4
# URL:
# Keywords:
# Compatibility:
#
#
# Commentary:
#
#
#
#
# Change log:
#
#
#
#
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public License as
# published by the Free Software Foundation; either version 3, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; see the file COPYING. If not, write to
# the Free Software Foundation, Inc., 51 Franklin Street, Fifth
# Floor, Boston, MA 02110-1301, USA.
#
#
# Code:
import numpy as np
import sys
import matplotlib.pyplot as plt
import moose
simtime = 1.0
def example():
demo = moose.Neutral('/model')
function = moose.Function('/model/function')
function.c['c0'] = 1.0
function.c['c1'] = 2.0
function.x.num = 1
function.expr = 'c0 * exp(c1*x0) * cos(y0) + sin(t)'
# mode 0 - evaluate function value, derivative and rate
# mode 1 - just evaluate function value,
# mode 2 - evaluate derivative,
# mode 3 - evaluate rate
function.mode = 0
function.independent = 'y0'
nsteps = 1000
xarr = np.linspace(0.0, 1.0, nsteps)
# Stimulus tables allow you to store sequences of numbers which
# are delivered via the 'output' message at each time step. This
# is a placeholder and in real scenario you will be using any
# sourceFinfo that sends out a double value.
input_x = moose.StimulusTable('/xtab')
input_x.vector = xarr
input_x.startTime = 0.0
input_x.stepPosition = xarr[0]
input_x.stopTime = simtime
moose.connect(input_x, 'output', function.x[0], 'input')
yarr = np.linspace(-np.pi, np.pi, nsteps)
input_y = moose.StimulusTable('/ytab')
input_y.vector = yarr
input_y.startTime = 0.0
input_y.stepPosition = yarr[0]
input_y.stopTime = simtime
moose.connect(function, 'requestOut', input_y, 'getOutputValue')
# data recording
result = moose.Table('/ztab')
moose.connect(result, 'requestOut', function, 'getValue')
derivative = moose.Table('/zprime')
moose.connect(derivative, 'requestOut', function, 'getDerivative')
rate = moose.Table('/dz_by_dt')
moose.connect(rate, 'requestOut', function, 'getRate')
x_rec = moose.Table('/xrec')
moose.connect(x_rec, 'requestOut', input_x, 'getOutputValue')
y_rec = moose.Table('/yrec')
moose.connect(y_rec, 'requestOut', input_y, 'getOutputValue')
dt = simtime/nsteps
for ii in range(32):
moose.setClock(ii, dt)
moose.reinit()
moose.start(simtime)
# Uncomment the following lines and the import matplotlib.pyplot as plt on top
# of this file to display the plot.
plt.subplot(3,1,1)
plt.plot(x_rec.vector, result.vector, 'r-', label='z = {}'.format(function.expr))
z = function.c['c0'] * np.exp(function.c['c1'] * xarr) * np.cos(yarr) + np.sin(np.arange(len(xarr)) * dt)
plt.plot(xarr, z, 'b--', label='numpy computed')
plt.xlabel('x')
plt.ylabel('z')
plt.legend()
plt.subplot(3,1,2)
plt.plot(y_rec.vector, derivative.vector, 'r-', label='dz/dy0')
# derivatives computed by putting x values in the analytical formula
dzdy = function.c['c0'] * np.exp(function.c['c1'] * xarr) * (- np.sin(yarr))
plt.plot(yarr, dzdy, 'b--', label='numpy computed')
plt.xlabel('y')
plt.ylabel('dz/dy')
plt.legend()
plt.subplot(3,1,3)
# *** BEWARE *** The first two entries are spurious. Entry 0 is
# *** from reinit sending out the defaults. Entry 2 is because
# *** there is no lastValue for computing real forward difference.
plt.plot(np.arange(2, len(rate.vector), 1) * dt, rate.vector[2:], 'r-', label='dz/dt')
dzdt = np.diff(z)/dt
plt.plot(np.arange(0, len(dzdt), 1.0) * dt, dzdt, 'b--', label='numpy computed')
plt.xlabel('t')
plt.ylabel('dz/dt')
plt.legend()
plt.tight_layout()
plt.show()
def main():
"""
Function objects can be used to evaluate expressions with arbitrary
number of variables and constants. We can assign expression of the
form::
f(c_, c1, ..., cM, x0, x1, ..., xN, y0,..., yP )
where `ci`'s are constants and `xi`'s and `yi`'s are variables.
The constants must be defined before setting the expression and
variables are connected via messages. The constants can have any
name, but the variable names must be of the form x{i} or y{i}
where i is increasing integer starting from 0.
The `xi`'s are field elements and you have to set their number
first (function.x.num = N). Then you can connect any source field
sending out double to the 'input' destination field of the
`x[i]`.
The `yi`'s are useful when the required variable is a value field
and is not available as a source field. In that case you connect
the `requestOut` source field of the function element to the
`get{Field}` destination field on the target element. The `yi`'s
are automatically added on connecting. Thus, if you call::
moose.connect(function, 'requestOut', a, 'getSomeField')
moose.connect(function, 'requestOut', b, 'getSomeField')
then ``a.someField`` will be assigned to ``y0`` and
``b.someField`` will be assigned to ``y1``.
In this example we evaluate the expression: ``z = c0 * exp(c1 *
x0) * cos(y0)``
with x0 ranging from -1 to +1 and y0 ranging from -pi to
+pi. These values are stored in two stimulus tables called xtab
and ytab respectively, so that at each timestep the next values of
x0 and y0 are assigned to the function.
Along with the value of the expression itself we also compute its
derivative with respect to y0 and its derivative with respect to
time (rate). The former uses a five-point stencil for the
numerical differentiation and has a glitch at y=0. The latter uses
backward difference divided by dt.
Unlike Func class, the number of variables and constants are
unlimited in Function and you can set all the variables via
messages.
"""
example()
if __name__ == '__main__':
main()
#
# function.py ends here
# intfire.py ---
#
# Filename: intfire.py
# Description:
# Author:Subhasis Ray
# Maintainer:
# Created: Thu Jun 21 16:40:25 2012 (+0530)
# Version:
# Last-Updated: Sat Jun 23 13:44:10 2012 (+0530)
# By: subha
# Update #: 35
# URL:
# Keywords:
# Compatibility:
#
#
# Commentary:
#
# Code snippet to show some operations on IntFire.
#
#
# Change log:
#
#
#
#
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public License as
# published by the Free Software Foundation; either version 3, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; see the file COPYING. If not, write to
# the Free Software Foundation, Inc., 51 Franklin Street, Fifth
# Floor, Boston, MA 02110-1301, USA.
#
#
# Code:
import moose
def connect_two_intfires():
"""
Connect two IntFire neurons so that spike events in one gets
transmitted to synapse of the other.
"""
if1 = moose.IntFire('if1')
if2 = moose.IntFire('if2')
sf1 = moose.SimpleSynHandler( 'if1/sh' )
moose.connect( sf1, 'activationOut', if1, 'activation' )
sf1.synapse.num = 1
syn1 = moose.element(sf1.synapse)
# Connect the spike message of if2 to the first synapse on if1
moose.connect(if2, 'spikeOut', syn1, 'addSpike')
def connect_spikegen():
"""
Connect a SpikeGen object to an IntFire neuron such that spike
events in spikegen get transmitted to the synapse of the IntFire
neuron.
"""
if3 = moose.IntFire('if3')
sf3 = moose.SimpleSynHandler( 'if3/sh' )
moose.connect( sf3, 'activationOut', if3, 'activation' )
sf3.synapse.num = 1
sg = moose.SpikeGen('sg')
syn = moose.element(sf3.synapse)
moose.connect(sg, 'spikeOut', syn, 'addSpike')
def setup_synapse():
"""
Create an intfire object and create two synapses on it.
"""
if4 = moose.IntFire('if4')
sf4 = moose.SimpleSynHandler( 'if4/sh' )
sg1 = moose.SpikeGen('sg1')
sg2 = moose.SpikeGen('sg2')
sf4.synapse.num = 2 # set synapse count to 2
sf4.synapse[0].weight = 0.5
sf4.synapse[0].delay = 1e-3
sf4.synapse[1].weight = 2.0
sf4.synapse[1].delay = 2e-3
moose.connect(sg1, 'spikeOut', sf4.synapse[0], 'addSpike')
moose.connect(sg2, 'spikeOut', sf4.synapse[1], 'addSpike')
def main():
"""
Demonstrates connection between 2 IntFire neurons to observe
spike generation.
"""
connect_two_intfires()
connect_spikegen()
setup_synapse()
if __name__ == '__main__':
main()
#
# intfire.py ends here
# loadKineticModel.py ---
#
# Filename: loadKineticModel.py
# Description:
# Author: Upi Bhalla
# Maintainer:
# Created: Sat Oct 04 12:14:15 2014 (+0530)
# Version:
# Last-Updated:
# By:
# Update #: 0
# URL:
# Keywords:
# Compatibility:
#
#
# Commentary:
#
#
#
#
# Change log:
#
#
#
#
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public License as
# published by the Free Software Foundation; either version 3, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; see the file COPYING. If not, write to
# the Free Software Foundation, Inc., 51 Franklin Street, Fifth
# Floor, Boston, MA 02110-1301, USA.
#
#
# Code:
import moose
import pylab
import numpy
import sys
def main():
"""
This example illustrates loading, running, and saving a kinetic
model defined in kkit format. It uses a default kkit model but
you can specify another using the command line
``python filename runtime solver``.
We use the gsl solver here.
The model already defines a couple of plots and sets the runtime 20 secs.
"""
solver = "gsl" # Pick any of gsl, gssa, ee..
mfile = '../genesis/kkit_objects_example.g'
runtime = 20.0
if ( len( sys.argv ) >= 3 ):
if sys.argv[1][0] == '/':
mfile = sys.argv[1]
else:
mfile = '../genesis/' + sys.argv[1]
runtime = float( sys.argv[2] )
if ( len( sys.argv ) == 4 ):
solver = sys.argv[3]
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
moose.reinit()
moose.start( runtime )
# Display all plots.
for x in moose.wildcardFind( '/model/#graphs/conc#/#' ):
t = numpy.arange( 0, x.vector.size, 1 ) * x.dt
pylab.plot( t, x.vector, label=x.name )
pylab.legend()
pylab.show()
quit()
# Run the 'main' if this script is executed standalone.
if __name__ == '__main__':
main()
import moogli
import moose
from moose import neuroml
from PyQt4 import Qt, QtCore, QtGui
import sys
import os
frameRunTime = 0.0002
runtime = 1.0
inject = 5e-10
simdt = 5e-5
def main():
app = QtGui.QApplication(sys.argv)
filename = 'barrionuevo_cell1zr.CNG.swc'
moose.Neutral( '/library' )
moose.Neutral( '/model' )
cell = moose.loadModel( filename, '/model/testSwc' )
for i in range( 8 ):
moose.setClock( i, simdt )
hsolve = moose.HSolve( '/model/testSwc/hsolve' )
hsolve.dt = simdt
hsolve.target = '/model/testSwc/soma'
moose.le( cell )
moose.reinit()
# Now we set up the display
compts = moose.wildcardFind( "/model/testSwc/#[ISA=CompartmentBase]" )
compts[0].inject = inject
ecomptPath = [x.path for x in compts]
morphology = moogli.read_morphology_from_moose(name = "", path = "/model/testSwc")
morphology.create_group( "group_all", ecomptPath, -0.08, 0.02, \
[0.0, 0.5, 1.0, 1.0], [1.0, 0.0, 0.0, 0.9] )
viewer = moogli.DynamicMorphologyViewerWidget(morphology)
def callback( morphology, viewer ):
moose.start( frameRunTime )
Vm = [moose.element( x ).Vm for x in compts]
morphology.set_color( "group_all", Vm )
currTime = moose.element( '/clock' ).currentTime
#print currTime, compts[0].Vm
if ( currTime < runtime ):
return True
return False
viewer.set_callback( callback, idletime = 0 )
viewer.showMaximized()
viewer.show()
app.exec_()
if __name__ == '__main__':
main()
# vclamp.py ---
#
# Filename: vclamp.py
# Description:
# Author:Subhasis Ray
# Maintainer:
# Created: Sat Feb 2 19:16:54 2013 (+0530)
# Version:
# Last-Updated: Tue Jun 11 17:35:20 2013 (+0530)
# By: subha
# Update #: 178
# URL:
# Keywords:
# Compatibility:
#
#
# Commentary:
#
#
#
#
# Change log:
#
#
#
#
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public License as
# published by the Free Software Foundation; either version 3, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; see the file COPYING. If not, write to
# the Free Software Foundation, Inc., 51 Franklin Street, Fifth
# Floor, Boston, MA 02110-1301, USA.
#
#
# Code:
import sys
sys.path.append('../../python')
import moose
sys.path.append('../squid')
from squid import SquidAxon
from pylab import *
def vclamp_demo(simtime=50.0, dt=1e-2):
## It is good practice to modularize test elements inside a
## container
container = moose.Neutral('/vClampDemo')
## Create a compartment with properties of a squid giant axon
comp = SquidAxon('/vClampDemo/axon')
# Create and setup the voltage clamp object
clamp = moose.VClamp('/vClampDemo/vclamp')
## The defaults should work fine
# clamp.mode = 2
# clamp.tau = 10*dt
# clamp.ti = dt
# clamp.td = 0
# clamp.gain = comp.Cm / dt
## Setup command voltage time course
command = moose.PulseGen('/vClampDemo/command')
command.delay[0] = 10.0
command.width[0] = 20.0
command.level[0] = 50.0
command.delay[1] = 1e9
moose.connect(command, 'output', clamp, 'commandIn')
## Connect the Voltage Clamp to the compartemnt
moose.connect(clamp, 'currentOut', comp, 'injectMsg')
moose.connect(comp, 'VmOut', clamp, 'sensedIn')
## setup stimulus recroding - this is the command pulse
stimtab = moose.Table('/vClampDemo/vclamp_command')
moose.connect(stimtab, 'requestOut', command, 'getOutputValue')
## Set up Vm recording
vmtab = moose.Table('/vClampDemo/vclamp_Vm')
moose.connect(vmtab, 'requestOut', comp, 'getVm')
## setup command potential recording - this is the filtered input
## to PID controller
commandtab = moose.Table('/vClampDemo/vclamp_filteredcommand')
moose.connect(commandtab, 'requestOut', clamp, 'getCommand')
## setup current recording
Imtab = moose.Table('/vClampDemo/vclamp_inject')
moose.connect(Imtab, 'requestOut', clamp, 'getCurrent')
# Scheduling
moose.setClock(0, dt)
moose.setClock(1, dt)
moose.setClock(2, dt)
moose.setClock(3, dt)
moose.useClock(0, '%s/##[TYPE=Compartment]' % (container.path), 'init')
moose.useClock(0, '%s/##[TYPE=PulseGen]' % (container.path), 'process')
moose.useClock(1, '%s/##[TYPE=Compartment]' % (container.path), 'process')
moose.useClock(2, '%s/##[TYPE=HHChannel]' % (container.path), 'process')
moose.useClock(2, '%s/##[TYPE=VClamp]' % (container.path), 'process')
moose.useClock(3, '%s/##[TYPE=Table]' % (container.path), 'process')
moose.reinit()
print(('RC filter in VClamp:: tau:', clamp.tau))
print(('PID controller in VClamp:: ti:', clamp.ti, 'td:', clamp.td, 'gain:', clamp.gain))
moose.start(simtime)
print(('Finished simulation for %g seconds' % (simtime)))
tseries = linspace(0, simtime, len(vmtab.vector))
subplot(211)
title('Membrane potential and clamp voltage')
plot(tseries, vmtab.vector, 'g-', label='Vm (mV)')
plot(tseries, commandtab.vector, 'b-', label='Filtered command (mV)')
plot(tseries, stimtab.vector, 'r-', label='Command (mV)')
xlabel('Time (ms)')
ylabel('Voltage (mV)')
legend()
# print len(commandtab.vector)
subplot(212)
title('Current through clamp circuit')
# plot(tseries, stimtab.vector, label='stimulus (uA)')
plot(tseries, Imtab.vector, label='injected current (uA)')
xlabel('Time (ms)')
ylabel('Current (uA)')
legend()
show()
def main():
"""
This snippet is to demonstrate modelling of voltage clamping.
"""
vclamp_demo()
if __name__ == '__main__':
main()
#
# vclamp.py ends here
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please to comment