def loadtree(hdfnode, moosenode):
    """Load the element tree saved under the group `hdfnode` into `moosenode`"""    
    pathclass = {}    
    basepath = hdfnode.attr['path']
    if basepath != '/':
        basepath = basepath + '/'
    emdata = hdfnode['ematrix'][:]
    sorted_paths = sorted(emdata['path'], key=lambda x: x.count('/'))
    dims = dict(emdata['path', 'dims'])
    classes = dict(emdata['path', 'class'])
    current = moose.getCwe()
    moose.setCwe(moosenode)
    # First create all the ematrices
    for path in sorted_paths:
        rpath = path[len(basepath):]
        classname = classes[path]
        shape = dims[path]
        em = moose.ematrix(rpath, shape, classname)
    wfields = {}
    for cinfo in moose__.element('/classes').children:
        cname = cinfo[0].name
        wfields[cname] = [f[len('set_'):] for f in moose__.getFieldNames(cname, 'destFinfo') 
                          if f.startswith('set_') and f not in ['set_this', 'set_name', 'set_lastDimension', 'set_runTime'] and not f.startswith('set_num_')]
        for key in hdfnode['/elements']:
            dset = hdfnode['/elements/'][key][:]
            fieldnames = dset.dtype.names
            for ii in range(len(dset)):
                obj = moose__.element(dset['path'][ii][len(basepath):])
                for f in wfields[obj.class_]:
                    obj.setField(f, dset[f][ii])

def make_Na( name ):
    if moose.exists( '/library/' + name ):
        return
    Na = moose.HHChannel( '/library/' + name )
    Na.Ek = 0.055             #    V
    Na.Gbar = 300 * SOMA_A    #    S
    Na.Gk = 0                 #    S
    Na.Xpower = 2
    Na.Ypower = 1
    Na.Zpower = 0

    xgate = moose.element( Na.path + '/gateX' )
    xA = np.array( [ 320e3 * (0.0131 + EREST_ACT),
        -320e3, -1.0, -1.0 * (0.0131 + EREST_ACT), -0.004, 
        -280e3 * (0.0401 + EREST_ACT), 280e3, -1.0, 
        -1.0 * (0.0401 + EREST_ACT), 5.0e-3, 
        3000, -0.1, 0.05 ] )
    xgate.alphaParms = xA

    ygate = moose.element( Na.path + '/gateY' )
    yA = np.array( [ 128.0, 0.0, 0.0, -1.0 * (0.017 + EREST_ACT), 0.018,
        4.0e3, 0.0, 1.0, -1.0 * (0.040 + EREST_ACT), -5.0e-3, 
        3000, -0.1, 0.05 ] )
    ygate.alphaParms = yA
    return Na

def make_LCa( name = 'LCa', parent = '/library' ):
        if moose.exists( parent + '/' + name ):
                return
        Ca = moose.HHChannel( parent + '/' + name )
        Ca.Ek = ECA
        Ca.Gbar = 0
        Ca.Gk = 0
        Ca.Xpower = 2
        Ca.Ypower = 1
        Ca.Zpower = 0

        xgate = moose.element( parent + '/' + name + '/gateX' )
        xA = np.array( [ 1.6e3, 0, 1.0, -1.0 * (0.065 + EREST_ACT), -0.01389, -20e3 * (0.0511 + EREST_ACT), 20e3, -1.0, -1.0 * (0.0511 + EREST_ACT), 5.0e-3, 3000, -0.1, 0.05 ] )
        xgate.alphaParms = xA
        ygate = moose.element( parent + '/' + name + '/gateY' )
        ygate.min = -0.1
        ygate.max = 0.05
        ygate.divs = 3000
        yA = np.zeros( (ygate.divs + 1), dtype=float)
        yB = np.zeros( (ygate.divs + 1), dtype=float)


#Fill the Y_A table with alpha values and the Y_B table with (alpha+beta)
        dx = (ygate.max - ygate.min)/ygate.divs
        x = ygate.min
        for i in range( ygate.divs + 1 ):
                if ( x > EREST_ACT):
                        yA[i] = 5.0 * math.exp( -50 * (x - EREST_ACT) )
                else:
                        yA[i] = 5.0
                yB[i] = 5.0
                x += dx
        ygate.tableA = yA
        ygate.tableB = yB
        return Ca

