Source code for crossComptNeuroMesh

#########################################################################
# crossComptNeuroMesh.py ---
#
# Filename:  crossComptNeuroMesh.py
# Author: Upinder S. Bhalla
# Maintainer:
# Created: Oct  12 16:26:05 2014 (+0530)
# Version:
# Last-Updated: May 15 2017
#           By:
#     Update #:
# URL:
# Keywords:
# Compatibility:
#
#

# Commentary:
#
#
#
#
# Change log: Indentation clean up
#
## 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 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' )

    try:
        ksolve0.numThread = 10
        ksolve1.numThread = 10
        ksolve2.numThread = 10
    except Exception as e:
        pass

    # 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' );


[docs]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. Concentrations 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()