# multiComptSigNeur.py ---
# Upi Bhalla NCBS Bangalore 2013.
# Commentary:
#
# 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
import os
import math
import pylab,numpy
import moose
import numpy
import matplotlib.pyplot as plt
EREST_ACT = -70e-3
# Gate equations have the form:
#
# y(x) = (A + B * x) / (C + exp((x + D) / F))
#
# where x is membrane voltage and y is the rate constant for gate
# closing or opening
Na_m_params = [1e5 * (25e-3 + EREST_ACT), # 'A_A':
-1e5, # 'A_B':
-1.0, # 'A_C':
-25e-3 - EREST_ACT, # 'A_D':
-10e-3, # 'A_F':
4e3, # 'B_A':
0.0, # 'B_B':
0.0, # 'B_C':
0.0 - EREST_ACT, # 'B_D':
18e-3 # 'B_F':
]
Na_h_params = [ 70.0, # 'A_A':
0.0, # 'A_B':
0.0, # 'A_C':
0.0 - EREST_ACT, # 'A_D':
0.02, # 'A_F':
1000.0, # 'B_A':
0.0, # 'B_B':
1.0, # 'B_C':
-30e-3 - EREST_ACT, # 'B_D':
-0.01 # 'B_F':
]
K_n_params = [ 1e4 * (10e-3 + EREST_ACT), # 'A_A':
-1e4, # 'A_B':
-1.0, # 'A_C':
-10e-3 - EREST_ACT, # 'A_D':
-10e-3, # 'A_F':
0.125e3, # 'B_A':
0.0, # 'B_B':
0.0, # 'B_C':
0.0 - EREST_ACT, # 'B_D':
80e-3 # 'B_F':
]
VMIN = -30e-3 + EREST_ACT
VMAX = 120e-3 + EREST_ACT
VDIVS = 3000
[docs]def createSquid():
"""Create a single compartment squid model."""
parent = moose.Neutral ('/n' )
compt = moose.SymCompartment( '/n/compt' )
Em = EREST_ACT + 10.613e-3
compt.Em = Em
compt.initVm = EREST_ACT
compt.Cm = 7.85e-9 * 0.5
compt.Rm = 4.2e5 * 5.0
compt.Ra = 7639.44e3
nachan = moose.HHChannel( '/n/compt/Na' )
nachan.Xpower = 3
xGate = moose.HHGate(nachan.path + '/gateX')
xGate.setupAlpha(Na_m_params + [VDIVS, VMIN, VMAX])
xGate.useInterpolation = 1
nachan.Ypower = 1
yGate = moose.HHGate(nachan.path + '/gateY')
yGate.setupAlpha(Na_h_params + [VDIVS, VMIN, VMAX])
yGate.useInterpolation = 1
nachan.Gbar = 0.942e-3
nachan.Ek = 115e-3+EREST_ACT
moose.connect(nachan, 'channel', compt, 'channel', 'OneToOne')
kchan = moose.HHChannel( '/n/compt/K' )
kchan.Xpower = 4.0
xGate = moose.HHGate(kchan.path + '/gateX')
xGate.setupAlpha(K_n_params + [VDIVS, VMIN, VMAX])
xGate.useInterpolation = 1
kchan.Gbar = 0.2836e-3
kchan.Ek = -12e-3+EREST_ACT
moose.connect(kchan, 'channel', compt, 'channel', 'OneToOne')
return compt
[docs]def createSpine( parentCompt, parentObj, index, frac, length, dia, theta ):
"""Create spine of specified dimensions and index"""
RA = 1.0
RM = 1.0
CM = 0.01
shaftDia = dia / 5.0
sname = 'shaft' + str(index)
hname = 'head' + str(index)
shaft = moose.SymCompartment( parentObj.path + '/' + sname )
moose.connect( parentCompt, 'cylinder', shaft, 'proximalOnly','Single' )
x = parentCompt.x0 + frac * ( parentCompt.x - parentCompt.x0 )
y = parentCompt.y0 + frac * ( parentCompt.y - parentCompt.y0 )
z = parentCompt.z0 + frac * ( parentCompt.z - parentCompt.z0 )
shaft.x0 = x
shaft.y0 = y
shaft.z0 = z
sy = y + length * math.cos( theta * math.pi / 180.0 )
sz = z + length * math.sin( theta * math.pi / 180.0 )
shaft.x = x
shaft.y = sy
shaft.z = sz
shaft.diameter = dia / 2.0
shaft.length = length
xa = math.pi * shaftDia * shaftDia / 4
circumference = math.pi * shaftDia
shaft.Ra = RA * length / xa
shaft.Rm = RM / ( length * circumference )
shaft.Cm = CM * length * circumference
shaft.Em = EREST_ACT
shaft.initVm = EREST_ACT
head = moose.SymCompartment( parentObj.path + '/' + hname )
moose.connect( shaft, 'distal', head, 'proximal', 'Single' )
head.x0 = x
head.y0 = sy
head.z0 = sz
hy = sy + length * math.cos ( theta * math.pi / 180.0 )
hz = sz + length * math.sin ( theta * math.pi / 180.0 )
head.x = x
head.y = hy
head.z = hz
head.diameter = dia
head.length = length
xa = math.pi * dia * dia / 4.0
circumference = math.pi * dia
head.Ra = RA * length / xa
head.Rm = RM / ( length * circumference )
head.Cm = CM * length * circumference
head.Em = EREST_ACT
head.initVm = EREST_ACT
#print head.Rm, head.Ra, head.Cm, head.diameter, head.length
#print shaft.Rm, shaft.Ra, shaft.Cm, shaft.diameter, shaft.length
return head
def createSpineWithReceptor( compt, cell, index, frac ):
FaradayConst = 96485.3415 # s A / mol
spineLength = 5.0e-6
spineDia = 4.0e-6
head = createSpine( compt, cell, index, frac, spineLength, spineDia, 0.0 )
gluR = moose.SynChan( head.path + '/gluR' )
gluR.tau1 = 4e-3
gluR.tau2 = 4e-3
gluR.Gbar = 1e-6
gluR.Ek= 10.0e-3
moose.connect( head, 'channel', gluR, 'channel', 'Single' )
synh = moose.SimpleSynHandler( gluR.path + '/synh' )
moose.connect( synh, 'activationOut', gluR, 'activation' )
caPool = moose.CaConc( head.path + '/ca' )
caPool.CaBasal = 1e-4 # 0.1 micromolar
caPool.tau = 0.01
B = 1.0 / ( FaradayConst * spineLength * spineDia * spineDia *math.pi/4)
B = B / 20.0 # scaling factor for Ca buffering
caPool.B = B
moose.connect( gluR, 'IkOut', caPool, 'current', 'Single' )
return synh
def addPlot( objpath, field, plot ):
assert moose.exists( objpath )
tab = moose.Table( '/graphs/' + plot )
obj = moose.element( objpath )
moose.connect( tab, 'requestOut', obj, field )
return tab
def makeElecPlots():
graphs = moose.Neutral( '/graphs' )
elec = moose.Neutral( '/graphs/elec' )
addPlot( '/model/elec/compt', 'getVm', 'elec/dendVm' )
#addPlot( '/n/compt/Na', 'getGbar', 'elec/NaGbar' )
addPlot( '/model/elec/compt/K', 'getGbar', 'elec/KGbar' )
#addPlot( '/n/compt/Na', 'getIk', 'elec/NaIk' )
#addPlot( '/n/compt/K', 'getIk', 'elec/KIk' )
#addPlot( '/n/compt/Na', 'getEk', 'elec/NaEk' )
#addPlot( '/n/compt/K', 'getEk', 'elec/KEk' )
addPlot( '/model/elec/head0', 'getVm', 'elec/head0Vm' )
addPlot( '/model/elec/head2', 'getVm', 'elec/head2Vm' )
#addPlot( '/n/head2', 'getIm', 'elec/head2Im' )
addPlot( '/model/elec/head0/ca', 'getCa', 'elec/head0Ca' )
addPlot( '/model/elec/head1/ca', 'getCa', 'elec/head1Ca' )
addPlot( '/model/elec/head2/ca', 'getCa', 'elec/head2Ca' )
addPlot( '/model/elec/head0/gluR', 'getIk', 'elec/head0Ik' )
addPlot( '/model/elec/head1/gluR', 'getIk', 'elec/head1Ik' )
addPlot( '/model/elec/head2/gluR', 'getIk', 'elec/head2Ik' )
addPlot( '/model/elec/head1/gluR', 'getGbar', 'elec/head1Gbar' )
addPlot( '/model/elec/head2/gluR', 'getGbar', 'elec/head2Gbar' )
#addPlot( '/n/head0/gluR', 'getGk', 'elec/head0Gk' )
#addPlot( '/n/head2/gluR', 'getGk', 'elec/head2Gk' )
def dumpPlots( fname,runtime ):
if ( os.path.exists( fname ) ):
os.remove( fname )
for x in moose.wildcardFind( '/graphs/##[ISA=Table]' ):
#x.xplot( fname, x.name )
t = numpy.arange( 0, x.vector.size, 1 ) * x.dt
pylab.plot( t, x.vector, label=x.name )
pylab.legend()
pylab.show()
def makeSpinyCompt():
comptLength = 30e-6
comptDia = 6e-6
numSpines = 5
compt = createSquid()
compt.inject = 0
compt.x0 = 0
compt.y0 = 0
compt.z0 = 0
compt.x = comptLength
compt.y = 0
compt.z = 0
compt.length = comptLength
compt.diameter = comptDia
#kchan = moose.element( '/n/compt/K' )
#kchan.Gbar = 0.2e-3
synInput = moose.SpikeGen( '/n/compt/synInput' )
synInput.refractT = 47e-3
synInput.threshold = -1.0
synInput.edgeTriggered = 0
synInput.Vm( 0 )
cell = moose.element( '/n' )
for i in range( numSpines ):
r = createSpineWithReceptor( compt, cell, i, i/float(numSpines) )
r.synapse.num = 1
syn = moose.element( r.path + '/synapse' )
moose.connect( synInput, 'spikeOut', syn, 'addSpike', 'Single' )
syn.weight = 0.2 * i * ( 4 - i )
syn.delay = i * 1.0e-3
def createPool( compt, name, concInit ):
meshEntries = moose.element( compt.path + '/mesh' )
pool = moose.Pool( compt.path + '/' + name )
pool.concInit = concInit
pool.diffConst = 1e-11
return pool
def createChemModel( neuroCompt, spineCompt, psdCompt ):
# Stuff in spine + psd
#psdCa = createPool( psdCompt, 'Ca', 0.0001 )
psdGluR = createPool( psdCompt, 'psdGluR', 1 )
headCa = createPool( spineCompt, 'Ca', 1e-4 )
headGluR = createPool( spineCompt, 'headGluR', 2 )
toPsd = createPool( spineCompt, 'toPsd', 0 )
toPsdInact = createPool( spineCompt, 'toPsdInact', 1e-3 )
turnOnPsd = moose.Reac( spineCompt.path + '/turnOnPsd' )
moose.connect( turnOnPsd, 'sub', headCa, 'reac', 'OneToOne' )
moose.connect( turnOnPsd, 'sub', toPsdInact, 'reac', 'OneToOne' )
moose.connect( turnOnPsd, 'prd', toPsd, 'reac', 'OneToOne' )
turnOnPsd.Kf = 1e3
turnOnPsd.Kb = 1
toPsdEnz = moose.Enz( toPsd.path + '/enz' )
toPsdEnzCplx = moose.Pool( toPsdEnz.path + '/cplx' )
toPsdEnzCplx.concInit = 0
moose.connect( toPsdEnz, 'enz', toPsd, 'reac', 'OneToOne' )
moose.connect( toPsdEnz, 'sub', headGluR, 'reac', 'OneToOne' )
moose.connect( toPsdEnz, 'prd', psdGluR, 'reac', 'OneToOne' )
moose.connect( toPsdEnz, 'cplx', toPsdEnzCplx, 'reac', 'OneToOne' )
toPsdEnz.Km = 1.0e-3
toPsdEnz.kcat = 10.0
fromPsd = moose.Reac( psdCompt.path + '/fromPsd' )
moose.connect( fromPsd, 'sub', psdGluR, 'reac', 'OneToOne' )
moose.connect( fromPsd, 'prd', headGluR, 'reac', 'OneToOne' )
fromPsd.Kf = 0.5
fromPsd.Kb = 0.0
# Stuff in dendrite
dendCa = createPool( neuroCompt, 'Ca', 1e-4 )
bufCa = moose.Pool( neuroCompt.path + '/bufCa' )
bufCa.concInit = 1e-4
pumpCa = moose.Reac( neuroCompt.path + '/pumpCa' )
moose.connect( pumpCa, 'sub', dendCa, 'reac', 'OneToOne' )
moose.connect( pumpCa, 'prd', bufCa, 'reac', 'OneToOne' )
pumpCa.Kf = 1
pumpCa.Kb = 1
dendKinaseInact = createPool( neuroCompt, 'inact_kinase', 1e-4 )
dendKinase = createPool( neuroCompt, 'Ca.kinase', 0.0 )
dendTurnOnKinase = moose.Reac( neuroCompt.path + '/turnOnKinase' )
moose.connect( dendTurnOnKinase, 'sub', dendCa, 'reac' )
moose.connect( dendTurnOnKinase, 'sub', dendKinaseInact, 'reac' )
moose.connect( dendTurnOnKinase, 'prd', dendKinase, 'reac' )
dendTurnOnKinase.Kf = 50000
dendTurnOnKinase.Kb = 1
dendKinaseEnz = moose.Enz( dendKinase.path + '/enz' )
dendKinaseEnzCplx = moose.Pool( dendKinase.path + '/enz/cplx' )
kChan = createPool( neuroCompt, 'kChan', 1e-3 )
kChan_p = createPool( neuroCompt, 'kChan_p', 0.0 )
moose.connect( dendKinaseEnz, 'enz', dendKinase, 'reac', 'OneToOne' )
moose.connect( dendKinaseEnz, 'sub', kChan, 'reac', 'OneToOne' )
moose.connect( dendKinaseEnz, 'prd', kChan_p, 'reac', 'OneToOne' )
moose.connect( dendKinaseEnz, 'cplx', dendKinaseEnzCplx, 'reac', 'OneToOne' )
dendKinaseEnz.Km = 1e-4
dendKinaseEnz.kcat = 20
dendPhosphatase = moose.Reac( neuroCompt.path + '/phosphatase' )
moose.connect( dendPhosphatase, 'sub', kChan_p, 'reac' )
moose.connect( dendPhosphatase, 'prd', kChan, 'reac' )
dendPhosphatase.Kf = 1
dendPhosphatase.Kb = 0.0
# Just for printf debugging
def printMolVecs( title ):
print(title)
"""
nCa = moose.vec( '/model/chem/neuroMesh/Ca' )
sCa = moose.vec( '/model/chem/spineMesh/Ca' )
sR = moose.vec( '/model/chem/spineMesh/headGluR' )
pR = moose.vec( '/model/chem/psdMesh/psdGluR' )
print 'sizes: nCa, sCa, sR, pR = ', len(nCa), len(sCa), len(sR), len(pR)
#print 'nCa=', nCa.conc, ', sCa=', sCa.conc, ', sR=', sR.n, ', pR=', pR.n
print 'nCaConcInit=', nCa.concInit, ', sCa=', sCa.concInit
#print 'sRnInit=', sR.nInit, ', pR=', pR.nInit
print 'sRconcInit=', sR.concInit, ', pR=', pR.concInit
#print 'nCaSize=', nCa.volume, ', sCa=', sCa.volume, ', sR=', sR.n, ', pR=', pR.n
"""
def makeChemInCubeMesh():
dendSide = 10.8e-6
spineSide = 6.8e-6
psdSide = 8.565e-7
parent = moose.Neutral ('/model/chem' )
neuroMesh = moose.CubeMesh( '/model/chem/neuroMesh' )
spineMesh = moose.CubeMesh( '/model/chem/spineMesh' )
psdMesh = moose.CubeMesh( '/model/chem/psdMesh' )
coords = [dendSide] * 9
coords[0] = 0
coords[1] = 0
coords[2] = 0
neuroMesh.coords = coords
neuroMesh.preserveNumEntries = 1
coords = [spineSide] * 9
coords[0] = dendSide
coords[1] = 0
coords[2] = 0
coords[3] = spineSide + dendSide
spineMesh.coords = coords
spineMesh.preserveNumEntries = 1
coords = [psdSide] * 9
coords[0] = dendSide + spineSide
coords[1] = 0
coords[2] = 0
coords[3] = psdSide + spineSide + dendSide
psdMesh.coords = coords
psdMesh.preserveNumEntries = 1
createChemModel( neuroMesh, spineMesh, psdMesh )
dendCa = moose.element( '/model/chem/neuroMesh/Ca' )
assert dendCa.volume == dendSide * dendSide * dendSide
spineCa = moose.element( '/model/chem/spineMesh/Ca' )
assert spineCa.volume == spineSide * spineSide * spineSide
psdGluR = moose.element( '/model/chem/psdMesh/psdGluR' )
assert psdGluR.volume == psdSide * psdSide * psdSide
dendKinaseEnzCplx = moose.element( '/model/chem/neuroMesh/Ca.kinase/enz/cplx' )
assert dendKinaseEnzCplx.volume == dendSide * dendSide * dendSide
def makeSolvers( elecDt ):
# Put in the solvers, see how they fare.
# Here we kludge in a single chem solver for the whole system.
ksolve = moose.Ksolve( '/model/ksolve' )
stoich = moose.Stoich( '/model/stoich' )
stoich.compartment = moose.element( '/model/chem/neuroMesh' )
stoich.ksolve = ksolve
stoich.path = '/model/chem/##'
#stoich.method = 'rk5'
moose.useClock( 5, '/model/ksolve', 'init' )
moose.useClock( 6, '/model/ksolve', 'process' )
# Here is the elec solver
hsolve = moose.HSolve( '/model/hsolve' )
moose.useClock( 1, '/model/hsolve', 'process' )
hsolve.dt = elecDt
hsolve.target = '/model/elec/compt'
def makeCubeMultiscale():
makeSpinyCompt()
model = moose.Neutral( '/model' )
elec = moose.element( '/n' )
elec.name = 'elec'
moose.move( elec, model )
synInput = moose.element( '/model/elec/compt/synInput' )
synInput.refractT = 47e-3
makeChemInCubeMesh()
# set up a reaction to fake diffusion between compts.
headCa = moose.element( '/model/chem/spineMesh/Ca' )
dendCa = moose.element( '/model/chem/neuroMesh/Ca' )
diffReac = moose.Reac( '/model/chem/spineMesh/diff' )
moose.connect( diffReac, 'sub', headCa, 'reac' )
moose.connect( diffReac, 'prd', dendCa, 'reac' )
diffReac.Kf = 1
diffReac.Kb = headCa.volume / dendCa.volume
# set up adaptors
headCa = moose.element( '/model/chem/spineMesh/Ca' )
dendCa = moose.element( '/model/chem/neuroMesh/Ca' )
adaptCa = moose.Adaptor( '/model/chem/adaptCa' )
elecCa = moose.element( '/model/elec/head2/ca' )
# There are 5 spine heads in the electrical model. Average their input.
for i in range( 5 ):
path = '/model/elec/head' + str( i ) + '/ca'
elecCa = moose.element( path )
moose.connect( elecCa, 'concOut', adaptCa, 'input', 'Single' )
moose.connect( adaptCa, 'output', headCa, 'setConc' )
adaptCa.outputOffset = 0.0001 # 100 nM offset in chem.
adaptCa.scale = 0.05 # 0.06 to 0.003 mM
adaptGluR = moose.Adaptor( '/model/chem/psdMesh/adaptGluR' )
chemR = moose.element( '/model/chem/psdMesh/psdGluR' )
# Here we connect up the chem adaptors to only 3 of the spine
# heads in the elec model, just to make it interesting.
elec1R = moose.element( '/model/elec/head1/gluR' )
elec2R = moose.element( '/model/elec/head2/gluR' )
elec3R = moose.element( '/model/elec/head3/gluR' )
moose.connect( adaptGluR, 'requestOut', chemR, 'getN', 'OneToAll' )
moose.connect( adaptGluR, 'output', elec1R, 'setGbar', 'OneToAll' )
moose.connect( adaptGluR, 'output', elec2R, 'setGbar', 'OneToAll' )
moose.connect( adaptGluR, 'output', elec3R, 'setGbar', 'OneToAll' )
adaptGluR.outputOffset = 1e-9 # pS
adaptGluR.scale = 1e-8 / 100 # from n to pS
adaptK = moose.Adaptor( '/model/chem/neuroMesh/adaptK' )
chemK = moose.element( '/model/chem/neuroMesh/kChan' )
elecK = moose.element( '/model/elec/compt/K' )
moose.connect( adaptK, 'requestOut', chemK, 'getConc', 'OneToAll' )
moose.connect( adaptK, 'output', elecK, 'setGbar', 'OneToAll' )
adaptK.scale = 0.3 # from mM to Siemens
def makeChemPlots():
graphs = moose.Neutral( '/graphs' )
addPlot( '/model/chem/psdMesh/psdGluR', 'getN', 'psd0R' )
addPlot( '/model/chem/spineMesh/Ca', 'getConc', 'spine0Ca' )
addPlot( '/model/chem/neuroMesh/Ca', 'getConc', 'dend0Ca' )
addPlot( '/model/chem/neuroMesh/kChan_p', 'getConc', 'kChan_p' )
addPlot( '/model/chem/neuroMesh/kChan', 'getConc', 'kChan' )
addPlot( '/model/chem/neuroMesh/Ca.kinase', 'getConc', 'dendKinase' )
addPlot( '/model/chem/spineMesh/toPsd', 'getConc', 'toPsd0' )
#addPlot( '/n/neuroMesh/Ca', 'getConc', 'dendCa' )
#addPlot( '/n/neuroMesh/inact_kinase', 'getConc', 'inactDendKinase' )
#addPlot( '/n/psdMesh/psdGluR', 'getN', 'psdGluR' )
def testCubeMultiscale( useSolver ):
elecDt = 10e-6
chemDt = 1e-4
plotDt = 5e-4
plotName = 'mc.plot'
if ( useSolver ):
elecDt = 50e-6
chemDt = 2e-3
plotName = 'mcs.plot'
makeCubeMultiscale()
makeChemPlots()
makeElecPlots()
moose.setClock( 0, elecDt )
moose.setClock( 1, elecDt )
moose.setClock( 2, elecDt )
moose.setClock( 5, chemDt )
moose.setClock( 6, chemDt )
moose.setClock( 7, plotDt )
moose.setClock( 8, plotDt )
moose.useClock( 0, '/model/elec/##[ISA=Compartment]', 'init' )
moose.useClock( 1, '/model/elec/##[ISA=Compartment],/model/elec/##[ISA=SpikeGen]', 'process' )
moose.useClock( 2, '/model/elec/##[ISA=SynBase],/model/elec/##[ISA=ChanBase],/model/elec/##[ISA=CaConc]','process')
moose.useClock( 5, '/model/##[ISA=ReacBase],/model/##[ISA=EnzBase]', 'process' )
moose.useClock( 6, '/model/##[ISA=PoolBase],/model/chem/##[ISA=Adaptor]', 'process' )
moose.useClock( 7, '/graphs/#', 'process' )
moose.useClock( 8, '/graphs/elec/#', 'process' )
if ( useSolver ):
makeSolvers( elecDt )
moose.reinit()
moose.start( 1.0 )
runtime = 1.0
dumpPlots( plotName,runtime )
[docs]def main():
"""
A toy compartmental neuronal + chemical model. The neuronal model is in
a dendrite and five dendritic spines. The chemical model is in three
compartments: one for the dendrite, one for the spine head, and one for the postsynaptic density.
However, the spatial geometry of the neuronal model is ignored and the chemical
model just has three cubic volumes for each compartment. So there
is a functional mapping but spatial considerations are lost.
The electrical model contributes the incoming calcium flux to the
chemical model. This comes from the synaptic channels.
The signalling here does two things to the electrical model.
- First, the amount of receptor in the chemical model controls the amount of glutamate
receptor in the PSD.
- Second, there is a small kinase reaction that phosphorylates and inactivates the dendritic potassium channel.
"""
testCubeMultiscale( 1 )
if __name__ == '__main__':
main()
#
# HsolveInstability.py ends here