def makeAdaptors():
    ##################################################################
    # set up adaptor for elec model Ca -> chem model Ca
    # Here it is easy because we don't have to deal with different
    # sizes of electrical and chemical compartments. 
    adaptCa = moose.Adaptor( '/model/chem/kinetics/adaptCa' )
    chemCa = moose.element( '/model/chem/kinetics/Ca' )
    elecCa = moose.element( '/model/elec/soma/Ca_conc' )
    moose.connect( elecCa, 'concOut', adaptCa, 'input' )
    moose.connect( adaptCa, 'output', chemCa, 'setConc' )
    adaptCa.inputOffset = 0.0    # 
    adaptCa.outputOffset = 0.00008    # 80 nM offset in chem.
    adaptCa.scale = 0.0008   

    # set up adaptor for chem model chan -> elec model chan.
    adaptChan = moose.Adaptor( '/model/chem/kinetics/adaptChan' )
    chemChan = moose.element( '/model/chem/kinetics/chan' )
    elecChan = moose.element( '/model/elec/soma/K_A' )
    # The Adaptor has to request the output conc of the chemical pool,
    # since there isn't an output message to deliver this value.
    moose.connect( adaptChan, 'requestOut', chemChan, 'getConc' )
    moose.connect( adaptChan, 'output', elecChan, 'setGbar' )
    adaptChan.inputOffset = 0.0    # 
    adaptChan.outputOffset = 0.0
    adaptChan.scale = 1e-5    #

def create_hhcomp(parent='/library', name='hhcomp', diameter=-30e-6, length=0.0):
    """Create a compartment with Hodgkin-Huxley type ion channels (Na and
    K).

    Returns a 3-tuple: (compartment, nachannel, kchannel)

    """
    comp, sarea = create_passive_comp(parent, name, diameter, length)
    if moose.exists('/library/na'):
        moose.copy('/library/na', comp.path, 'na')
    else:
        create_na_chan(parent=comp.path)
    na = moose.element('%s/na' % (comp.path))
    # Na-conductance 120 mS/cm^2
    na.Gbar = 120e-3 * sarea * 1e4 
    na.Ek = 115e-3 + EREST_ACT
    moose.connect(comp, 'channel', na, 'channel')
    if moose.exists('/library/k'):
        moose.copy('/library/k', comp.path, 'k')
    else:
        create_k_chan(parent=comp.path)
    k = moose.element('%s/k' % (comp.path))
    # K-conductance 36 mS/cm^2
    k.Gbar = 36e-3 * sarea * 1e4 
    k.Ek = -12e-3 + EREST_ACT
    moose.connect(comp, 'channel', k, 'channel')
    return comp, na, k

def deliverStim(currTime):
	global injectionCurrent	
	global spineVm
	global somaVm
	if numpy.fabs( currTime - baselineTime ) < frameRunTime/2.0 :
		#start
		eList = moose.wildcardFind( '/model/elec/#soma#' )
		assert( len(eList) > 0 )
		eList[0].inject = injectionCurrent
		#print "1. injected current = ", injectionCurrent
		injectionCurrent += deltaCurrent
		#print "del stim first ", moose.element('/clock').currentTime
	if numpy.fabs( currTime - baselineTime - currPulseTime) < frameRunTime/2.0 :
		#end
		eList = moose.wildcardFind( '/model/elec/#soma#' )
		assert( len(eList) > 0 )
		eList[0].inject = 0.0
		#print "2. injected current = ", injectionCurrent
		#print "del stim second ", moose.element('/clock').currentTime
	if runtime - currTime < frameRunTime * 2.0 :
		#print "3. reinit-ing"
		somaVm.append( moose.element( '/graphs/VmTab' ).vector )
		spineVm.append( moose.element( '/graphs/eSpineVmTab' ).vector )
		iList.append(injectionCurrent)
		if injectionCurrent < maxCurrent :
			moose.reinit()	

def make_Na():
	if moose.exists( 'Na' ):
		return
	Na = moose.HHChannel( 'Na' )
	Na.Ek = ENA				#	V
	Na.Gbar = 300 * SOMA_A	#	S
	Na.Gk = 0				#	S
	Na.Xpower = 2
	Na.Ypower = 1
	Na.Zpower = 0

	xgate = moose.element( 'Na/gateX' )
	xA = numpy.array( [ 320e3 * (0.0131 + EREST_ACT),
		-320e3, -1.0, -1.0 * (0.0131 + EREST_ACT), -0.004, 
		-280e3 * (0.0401 + EREST_ACT), 280e3, -1.0, 
		-1.0 * (0.0401 + EREST_ACT), 5.0e-3, 
		3000, -0.1, 0.05 ] )
	xgate.alphaParms = xA


	#xgate.alpha( 320e3 * (0.0131 + EREST_ACT), -320e3, -1.0, -1.0 * (0.0131 + EREST_ACT), -0.004 )
	#xgate.beta( -280e3 * (0.0401 + EREST_ACT), 280e3, -1.0, -1.0 * (0.0401 + EREST_ACT), 5.0e-3 )

	ygate = moose.element( 'Na/gateY' )
	yA = numpy.array( [ 128.0, 0.0, 0.0, -1.0 * (0.017 + EREST_ACT), 0.018,
		4.0e3, 0.0, 1.0, -1.0 * (0.040 + EREST_ACT), -5.0e-3, 
		3000, -0.1, 0.05 ] )
	ygate.alphaParms = yA

def spinetabs(model,neuron,comps='all'):
    if not moose.exists(DATA_NAME):
        moose.Neutral(DATA_NAME)
    #creates tables of calcium and vm for spines
    spcatab = defaultdict(list)
    spvmtab = defaultdict(list)
    for typenum,neurtype in enumerate(neuron.keys()):
        if type(comps)==str and comps in {'*', 'all'}:
            spineHeads=[moose.wildcardFind(neurtype+'/##/#head#[ISA=CompartmentBase]')]
        else:
            spineHeads=[moose.wildcardFind(neurtype+'/'+c+'/#head#[ISA=CompartmentBase]') for c in comps]
        for spinelist in spineHeads:
            for spinenum,spine in enumerate(spinelist):
                compname = spine.parent.name
                sp_num=spine.name.split(NAME_HEAD)[0]
                spvmtab[typenum].append(moose.Table(vm_table_path(neurtype, spine=sp_num, comp=compname)))
                log.debug('{} {} {}',spinenum, spine.path, spvmtab[typenum][-1].path)
                moose.connect(spvmtab[typenum][-1], 'requestOut', spine, 'getVm')
                if model.calYN:
                    for child in spine.children:
                        if child.className == "CaConc" or  child.className == "ZombieCaConc" :
                            spcatab[typenum].append(moose.Table(DATA_NAME+'/%s_%s%s'% (neurtype,sp_num,compname)+child.name))
                            spcal = moose.element(spine.path+'/'+child.name)
                            moose.connect(spcatab[typenum][-1], 'requestOut', spcal, 'getCa')
                        elif child.className == 'DifShell':
                            spcatab[typenum].append(moose.Table(DATA_NAME+'/%s_%s%s'% (neurtype,sp_num,compname)+child.name))
                            spcal = moose.element(spine.path+'/'+child.name)
                            moose.connect(spcatab[typenum][-1], 'requestOut', spcal, 'getC')
    return spcatab,spvmtab

def storePlotMsgs( tgraphs,f):
        s = ""
        if tgraphs:
                for graph in tgraphs:
                        slash = graph.path.find('graphs')
                        if not slash > -1:
                                slash = graph.path.find('graph_0')
                        if slash > -1:
                                conc = graph.path.find('conc')
                                if conc > -1 :
                                        tabPath = graph.path[slash:len(graph.path)]
                                else:
                                        slash1 = graph.path.find('/',slash)
                                        tabPath = "/graphs/conc1" +graph.path[slash1:len(graph.path)]

                                if len(element(graph).msgOut):
                                        poolPath = (element(graph).msgOut)[0].e2.path
                                        poolEle = element(poolPath)
                                        poolName = poolEle.name
                                        bgPath = (poolEle.path+'/info')
                                        bg = Annotator(bgPath).color
                                        bg = getColorCheck(bg,GENESIS_COLOR_SEQUENCE)
                                        tabPath = re.sub("\[[0-9]+\]", "", tabPath)
                                        s = s+"addmsg /kinetics/" + trimPath( poolEle ) + " " + tabPath + \
                                                " PLOT Co *" + poolName + " *" + bg +"\n";
        f.write(s)

def dumpPlots( fname ):
    if ( os.path.exists( fname ) ):
        os.remove( fname )
    for x in moose.wildcardFind( '/graphs/#[ISA=Table]' ):
        moose.element( x[0] ).xplot( fname, x[0].name )
    for x in moose.wildcardFind( '/graphs/elec/#[ISA=Table]' ):
        moose.element( x[0] ).xplot( fname, x[0].name + '_e' )

def setup_hdf5_output(model, neuron, filename=None, compartments=DEFAULT_HDF5_COMPARTMENTS):
    # Make sure /hdf5 exists
    if not moose.exists(HDF5WRITER_NAME):
        print('creating', HDF5WRITER_NAME)
        writer = moose.HDF5DataWriter(HDF5WRITER_NAME)
        #writer = moose.NSDFWriter(HDF5WRITER_NAME)
        writer.mode = 2 # Truncate existing file
        if filename is not None:
            writer.filename = filename
        moose.useClock(8, HDF5WRITER_NAME, 'process')
    else:
        print('using', HDF5WRITER_NAME)
        writer = moose.element(HDF5WRITER_NAME)

    for typenum,neur_type in enumerate(neuron.keys()):
        for ii,compname in enumerate(compartments):  #neur_comps):
            comp=moose.element(neur_type+'/'+compname)
            moose.connect(writer, 'requestOut', comp, 'getVm')

    if model.calYN:
        for typenum,neur_type in enumerate(neuron.keys()):
            for ii,compname in enumerate(compartments):  #neur_comps):
                comp=moose.element(neur_type+'/'+compname)
                for child in comp.children:
                    if child.className in {"CaConc", "ZombieCaConc"}:
                        cal = moose.element(comp.path+'/'+child.name)
                        moose.connect(writer, 'requestOut', cal, 'getCa')
                    elif  child.className == 'DifShell':
                        cal = moose.element(comp.path+'/'+child.name)
                        moose.connect(writer, 'requestOut', cal, 'getC')
    return writer

    def installCellFromProtos( self ):
        if self.stealCellFromLibrary:
            moose.move( self.elecid, self.model )
            if self.elecid.name != 'elec':
                self.elecid.name = 'elec'
        else:
            moose.copy( self.elecid, self.model, 'elec' )
            self.elecid = moose.element( self.model.path + '/elec' )
            self.elecid.buildSegmentTree() # rebuild: copy has happened.
        if hasattr( self, 'chemid' ):
            self.validateChem()
            if self.stealCellFromLibrary:
                moose.move( self.chemid, self.model )
                if self.chemid.name != 'chem':
                    self.chemid.name = 'chem'
            else:
                moose.copy( self.chemid, self.model, 'chem' )
                self.chemid = moose.element( self.model.path + '/chem' )

        ep = self.elecid.path
        somaList = moose.wildcardFind( ep + '/#oma#[ISA=CompartmentBase]' )
        if len( somaList ) == 0:
            somaList = moose.wildcardFind( ep + '/#[ISA=CompartmentBase]' )
        if len( somaList ) == 0:
            raise BuildError( "installCellFromProto: No soma found" )
        maxdia = 0.0
        for i in somaList:
            if ( i.diameter > maxdia ):
                self.soma = i

def test_crossing_single():
    """This function creates an ematrix of two PulseGen elements and
    another ematrix of two Table elements.

    The two pulsegen elements have same amplitude but opposite phase.

    Table[0] is connected to PulseGen[1] and Table[1] to Pulsegen[0].

    In the plot you should see two square pulses of opposite phase.

    """
    size = 2
    pg = moose.PulseGen('pulsegen', size)
    for ix, ii in enumerate(pg.vec):
        pulse = moose.element(ii)
        pulse.delay[0] = 1.0
        pulse.width[0] = 2.0
        pulse.level[0] = (-1)**ix
    tab = moose.Table('table', size)
    moose.connect(tab.vec[0], 'requestOut', pg.vec[1], 'getOutputValue', 'Single')
    moose.connect(tab.vec[1], 'requestOut', pg.vec[0], 'getOutputValue', 'Single')
    print 'Neighbors:'
    for t in tab.vec:
        print t.path
        for n in moose.element(t).neighbors['requestOut']:
            print 'requestOut <-', n.path
    moose.setClock(0, 0.1)
    moose.useClock(0, '/##', 'process')
    moose.start(5)
    for ii in tab.vec:
        t = moose.Table(ii).vector
        print len(t)
        pylab.plot(t)
    pylab.show()

 def getChannels(self, root='/library'):
     if isinstance(root, str):
         root = moose.element(root)
     for chan in moose.wildcardFind('%s/#[ISA=HHChannel]' % (root.path)):
         channel = moose.element(chan)
         self.channels[channel.name] = channel
     return self.channels

def main():
        """
        This example illustrates loading, and running a kinetic model 
        for a bistable positive feedback system, defined in kkit format. 
        This is based on Bhalla, Ram and Iyengar, Science 2002.

        The core of this model is a positive feedback loop comprising of
        the MAPK cascade, PLA2, and PKC. It receives PDGF and Ca2+ as 
        inputs.

        This model is quite a large one and due to some stiffness in its
        equations, it runs somewhat slowly. 

        The simulation illustrated here shows how the model starts out in
        a state of low activity. It is induced to 'turn on' when a 
        a PDGF stimulus is given for 400 seconds. 
        After it has settled to the new 'on' state, model is made to 
        'turn off'
        by setting the system calcium levels to zero for a while. This
        is a somewhat unphysiological manipulation!
        """
        solver = "gsl"  # Pick any of gsl, gssa, ee..
        #solver = "gssa"  # Pick any of gsl, gssa, ee..
	mfile = '../../genesis/acc35.g'
	runtime = 2000.0
	if ( len( sys.argv ) == 2 ):
                solver = sys.argv[1]
	modelId = moose.loadModel( mfile, 'model', solver )
        # Increase volume so that the stochastic solver gssa 
        # gives an interesting output
        compt = moose.element( '/model/kinetics' )
        compt.volume = 5e-19 

	moose.reinit()
	moose.start( 500 ) 
        moose.element( '/model/kinetics/PDGFR/PDGF' ).concInit = 0.0001
	moose.start( 400 ) 
        moose.element( '/model/kinetics/PDGFR/PDGF' ).concInit = 0.0
	moose.start( 2000 ) 
        moose.element( '/model/kinetics/Ca' ).concInit = 0.0
	moose.start( 500 ) 
        moose.element( '/model/kinetics/Ca' ).concInit = 0.00008
	moose.start( 2000 ) 

	# Display all plots.
        img = mpimg.imread( 'mapkFB.png' )
        fig = plt.figure( figsize=(12, 10 ) )
        png = fig.add_subplot( 211 )
        imgplot = plt.imshow( img )
        ax = fig.add_subplot( 212 )
	x = moose.wildcardFind( '/model/#graphs/conc#/#' )
        t = numpy.arange( 0, x[0].vector.size, 1 ) * x[0].dt
        ax.plot( t, x[0].vector, 'b-', label=x[0].name )
        ax.plot( t, x[1].vector, 'c-', label=x[1].name )
        ax.plot( t, x[2].vector, 'r-', label=x[2].name )
        ax.plot( t, x[3].vector, 'm-', label=x[3].name )
        plt.ylabel( 'Conc (mM)' )
        plt.xlabel( 'Time (seconds)' )
        pylab.legend()
        pylab.show()

def creaetHHComp(parent='/library', name='hhcomp', diameter=1e-6, length=1e-6):
    """Create a compartment with Hodgkin-Huxley type ion channels (Na and
    K).

    Returns a 3-tuple: (compartment, nachannel, kchannel)

    """
    compPath = '{}/{}'.format(parent, name)
    mc = MooseCompartment( compPath, length, diameter, {})
    c = mc.mc_
    sarea = mc.surfaceArea

    if moose.exists('/library/na'):
        moose.copy('/library/na', c.path, 'na')
    else:
        create_na_chan(parent = c.path)

    na = moose.element('%s/na' % (c.path))

    # Na-conductance 120 mS/cm^2
    na.Gbar = 120e-3 * sarea * 1e4
    na.Ek = 115e-3 + EREST_ACT

    moose.connect(c, 'channel', na, 'channel')
    if moose.exists('/library/k'):
        moose.copy('/library/k', c.path, 'k')
    else:
        create_k_chan(parent = c.path)

    k = moose.element('%s/k' % (c.path))
    # K-conductance 36 mS/cm^2
    k.Gbar = 36e-3 * sarea * 1e4
    k.Ek = -12e-3 + EREST_ACT
    moose.connect(c, 'channel', k, 'channel')
    return (c, na, k)

def main():
    """ This example illustrates loading, running of an SBML model defined in XML format.\n
	The model 00001-sbml-l3v1.xml is taken from l3v1 SBML testcase.\n
	Plots are setup.\n
	Model is run for 20sec.\n
	As a general rule we created model under '/path/model' and plots under '/path/graphs'.\n
    """

    mfile =  os.path.join( script_dir, 'chem_models/00001-sbml-l3v1.xml')
    runtime = 20.0
        
    # Loading the sbml file into MOOSE, models are loaded in path/model
    sbmlId = moose.readSBML(mfile,'sbml')
    

    s1 = moose.element('/sbml/model/compartment/S1')
    s2= moose.element('/sbml/model/compartment/S2')
                      
    # Creating MOOSE Table, Table2 is for the chemical model
    graphs = moose.Neutral( '/sbml/graphs' )
    outputs1 = moose.Table2 ( '/sbml/graphs/concS1')
    outputs2 = moose.Table2 ( '/sbml/graphs/concS2')

    # connect up the tables
    moose.connect( outputs1,'requestOut', s1, 'getConc' );
    moose.connect( outputs2,'requestOut', s2, 'getConc' );

        
    # Reset and Run
    moose.reinit()
    moose.start(runtime)

def create_compartment(path):
    comp = moose.Compartment(path)
    comp.diameter = soma_dia
    comp.Em = EREST_ACT + 10.613e-3
    comp.initVm = EREST_ACT
    sarea = np.pi * soma_dia * soma_dia
    comp.Rm = 1/(0.3e-3 * 1e4 * sarea)
    comp.Cm = 1e-6 * 1e4 * sarea
    if moose.exists('/library/na'):
        nachan = moose.element(moose.copy('/library/na', comp, 'na'))
    else:
        nachan = create_na_chan(comp.path)
    nachan.Gbar = 120e-3 * sarea * 1e4
    moose.showfield(nachan)
    moose.connect(nachan, 'channel', comp, 'channel')
    if moose.exists('/library/k'):
        kchan = moose.element(moose.copy('/library/k', comp, 'k'))
    else:
        kchan = create_k_chan(comp.path)
    kchan.Gbar = 36e-3 * sarea * 1e4
    moose.connect(kchan, 'channel', comp, 'channel')
    synchan = moose.SynChan(comp.path + '/synchan')
    synchan.Gbar = 1e-8
    synchan.tau1 = 2e-3
    synchan.tau2 = 2e-3
    synchan.Ek = 0.0
    m = moose.connect(comp, 'channel', synchan, 'channel')
    spikegen = moose.SpikeGen(comp.path + '/spikegen')
    spikegen.threshold = 0.0
    m = moose.connect(comp, 'VmOut', spikegen, 'Vm')
    return comp

def make_prototype(passive=True):
    path = '%s/rc19' % (library.path)
    pfile = 'rc19.p'
    try:
        return moose.element(path)
    except ValueError:
        pass
    if not passive:
        make_na()
        make_kv()
        make_km()
        make_kca()
        make_cat()
        make_cahva()
        make_h()
    try:
        proto = moose.element(path)
    except ValueError:
        print('Loading model %s to %s' % (pfile, path))
        proto = moose.loadModel(pfile, path, 'ee') # hsolve is not functional yet
        for comp in proto[0].children:
            comp.initVm = -75e-3
        for chan in moose.wildcardFind('%s/##[ISA=HHChannel]'):
            chan.Gbar *= tadj
    return proto

def create_population(container, size):
    """Create a population of `size` single compartmental neurons with
    Na and K channels. Also create SpikeGen objects and SynChan
    objects connected to these which can act as plug points for
    setting up synapses later."""
    path = container.path
    # Contrast this with
    # comps = moose.vec(path+'/soma', size, 'Compartment')
    comps = [create_compartment(path+'/soma_%d' % (ii)) for ii in range(size)]
    spikegens = []
    synchans = []
    Em = EREST_ACT + 10.613e-3
    initVm_array = [EREST_ACT] * size
    Em_array = [Em] * size
    # initVm_array = np.random.normal(EREST_ACT, np.abs(EREST_ACT) * 0.1, size)
    # Em_array = np.random.normal(Em, np.abs(Em) * 0.1, size)
    for comp, initVm, Em in zip(comps, initVm_array, Em_array):
        comp.Em = Em
        comp.initVm = initVm
        synchan = moose.element(comp.path + '/synchan')
        synchans.append(synchan)
        spikegen = moose.element(comp.path + '/spikegen')
        spikegens.append(spikegen)

    return {'compartment': comps,
            'spikegen': spikegens,
            'synchan': synchans}