diff --git a/.github/workflows/ebrains-mirror.yml b/.github/workflows/ebrains-mirror.yml new file mode 100644 index 000000000..e4d84b622 --- /dev/null +++ b/.github/workflows/ebrains-mirror.yml @@ -0,0 +1,33 @@ +name: Mirror to Ebrains + +# Configure the events that are going to trigger tha automated update of the mirror +on: + push: + branches: [ master ] + +# Configure what will be updated +jobs: + # set the job name + to_ebrains: + runs-on: ubuntu-latest + steps: + # this task will push the master branch of the source_repo (github) to the + # destination_repo (ebrains gitlab) + - name: syncmaster + uses: wei/git-sync@v3 + # component owners need to set their own variables + # the destination_repo format is + # https://gitlab_service_account_name:${{ secrets.EBRAINS_GITLAB_ACCESS_TOKEN }}@gitlab.ebrains.eu/name_of_mirror.git + with: + source_repo: "suny-downstate-medical-center/netpyne" + source_branch: "master" + destination_repo: "https://github-pusher:${{ secrets.EBRAINS_GITLAB_ACCESS_TOKEN }}@gitlab.ebrains.eu/vbragin/netpyne.git" + destination_branch: "master" + # this task will push all tags from the source_repo to the destination_repo + - name: synctags + uses: wei/git-sync@v3 + with: + source_repo: "suny-downstate-medical-center/netpyne" + source_branch: "refs/tags/*" + destination_repo: "https://github-pusher:${{ secrets.EBRAINS_GITLAB_ACCESS_TOKEN }}@gitlab.ebrains.eu/vbragin/netpyne.git" + destination_branch: "refs/tags/*" diff --git a/.gitignore b/.gitignore index f7a7de89d..a95b201e3 100644 --- a/.gitignore +++ b/.gitignore @@ -1,4 +1,5 @@ -/doc/source/*.rst +/doc/source/netpyne.*.rst +doc/source/package_reference.rst /examples/sandbox/LEMS_sandbox.xml /examples/M1/LEMS_M1.xml diff --git a/CHANGES.md b/CHANGES.md index fd15f67a8..c9871d8b1 100644 --- a/CHANGES.md +++ b/CHANGES.md @@ -1,3 +1,63 @@ +# Version in development + +**New features** + +- Raster plot colored by phase + +- Examples based on the CA3 model using the 'colorbyPhase' option in the plotRaster + +**Bug fixes** + +- Fixed loading point cell params from legacy models (issue 607) + +- Fix voltage movie tutorial + +- Fix to automatically include netstims in the sim.allSimData object when plotRaster 'include' selects 'all' + +- Fixed loading netParams in some scenarios (bug caused by srting functions pre-processing) + +# Version 1.0.5 + +**New features** + +- Time series and PSD plots for current source density (CSD) + +- Added colorbar to CSD plot + +- Support for Sun Grid Engine HPC + +- Extended sim.gatherData() with more optional arguments for flexibility + +- Specify linear/log scale in `plotRatePSD()` + +- Print more info about exceptions in plotting functions + +- Check RxD specification for potential syntax issues + +- Prevent zero population size in scaled-down models + +- Better messages for validation errors + +**Bug fixes** + +- Fixed errors in plotting with sim.cfg.gatherOnlySimData=True + +- Support saving and loading data in .mat and .hdf5 formats + +- Multiple fixes in `plotRatePSD()` - popColors, fft, minFreq + +- Fixed sync lines in `plotRaster()` + +- Fixed performance issue in `plotConn()` with `groupBy='pop'` (default) + +- Fixed `recordTraces` to not require more presicion than segment size (for stim and synMech) + +# Version 1.0.4.2 + +**Bug fixes** + +- Unpin matplotlib version (params copying issue fixed) + # Version 1.0.4.1 **New features** @@ -52,28 +112,6 @@ - Added ability to load PointCell from saved network -- Added MultiPlotter class to allow plotting line data on multiple axes - -- Added option to use separate axis for each PSD trace (set axis='multi') - -- Added new Batch method named SBI (Simulation Based Inference) with example folder (sbiOptim) - -- Added support for string functions in properties of cell mechanism, in cell geometry (in netParams.cellParams) - -- Added support for string functions in synMech parameters - -- Massive update of schemas (validator.py and setup.py) - -- More control over POINTER variables through synMechParams (e.g. for gap junctions) - -- Introduced cell variables in cellParams - -- Added plotRateSpectrogram to utils.py - -- Functions prepareCSD() and plotCSD() are now available - -- Updated documentation on install and about - **Bug fixes** - Fixed bug with loading CompartCell with custom mechanisms from saved network @@ -92,9 +130,9 @@ - Fixed some misinformation in reference.rst about the subconn -- Fixed bug in dipole calculation units - changed from mA to uA +- Fixed bug in dipole calculation units - changed from mA to uA -- Fixed bug in conditional logic when gathering LFP / dipoles +- Fixed bug in conditional logic when gathering LFP / dipoles - Allow tuples to specify population's cells in 'include' for plotSpike diff --git a/README.md b/README.md index df6555c8e..110eb3386 100644 --- a/README.md +++ b/README.md @@ -7,3 +7,14 @@ For more details, installation instructions, documentation, tutorials, forums, v This package is developed and maintained by Dura-Bernal lab (http://dura-bernal.org) and the Neurosim lab (www.neurosimlab.org). [![Build](https://github.com/Neurosim-lab/netpyne/actions/workflows/tests.yml/badge.svg)](https://github.com/Neurosim-lab/netpyne/actions/workflows/tests.yml) + +## Acknowledgements +This work was funded by grants from the NIH, NSF, NYS SCIRB, UK Welcome Trust and Australian Research Council. We are thankful to all the contributors that have collaborated in the development of this open source tool via GitHub. + +This project was supported by: +- HBP Brain Simulation Platform funded from the European Union’s Horizon 2020 Framework Programme for Research and Innovation under the Specific Grant Agreement No. 785907 (Human Brain Project SGA2). +- EBRAINS research infrastructure, funded from the European Union’s Horizon 2020 Framework Programme for Research and Innovation under the Specific Grant Agreement No. 945539 (Human Brain Project SGA3). + +This project has received a Voucher "Integration of NetPyNE into EBRAINS Platform (EBRnetpyne)" from the European Union’s Horizon 2020 Framework Programme for Research and Innovation under the Specific Grant Agreement: +- No. 785907 (Human Brain Project SGA2). +- No. 945539 (Human Brain Project SGA3). diff --git a/doc/build.py b/doc/build.py index ca3ee3e63..f9d89529a 100644 --- a/doc/build.py +++ b/doc/build.py @@ -74,7 +74,7 @@ shutil.rmtree('build', ignore_errors=True) # All .rst files but those listed here will be deleted during this process -keep = ['about.rst', 'advanced.rst', 'index.rst', 'install.rst', 'reference.rst', 'tutorial.rst', 'contrib.rst', 'modeling-specification-v1.0.rst'] +keep = ['about.rst', 'index.rst', 'install.rst', 'user_documentation.rst', 'tutorial.rst', 'contrib.rst', 'modeling-specification-v1.0.rst'] print('Deleting old .rst files.') for file in os.listdir('source'): @@ -94,12 +94,12 @@ # sphinx-apidoc produces a file called "netpyne" that we want to call "Package Index" print('Fixing Package Index file.') -os.system('mv source/netpyne.rst source/package_index.rst') -with open('source/package_index.rst') as f: +os.system('mv source/netpyne.rst source/package_reference.rst') +with open('source/package_reference.rst') as f: lines = f.readlines() - lines[0] = 'Package Index\n' + lines[0] = 'Package Reference\n' lines[1] = '=============\n' -with open('source/package_index.rst', 'w') as f: +with open('source/package_reference.rst', 'w') as f: f.writelines(lines) # Generate the html files diff --git a/doc/source/advanced.rst b/doc/source/advanced.rst deleted file mode 100644 index 14b7b207d..000000000 --- a/doc/source/advanced.rst +++ /dev/null @@ -1,780 +0,0 @@ -Advanced Features -======================================= - -.. _importing_cells: - -Importing externally-defined cell models ---------------------------------------- - -NetPyNE provides support for internally defining cell properties of for example Hodgkin-Huxley type cells with one or multiple compartments, or Izhikevich type cells (eg. see :ref:`tutorial`). However, it is also possible to import previously defined cells in external files eg. in hoc cell templates, or cell classes, using the ``importCellParams()`` method. This method will convert all the cell information into the required NetPyNE format. This way it is possible to make use of cells which have been implemented separately. - -The ``cellRule = netParams.importCellParams(label, conds, fileName, cellName, cellArgs={}, importSynMechs=False)`` method takes as arguments the label of the new cell rule, the name of the file where the cell is defined (either .py or .hoc files), and the name of the cell template (hoc) or class (python). Optionally, a set of arguments can be passed to the cell template/class (eg. ``{'type': 'RS'}``). If you wish to import the synaptic mechanisms parameters, you can set the ``importSynMechs=True``. The method returns the new cell rule so that it can be further modified. - - -NetPyNE contains NO built-in information about any of the cell models being imported. Importing is based on temporarily instantiating the external cell model and reading all the required information (geometry, topology, distributed mechanisms, point processes, etc.). - -Below we show example of importing 10 different cell models from external files. For each one we provide the required files as well as the NetPyNE code. Make sure you run ``nrnivmodl`` to compile the mod files for each example. The list of example cell models is: - -* :ref:`import_HH` -* :ref:`import_HH3D_hoc` -* :ref:`import_HH3D_swc` -* :ref:`import_Traub` -* :ref:`import_Mainen` -* :ref:`import_Friesen` -* :ref:`import_Izhi03a` -* :ref:`import_Izhi03b` -* :ref:`import_Izhi07a` -* :ref:`import_Izhi07b` - - -Additionally, we provide an example NetPyNE file (:download:`tut_import.py `) which imports all 10 cell models, creates a population of each type, provides background inputs, and randomly connects all cells. To run the example you also need to download all the files where cells models are defined and the mod files (see below). The resulting raster is shown below: - -.. image:: figs/tut_import_raster.png - :width: 50% - :align: center - -.. _import_HH: - -Hodgkin-Huxley model -^^^^^^^^^^^^^^^^^^^^ - -*Description:* A 2-compartment (soma and dendrite) cell with ``hh`` and ``pas`` mechanisms, and synaptic mechanisms. Defined as a Python class. - -*Required files:* -:download:`HHCellFile.py ` - -*NetPyNE Code* :: - - netParams.importCellParams( - label='PYR_HH_rule', - conds={'cellType': 'PYR', 'cellModel': 'HH'}, - fileName='HHCellFile.py', - cellName='HHCellClass', - importSynMechs=True) - - -.. _import_HH3D_hoc: - -Hodgkin-Huxley model with 3D geometry (from .hoc) -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -*Description:* A multi-compartment cell. Defined as a HOC cell template. Only the cell geometry is included. Example of importing only geometry, and then adding biophysics (``hh`` and ``pas`` channels) from NetPyNE. - -*Required files:* -:download:`geom.hoc ` - -*NetPyNE Code:* :: - - cellRule = netParams.importCellParams( - label='PYR_HH3D_hoc', - conds={'cellType': 'PYR', 'cellModel': 'HH3D'}, - fileName='geom.hoc', - cellName='E21', - importSynMechs=False) - - cellRule['secs']['soma']['mechs']['hh'] = {'gnabar': 0.12, 'gkbar': 0.036, 'gl': 0.003, 'el': -70} # soma hh mechanism - - for secName in cellRule['secs']: - cellRule['secs'][secName]['mechs']['pas'] = {'g': 0.0000357, 'e': -70} - cellRule['secs'][secName]['geom']['cm'] = 1 - -.. _import_HH3D_swc: - -Hodgkin-Huxley model with 3D geometry (from .swc) -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -*Description:* A multi-compartment cell, with imported morphology from an SWC file. Only the cell geometry is included. Example of importing only geometry, and then adding biophysics (``hh`` and ``pas`` channels) from NetPyNE. - -Importing a morphology into NetPyNE from an SWC file is simple, but NetPyNE does no testing or validation of morphologies, so you should ensure your morphology file is accurate and valid before using it in NetPyNE. - -*Required files:* -:download:`BS0284.swc ` - -*NetPyNE Code:* :: - - cellRule = netParams.importCellParams( - label='PYR_HH3D_swc', - conds={'cellType': 'PYR', 'cellModel': 'HH3D'}, - fileName='BS0284.swc', - cellName='swc_cell') - - netParams.renameCellParamsSec('PYR_HH3D_swc_rule', 'soma_0', 'soma') # rename imported section 'soma_0' to 'soma' - - for secName in cellRule['secs']: - cellRule['secs'][secName]['mechs']['pas'] = {'g': 0.0000357, 'e': -70} - cellRule['secs'][secName]['geom']['cm'] = 1 - if secName.startswith('soma'): - cellRule['secs'][secName]['mechs']['hh'] = {'gnabar': 0.12, 'gkbar': 0.036, 'gl': 0.003, 'el': -70} - - -.. _import_Traub: - -Traub model -^^^^^^^^^^^^ - -*Description:* Traub cell model defined as hoc cell template. Requires multiple mechanisms defined in mod files. Downloaded from ModelDB and modified to remove calls to figure plotting and others. The ``km`` mechanism was renamed ``km2`` to avoid collision with a different ``km`` mechanism required for the Traub cell model. Synapse added from NetPyNE. - -ModelDB link: http://senselab.med.yale.edu/ModelDB/showmodel.cshtml?model=20756 - -*Required files:* -:download:`pyr3_traub.hoc `, -:download:`ar.mod `, -:download:`cad.mod `, -:download:`cal.mod `, -:download:`cat.mod `, -:download:`k2.mod `, -:download:`ka.mod `, -:download:`kahp.mod `, -:download:`kc.mod `, -:download:`kdr.mod `, -:download:`km2.mod `, -:download:`naf.mod `, -:download:`nap.mod ` - -*NetPyNE Code:* :: - - cellRule = netParams.importCellParams( - label='PYR_Traub_rule', - conds= {'cellType': 'PYR', 'cellModel': 'Traub'}, - fileName='pyr3_traub.hoc', - cellName='pyr3') - - somaSec = cellRule['secLists']['Soma'][0] - - cellRule['secs'][somaSec]['spikeGenLoc'] = 0.5 - - -.. _import_Mainen: - -Mainen model -^^^^^^^^^^^^ - -*Description:* Mainen cell model defined as python class. Requires multiple mechanisms defined in mod files. Adapted to python from hoc ModelDB version. Synapse added from NetPyNE. - -ModelDB link: http://senselab.med.yale.edu/ModelDB/showModel.cshtml?model=2488 (old hoc version) - -*Required files:* -:download:`mainen.py `, -:download:`cadad.mod `, -:download:`kca.mod `, -:download:`km.mod `, -:download:`kv.mod `, -:download:`naz.mod `, -:download:`Nca.mod ` - -*NetPyNE Code:* :: - - netParams.importCellParams( - label='PYR_Mainen_rule', - conds={'cellType': 'PYR', 'cellModel': 'Mainen'}, - fileName='mainen.py', - cellName='PYR2') - - -.. _import_Friesen: - -Friesen model -^^^^^^^^^^^^^^ - -*Required files:* Friesen cell model defined as python class. Requires multiple mechanisms (including point processes) defined in mod files. Spike generation happens at the ``axon`` section (not the ``soma``). This is indicated in NetPyNE adding the ``spikeGenLoc`` item to the ``axon`` section entry, and specifying the section location (eg. 0.5). - -*Required files:* -:download:`friesen.py `, -:download:`A.mod `, -:download:`GABAa.mod `, -:download:`AMPA.mod `, -:download:`NMDA.mod `, -:download:`OFThpo.mod `, -:download:`OFThresh.mod `, -:download:`netcon.inc `, -:download:`ofc.inc ` - -*NetPyNE Code:* :: - - cellRule = netParams.importCellParams( - label='PYR_Friesen_rule', - conds={'cellType': 'PYR', 'cellModel': 'Friesen'}, - fileName='friesen.py', - cellName='MakeRSFCELL') - - cellRule['secs']['axon']['spikeGenLoc'] = 0.5 # spike generator location. - -.. _import_Izhi03a: - -Izhikevich 2003a model (independent voltage variable) -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -*Description:* Izhikevich, 2003 cell model defined as python class. Requires point process defined in mod file. This version is added to a section but does not employ the section voltage or synaptic mechanisms. Instead it uses its own internal voltage variable and synaptic mechanism. This is indicated in NetPyNE adding the ``vref`` item to the point process entry, and specifying the name of the internal voltage variable (``V``). - -Modeldb link: https://senselab.med.yale.edu/modeldb/showModel.cshtml?model=39948 - -*Required files:* -:download:`izhi2003Wrapper.py `, -:download:`izhi2003a.mod ` - -*NetPyNE Code:* :: - - cellRule = netParams.importCellParams( - label='PYR_Izhi03a_rule', - conds={'cellType': 'PYR', 'cellModel':'Izhi2003a'}, - fileName='izhi2003Wrapper.py', - cellName='IzhiCell', - cellArgs={'type':'tonic spiking', 'host':'dummy'}) - - cellRule['secs']['soma']['pointps']['Izhi2003a_0']['vref'] = 'V' # specify that uses its own voltage V - - -.. _import_Izhi03b: - -Izhikevich 2003b model (uses section voltage) -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -*Description:* Izhikevich, 2003 cell model defined as python class. Requires point process defined in mod file. This version is added to a section and shares the section voltage and synaptic mechanisms. A synaptic mechanism is added from NetPyNE during the connection phase. - -Modeldb link: https://senselab.med.yale.edu/modeldb/showModel.cshtml?model=39948 - -*Required files:* -:download:`izhi2003Wrapper.py `, -:download:`izhi2003b.mod ` - -*NetPyNE Code:* :: - - netParams.importCellParams( - label='PYR_Izhi03b_rule', - conds={'cellType': 'PYR', 'cellModel':'Izhi2003b'}, - fileName='izhi2003Wrapper.py', - cellName='IzhiCell', - cellArgs={'type':'tonic spiking'}) - - -.. _import_Izhi07a: - -Izhikevich 2007a model (independent voltage variable) -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -*Description:* Izhikevich, 2007 cell model defined as python clas. Requires point process defined in mod file. This version is added to a section but does not employ the section voltage or synaptic mechanisms. Instead it uses its own internal voltage variable and synaptic mechanism. This is indicated in NetPyNE adding the ``vref`` item to the point process entry, and specifying the name of the internal voltage variable (``V``). The cell model includes several internal synaptic mechanisms, which can be specified as a list in NetPyNE by adding the ``synList`` item to the point process entry. - -Modeldb link: https://senselab.med.yale.edu/modeldb/showModel.cshtml?model=39948 - -*Required files:* -:download:`izhi2007Wrapper.py `, -:download:`izhi2007a.mod ` - -*NetPyNE Code:* :: - - cellRule = netParams.importCellParams( - label='PYR_Izhi07a_rule', - conds={'cellType': 'PYR', 'cellModel':'Izhi2007a'}, - fileName='izhi2007Wrapper.py', - cellName='IzhiCell', - cellArgs={'type':'RS', 'host':'dummy'}) - - cellRule['secs']['soma']['pointps']['Izhi2007a_0']['vref'] = 'V' # specify that uses its own voltage V - - cellRule['secs']['soma']['pointps']['Izhi2007a_0']['synList'] = ['AMPA', 'NMDA', 'GABAA', 'GABAB'] # specify its own synapses - - -.. _import_Izhi07b: - -Izhikevich 2007b model (uses section voltage) -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -*Description:* Izhikevich, 2007 cell model defined as python class. Requires point process defined in mod file. This version is added to a section and shares the section voltage and synaptic mechanisms. - -Modeldb link: https://senselab.med.yale.edu/modeldb/showModel.cshtml?model=39948 - -*Required files:* -:download:`izhi2007Wrapper.py `, -:download:`izhi2007b.mod ` - -*NetPyNE Code:* :: - - netParams.importCellParams( - label='PYR_Izhi07b_rule', - conds={'cellType': 'PYR', 'cellModel':'Izhi2007b'}, - fileName='izhi2007Wrapper.py', - cellName='IzhiCell', - cellArgs={'type':'RS'}) - - -The full code to import all cell models above and create a network with them is available here: :download:`tut_import.py `. - - -Parameter Optimization of a Simple Neural Network Using An Evolutionary Algorithm ---------------------------------------------------------------------------------- - -This tutorial provides an example of how to use -\ `inspyred `__\ , -an evolutionary algorithm toolkit, to optimize parameters in our prior -\ `tut2.py `__\ \*\* -neural network--modified to remove any code relating to initiating -network simulation and output display--, such that it achieves a target -average firing rate around (~) 17 Hz. - -\*\*Some modification is required near the end of the tut2.py code, to -remove any code relating to initiating network simulation and output -display, all of which has now been handled in the new top level code -(:download:`tut_optimization.py `): - -.. code-block:: python - - # Create network and run simulation - # sim.createSimulateAnalyze(netParams = netParams, simConfig = simConfig) # line commented out - - # import pylab; pylab.show()  # if figures appear empty # line commented out - -excerpt from tut2.py - -Additional Background Reading -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -`A description of the -algorithm `__\  methodology -that will be used to optimize the simple neural network in this example. - -Introduction -^^^^^^^^^^^^^ -Using the inspyred python package to find neural network parameters so -that some property of the network (e.g. firing rate) matches a desired -target can be broken down into 3 steps. First, 1) defining a desired -target model (in this case, some measurable value) and fitness -function--fitness defined here as a calculable value that represents how -close a neural network with a given parameters matches the target. -Subsequently, it is necessary to 2) determine the appropriate neural -network parameters to modify to achieve that model/value. Finally, -3) appropriate parameters for the evolutionary algorithm are defined. -Ultimately, If the inputs to the evolutionary algorithm are appropriate, -then over successive iterations, the parameters determined by the -evolutionary algorithm should generate models closer to the target. - -Particularizing these 3 steps to our example we get: - -.. image:: figs/tut_optimization_diagram.png - :width: 80% - :align: center - -1. Defining a desired target model and fitness function. - -Defining a desired target model is largely arbitrary, some constraints -being that there must be a way to adjust parameters such that the -results are closer to the target model than before (or that fitness is -improved), and that there must be a way to evaluate the fitness of a -model with given parameters. In this case, our target model is a neural -network that achieves an average firing rate of 17 Hz. The fitness for -such a model can be defined as the difference between the average firing -rate of a certain model and the target firing rate of 17 Hz. - -2. Selecting the model parameters to be optimized. - -If a parameter can in some way alter the fitness of the final model, it -may be an appropriate candidate for optimization, depending on what the -model is seeking to achieve. As well as a host of other parameters, -altering the probability, weight or delay of the synaptic connections in -the neural network can affect the average firing rate. In this example, -we will optimize the values of the probability, weight and delay of -connections from the sensory to the motor population. - -3. Selecting appropriate parameters for the evolutionary algorithm. - -inspyred allows customization of the various components of the -evolutionary algorithm, including: - --  a selector that determines which sets of parameter values become - parents and thus which parameter values will be used to form the next - generation in the evolutionary iteration, -- a variator that determines how each current iteration of parameter - sets is formed from the previous iteration, -- a replacer which determines whether previous sets of parameter values - are brought into the next iteration, -- a terminator which defines when to end evolutionary iterations, -- an observer which allows for tracking of parameter values through - each evolutionary iteration. - -         - -Using inspyred -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -The evolutionary algorithm is implemented the ec module from the -inspyred package: - -.. code-block:: python - - from inspyred import ec # import evolutionary computation from inspyred - -excerpt from tut\_optimization.py - -ec includes a class for the evolutionary computation algorithm: -ec.EvolutionaryComputation(), which allows entering parameters to -customize the algorithm. The evolutionary algorithm involves random -processes (e.g. randomly mutating genes) and so requires random number -generator. In this case we will use python's Random() method, which we -initialize using a specific seed value so that we can reproduce the -results in the future: - -.. code-block:: python - - # create random seed for evolutionary computation algorithm - rand = Random() - rand.seed(1) - - # instantiate evolutionary computation algorithm - my_ec = ec.EvolutionaryComputation(rand) - - -excerpt from tut\_optimization.py - -Parameters for the evolutionary algorithm are then established for our -ec evolutionary computation instance by assigning various variator, -replacer, terminator and observer elements--essentially toggling -specific components of the algorithm-- to ec.selectors, ec.variators, -ec.replacers, ec.terminators, ec.observers: - -.. code-block:: python - - #toggle variators - my_ec.variator = [ec.variators.uniform_crossover, # implement uniform crossover & gaussian replacement - ec.variators.gaussian_mutation]    - my_ec.replacer = ec.replacers.generational_replacement   # implement generational replacement - - my_ec.terminator = ec.terminators.evaluation_termination # termination dictated by no. evaluations - - #toggle observers - my_ec.observer = [ec.observers.stats_observer,  # print evolutionary computation statistics -                 ec.observers.plot_observer,   # plot output of the evolutionary computation as graph -                 ec.observers.best_observer]   # print the best individual in the population to screen - -excerpt from ex_optimization.py - -where: - -+----------------------------------------+--------------------------------------+ -| ec.variators.uniform\_crossover | variator where coin flip to | -| | determine whether 'mom' or 'dad' | -| | element is inherited by offspring | -+----------------------------------------+--------------------------------------+ -| ec.variators.gaussian\_mutation | variator implements gaussian | -| | mutation which makes use of bounder | -| | function as specified | -| | in: my\_ec.evolve(...,bounder=ec.Bou | -| | nder(minParamValues, maxParamValues) | -| | ,...) | -| | | -+----------------------------------------+--------------------------------------+ -| ec.replacers.generational\_replacement | replacer implements generational | -| | replacement with elitism (as | -| | specified in | -| | my\_ec.evolve(...,num\_elites=1,...) | -| | , | -| | where the existing generation is | -| | replaced by offspring, and | -| | existing individuals | -| | will survive if they have better | -| | fitness than the offspring | -+----------------------------------------+--------------------------------------+ -| ec.terminators.evaluation\_termination | terminator runs based on the number | -| | of evaluations that have occurred | -+----------------------------------------+--------------------------------------+ -| ec.observers.stats\_observer | indicates how many of the generated | -| | individuals (parameter sets) will be | -| | selected for the next evolutionary | -| | iteration. | -+----------------------------------------+--------------------------------------+ -| ec.observers.plot\_observer | indicates the rate of mutation, or | -| | the rate at which values for each | -| | parameter (probability, weight and | -| | delay) taken from a prior generation | -| | are altered in the next generation | -+----------------------------------------+--------------------------------------+ -| ec.observers.best\_observer | sets the number of parameters that | -| | will be optimized to 3, | -| | corresponding to the length of | -| | [probability, weight, delay]. | -+----------------------------------------+--------------------------------------+ - -These predefined selector, variator, replacer, terminator and observer -elements as well as other options can be found in the \ `inspyred -documentation `__\ . - -FInally, the evolutionary computation algorithm instance includes a -method: my\_ec.evolve() , which will move through successive -evolutionary iterations evaluating different parameter sets until the -terminating condition is achieved. This function comes with multiple -arguments, with two significant arguments being the generator and -evaluator functions. A function call for  my\_ec.evolve() will look -similar to the following: - -.. code-block:: python - - # call evolution iterator - - final_pop = my_ec.evolve(generator=generate_netparams, # assign model parameter generator to iterator generator -                       evaluator=evaluate_netparams, # assign fitness function to iterator evaluator -                     pop_size=10, -                     maximize=False,                    - bounder=ec.Bounder(minParamValues, maxParamValues), - max_evaluations=50, - num_selected=10, - mutation_rate=0.2, - num_inputs=3, - num_elites=1) - - -excerpt from tut\_optimization.py - -where: - -+--------------------------------------+--------------------------------------+ -| pop\_size=10 | means that each generation of | -| | parameter sets will consist of 10 | -| | individuals | -+--------------------------------------+--------------------------------------+ -| maximize=False | means that we are taking higher | -| | fitness to correspond to minimal | -| | values in terms of difference | -| | between model firing frequency and | -| | 17 Hz | -+--------------------------------------+--------------------------------------+ -| bounder=ec.Bounder(minParamValues, | defines boundaries for each of the | -| maxParamValues) | parameters. The format to describe | -|           | the minimum and maximum values for | -| | the parameters we are seeking to | -|       | optimize: minParamValues is an array | -| | of minimum of values corresponding | -| | to [probability, weight, delay], and | -| | maxParamValues is the array of | -| | maximum values. | -+--------------------------------------+--------------------------------------+ -| max\_evaluations=50 | indicates how many parameter sets | -| | are evaluated prior termination of | -| | the evolutionary iterations | -+--------------------------------------+--------------------------------------+ -| num\_selected=10 | indicates how many of the generated | -| | individuals (parameter sets) will be | -| | selected for the next evolutionary | -| | iteration. | -+--------------------------------------+--------------------------------------+ -| mutation\_rate=0.2 | indicates the rate of mutation, or | -| | the rate at which values for each | -| | parameter (probability, weight and | -| | delay) taken from a prior generation | -| | are altered in the next generation | -+--------------------------------------+--------------------------------------+ -| num\_inputs=3 | sets the number of parameters that | -| | will be optimized to 3, | -| | corresponding to the length of | -| | [probability, weight, delay]. | -+--------------------------------------+--------------------------------------+ -| num\_elites=1 | sets the number of elites to 1. That | -| | is, one individual from the existing | -| | generation may be retained (as | -| | opposed to a complete generational | -| | replacement) if it has better | -| | fitness than an individual selected | -| | from the offspring. | -+--------------------------------------+--------------------------------------+ - -The generator and evaluator arguments expect user defined functions as -inputs, with generator used to define a population of initial parameter -value sets for the very first iteration, and evaluator being the fitness -function that will be used to evaluate each model for how close it is to -the target. In this example, the generator is a fairly straightforward -function which creates an initial set of parameter values (i.e. -[probability, weight, delay] ) by drawing from a parameterized uniform -distribution: - -.. code-block:: python - - # return a set of initialParams which contains a [probability, weight, delay] - - def generate_netparams(random, args): - -     size = args.get('num_inputs') -   initialParams = [random.uniform(minParamValues[i], maxParamValues[i]) for i in range(size)] - - return initialParams - -excerpt from tut\_optimization.py - -The fitness function involves taking a list of sets of parameter values, -i.e. : [ [ a0, b0, c0], [a1, b1, c1], [a2, b2, c2], ... , [an, bn, cn ] -] where a, b, c represent the parameter values and 1 through n -representing the individual number within the population, and -calculating a fitness score for each element of the list, which is then -returned as a list of fitness values (i.e. : [ f0, f1, f2, ... , fn ] -) corresponding to the initial sets of parameter values. It follows the -general template: - -.. code-block:: python - - def evaluate_fitness(candidates, args): -    fitness = [] -    for candidate in candidates: -        fit = some_fitness_function(candidate) -        fitness.append(fit) -    return fitness - -excerpt from tut\_optimization.py - -The actual code that is used to serve as     - some\_fitness\_function(candidate)    is described below: - -  - -Overview of the Fitness Function -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -The fitness function in this case involves 1) creating a neural network -with the given parameters, 2) simulating it to find the average firing -rate, then 3) comparing this firing rate to a target firing rate. - -1. Creating a neural network with the parameters to evaluate - -We will employ the NetPyNE defined network in tut2.py, and modify -the [probability, weight, delay] parameters. This  involves redefining -specific values found in tut2.py found within the connectivity rule -between the S and M populations:    netParams.connParams['S->M']    - -.. code-block:: python - - ## Cell connectivity rules - netParams.connParams['S->M'] = {      #  S -> M label -       'preConds': {'popLabel': 'S'},  # conditions of presyn cells -       'postConds': {'popLabel': 'M'}, # conditions of postsyn cells -       'probability': 0.5,             # probability of connection <-- to be optimized by evolutionary algorithm -       'weight': 0.01,                 # synaptic weight           <-- to be optimized by evolutionary algorithm -       'delay': 5,                     # transmission delay (ms) <-- to be optimized by evolutionary algorithm -       'synMech': 'exc'}               # synaptic mechanism - -excerpt from tut2.py - -these values are replaced in the fitness function with the parameter -values generated by the evolutionary algorithm. As the fitness function -resides within a for loop iterating through the list of candidates (     -for icand,cand in enumerate(candidates):    ), the individual parameters -can be accessed as cand[0], cand[1], and cand[2]. Reassigning values to -the parameters in tut\_optimization.py can be done via the following -line: - -.. code-block:: python - - tut2.netParams.connParams['S->M'][''] =  - -2. Simulating the created neural network and finding the average firing - rate - -Once the network parameters have been modified we can call the -sim.createSimulate() NetPyNE function to run the simulation. We will -pass as arguments the tut2 netParams and simConfig objects that we just -modified. Once the simulation has ran we will have access to the -simulation output via sim.simData.  - -:: - - # create network - sim.createSimulate(netParams=tut2.netParams, simConfig=tut2.simConfig) - -excerpt from tut\_optimization.py - -3. Comparing the average firing rate to a target average firing rate - -To calculate the average firing rate (in spikes/sec = Hz) of the -network, we divide the spikes that have occurred during the simulation, -by the number of neurons and the duration. A list of spike times and a -list of neurons can be accessed via the NetPyNE sim module: - sim.simData['spkt'] and   sim.net.cells   . These are populated after -running   sim.createSimulate()  . From these lists, getting the number -of spike times and neurons is done by using python’s   len()   function. -The duration of the simulation can be accessed in the -tut\_optimization.py code  via        tut2.simConfig.duration    .  The -calculation for average firing rate is thus as follows: - -.. code-block:: python - - # calculate firing rate - numSpikes = float(len(sim.simData['spkt'])) - numCells = float(len(sim.net.cells)) - duration = tut2.simConfig.duration/1000.0 - netFiring = numSpikes/numCells/duration - -excerpt from tut\_optimization.py - -Finally, the average firing rate of the model is compared to the target -firing rate as follows: - -.. code-block:: python - - # calculate fitness for this candidate - fitness = abs(targetFiring - netFiring)  # minimize absolute difference in firing rate - -excerpt from tut\_optimization.py - -Displaying Findings -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -The results of the evolutionary algorithm are displayed on the standard -output (terminal) as well as plotted using the matplotlib package. The -following lines are relevant to showing results of the various -candidates within the iterator: - -.. code-block:: python - - for icand,cand in enumerate(candidates): -       ... -       print '\n CHILD/CANDIDATE %d: Network with prob:%.2f, weight:%.2f, delay:%.1f \n  firing rate: %.1f, FITNESS = %.2f \n'\ -       %(icand, cand[0], cand[1], cand[2], netFiring, fitness) - -excerpt from tut\_optimization.py - -The first line:  for icand,cand in enumerate(candidates): is analogous -to the the iterator  for candidate in candidates:  used in the -pseudocode example above, except that the  enumerate() function will -also return an index--starting from 0-- for each element in the list, -and is used in the subsequent print statement. - -This example also displays the generated candidate with average -frequency closest to 17 Hz. This candidate will exist in the final -generation, and possess the best fitness score (corresponding to a -minimum difference). Since   num\_elites=1   there is no risk that a -prior generation will have a candidate with a better fitness. - -After the evolution finishes, to access the candidate with the best -fitness score, the final generation of candidates, which is returned by -the  my\_ec.evolve()   function is then sorted in reverse (least to -greatest), placing the candidate that achieves an average firing rate -closest to 17 Hz (and therefore has the minimum difference) at the start -of the list (or at position 0). We will use NetPyNE to visualize the -output of this network, by setting the optimized parameters, simulating -the network and plotting a raster plot. The code that performs this task -is isolated below: - -.. code-block:: python - - final_pop = my_ec.evolve(...) - ... - # plot raster of top solutions - final_pop.sort(reverse=True)        # sort final population so best fitness (minimum difference) is first in list - bestCand = final_pop[0].candidate   # bestCand <-- candidate in first position of list - tut2.simConfig.analysis['plotRaster'] = True                      # plotting - tut2.netParams.connParams['S->M']['probability'] = bestCand[0]    # set tut2 values to corresponding - tut2.netParams.connParams['S->M']['weight'] = bestCand[1]         # best candidate values - tut2.netParams.connParams['S->M']['delay'] = bestCand[2] - sim.createSimulateAnalyze(netParams=tut2.netParams, simConfig=tut2.simConfig) # run simulation - -excerpt from tut\_optimization.py - -The code for neural network optimization through evolutionary algorithm used in this tutorial can be found here: :download:`tut_optimization.py `. - - -.. Cell density and connectivity as a function of cell location -.. ------------------------------------------------------------ - - -.. Create population as list of individual cells -.. ------------------------------------------------ -.. (eg. measured experimentally) - - -.. Adding connectivity functions -.. ------------------------------ - - -.. Adding cell classes -.. -------------------- - diff --git a/doc/source/code/tut1.py b/doc/source/code/tut1.py index 1ec299cd3..8ca5a43f8 100644 --- a/doc/source/code/tut1.py +++ b/doc/source/code/tut1.py @@ -1,4 +1,3 @@ +import HHTut from netpyne import sim -from netParams import netParams -from cfg import cfg -sim.createSimulateAnalyze(netParams = netParams, simConfig = cfg) +sim.createSimulateAnalyze(netParams = HHTut.netParams, simConfig = HHTut.simConfig) diff --git a/doc/source/conf.py b/doc/source/conf.py index e1521168f..e5f737a02 100644 --- a/doc/source/conf.py +++ b/doc/source/conf.py @@ -67,9 +67,9 @@ # built documents. # # The short X.Y version. -version = '1.0.4.1' +version = '1.0.5' # The full version, including alpha/beta/rc tags. -release = '1.0.4.1' +release = '1.0.5' # The language for content autogenerated by Sphinx. Refer to documentation # for a list of supported languages. diff --git a/doc/source/index.rst b/doc/source/index.rst index b2a492b3d..116f65703 100644 --- a/doc/source/index.rst +++ b/doc/source/index.rst @@ -37,9 +37,8 @@ Table of Contents about install tutorial - reference - advanced - package_index + user_documentation + package_reference modeling-specification-v1.0 diff --git a/doc/source/install.rst b/doc/source/install.rst index cf7a1920d..5f4606123 100644 --- a/doc/source/install.rst +++ b/doc/source/install.rst @@ -9,31 +9,47 @@ Requirements Before installing NetPyNE, please ensure you have the following installed: -1) Python 2 or 3 (2.7, 3.6 and 3.7 are supported). Download from the `official Python website `_. Alternatively, you can download the `Anaconda Distribution `_ which also includes several data science and visualization packages. If you are using Windows, it is highly recommended that you download the Anaconda distribution. +1) Python 2 or 3 (2.7, and versions > 3.6 are supported). We recommend the package/environment manager `Anaconda `_ which comes with Python and many necessary packages already installed. Alternatively, installing Python can be handled through the `official Python website `_. If you are using Windows, see the quickstart guides for using `Windows Subsystem Linux `_ (preferred) or using Python through the `Anaconda distribution on Windows `_. 2) The ``pip`` tool for installing Python packages. See `pip installation here `_. -3) The NEURON simulator. See NEURON's `installation instructions `_. If you would like to run parallelized simulations, please ensure you install NEURON with MPI support (see also `Quick start guide `_). Note: the latest NEURON version can be installed simply via: ``pip install neuron`` +3) The NEURON simulator. See NEURON's `installation instructions `_. If you would like to run parallelized simulations, please ensure you install NEURON with MPI support (see also `Quick start guide `_). Note for Linux or Mac: the latest NEURON version on Linux or Mac can be installed simply via: ``pip install neuron``. -N.B. Windows users can make use of `Windows Subsystem Linux `_. Getting started is `relatively simple `_. +N.B. Windows users can make use of `Windows Subsystem Linux `_ to recreate a Linux environment. -N.B. Please make sure that all path definitions are included with the installation +F.A.Q. Historically, path definitions have caused issues with running Python, NEURON or NetPyNE. Using a path manager helps with resolving this (Anaconda, Miniconda, etc.). When installing be careful that all *recommended* path definitions are included with the installation. Install the latest released version of NetPyNE via pip (Recommended) ------------------------------------------------------ -Linux or Mac OS (from a terminal): ``pip install netpyne`` +Linux (or Windows through WSL) or Mac OS: -Windows (from a terminal): ``python -m pip install netpyne`` +From the terminal: + +``pip install netpyne`` + +Windows + +From Anaconda PowerShell or a user preconfigured terminal: + +``pip install netpyne`` Upgrade to the latest released version of NetPyNE via pip ------------------------------------------------------ Use this option if you already have NetPyNE installed and just want to update to the latest version. -Linux or Mac OS (from a terminal): ``pip install netpyne -U`` +Linux (or Windows through WSL) or Mac OS: + +From the terminal: + +``pip install netpyne -U`` + +Windows: + +From Anaconda PowerShell or a user preconfigured terminal: -Windows (from a terminal): ``python -m pip install -U netpyne`` +``pip install -U netpyne`` Install the development version of NetPyNE via GitHub and pip @@ -41,7 +57,7 @@ Install the development version of NetPyNE via GitHub and pip The NetPyNE package source files, as well as example models, are available via GitHub at: https://github.com/Neurosim-lab/netpyne. The following instructions will install the version in the GitHub "development" branch -- it will include some of the latest enhancements and bug fixes, but could also include temporary bugs: -1) ``git clone https://github.com/Neurosim-lab/netpyne.git`` +1) ``git clone https://github.com/suny-downstate-medical-center/netpyne.git`` 2) ``cd netpyne`` 3) ``git checkout development`` 4) ``pip install -e .`` @@ -55,7 +71,7 @@ This version can also be used by developers interested in extending the package. Use a browser-based online version of NetPyNE GUI (beta version) ------------------------------------------------------ -The NetPyNE GUI is available online at: `gui.netpyne.org `_. There is a maximum number of simultaneous users for this online version, so if you can't log in, please try again later. +The NetPyNE GUI is available online at: `v2.opensourcebrain.org `_. There is a maximum number of simultaneous users for this online version, so if you can't log in, please try again later. Note: the GUI also includes an interactive Python Jupyter Notebook (click "Python" icon at bottom-left) that you can use to directly run NetPyNE code/models (i.e. without using the actual graphical interface). diff --git a/doc/source/modeling-specification-v1.0.rst b/doc/source/modeling-specification-v1.0.rst index fe12c497c..6e3409dba 100644 --- a/doc/source/modeling-specification-v1.0.rst +++ b/doc/source/modeling-specification-v1.0.rst @@ -34,7 +34,7 @@ Authors and contributors ------------------------------- The original author of the specification is Salvador Dura-Bernal (salvadordura@gmail.com). Contributors to the specification include: -Padraig Gleeson, Joe W Graham, Eugenio Urdapilleta, Valery Bragin, William W Lytton, Michael L Hines, Ben Suter, Greg Patrick, +Padraig Gleeson, Joe W Graham, Eugenio Urdapilleta, Valery Bragin, James Chen, William W Lytton, Michael L Hines, Ben Suter, Greg Patrick, Matteo Cantarelli, Dario del Piano, Filippo Ledda, Facundo Rodriguez, Nicolas Gomez, Afonso Pinto, Craig Kelley, Adrian Quintana, Siddartha Mitra, Adam Newton, Joao Vitor, Cliff Kerr and David Kedziora. diff --git a/doc/source/tutorial.rst b/doc/source/tutorial.rst index 01d1a4f61..712a890f9 100644 --- a/doc/source/tutorial.rst +++ b/doc/source/tutorial.rst @@ -17,7 +17,7 @@ Tutorial 1: Quick and easy example ----------------------------------------- To start in an encouraging way, we will implement the simplest example possible consisting of just three lines! This will create a simple network (200 randomly connected cells), run a one-second simulation, and plot the network spiking raster plot and the voltage trace of a cell. -You will need to **download** the ``HHTut.py`` example parameter file (`download here `_ - right click and "Save as" to same folder you are working in). +You will need to **download** the ``HHTut.py`` example parameter file (`download here `_ - right click and "Save as" to same folder you are working in). The code looks like this (available here :download:`tut1.py `):: @@ -910,7 +910,7 @@ Notice how the rate initially increases as a function of connection weight, but Tutorial 9: Recording and plotting LFPs --------------------------------------- -Examples of how to record and analyze local field potentials (LFP) in single cells and networks are included in the \examples folder: `LFP recording example `_ . LFP recording also works with parallel simulations. +Examples of how to record and analyze local field potentials (LFP) in single cells and networks are included in the \examples folder: `LFP recording example `_ . LFP recording also works with parallel simulations. To record LFP just set the list of 3D locations of the LFP electrodes in the `simConfig` attribute `recordLFP` e.g. ``simConfig.recordLFP = e.g. [[50, 100, 50], [50, 200, 50]]`` (note the y coordinate represents depth, so will be represented as a negative value when plotted). The LFP signal in each electrode is obtained by summing the extracellular potential contributed by each segment of each neuron. Extracellular potentials are calculated using the "line source approximation" and assuming an Ohmic medium with conductivity sigma = 0.3 mS/mm. For more information on modeling LFPs see `Scholarpedia `_ or `this article `_ . @@ -918,13 +918,13 @@ The recorded LFP signal will be stored in ``sim.allSimData['LFP']`` as a 2D list To plot the LFP use the ``sim.analysis.plotLFP()`` method. This allows to plot for each electrode: 1) the time-resolved LFP signal ('timeSeries'), 2) the power spectral density ('PSD'), 3) the spectrogram / time-frequency profile ('spectrogram'), 4) and the 3D locations of the electrodes overlaid over the model neurons ('locations'). See :ref:`analysis_functions` for the full list of ``plotLFP()`` arguments. -The first example ( :download:`cell_lfp.py <../../examples/LFPrecording/cell_lfp.py>`) shows LFP recording for a single cell -- M1 corticospinal neuron -- with 700+ compartments (segments) and multiple somatic and dendritic ionic channels. The cell parameters are loaded from a .json file. The cell receives NetStim input to its soma via an excitatory synapse. Ten LFP electrodes are placed at both sides of the neuron at 5 different cortical depths. The soma voltage, LFP time-resolved signal and the 3D locations of the electrodes are plotted: +The `first example `_ shows LFP recording for a single cell -- M1 corticospinal neuron -- with 700+ compartments (segments) and multiple somatic and dendritic ionic channels. The cell parameters are loaded from a .json file. The cell receives NetStim input to its soma via an excitatory synapse. Ten LFP electrodes are placed at both sides of the neuron at 5 different cortical depths. The soma voltage, LFP time-resolved signal and the 3D locations of the electrodes are plotted: .. image:: figs/lfp_cell.png :width: 60% :align: center -The second example ( :download:`net_lfp.py <../../examples/LFPrecording/net_lfp.py>`) shows LFP recording for a network very similar to that shown in Tutorial 5. However, in this case, the cells have been replaced with a more realistic model: a 6-compartment M1 corticostriatal neuron with multiple ionic channels. The cell parameters are loaded from a .json file. Cell receive NetStim inputs and include excitatory and inhibitory connections. Four LFP electrodes are placed at different cortical depths. The raster plot and LFP time-resolved signal, PSD, spectrogram and 3D locations of the electrodes are plotted: +The `second example `_ shows LFP recording for a network very similar to that shown in Tutorial 5. However, in this case, the cells have been replaced with a more realistic model: a 6-compartment M1 corticostriatal neuron with multiple ionic channels. The cell parameters are loaded from a .json file. Cell receive NetStim inputs and include excitatory and inhibitory connections. Four LFP electrodes are placed at different cortical depths. The raster plot and LFP time-resolved signal, PSD, spectrogram and 3D locations of the electrodes are plotted: .. image:: figs/lfp_net.png :width: 90% diff --git a/doc/source/reference.rst b/doc/source/user_documentation.rst similarity index 71% rename from doc/source/reference.rst rename to doc/source/user_documentation.rst index 307ee7c4f..8843f4c64 100644 --- a/doc/source/reference.rst +++ b/doc/source/user_documentation.rst @@ -1,6 +1,6 @@ .. _package_reference: -Package Reference +User Documentation ======================================= Model components and structure @@ -305,7 +305,7 @@ Currently, 'rhythmic', 'evoked', 'poisson' and 'gauss' spike pattern generators * **evoked** - creates the ongoing external inputs (rhythmic) - * **start** - time of first spike. if -1, uniform distribution between startMin and startMax (ms) + * **start** - time of first spike * **inc** - increase in time of first spike; from cfg.inc_evinput (ms) @@ -564,7 +564,7 @@ Some of parameters of cells, synapses and connectivity rules can be provided usi * Contextual variables like cell location, segment position in cell, distance between connected cells etc. Such variables' names are pre-defined and specific to where the function is used (see below). -Functios as strings can be used as parameters in the entries of the following model components (in each case, there is a specific set of variables that can be contained in the function, in addition to elements listed above). +Functions as strings can be used as parameters in the entries of the following model components (in each case, there is a specific set of variables that can be contained in the function, in addition to elements listed above). * ``weight``, ``delay``, ``probability``, ``convergence`` and ``divergence`` in connectivity rules. Available variables are: @@ -1767,3 +1767,783 @@ Data saved to file * netParams * net * simData + + + +.. _importing_cells: + +Importing externally-defined cell models +--------------------------------------- + +NetPyNE provides support for internally defining cell properties of for example Hodgkin-Huxley type cells with one or multiple compartments, or Izhikevich type cells (eg. see :ref:`tutorial`). However, it is also possible to import previously defined cells in external files eg. in hoc cell templates, or cell classes, using the ``importCellParams()`` method. This method will convert all the cell information into the required NetPyNE format. This way it is possible to make use of cells which have been implemented separately. + +The ``cellRule = netParams.importCellParams(label, conds, fileName, cellName, cellArgs={}, importSynMechs=False)`` method takes as arguments the label of the new cell rule, the name of the file where the cell is defined (either .py or .hoc files), and the name of the cell template (hoc) or class (python). Optionally, a set of arguments can be passed to the cell template/class (eg. ``{'type': 'RS'}``). If you wish to import the synaptic mechanisms parameters, you can set the ``importSynMechs=True``. The method returns the new cell rule so that it can be further modified. + + +NetPyNE contains NO built-in information about any of the cell models being imported. Importing is based on temporarily instantiating the external cell model and reading all the required information (geometry, topology, distributed mechanisms, point processes, etc.). + +Below we show example of importing 10 different cell models from external files. For each one we provide the required files as well as the NetPyNE code. Make sure you run ``nrnivmodl`` to compile the mod files for each example. The list of example cell models is: + +* :ref:`import_HH` +* :ref:`import_HH3D_hoc` +* :ref:`import_HH3D_swc` +* :ref:`import_Traub` +* :ref:`import_Mainen` +* :ref:`import_Friesen` +* :ref:`import_Izhi03a` +* :ref:`import_Izhi03b` +* :ref:`import_Izhi07a` +* :ref:`import_Izhi07b` + + +Additionally, we provide an example NetPyNE file (:download:`tut_import.py `) which imports all 10 cell models, creates a population of each type, provides background inputs, and randomly connects all cells. To run the example you also need to download all the files where cells models are defined and the mod files (see below). The resulting raster is shown below: + +.. image:: figs/tut_import_raster.png + :width: 50% + :align: center + +.. _import_HH: + +Hodgkin-Huxley model +^^^^^^^^^^^^^^^^^^^^ + +*Description:* A 2-compartment (soma and dendrite) cell with ``hh`` and ``pas`` mechanisms, and synaptic mechanisms. Defined as a Python class. + +*Required files:* +:download:`HHCellFile.py ` + +*NetPyNE Code* :: + + netParams.importCellParams( + label='PYR_HH_rule', + conds={'cellType': 'PYR', 'cellModel': 'HH'}, + fileName='HHCellFile.py', + cellName='HHCellClass', + importSynMechs=True) + + +.. _import_HH3D_hoc: + +Hodgkin-Huxley model with 3D geometry (from .hoc) +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +*Description:* A multi-compartment cell. Defined as a HOC cell template. Only the cell geometry is included. Example of importing only geometry, and then adding biophysics (``hh`` and ``pas`` channels) from NetPyNE. + +*Required files:* +:download:`geom.hoc ` + +*NetPyNE Code:* :: + + cellRule = netParams.importCellParams( + label='PYR_HH3D_hoc', + conds={'cellType': 'PYR', 'cellModel': 'HH3D'}, + fileName='geom.hoc', + cellName='E21', + importSynMechs=False) + + cellRule['secs']['soma']['mechs']['hh'] = {'gnabar': 0.12, 'gkbar': 0.036, 'gl': 0.003, 'el': -70} # soma hh mechanism + + for secName in cellRule['secs']: + cellRule['secs'][secName]['mechs']['pas'] = {'g': 0.0000357, 'e': -70} + cellRule['secs'][secName]['geom']['cm'] = 1 + +.. _import_HH3D_swc: + +Hodgkin-Huxley model with 3D geometry (from .swc) +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +*Description:* A multi-compartment cell, with imported morphology from an SWC file. Only the cell geometry is included. Example of importing only geometry, and then adding biophysics (``hh`` and ``pas`` channels) from NetPyNE. + +Importing a morphology into NetPyNE from an SWC file is simple, but NetPyNE does no testing or validation of morphologies, so you should ensure your morphology file is accurate and valid before using it in NetPyNE. + +*Required files:* +:download:`BS0284.swc ` + +*NetPyNE Code:* :: + + cellRule = netParams.importCellParams( + label='PYR_HH3D_swc', + conds={'cellType': 'PYR', 'cellModel': 'HH3D'}, + fileName='BS0284.swc', + cellName='swc_cell') + + netParams.renameCellParamsSec('PYR_HH3D_swc_rule', 'soma_0', 'soma') # rename imported section 'soma_0' to 'soma' + + for secName in cellRule['secs']: + cellRule['secs'][secName]['mechs']['pas'] = {'g': 0.0000357, 'e': -70} + cellRule['secs'][secName]['geom']['cm'] = 1 + if secName.startswith('soma'): + cellRule['secs'][secName]['mechs']['hh'] = {'gnabar': 0.12, 'gkbar': 0.036, 'gl': 0.003, 'el': -70} + + +.. _import_Traub: + +Traub model +^^^^^^^^^^^^ + +*Description:* Traub cell model defined as hoc cell template. Requires multiple mechanisms defined in mod files. Downloaded from ModelDB and modified to remove calls to figure plotting and others. The ``km`` mechanism was renamed ``km2`` to avoid collision with a different ``km`` mechanism required for the Traub cell model. Synapse added from NetPyNE. + +ModelDB link: http://senselab.med.yale.edu/ModelDB/showmodel.cshtml?model=20756 + +*Required files:* +:download:`pyr3_traub.hoc `, +:download:`ar.mod `, +:download:`cad.mod `, +:download:`cal.mod `, +:download:`cat.mod `, +:download:`k2.mod `, +:download:`ka.mod `, +:download:`kahp.mod `, +:download:`kc.mod `, +:download:`kdr.mod `, +:download:`km2.mod `, +:download:`naf.mod `, +:download:`nap.mod ` + +*NetPyNE Code:* :: + + cellRule = netParams.importCellParams( + label='PYR_Traub_rule', + conds= {'cellType': 'PYR', 'cellModel': 'Traub'}, + fileName='pyr3_traub.hoc', + cellName='pyr3') + + somaSec = cellRule['secLists']['Soma'][0] + + cellRule['secs'][somaSec]['spikeGenLoc'] = 0.5 + + +.. _import_Mainen: + +Mainen model +^^^^^^^^^^^^ + +*Description:* Mainen cell model defined as python class. Requires multiple mechanisms defined in mod files. Adapted to python from hoc ModelDB version. Synapse added from NetPyNE. + +ModelDB link: http://senselab.med.yale.edu/ModelDB/showModel.cshtml?model=2488 (old hoc version) + +*Required files:* +:download:`mainen.py `, +:download:`cadad.mod `, +:download:`kca.mod `, +:download:`km.mod `, +:download:`kv.mod `, +:download:`naz.mod `, +:download:`Nca.mod ` + +*NetPyNE Code:* :: + + netParams.importCellParams( + label='PYR_Mainen_rule', + conds={'cellType': 'PYR', 'cellModel': 'Mainen'}, + fileName='mainen.py', + cellName='PYR2') + + +.. _import_Friesen: + +Friesen model +^^^^^^^^^^^^^^ + +*Required files:* Friesen cell model defined as python class. Requires multiple mechanisms (including point processes) defined in mod files. Spike generation happens at the ``axon`` section (not the ``soma``). This is indicated in NetPyNE adding the ``spikeGenLoc`` item to the ``axon`` section entry, and specifying the section location (eg. 0.5). + +*Required files:* +:download:`friesen.py `, +:download:`A.mod `, +:download:`GABAa.mod `, +:download:`AMPA.mod `, +:download:`NMDA.mod `, +:download:`OFThpo.mod `, +:download:`OFThresh.mod `, +:download:`netcon.inc `, +:download:`ofc.inc ` + +*NetPyNE Code:* :: + + cellRule = netParams.importCellParams( + label='PYR_Friesen_rule', + conds={'cellType': 'PYR', 'cellModel': 'Friesen'}, + fileName='friesen.py', + cellName='MakeRSFCELL') + + cellRule['secs']['axon']['spikeGenLoc'] = 0.5 # spike generator location. + +.. _import_Izhi03a: + +Izhikevich 2003a model (independent voltage variable) +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +*Description:* Izhikevich, 2003 cell model defined as python class. Requires point process defined in mod file. This version is added to a section but does not employ the section voltage or synaptic mechanisms. Instead it uses its own internal voltage variable and synaptic mechanism. This is indicated in NetPyNE adding the ``vref`` item to the point process entry, and specifying the name of the internal voltage variable (``V``). + +Modeldb link: https://senselab.med.yale.edu/modeldb/showModel.cshtml?model=39948 + +*Required files:* +:download:`izhi2003Wrapper.py `, +:download:`izhi2003a.mod ` + +*NetPyNE Code:* :: + + cellRule = netParams.importCellParams( + label='PYR_Izhi03a_rule', + conds={'cellType': 'PYR', 'cellModel':'Izhi2003a'}, + fileName='izhi2003Wrapper.py', + cellName='IzhiCell', + cellArgs={'type':'tonic spiking', 'host':'dummy'}) + + cellRule['secs']['soma']['pointps']['Izhi2003a_0']['vref'] = 'V' # specify that uses its own voltage V + + +.. _import_Izhi03b: + +Izhikevich 2003b model (uses section voltage) +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +*Description:* Izhikevich, 2003 cell model defined as python class. Requires point process defined in mod file. This version is added to a section and shares the section voltage and synaptic mechanisms. A synaptic mechanism is added from NetPyNE during the connection phase. + +Modeldb link: https://senselab.med.yale.edu/modeldb/showModel.cshtml?model=39948 + +*Required files:* +:download:`izhi2003Wrapper.py `, +:download:`izhi2003b.mod ` + +*NetPyNE Code:* :: + + netParams.importCellParams( + label='PYR_Izhi03b_rule', + conds={'cellType': 'PYR', 'cellModel':'Izhi2003b'}, + fileName='izhi2003Wrapper.py', + cellName='IzhiCell', + cellArgs={'type':'tonic spiking'}) + + +.. _import_Izhi07a: + +Izhikevich 2007a model (independent voltage variable) +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +*Description:* Izhikevich, 2007 cell model defined as python clas. Requires point process defined in mod file. This version is added to a section but does not employ the section voltage or synaptic mechanisms. Instead it uses its own internal voltage variable and synaptic mechanism. This is indicated in NetPyNE adding the ``vref`` item to the point process entry, and specifying the name of the internal voltage variable (``V``). The cell model includes several internal synaptic mechanisms, which can be specified as a list in NetPyNE by adding the ``synList`` item to the point process entry. + +Modeldb link: https://senselab.med.yale.edu/modeldb/showModel.cshtml?model=39948 + +*Required files:* +:download:`izhi2007Wrapper.py `, +:download:`izhi2007a.mod ` + +*NetPyNE Code:* :: + + cellRule = netParams.importCellParams( + label='PYR_Izhi07a_rule', + conds={'cellType': 'PYR', 'cellModel':'Izhi2007a'}, + fileName='izhi2007Wrapper.py', + cellName='IzhiCell', + cellArgs={'type':'RS', 'host':'dummy'}) + + cellRule['secs']['soma']['pointps']['Izhi2007a_0']['vref'] = 'V' # specify that uses its own voltage V + + cellRule['secs']['soma']['pointps']['Izhi2007a_0']['synList'] = ['AMPA', 'NMDA', 'GABAA', 'GABAB'] # specify its own synapses + + +.. _import_Izhi07b: + +Izhikevich 2007b model (uses section voltage) +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +*Description:* Izhikevich, 2007 cell model defined as python class. Requires point process defined in mod file. This version is added to a section and shares the section voltage and synaptic mechanisms. + +Modeldb link: https://senselab.med.yale.edu/modeldb/showModel.cshtml?model=39948 + +*Required files:* +:download:`izhi2007Wrapper.py `, +:download:`izhi2007b.mod ` + +*NetPyNE Code:* :: + + netParams.importCellParams( + label='PYR_Izhi07b_rule', + conds={'cellType': 'PYR', 'cellModel':'Izhi2007b'}, + fileName='izhi2007Wrapper.py', + cellName='IzhiCell', + cellArgs={'type':'RS'}) + + +The full code to import all cell models above and create a network with them is available here: :download:`tut_import.py `. + + +Parameter Optimization of a Simple Neural Network Using An Evolutionary Algorithm +--------------------------------------------------------------------------------- + +This tutorial provides an example of how to use +\ `inspyred `__\ , +an evolutionary algorithm toolkit, to optimize parameters in our prior +\ `tut2.py `__\ \*\* +neural network--modified to remove any code relating to initiating +network simulation and output display--, such that it achieves a target +average firing rate around (~) 17 Hz. + +\*\*Some modification is required near the end of the tut2.py code, to +remove any code relating to initiating network simulation and output +display, all of which has now been handled in the new top level code +(:download:`tut_optimization.py `): + +.. code-block:: python + + # Create network and run simulation + # sim.createSimulateAnalyze(netParams = netParams, simConfig = simConfig) # line commented out + + # import pylab; pylab.show()  # if figures appear empty # line commented out + +excerpt from tut2.py + +Additional Background Reading +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +`A description of the +algorithm `__\  methodology +that will be used to optimize the simple neural network in this example. + +Introduction +^^^^^^^^^^^^^ +Using the inspyred python package to find neural network parameters so +that some property of the network (e.g. firing rate) matches a desired +target can be broken down into 3 steps. First, 1) defining a desired +target model (in this case, some measurable value) and fitness +function--fitness defined here as a calculable value that represents how +close a neural network with a given parameters matches the target. +Subsequently, it is necessary to 2) determine the appropriate neural +network parameters to modify to achieve that model/value. Finally, +3) appropriate parameters for the evolutionary algorithm are defined. +Ultimately, If the inputs to the evolutionary algorithm are appropriate, +then over successive iterations, the parameters determined by the +evolutionary algorithm should generate models closer to the target. + +Particularizing these 3 steps to our example we get: + +.. image:: figs/tut_optimization_diagram.png + :width: 80% + :align: center + +1. Defining a desired target model and fitness function. + +Defining a desired target model is largely arbitrary, some constraints +being that there must be a way to adjust parameters such that the +results are closer to the target model than before (or that fitness is +improved), and that there must be a way to evaluate the fitness of a +model with given parameters. In this case, our target model is a neural +network that achieves an average firing rate of 17 Hz. The fitness for +such a model can be defined as the difference between the average firing +rate of a certain model and the target firing rate of 17 Hz. + +2. Selecting the model parameters to be optimized. + +If a parameter can in some way alter the fitness of the final model, it +may be an appropriate candidate for optimization, depending on what the +model is seeking to achieve. As well as a host of other parameters, +altering the probability, weight or delay of the synaptic connections in +the neural network can affect the average firing rate. In this example, +we will optimize the values of the probability, weight and delay of +connections from the sensory to the motor population. + +3. Selecting appropriate parameters for the evolutionary algorithm. + +inspyred allows customization of the various components of the +evolutionary algorithm, including: + +-  a selector that determines which sets of parameter values become + parents and thus which parameter values will be used to form the next + generation in the evolutionary iteration, +- a variator that determines how each current iteration of parameter + sets is formed from the previous iteration, +- a replacer which determines whether previous sets of parameter values + are brought into the next iteration, +- a terminator which defines when to end evolutionary iterations, +- an observer which allows for tracking of parameter values through + each evolutionary iteration. + +         + +Using inspyred +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +The evolutionary algorithm is implemented the ec module from the +inspyred package: + +.. code-block:: python + + from inspyred import ec # import evolutionary computation from inspyred + +excerpt from tut\_optimization.py + +ec includes a class for the evolutionary computation algorithm: +ec.EvolutionaryComputation(), which allows entering parameters to +customize the algorithm. The evolutionary algorithm involves random +processes (e.g. randomly mutating genes) and so requires random number +generator. In this case we will use python's Random() method, which we +initialize using a specific seed value so that we can reproduce the +results in the future: + +.. code-block:: python + + # create random seed for evolutionary computation algorithm + rand = Random() + rand.seed(1) + + # instantiate evolutionary computation algorithm + my_ec = ec.EvolutionaryComputation(rand) + + +excerpt from tut\_optimization.py + +Parameters for the evolutionary algorithm are then established for our +ec evolutionary computation instance by assigning various variator, +replacer, terminator and observer elements--essentially toggling +specific components of the algorithm-- to ec.selectors, ec.variators, +ec.replacers, ec.terminators, ec.observers: + +.. code-block:: python + + #toggle variators + my_ec.variator = [ec.variators.uniform_crossover, # implement uniform crossover & gaussian replacement + ec.variators.gaussian_mutation]    + my_ec.replacer = ec.replacers.generational_replacement   # implement generational replacement + + my_ec.terminator = ec.terminators.evaluation_termination # termination dictated by no. evaluations + + #toggle observers + my_ec.observer = [ec.observers.stats_observer,  # print evolutionary computation statistics +                 ec.observers.plot_observer,   # plot output of the evolutionary computation as graph +                 ec.observers.best_observer]   # print the best individual in the population to screen + +excerpt from ex_optimization.py + +where: + ++----------------------------------------+--------------------------------------+ +| ec.variators.uniform\_crossover | variator where coin flip to | +| | determine whether 'mom' or 'dad' | +| | element is inherited by offspring | ++----------------------------------------+--------------------------------------+ +| ec.variators.gaussian\_mutation | variator implements gaussian | +| | mutation which makes use of bounder | +| | function as specified | +| | in: my\_ec.evolve(...,bounder=ec.Bou | +| | nder(minParamValues, maxParamValues) | +| | ,...) | +| | | ++----------------------------------------+--------------------------------------+ +| ec.replacers.generational\_replacement | replacer implements generational | +| | replacement with elitism (as | +| | specified in | +| | my\_ec.evolve(...,num\_elites=1,...) | +| | , | +| | where the existing generation is | +| | replaced by offspring, and | +| | existing individuals | +| | will survive if they have better | +| | fitness than the offspring | ++----------------------------------------+--------------------------------------+ +| ec.terminators.evaluation\_termination | terminator runs based on the number | +| | of evaluations that have occurred | ++----------------------------------------+--------------------------------------+ +| ec.observers.stats\_observer | indicates how many of the generated | +| | individuals (parameter sets) will be | +| | selected for the next evolutionary | +| | iteration. | ++----------------------------------------+--------------------------------------+ +| ec.observers.plot\_observer | indicates the rate of mutation, or | +| | the rate at which values for each | +| | parameter (probability, weight and | +| | delay) taken from a prior generation | +| | are altered in the next generation | ++----------------------------------------+--------------------------------------+ +| ec.observers.best\_observer | sets the number of parameters that | +| | will be optimized to 3, | +| | corresponding to the length of | +| | [probability, weight, delay]. | ++----------------------------------------+--------------------------------------+ + +These predefined selector, variator, replacer, terminator and observer +elements as well as other options can be found in the \ `inspyred +documentation `__\ . + +FInally, the evolutionary computation algorithm instance includes a +method: my\_ec.evolve() , which will move through successive +evolutionary iterations evaluating different parameter sets until the +terminating condition is achieved. This function comes with multiple +arguments, with two significant arguments being the generator and +evaluator functions. A function call for  my\_ec.evolve() will look +similar to the following: + +.. code-block:: python + + # call evolution iterator + + final_pop = my_ec.evolve(generator=generate_netparams, # assign model parameter generator to iterator generator +                       evaluator=evaluate_netparams, # assign fitness function to iterator evaluator +                     pop_size=10, +                     maximize=False,                    + bounder=ec.Bounder(minParamValues, maxParamValues), + max_evaluations=50, + num_selected=10, + mutation_rate=0.2, + num_inputs=3, + num_elites=1) + + +excerpt from tut\_optimization.py + +where: + ++--------------------------------------+--------------------------------------+ +| pop\_size=10 | means that each generation of | +| | parameter sets will consist of 10 | +| | individuals | ++--------------------------------------+--------------------------------------+ +| maximize=False | means that we are taking higher | +| | fitness to correspond to minimal | +| | values in terms of difference | +| | between model firing frequency and | +| | 17 Hz | ++--------------------------------------+--------------------------------------+ +| bounder=ec.Bounder(minParamValues, | defines boundaries for each of the | +| maxParamValues) | parameters. The format to describe | +|           | the minimum and maximum values for | +| | the parameters we are seeking to | +|       | optimize: minParamValues is an array | +| | of minimum of values corresponding | +| | to [probability, weight, delay], and | +| | maxParamValues is the array of | +| | maximum values. | ++--------------------------------------+--------------------------------------+ +| max\_evaluations=50 | indicates how many parameter sets | +| | are evaluated prior termination of | +| | the evolutionary iterations | ++--------------------------------------+--------------------------------------+ +| num\_selected=10 | indicates how many of the generated | +| | individuals (parameter sets) will be | +| | selected for the next evolutionary | +| | iteration. | ++--------------------------------------+--------------------------------------+ +| mutation\_rate=0.2 | indicates the rate of mutation, or | +| | the rate at which values for each | +| | parameter (probability, weight and | +| | delay) taken from a prior generation | +| | are altered in the next generation | ++--------------------------------------+--------------------------------------+ +| num\_inputs=3 | sets the number of parameters that | +| | will be optimized to 3, | +| | corresponding to the length of | +| | [probability, weight, delay]. | ++--------------------------------------+--------------------------------------+ +| num\_elites=1 | sets the number of elites to 1. That | +| | is, one individual from the existing | +| | generation may be retained (as | +| | opposed to a complete generational | +| | replacement) if it has better | +| | fitness than an individual selected | +| | from the offspring. | ++--------------------------------------+--------------------------------------+ + +The generator and evaluator arguments expect user defined functions as +inputs, with generator used to define a population of initial parameter +value sets for the very first iteration, and evaluator being the fitness +function that will be used to evaluate each model for how close it is to +the target. In this example, the generator is a fairly straightforward +function which creates an initial set of parameter values (i.e. +[probability, weight, delay] ) by drawing from a parameterized uniform +distribution: + +.. code-block:: python + + # return a set of initialParams which contains a [probability, weight, delay] + + def generate_netparams(random, args): + +     size = args.get('num_inputs') +   initialParams = [random.uniform(minParamValues[i], maxParamValues[i]) for i in range(size)] + + return initialParams + +excerpt from tut\_optimization.py + +The fitness function involves taking a list of sets of parameter values, +i.e. : [ [ a0, b0, c0], [a1, b1, c1], [a2, b2, c2], ... , [an, bn, cn ] +] where a, b, c represent the parameter values and 1 through n +representing the individual number within the population, and +calculating a fitness score for each element of the list, which is then +returned as a list of fitness values (i.e. : [ f0, f1, f2, ... , fn ] +) corresponding to the initial sets of parameter values. It follows the +general template: + +.. code-block:: python + + def evaluate_fitness(candidates, args): +    fitness = [] +    for candidate in candidates: +        fit = some_fitness_function(candidate) +        fitness.append(fit) +    return fitness + +excerpt from tut\_optimization.py + +The actual code that is used to serve as     + some\_fitness\_function(candidate)    is described below: + +  + +Overview of the Fitness Function +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +The fitness function in this case involves 1) creating a neural network +with the given parameters, 2) simulating it to find the average firing +rate, then 3) comparing this firing rate to a target firing rate. + +1. Creating a neural network with the parameters to evaluate + +We will employ the NetPyNE defined network in tut2.py, and modify +the [probability, weight, delay] parameters. This  involves redefining +specific values found in tut2.py found within the connectivity rule +between the S and M populations:    netParams.connParams['S->M']    + +.. code-block:: python + + ## Cell connectivity rules + netParams.connParams['S->M'] = {      #  S -> M label +       'preConds': {'popLabel': 'S'},  # conditions of presyn cells +       'postConds': {'popLabel': 'M'}, # conditions of postsyn cells +       'probability': 0.5,             # probability of connection <-- to be optimized by evolutionary algorithm +       'weight': 0.01,                 # synaptic weight           <-- to be optimized by evolutionary algorithm +       'delay': 5,                     # transmission delay (ms) <-- to be optimized by evolutionary algorithm +       'synMech': 'exc'}               # synaptic mechanism + +excerpt from tut2.py + +these values are replaced in the fitness function with the parameter +values generated by the evolutionary algorithm. As the fitness function +resides within a for loop iterating through the list of candidates (     +for icand,cand in enumerate(candidates):    ), the individual parameters +can be accessed as cand[0], cand[1], and cand[2]. Reassigning values to +the parameters in tut\_optimization.py can be done via the following +line: + +.. code-block:: python + + tut2.netParams.connParams['S->M'][''] =  + +2. Simulating the created neural network and finding the average firing + rate + +Once the network parameters have been modified we can call the +sim.createSimulate() NetPyNE function to run the simulation. We will +pass as arguments the tut2 netParams and simConfig objects that we just +modified. Once the simulation has ran we will have access to the +simulation output via sim.simData.  + +:: + + # create network + sim.createSimulate(netParams=tut2.netParams, simConfig=tut2.simConfig) + +excerpt from tut\_optimization.py + +3. Comparing the average firing rate to a target average firing rate + +To calculate the average firing rate (in spikes/sec = Hz) of the +network, we divide the spikes that have occurred during the simulation, +by the number of neurons and the duration. A list of spike times and a +list of neurons can be accessed via the NetPyNE sim module: + sim.simData['spkt'] and   sim.net.cells   . These are populated after +running   sim.createSimulate()  . From these lists, getting the number +of spike times and neurons is done by using python’s   len()   function. +The duration of the simulation can be accessed in the +tut\_optimization.py code  via        tut2.simConfig.duration    .  The +calculation for average firing rate is thus as follows: + +.. code-block:: python + + # calculate firing rate + numSpikes = float(len(sim.simData['spkt'])) + numCells = float(len(sim.net.cells)) + duration = tut2.simConfig.duration/1000.0 + netFiring = numSpikes/numCells/duration + +excerpt from tut\_optimization.py + +Finally, the average firing rate of the model is compared to the target +firing rate as follows: + +.. code-block:: python + + # calculate fitness for this candidate + fitness = abs(targetFiring - netFiring)  # minimize absolute difference in firing rate + +excerpt from tut\_optimization.py + +Displaying Findings +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +The results of the evolutionary algorithm are displayed on the standard +output (terminal) as well as plotted using the matplotlib package. The +following lines are relevant to showing results of the various +candidates within the iterator: + +.. code-block:: python + + for icand,cand in enumerate(candidates): +       ... +       print '\n CHILD/CANDIDATE %d: Network with prob:%.2f, weight:%.2f, delay:%.1f \n  firing rate: %.1f, FITNESS = %.2f \n'\ +       %(icand, cand[0], cand[1], cand[2], netFiring, fitness) + +excerpt from tut\_optimization.py + +The first line:  for icand,cand in enumerate(candidates): is analogous +to the the iterator  for candidate in candidates:  used in the +pseudocode example above, except that the  enumerate() function will +also return an index--starting from 0-- for each element in the list, +and is used in the subsequent print statement. + +This example also displays the generated candidate with average +frequency closest to 17 Hz. This candidate will exist in the final +generation, and possess the best fitness score (corresponding to a +minimum difference). Since   num\_elites=1   there is no risk that a +prior generation will have a candidate with a better fitness. + +After the evolution finishes, to access the candidate with the best +fitness score, the final generation of candidates, which is returned by +the  my\_ec.evolve()   function is then sorted in reverse (least to +greatest), placing the candidate that achieves an average firing rate +closest to 17 Hz (and therefore has the minimum difference) at the start +of the list (or at position 0). We will use NetPyNE to visualize the +output of this network, by setting the optimized parameters, simulating +the network and plotting a raster plot. The code that performs this task +is isolated below: + +.. code-block:: python + + final_pop = my_ec.evolve(...) + ... + # plot raster of top solutions + final_pop.sort(reverse=True)        # sort final population so best fitness (minimum difference) is first in list + bestCand = final_pop[0].candidate   # bestCand <-- candidate in first position of list + tut2.simConfig.analysis['plotRaster'] = True                      # plotting + tut2.netParams.connParams['S->M']['probability'] = bestCand[0]    # set tut2 values to corresponding + tut2.netParams.connParams['S->M']['weight'] = bestCand[1]         # best candidate values + tut2.netParams.connParams['S->M']['delay'] = bestCand[2] + sim.createSimulateAnalyze(netParams=tut2.netParams, simConfig=tut2.simConfig) # run simulation + +excerpt from tut\_optimization.py + +The code for neural network optimization through evolutionary algorithm used in this tutorial can be found here: :download:`tut_optimization.py `. + + +.. Cell density and connectivity as a function of cell location +.. ------------------------------------------------------------ + + +.. Create population as list of individual cells +.. ------------------------------------------------ +.. (eg. measured experimentally) + + +.. Adding connectivity functions +.. ------------------------------ + + +.. Adding cell classes +.. -------------------- + diff --git a/examples/CA3model_3pops/src/netParams.py b/examples/CA3model_3pops/src/netParams.py index ed795439c..083a02dde 100644 --- a/examples/CA3model_3pops/src/netParams.py +++ b/examples/CA3model_3pops/src/netParams.py @@ -1,4 +1,3 @@ -import imp from netpyne import specs try: from __main__ import cfg diff --git a/examples/CA3model_RasterColoredbyPhase/index.npjson b/examples/CA3model_RasterColoredbyPhase/index.npjson new file mode 100644 index 000000000..20c14b4a6 --- /dev/null +++ b/examples/CA3model_RasterColoredbyPhase/index.npjson @@ -0,0 +1,6 @@ +{ + "mod_folder": "mod", + "simConfig": "src/cfg.py", + "python_run": "src/init.py", + "netParams": "src/netParams.py" +} \ No newline at end of file diff --git a/examples/CA3model_RasterColoredbyPhase/mod/CA1ih.mod b/examples/CA3model_RasterColoredbyPhase/mod/CA1ih.mod new file mode 100644 index 000000000..93d435e30 --- /dev/null +++ b/examples/CA3model_RasterColoredbyPhase/mod/CA1ih.mod @@ -0,0 +1,64 @@ +: $Id: CA1ih.mod,v 1.4 2010/12/13 21:35:47 samn Exp $ +TITLE Ih CA3 + +UNITS { + (mA) = (milliamp) + (mV) = (millivolt) +} + +NEURON { + SUFFIX hcurrent + NONSPECIFIC_CURRENT ih + RANGE g, e, v50, htau, hinf + RANGE gfactor +} + +PARAMETER { + celsius (degC) + g= 0.0001 (mho/cm2) + e= -30 (mV) + v50=-82 (mV) + gfactor = 1 +} + +STATE { + h +} + +ASSIGNED { + ih (mA/cm2) + hinf + htau (ms) + v (mV) +} + +PROCEDURE iassign () { ih=g*h*(v-e)*gfactor } + +BREAKPOINT { + SOLVE states METHOD cnexp + iassign() +} + +DERIVATIVE states { + rates(v) + h'= (hinf- h)/ htau +} + +INITIAL { + rates(v) + h = hinf + iassign() +} + +PROCEDURE rates(v (mV)) { + UNITSOFF + : HCN1 + :hinf = 1/(1+exp(0.151*(v-v50))) + :htau = exp((0.033*(v+75)))/(0.011*(1+exp(0.083*(v+75)))) + + : HCN2 + hinf = 1/(1+exp((v-v50)/10.5)) + htau = (1/(exp(-14.59-0.086*v)+exp(-1.87+0.0701*v))) + UNITSON +} + diff --git a/examples/CA3model_RasterColoredbyPhase/mod/CA1ika.mod b/examples/CA3model_RasterColoredbyPhase/mod/CA1ika.mod new file mode 100644 index 000000000..9e4fe6922 --- /dev/null +++ b/examples/CA3model_RasterColoredbyPhase/mod/CA1ika.mod @@ -0,0 +1,85 @@ +: $Id: CA1ika.mod,v 1.2 2010/12/01 05:06:07 samn Exp $ +TITLE Ika CA1 + +UNITS { + (mA) = (milliamp) + (mV) = (millivolt) +} + +NEURON { + SUFFIX kacurrent + NONSPECIFIC_CURRENT ika, ikad + RANGE g, gd, e, ninf, ntau, ndinf, ndtau, linf, ltau +} + +PARAMETER { + celsius (degC) + g= 0.048 (mho/cm2) + gd= 0 (mho/cm2) + e= -90 (mV) +} + +STATE { + n + nd : distal + l +} + +ASSIGNED { + v (mV) + ika (mA/cm2) + ikad (mA/cm2) + ninf + ntau (ms) + ndinf + ndtau (ms) + linf + ltau (ms) +} + +PROCEDURE iassign () { + ika=g*n*l*(v-e) + ikad=gd*nd*l*(v-e) +} + +BREAKPOINT { + SOLVE states METHOD cnexp + iassign() +} + +DERIVATIVE states { + rates(v) + n'= (ninf- n)/ ntau + l'= (linf- l)/ ltau + nd'= (ndinf-nd)/ndtau +} + +INITIAL { + rates(v) + n = ninf + l = linf + iassign() +} + +PROCEDURE rates(v (mV)) { + LOCAL a, b + UNITSOFF + a = exp(-0.038*(1.5+1/(1+exp(v+40)/5))*(v-11)) + b = exp(-0.038*(0.825+1/(1+exp(v+40)/5))*(v-11)) + ntau=4*b/(1+a) + if (ntau<0.1) {ntau=0.1} + ninf=1/(1+a) + + a=exp(-0.038*(1.8+1/(1+exp(v+40)/5))*(v+1)) + b=exp(-0.038*(0.7+1/(1+exp(v+40)/5))*(v+1)) + ndtau=2*b/(1+a) + if (ndtau<0.1) {ndtau=0.1} + ndinf=1/(1+a) + + a = exp(0.11*(v+56)) + ltau=0.26*(v+50) + if (ltau<2) {ltau=2} + linf=1/(1+a) + UNITSON +} + diff --git a/examples/CA3model_RasterColoredbyPhase/mod/CA1ikdr.mod b/examples/CA3model_RasterColoredbyPhase/mod/CA1ikdr.mod new file mode 100644 index 000000000..4c5236362 --- /dev/null +++ b/examples/CA3model_RasterColoredbyPhase/mod/CA1ikdr.mod @@ -0,0 +1,60 @@ +: $Id: CA1ikdr.mod,v 1.2 2010/12/01 05:10:52 samn Exp $ +TITLE IKDR CA1 + +UNITS { + (mA) = (milliamp) + (mV) = (millivolt) +} + +NEURON { + SUFFIX kdrcurrent + NONSPECIFIC_CURRENT ik + RANGE g, e, ninf, ntau +} + +PARAMETER { + celsius (degC) + g = 0.010 (mho/cm2) + e = -90 (mV) +} + +STATE { + n +} + +ASSIGNED { + v (mV) + ik (mA/cm2) + ninf + ntau (ms) +} + +PROCEDURE iassign () { ik=g*n*(v-e) } + +BREAKPOINT { + SOLVE states METHOD cnexp + iassign() +} + +DERIVATIVE states { + rates(v) + n'= (ninf- n)/ ntau +} + +INITIAL { + rates(v) + n = ninf + iassign() +} + +PROCEDURE rates(v (mV)) { + LOCAL a, b + UNITSOFF + a = exp(-0.11*(v-13)) + b = exp(-0.08*(v-13)) + ntau=50*b/(1+a) + if (ntau<2) {ntau=2} + ninf=1/(1+a) + UNITSON +} + diff --git a/examples/CA3model_RasterColoredbyPhase/mod/CA1ina.mod b/examples/CA3model_RasterColoredbyPhase/mod/CA1ina.mod new file mode 100644 index 000000000..d33ab9739 --- /dev/null +++ b/examples/CA3model_RasterColoredbyPhase/mod/CA1ina.mod @@ -0,0 +1,89 @@ +: $Id: CA1ina.mod,v 1.4 2010/11/30 19:50:00 samn Exp $ +TITLE INa CA1 + +UNITS { + (mA) = (milliamp) + (mV) = (millivolt) +} + +NEURON { + SUFFIX nacurrent + NONSPECIFIC_CURRENT ina + RANGE g, e, vi, ki + RANGE minf,hinf,iinf,mtau,htau,itau : testing +} + +PARAMETER { + : v (mV) + celsius (degC) + g = 0.032 (mho/cm2) + e = 55 (mV) + vi = -60 (mV) + ki = 0.8 +} + +STATE { + m + h + I : i +} + +ASSIGNED { + i (mA/cm2) + ina (mA/cm2) + minf + mtau (ms) + hinf + htau (ms) + iinf + itau (ms) + v (mV) : testing +} + +: PROCEDURE iassign () { ina=g*m*m*m*h*i*(v-e) } +PROCEDURE iassign () { i=g*m*m*m*h*I*(v-e) ina=i} + +BREAKPOINT { + SOLVE states METHOD cnexp + iassign() +} + +DERIVATIVE states { + rates(v) + m' = (minf - m) / mtau + h' = (hinf - h) / htau + : i' = (iinf - i) / itau + I' = (iinf - I) / itau +} + +INITIAL { + rates(v) + h = hinf + m = minf + : i = iinf + I = iinf + iassign() : testing +} + + +PROCEDURE rates(v (mV)) { + LOCAL a, b + UNITSOFF + a = 0.4*(v+30)/(1-exp(-(v+30)/7.2)) + b = 0.124*(v+30)/(exp((v+30)/7.2)-1) + mtau=0.5/(a+b) + if (mtau<0.02) {mtau=0.02} + minf=a/(a+b) + a = 0.03*(v+45)/(1-exp(-(v+45)/1.5)) + b = 0.01*(v+45)/(exp((v+45)/1.5)-1) + htau=0.5/(a+b) + if (htau<0.5) {htau=0.5} + hinf=1/(1+exp((v+50)/4)) + a = exp(0.45*(v+66)) + b = exp(0.09*(v+66)) + itau=3000*b/(1+a) + if (itau<10) {itau=10} + iinf=(1+ki*exp((v-vi)/2))/(1+exp((v-vi)/2)) + UNITSON +} + diff --git a/examples/CA3model_RasterColoredbyPhase/mod/MyExp2SynBB.mod b/examples/CA3model_RasterColoredbyPhase/mod/MyExp2SynBB.mod new file mode 100644 index 000000000..9a68baef1 --- /dev/null +++ b/examples/CA3model_RasterColoredbyPhase/mod/MyExp2SynBB.mod @@ -0,0 +1,67 @@ +: $Id: MyExp2SynBB.mod,v 1.4 2010/12/13 21:27:51 samn Exp $ +NEURON { +: THREADSAFE + POINT_PROCESS MyExp2SynBB + RANGE tau1, tau2, e, i, g, Vwt, gmax + NONSPECIFIC_CURRENT i +} + +UNITS { + (nA) = (nanoamp) + (mV) = (millivolt) + (uS) = (microsiemens) +} + +PARAMETER { + tau1=.1 (ms) <1e-9,1e9> + tau2 = 10 (ms) <1e-9,1e9> + e=0 (mV) + gmax = 1e9 (uS) + Vwt = 0 : weight for inputs coming in from vector +} + +ASSIGNED { + v (mV) + i (nA) + g (uS) + factor + etime (ms) +} + +STATE { + A (uS) + B (uS) +} + +INITIAL { + LOCAL tp + + Vwt = 0 : testing + + if (tau1/tau2 > .9999) { + tau1 = .9999*tau2 + } + A = 0 + B = 0 + tp = (tau1*tau2)/(tau2 - tau1) * log(tau2/tau1) + factor = -exp(-tp/tau1) + exp(-tp/tau2) + factor = 1/factor +} + +BREAKPOINT { + SOLVE state METHOD cnexp + g = B - A + if (g>gmax) {g=gmax}: saturation + i = g*(v - e) +} + +DERIVATIVE state { + A' = -A/tau1 + B' = -B/tau2 +} + +NET_RECEIVE(w (uS)) {LOCAL ww + ww=w + A = A + ww*factor + B = B + ww*factor +} diff --git a/examples/CA3model_RasterColoredbyPhase/mod/MyExp2SynNMDABB.mod b/examples/CA3model_RasterColoredbyPhase/mod/MyExp2SynNMDABB.mod new file mode 100644 index 000000000..01291643a --- /dev/null +++ b/examples/CA3model_RasterColoredbyPhase/mod/MyExp2SynNMDABB.mod @@ -0,0 +1,108 @@ +: $Id: MyExp2SynNMDABB.mod,v 1.4 2010/12/13 21:28:02 samn Exp $ +NEURON { +: THREADSAFE + POINT_PROCESS MyExp2SynNMDABB + RANGE tau1, tau2, e, i, iNMDA, s, sNMDA, r, tau1NMDA, tau2NMDA, Vwt, smax, sNMDAmax + NONSPECIFIC_CURRENT i, iNMDA +} + +UNITS { + (nA) = (nanoamp) + (mV) = (millivolt) + (uS) = (microsiemens) +} + +PARAMETER { + tau1 = 0.1 (ms) <1e-9,1e9> + tau2 = 10 (ms) <1e-9,1e9> + tau1NMDA = 15 (ms) + tau2NMDA = 150 (ms) + e = 0 (mV) + mg = 1 + r = 1 + smax = 1e9 (1) + sNMDAmax = 1e9 (1) + + Vwt = 0 : weight for inputs coming in from vector +} + +ASSIGNED { + v (mV) + i (nA) + iNMDA (nA) + s (1) + sNMDA (1) + mgblock (1) + factor (1) + factor2 (1) + + etime (ms) +} + +STATE { + A (1) + B (1) + A2 (1) + B2 (1) +} + +INITIAL { + + LOCAL tp + + Vwt = 0 : testing + + if (tau1/tau2 > .9999) { + tau1 = .9999*tau2 + } + A = 0 + B = 0 + tp = (tau1*tau2)/(tau2 - tau1) * log(tau2/tau1) + factor = -exp(-tp/tau1) + exp(-tp/tau2) + factor = 1/factor + + if (tau1NMDA/tau2NMDA > .9999) { + tau1NMDA = .9999*tau2NMDA + } + A2 = 0 + B2 = 0 + tp = (tau1NMDA*tau2NMDA)/(tau2NMDA - tau1NMDA) * log(tau2NMDA/tau1NMDA) + factor2 = -exp(-tp/tau1NMDA) + exp(-tp/tau2NMDA) + factor2 = 1/factor2 +} + +BREAKPOINT { + SOLVE state METHOD cnexp + : Jahr Stevens 1990 J. Neurosci + mgblock = 1.0 / (1.0 + 0.28 * exp(-0.062(/mV) * v) ) + s = B - A + sNMDA = B2 - A2 + if (s >smax) {s =smax }: saturation + if (sNMDA>sNMDAmax) {sNMDA=sNMDAmax}: saturation + i = s * (v - e) + iNMDA = sNMDA * (v - e) * mgblock +} + +DERIVATIVE state { + A' = -A/tau1 + B' = -B/tau2 + A2' = -A2/tau1NMDA + B2' = -B2/tau2NMDA +} + +NET_RECEIVE(w (uS)) {LOCAL ww + ww=w + :printf("NMDA Spike: %g\n", t) + if(r>=0){ : if r>=0, g = AMPA + NMDA*r + A = A + factor *ww + B = B + factor *ww + A2 = A2 + factor2*ww*r + B2 = B2 + factor2*ww*r + }else{ + if(r>-1000){ : if r>-1, g = NMDA*r + A2 = A2 - factor2*ww*r + B2 = B2 - factor2*ww*r + } + : if r<0 and r<>-1, g = 0 + } +} diff --git a/examples/CA3model_RasterColoredbyPhase/mod/aux_fun.inc b/examples/CA3model_RasterColoredbyPhase/mod/aux_fun.inc new file mode 100644 index 000000000..ccb579afb --- /dev/null +++ b/examples/CA3model_RasterColoredbyPhase/mod/aux_fun.inc @@ -0,0 +1,43 @@ +: $Id: aux_fun.inc,v 1.1 2009/11/04 01:24:52 samn Exp $ +COMMENT + +//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +// +// NOTICE OF COPYRIGHT AND OWNERSHIP OF SOFTWARE +// +// Copyright 2007, The University Of Pennsylvania +// School of Engineering & Applied Science. +// All rights reserved. +// For research use only; commercial use prohibited. +// Distribution without permission of Maciej T. Lazarewicz not permitted. +// mlazarew@seas.upenn.edu +// +//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +ENDCOMMENT + + + +:------------------------------------------------------------------- +FUNCTION fun1(v(mV),V0(mV),A(/ms),B(mV))(/ms) { + + fun1 = A*exp((v-V0)/B) +} + +FUNCTION fun2(v(mV),V0(mV),A(/ms),B(mV))(/ms) { + + fun2 = A/(exp((v-V0)/B)+1) +} + +FUNCTION fun3(v(mV),V0(mV),A(/ms),B(mV))(/ms) { + + if(fabs((v-V0)/B)<1e-6) { + :if(v==V0) { + fun3 = A*B/1(mV) * (1- 0.5 * (v-V0)/B) + } else { + fun3 = A/1(mV)*(v-V0)/(exp((v-V0)/B)-1) + } +} + +FUNCTION min(x,y) { if (x<=y){ min = x }else{ min = y } } +FUNCTION max(x,y) { if (x>=y){ max = x }else{ max = y } } diff --git a/examples/CA3model_RasterColoredbyPhase/mod/caolmw.mod b/examples/CA3model_RasterColoredbyPhase/mod/caolmw.mod new file mode 100644 index 000000000..3ea21a7ef --- /dev/null +++ b/examples/CA3model_RasterColoredbyPhase/mod/caolmw.mod @@ -0,0 +1,47 @@ +: $Id: caolmw.mod,v 1.2 2010/11/30 16:40:09 samn Exp $ +COMMENT + +//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +// +// NOTICE OF COPYRIGHT AND OWNERSHIP OF SOFTWARE +// +// Copyright 2007, The University Of Pennsylvania +// School of Engineering & Applied Science. +// All rights reserved. +// For research use only; commercial use prohibited. +// Distribution without permission of Maciej T. Lazarewicz not permitted. +// mlazarew@seas.upenn.edu +// +//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +ENDCOMMENT + +UNITS { + (mollar) = (1/liter) + (M) = (mollar) + (mM) = (millimollar) + (mA) = (milliamp) + (mV) = (millivolt) + (mS) = (millisiemens) +} + +NEURON { + SUFFIX Caolmw + USEION ca READ ica, cai WRITE cai + RANGE alpha, tau +} + +PARAMETER { + alpha = 0.002 (cm2-M/mA-ms) + tau = 80 (ms) +} + +ASSIGNED { ica (mA/cm2) } + +INITIAL { cai = 0 } + +STATE { cai (mM) } + +BREAKPOINT { SOLVE states METHOD cnexp } + +DERIVATIVE states { cai' = -(1000) * alpha * ica - cai/tau } diff --git a/examples/CA3model_RasterColoredbyPhase/mod/icaolmw.mod b/examples/CA3model_RasterColoredbyPhase/mod/icaolmw.mod new file mode 100644 index 000000000..51112d099 --- /dev/null +++ b/examples/CA3model_RasterColoredbyPhase/mod/icaolmw.mod @@ -0,0 +1,51 @@ +: $Id: icaolmw.mod,v 1.2 2010/11/30 16:44:13 samn Exp $ +COMMENT + +//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +// +// NOTICE OF COPYRIGHT AND OWNERSHIP OF SOFTWARE +// +// Copyright 2007, The University Of Pennsylvania +// School of Engineering & Applied Science. +// All rights reserved. +// For research use only; commercial use prohibited. +// Distribution without permission of Maciej T. Lazarewicz not permitted. +// mlazarew@seas.upenn.edu +// +//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +ENDCOMMENT + +UNITS { + (mA) = (milliamp) + (mV) = (millivolt) + (mS) = (millisiemens) +} + +NEURON { + SUFFIX ICaolmw + USEION ca WRITE ica + RANGE gca,eca +} + +PARAMETER { + gca = 1 (mS/cm2) + eca = 120 (mV) +} + +ASSIGNED { + ica (mA/cm2) + v (mV) +} + +PROCEDURE iassign () { ica = (1e-3) * gca * mcainf(v)^2 * (v-eca) } + +INITIAL { + iassign() +} + +BREAKPOINT { iassign() } + +FUNCTION mcainf(v(mV)) { mcainf = fun2(v, -20, 1, -9)*1(ms) } + +INCLUDE "aux_fun.inc" diff --git a/examples/CA3model_RasterColoredbyPhase/mod/iholmw.mod b/examples/CA3model_RasterColoredbyPhase/mod/iholmw.mod new file mode 100644 index 000000000..ccd919202 --- /dev/null +++ b/examples/CA3model_RasterColoredbyPhase/mod/iholmw.mod @@ -0,0 +1,60 @@ +: $Id: iholmw.mod,v 1.2 2010/11/30 16:34:22 samn Exp $ +COMMENT + +//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +// +// NOTICE OF COPYRIGHT AND OWNERSHIP OF SOFTWARE +// +// Copyright 2007, The University Of Pennsylvania +// School of Engineering & Applied Science. +// All rights reserved. +// For research use only; commercial use prohibited. +// Distribution without permission of Maciej T. Lazarewicz not permitted. +// mlazarew@seas.upenn.edu +// +//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +ENDCOMMENT + +UNITS { + (mA) = (milliamp) + (mV) = (millivolt) + (mS) = (millisiemens) +} + +NEURON { + SUFFIX Iholmw + NONSPECIFIC_CURRENT i + RANGE gh,eh +} + +PARAMETER { + gh = 0.15 (mS/cm2) + eh = -40 (mV) +} + +ASSIGNED { + v (mV) + i (mA/cm2) +} + +STATE { q } + +PROCEDURE iassign () { i = (1e-3) * gh * q * (v-eh) } + +INITIAL { + q = qinf(v) + iassign() +} + +BREAKPOINT { + SOLVE states METHOD cnexp + iassign() +} + +DERIVATIVE states { q' = (qinf(v)-q)/qtau(v) } + +FUNCTION qinf(v(mV)) { qinf = fun2(v, -80, 1, 10)*1(ms) } +FUNCTION qtau(v(mV))(ms) { qtau = 200(ms)/(exp((v+70(mV))/20(mV))+exp(-(v+70(mV))/20(mV))) + 5(ms) } + +INCLUDE "aux_fun.inc" diff --git a/examples/CA3model_RasterColoredbyPhase/mod/kcaolmw.mod b/examples/CA3model_RasterColoredbyPhase/mod/kcaolmw.mod new file mode 100644 index 000000000..b2368787e --- /dev/null +++ b/examples/CA3model_RasterColoredbyPhase/mod/kcaolmw.mod @@ -0,0 +1,52 @@ +: $Id: kcaolmw.mod,v 1.2 2010/11/30 16:47:18 samn Exp $ +COMMENT + +//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +// +// NOTICE OF COPYRIGHT AND OWNERSHIP OF SOFTWARE +// +// Copyright 2007, The University Of Pennsylvania +// School of Engineering & Applied Science. +// All rights reserved. +// For research use only; commercial use prohibited. +// Distribution without permission of Maciej T. Lazarewicz not permitted. +// mlazarew@seas.upenn.edu +// +//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +ENDCOMMENT + +UNITS { + (mA) = (milliamp) + (mV) = (millivolt) + (mS) = (millisiemens) + (mollar) = (1/liter) + (mM) = (millimollar) +} + +NEURON { + SUFFIX KCaolmw + USEION k WRITE ik + USEION ca READ cai + RANGE gkca,ek,kd +} + +PARAMETER { + gkca = 10 (mS/cm2) + ek = -90 (mV) + kd = 30 (mM) +} + +ASSIGNED { + cai (mM) + v (mV) + ik (mA/cm2) +} + +PROCEDURE iassign () { ik = (1e-3) * gkca * cai/(cai+kd) * (v-ek) } + +INITIAL { + iassign() +} + +BREAKPOINT { iassign() } diff --git a/examples/CA3model_RasterColoredbyPhase/mod/kdrbwb.mod b/examples/CA3model_RasterColoredbyPhase/mod/kdrbwb.mod new file mode 100644 index 000000000..fc52ae534 --- /dev/null +++ b/examples/CA3model_RasterColoredbyPhase/mod/kdrbwb.mod @@ -0,0 +1,76 @@ +: $Id: kdrbwb.mod,v 1.4 2010/12/13 21:35:26 samn Exp $ +COMMENT + +//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +// +// NOTICE OF COPYRIGHT AND OWNERSHIP OF SOFTWARE +// +// Copyright 2007, The University Of Pennsylvania +// School of Engineering & Applied Science. +// All rights reserved. +// For research use only; commercial use prohibited. +// Distribution without permission of Maciej T. Lazarewicz not permitted. +// mlazarew@seas.upenn.edu +// +//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +ENDCOMMENT + +UNITS { + (mA) = (milliamp) + (mV) = (millivolt) + (mS) = (millisiemens) +} + +NEURON { + SUFFIX Kdrbwb + USEION k WRITE ik + RANGE phin,gkdr,ek + RANGE taon,ninf +} + +PARAMETER { + gkdr = 9 (mS/cm2) + ek = -90 (mV) + phin = 5 +} + +ASSIGNED { + v (mV) + ik (mA/cm2) + celsius (degC) + ninf (1) + taon (ms) +} + +STATE { n } + +PROCEDURE iassign () { ik = (1e-3) * gkdr * n^4 * (v-ek) } + +INITIAL { + rates(v) + n = ninf + iassign() +} + +BREAKPOINT { + SOLVE states METHOD cnexp + iassign() +} + +DERIVATIVE states { + rates(v) + n' = (ninf-n)/taon +} + +PROCEDURE rates(v(mV)) { LOCAL an, bn, q10 + q10 = phin:^((celsius-27.0(degC))/10.0(degC)) + + an = fun3(v, -34, -0.01, -10) + bn = fun1(v, -44, 0.125, -80) + + ninf = an/(an+bn) + taon = 1./((an+bn)*q10) +} + +INCLUDE "aux_fun.inc" diff --git a/examples/CA3model_RasterColoredbyPhase/mod/nafbwb.mod b/examples/CA3model_RasterColoredbyPhase/mod/nafbwb.mod new file mode 100644 index 000000000..37281dc94 --- /dev/null +++ b/examples/CA3model_RasterColoredbyPhase/mod/nafbwb.mod @@ -0,0 +1,81 @@ +: $Id: nafbwb.mod,v 1.4 2010/12/13 21:35:08 samn Exp $ +COMMENT + +//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +// +// NOTICE OF COPYRIGHT AND OWNERSHIP OF SOFTWARE +// +// Copyright 2007, The University Of Pennsylvania +// School of Engineering & Applied Science. +// All rights reserved. +// For research use only; commercial use prohibited. +// Distribution without permission of Maciej T. Lazarewicz not permitted. +// mlazarew@seas.upenn.edu +// +//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +ENDCOMMENT + +UNITS { + (mA) = (milliamp) + (mV) = (millivolt) + (mS) = (millisiemens) +} + +NEURON { + SUFFIX Nafbwb + USEION na WRITE ina + RANGE phih + RANGE gna, ena, taoh : testing +} + +PARAMETER { + gna = 35 (mS/cm2) + ena = 55 (mV) + phih = 5 +} + +ASSIGNED { + v (mV) + ina (mA/cm2) + minf (1) + hinf (1) + taoh (ms) + celsius (degC) +} + +STATE { h } + +PROCEDURE iassign () { ina = (1e-3) * gna * minf^3 * h * (v-ena) } + +INITIAL { + rates(v) + h = hinf + iassign() +} + +BREAKPOINT { + SOLVE states METHOD cnexp + iassign() +} + +DERIVATIVE states { + rates(v) + h' = (hinf-h)/taoh +} + +PROCEDURE rates(v(mV)) { LOCAL am, bm, ah, bh, q10 + + q10 = phih:^((celsius-27.0(degC))/10.0(degC)) + + am = fun3(v, -35, -0.1, -10) + bm = fun1(v, -60, 4, -18) + minf = am/(am+bm) + + ah = fun1(v, -58, 0.07, -20) + bh = fun2(v, -28, 1, -10) + hinf = ah/(ah+bh) + taoh = 1./((ah+bh)*q10) +} + +INCLUDE "aux_fun.inc" diff --git a/examples/CA3model_RasterColoredbyPhase/src/cfg.py b/examples/CA3model_RasterColoredbyPhase/src/cfg.py new file mode 100644 index 000000000..85483afd6 --- /dev/null +++ b/examples/CA3model_RasterColoredbyPhase/src/cfg.py @@ -0,0 +1,34 @@ +from netpyne import specs + +cfg = specs.SimConfig() + +cfg.duration = 750 +cfg.dt = 0.1 +cfg.hparams = {'v_init': -65.0} +cfg.verbose = False +cfg.recordTraces = {'V_soma':{'sec':'soma','loc':0.5,'var':'v'}} # Dict with traces to record +cfg.recordStim = False +cfg.recordStep = 0.1 # Step size in ms to save data (eg. V traces, LFP, etc) +cfg.filename = '00' # Set file output name +cfg.savePickle = False # Save params, network and sim output to pickle file +cfg.saveDat = False +cfg.printRunTime = 0.1 +cfg.recordLFP = [[50, 50, 50]] + +cfg.analysis['plotRaster']={ + 'colorbyPhase':{ +# 'signal': LFP_array, # when signal is provided as a np.array +# 'signal': 'LFP_signal.pkl', # when signal is provided as a list in a .pkl +# 'fs': 10000, # in both cases, sampling frequency would be neccessary + 'signal': 'LFP', # internal signal, from the computation of LFPs + 'electrode':1, + 'filtFreq':[4,8], + 'filtOrder': 3, + 'pop_background': True, + 'include_signal': True + }, + 'include':['allCells'], + 'saveFig':True, 'timeRange': [100,600]} # Plot a raster +cfg.analysis['plotTraces'] = {'include': [0, 1, 800, 801, 1000, 1001], 'saveFig': True} # Plot recorded traces for this list of cells +cfg.analysis['plotLFPTimeSeries'] = {'showFig': True, 'saveFig': True} +#cfg.analysis['plotShape'] = {'includePre': ['all'],'includePost': [0,800,1000],'cvar':'numSyns','dist':0.7, 'saveFig': True} diff --git a/examples/CA3model_RasterColoredbyPhase/src/init.py b/examples/CA3model_RasterColoredbyPhase/src/init.py new file mode 100644 index 000000000..5c371b5c4 --- /dev/null +++ b/examples/CA3model_RasterColoredbyPhase/src/init.py @@ -0,0 +1,10 @@ +""" +init.py + +Starting script to run NetPyNE-based CA3 model. +""" + +from netpyne import sim + +cfg, netParams = sim.loadFromIndexFile('index.npjson') +sim.createSimulateAnalyze(netParams, cfg) diff --git a/examples/CA3model_RasterColoredbyPhase/src/netParams.py b/examples/CA3model_RasterColoredbyPhase/src/netParams.py new file mode 100644 index 000000000..083a02dde --- /dev/null +++ b/examples/CA3model_RasterColoredbyPhase/src/netParams.py @@ -0,0 +1,322 @@ +from netpyne import specs +try: + from __main__ import cfg +except: + from cfg import cfg + +# Network parameters +netParams = specs.NetParams() # object of class NetParams to store the network parameters +netParams.defaultThreshold = 0.0 +netParams.defineCellShapes = True # sets 3d geometry aligned along the y-axis + + +############################################################################### +## Cell types +############################################################################### +# Basket cell +BasketCell = {'secs':{}} +BasketCell['secs']['soma'] = {'geom': {}, 'mechs': {}} +BasketCell['secs']['soma']['geom'] = {'diam': 100, 'L': 31.831, 'nseg': 1, 'cm': 1} +BasketCell['secs']['soma']['mechs'] = {'pas': {'g': 0.1e-3, 'e': -65}, 'Nafbwb': {}, 'Kdrbwb': {}} +netParams.cellParams['BasketCell'] = BasketCell + + +# OLM cell +OlmCell = {'secs':{}} +OlmCell['secs']['soma'] = {'geom': {}, 'mechs': {}} +OlmCell['secs']['soma']['geom'] = {'diam': 100, 'L': 31.831, 'nseg': 1, 'cm': 1} +OlmCell['secs']['soma']['mechs'] = { + 'pas': {'g': 0.1e-3, 'e': -65}, + 'Nafbwb': {}, + 'Kdrbwb': {}, + 'Iholmw': {}, + 'Caolmw': {}, + 'ICaolmw': {}, + 'KCaolmw': {}} +netParams.cellParams['OlmCell'] = OlmCell + + +# Pyramidal cell +PyrCell = {'secs':{}} +PyrCell['secs']['soma'] = {'geom': {}, 'mechs': {}} +PyrCell['secs']['soma']['geom'] = {'diam': 20, 'L': 20, 'cm': 1, 'Ra': 150} +PyrCell['secs']['soma']['mechs'] = { + 'pas': {'g': 0.0000357, 'e': -70}, + 'nacurrent': {}, + 'kacurrent': {}, + 'kdrcurrent': {}, + 'hcurrent': {}} +PyrCell['secs']['Bdend'] = {'geom': {}, 'mechs': {}} +PyrCell['secs']['Bdend']['geom'] = {'diam': 2, 'L': 200, 'cm': 1, 'Ra': 150} +PyrCell['secs']['Bdend']['topol'] = {'parentSec': 'soma', 'parentX': 0, 'childX': 0} +PyrCell['secs']['Bdend']['mechs'] = { + 'pas': {'g': 0.0000357, 'e': -70}, + 'nacurrent': {'ki': 1}, + 'kacurrent': {}, + 'kdrcurrent': {}, + 'hcurrent': {}} +PyrCell['secs']['Adend1'] = {'geom': {}, 'mechs': {}} +PyrCell['secs']['Adend1']['geom'] = {'diam': 2, 'L': 150, 'cm': 1, 'Ra': 150} +PyrCell['secs']['Adend1']['topol'] = {'parentSec': 'soma', 'parentX': 1.0, 'childX': 0} # here there is a change: connected to end soma(1) instead of soma(0.5) +PyrCell['secs']['Adend1']['mechs'] = { + 'pas': {'g': 0.0000357, 'e': -70}, + 'nacurrent': {'ki': 0.5}, + 'kacurrent': {'g': 0.072}, + 'kdrcurrent': {}, + 'hcurrent': {'v50': -82, 'g': 0.0002}} +PyrCell['secs']['Adend2'] = {'geom': {}, 'mechs': {}} +PyrCell['secs']['Adend2']['geom'] = {'diam': 2, 'L': 150, 'cm': 1, 'Ra': 150} +PyrCell['secs']['Adend2']['topol'] = {'parentSec': 'Adend1', 'parentX': 1, 'childX': 0} +PyrCell['secs']['Adend2']['mechs'] = { + 'pas': {'g': 0.0000357, 'e': -70}, + 'nacurrent': {'ki': 0.5}, + 'kacurrent': {'g': 0, 'gd': 0.120}, + 'kdrcurrent': {}, + 'hcurrent': {'v50': -90, 'g': 0.0004}} +PyrCell['secs']['Adend3'] = {'geom': {}, 'mechs': {}} +PyrCell['secs']['Adend3']['geom'] = {'diam': 2, 'L': 150, 'cm': 2, 'Ra': 150} +PyrCell['secs']['Adend3']['topol'] = {'parentSec': 'Adend2', 'parentX': 1, 'childX': 0} +PyrCell['secs']['Adend3']['mechs'] = { + 'pas': {'g': 0.0000714, 'e': -70}, + 'nacurrent': {'ki': 0.5}, + 'kacurrent': {'g': 0, 'gd': 0.200}, + 'kdrcurrent': {}, + 'hcurrent': {'v50': -90, 'g': 0.0007}} +netParams.cellParams['PyrCell'] = PyrCell + + +############################################################################### +## Synaptic mechs +############################################################################### + +netParams.synMechParams['AMPAf'] = {'mod': 'MyExp2SynBB', 'tau1': 0.05, 'tau2': 5.3, 'e': 0} +netParams.synMechParams['NMDA'] = {'mod': 'MyExp2SynNMDABB', 'tau1': 0.05, 'tau2': 5.3, 'tau1NMDA': 15, 'tau2NMDA': 150, 'r': 1, 'e': 0} +netParams.synMechParams['GABAf'] = {'mod': 'MyExp2SynBB', 'tau1': 0.07, 'tau2': 9.1, 'e': -80} +netParams.synMechParams['GABAs'] = {'mod': 'MyExp2SynBB', 'tau1': 0.2, 'tau2': 20, 'e': -80} +netParams.synMechParams['GABAss'] = {'mod': 'MyExp2SynBB', 'tau1': 20, 'tau2': 40, 'e': -80} + + +############################################################################### +## Populations +############################################################################### +netParams.popParams['PYR'] = {'cellType': 'PyrCell', 'numCells': 800} +netParams.popParams['BC'] = {'cellType': 'BasketCell', 'numCells': 200} +netParams.popParams['OLM'] = {'cellType': 'OlmCell', 'numCells': 200} + + +############################################################################### +# Current-clamp to cells +############################################################################### +netParams.stimSourceParams['IClamp_PYR'] = {'type': 'IClamp', 'del': 2*cfg.dt, 'dur': 1e9, 'amp': 50e-3} +netParams.stimSourceParams['IClamp_OLM'] = {'type': 'IClamp', 'del': 2*cfg.dt, 'dur': 1e9, 'amp': -25e-3} + +netParams.stimTargetParams['IClamp_PYR->PYR'] = { + 'source': 'IClamp_PYR', + 'sec': 'soma', + 'loc': 0.5, + 'conds': {'pop': 'PYR'}} + +netParams.stimTargetParams['IClamp_OLM->OLM'] = { + 'source': 'IClamp_OLM', + 'sec': 'soma', + 'loc': 0.5, + 'conds': {'pop': 'OLM'}} + + +############################################################################### +# Setting connections +############################################################################### + +# PYR -> X, NMDA +netParams.connParams['PYR->BC_NMDA'] = {'preConds': {'pop': 'PYR'}, 'postConds': {'pop': 'BC'}, + 'convergence': 100, + 'weight': 1.38e-3, + 'delay': 2, + 'sec': 'soma', + 'loc': 0.5, + 'synMech': 'NMDA'} + +netParams.connParams['PYR->OLM_NMDA'] = {'preConds': {'pop': 'PYR'}, 'postConds': {'pop': 'OLM'}, + 'convergence': 10, + 'weight': 0.7e-3, + 'delay': 2, + 'sec': 'soma', + 'loc': 0.5, + 'synMech': 'NMDA'} + +netParams.connParams['PYR->PYR_NMDA'] = {'preConds': {'pop': 'PYR'}, 'postConds': {'pop': 'PYR'}, + 'convergence': 25, + 'weight': 0.004e-3, + 'delay': 2, + 'sec': 'Bdend', + 'loc': 1.0, + 'synMech': 'NMDA'} + +# PYR -> X, AMPA +netParams.connParams['PYR->BC_AMPA'] = {'preConds': {'pop': 'PYR'}, 'postConds': {'pop': 'BC'}, + 'convergence': 100, + 'weight': 0.36e-3, + 'delay': 2, + 'sec': 'soma', + 'loc': 0.5, + 'synMech': 'AMPAf'} + +netParams.connParams['PYR->OLM_AMPA'] = {'preConds': {'pop': 'PYR'}, 'postConds': {'pop': 'OLM'}, + 'convergence': 10, + 'weight': 0.36e-3, + 'delay': 2, + 'sec': 'soma', + 'loc': 0.5, + 'synMech': 'AMPAf'} + +netParams.connParams['PYR->PYR_AMPA'] = {'preConds': {'pop': 'PYR'}, 'postConds': {'pop': 'PYR'}, + 'convergence': 25, + 'weight': 0.02e-3, + 'delay': 2, + 'sec': 'Bdend', + 'loc': 1.0, + 'synMech': 'AMPAf'} + +# BC -> X, GABA +netParams.connParams['BC->BC_GABA'] = {'preConds': {'pop': 'BC'}, 'postConds': {'pop': 'BC'}, + 'convergence': 60, + 'weight':4.5e-3, + 'delay': 2, + 'sec': 'soma', + 'loc': 0.5, + 'synMech': 'GABAf'} + +netParams.connParams['BC->PYR_GABA'] = {'preConds': {'pop': 'BC'}, 'postConds': {'pop': 'PYR'}, + 'convergence': 50, + 'weight': 0.72e-3, + 'delay': 2, + 'sec': 'soma', + 'loc': 0.5, + 'synMech': 'GABAf'} + + +# OLM -> PYR, GABA +netParams.connParams['OLM->PYR_GABA'] = {'preConds': {'pop': 'OLM'}, 'postConds': {'pop': 'PYR'}, + 'convergence': 20, + 'weight': 72e-3, + 'delay': 2, + 'sec': 'Adend2', + 'loc': 0.5, + 'synMech': 'GABAs'} + + +############################################################################### +# Setting NetStims +############################################################################### +# to PYR +netParams.stimSourceParams['NetStim_PYR_SOMA_AMPA'] = {'type': 'NetStim', 'interval': 1, 'number': 1000*cfg.duration, 'start': 0, 'noise': 1} +netParams.stimTargetParams['NetStim_PYR_SOMA_AMPA->PYR'] = { + 'source': 'NetStim_PYR_SOMA_AMPA', + 'conds': {'pop': 'PYR'}, + 'sec': 'soma', + 'loc': 0.5, + 'weight': 4*0.05e-3, # different from published value + 'delay': 2*cfg.dt, + 'synMech': 'AMPAf'} + +netParams.stimSourceParams['NetStim_PYR_ADEND3_AMPA'] = {'type': 'NetStim', 'interval': 1, 'number': 1000*cfg.duration, 'start': 0, 'noise': 1} +netParams.stimTargetParams['NetStim_PYR_ADEND3_AMPA->PYR'] = { + 'source': 'NetStim_PYR_ADEND3_AMPA', + 'conds': {'pop': 'PYR'}, + 'sec': 'Adend3', + 'loc': 0.5, + 'weight': 4*0.05e-3, # different from published value + 'delay': 2*cfg.dt, + 'synMech': 'AMPAf'} + +netParams.stimSourceParams['NetStim_PYR_SOMA_GABA'] = {'type': 'NetStim', 'interval': 1, 'number': 1000*cfg.duration, 'start': 0, 'noise': 1} +netParams.stimTargetParams['NetStim_PYR_SOMA_GABA->PYR'] = { + 'source': 'NetStim_PYR_SOMA_GABA', + 'conds': {'pop': 'PYR'}, + 'sec': 'soma', + 'loc': 0.5, + 'weight': 0.012e-3, + 'delay': 2*cfg.dt, + 'synMech': 'GABAf'} + +netParams.stimSourceParams['NetStim_PYR_ADEND3_GABA'] = {'type': 'NetStim', 'interval': 1, 'number': 1000*cfg.duration, 'start': 0, 'noise': 1} +netParams.stimTargetParams['NetStim_PYR_ADEND3_GABA->PYR'] = { + 'source': 'NetStim_PYR_ADEND3_GABA', + 'conds': {'pop': 'PYR'}, + 'sec': 'Adend3', + 'loc': 0.5, + 'weight': 0.012e-3, + 'delay': 2*cfg.dt, + 'synMech': 'GABAf'} + +netParams.stimSourceParams['NetStim_PYR_ADEND3_NMDA'] = {'type': 'NetStim', 'interval': 100, 'number': int((1000/100.0)*cfg.duration), 'start': 0, 'noise': 1} +netParams.stimTargetParams['NetStim_PYR_ADEND3_NMDA->PYR'] = { + 'source': 'NetStim_PYR_ADEND3_NMDA', + 'conds': {'pop': 'PYR'}, + 'sec': 'Adend3', + 'loc': 0.5, + 'weight': 6.5e-3, + 'delay': 2*cfg.dt, + 'synMech': 'NMDA'} + +# to BC +netParams.stimSourceParams['NetStim_BC_SOMA_AMPA'] = {'type': 'NetStim', 'interval': 1, 'number': 1000*cfg.duration, 'start': 0, 'noise': 1} +netParams.stimTargetParams['NetStim_BC_SOMA_AMPA->BC'] = { + 'source': 'NetStim_BC_SOMA_AMPA', + 'conds': {'pop': 'BC'}, + 'sec': 'soma', + 'loc': 0.5, + 'weight': 0.02e-3, + 'delay': 2*cfg.dt, + 'synMech': 'AMPAf'} + +netParams.stimSourceParams['NetStim_BC_SOMA_GABA'] = {'type': 'NetStim', 'interval': 1, 'number': 1000*cfg.duration, 'start': 0, 'noise': 1} +netParams.stimTargetParams['NetStim_BC_SOMA_GABA->BC'] = { + 'source': 'NetStim_BC_SOMA_GABA', + 'conds': {'pop': 'BC'}, + 'sec': 'soma', + 'loc': 0.5, + 'weight': 0.2e-3, + 'delay': 2*cfg.dt, + 'synMech': 'GABAf'} + +# to OLM +netParams.stimSourceParams['NetStim_OLM_SOMA_AMPA'] = {'type': 'NetStim', 'interval': 1, 'number': 1000*cfg.duration, 'start': 0, 'noise': 1} +netParams.stimTargetParams['NetStim_OLM_SOMA_AMPA->OLM'] = { + 'source': 'NetStim_OLM_SOMA_AMPA', + 'conds': {'pop': 'OLM'}, + 'sec': 'soma', + 'loc': 0.5, + 'weight': 0.0625e-3, + 'delay': 2*cfg.dt, + 'synMech': 'AMPAf'} + +netParams.stimSourceParams['NetStim_OLM_SOMA_GABA'] = {'type': 'NetStim', 'interval': 1, 'number': 1000*cfg.duration, 'start': 0, 'noise': 1} +netParams.stimTargetParams['NetStim_OLM_SOMA_GABA->OLM'] = { + 'source': 'NetStim_OLM_SOMA_GABA', + 'conds': {'pop': 'OLM'}, + 'sec': 'soma', + 'loc': 0.5, + 'weight': 0.2e-3, + 'delay': 2*cfg.dt, + 'synMech': 'GABAf'} + +# Medial Septal inputs to BC and OLM cells +netParams.stimSourceParams['Septal'] = {'type': 'NetStim', 'interval': 150, 'number': int((1000/150)*cfg.duration), 'start': 0, 'noise': 0} +netParams.stimTargetParams['Septal->BC'] = { + 'source': 'Septal', + 'conds': {'pop': 'BC'}, + 'sec': 'soma', + 'loc': 0.5, + 'weight': 1.6e-3, + 'delay': 2*cfg.dt, + 'synMech': 'GABAss'} + +netParams.stimTargetParams['Septal->OLM'] = { + 'source': 'Septal', + 'conds': {'pop': 'OLM'}, + 'sec': 'soma', + 'loc': 0.5, + 'weight': 1.6e-3, + 'delay': 2*cfg.dt, + 'synMech': 'GABAss'} diff --git a/examples/HHTut/index.npjson b/examples/HHTut/index.npjson index 073018a1b..2ea0fcda6 100644 --- a/examples/HHTut/index.npjson +++ b/examples/HHTut/index.npjson @@ -1,5 +1,6 @@ { - "simConfig": "src/cfg.py", + "simConfig": "src/HHTut.py", + "simConfig_variable": "simConfig", "python_run": "src/init.py", - "netParams": "src/netParams.py" + "netParams": "src/HHTut.py" } \ No newline at end of file diff --git a/examples/HHTut/src/HHTut.py b/examples/HHTut/src/HHTut.py new file mode 100644 index 000000000..9ff40e99a --- /dev/null +++ b/examples/HHTut/src/HHTut.py @@ -0,0 +1,83 @@ +""" +HHTut.py +""" + +from netpyne import specs, sim + +netParams = specs.NetParams() # object of class NetParams to store the network parameters +simConfig = specs.SimConfig() # object of class SimConfig to store the simulation configuration + + +############################################################################### +# +# HH TUTORIAL +# +############################################################################### + +############################################################################### +# NETWORK PARAMETERS +############################################################################### + +# Population parameters +netParams.popParams['PYR'] = {'cellType': 'PYR', 'numCells': 200} # add dict with params for this pop + + +# Cell parameters +## PYR cell properties +PYRcell = {'secs': {}} # cell rule dict +PYRcell['secs']['soma'] = {'geom': {}, 'mechs': {}} # soma params dict +PYRcell['secs']['soma']['geom'] = {'diam': 18.8, 'L': 18.8, 'Ra': 123.0} # soma geometry +PYRcell['secs']['soma']['mechs']['hh'] = {'gnabar': 0.12, 'gkbar': 0.036, 'gl': 0.003, 'el': -70} # soma hh mechanism +PYRcell['secs']['soma']['vinit'] = -71 # set initial membrane potential +netParams.cellParams['PYR'] = PYRcell # add dict to list of cell params + + +# Synaptic mechanism parameters +netParams.synMechParams['AMPA'] = {'mod': 'Exp2Syn', 'tau1': 0.1, 'tau2': 1.0, 'e': 0} + + +# Stimulation parameters +netParams.stimSourceParams['bkg'] = {'type': 'NetStim', 'rate': 10, 'noise': 0.5, 'start': 1} +netParams.stimTargetParams['bkg->PYR'] = {'source': 'bkg', 'conds': {'pop': 'PYR'}, 'weight': 0.1, 'delay': 'uniform(1,5)'} + + +# Connectivity parameters +netParams.connParams['PYR->PYR'] = { + 'preConds': {'pop': 'PYR'}, 'postConds': {'pop': 'PYR'}, + 'weight': 0.002, # weight of each connection + 'delay': '0.2+normal(13.0,1.4)', # delay min=0.2, mean=13.0, var = 1.4 + 'threshold': 10, # threshold + 'convergence': 'uniform(1,15)'} # convergence (num presyn targeting postsyn) is uniformly distributed between 1 and 15 + + + +############################################################################### +# SIMULATION PARAMETERS +############################################################################### + +# Simulation parameters +simConfig.duration = 1*1e3 # Duration of the simulation, in ms +simConfig.dt = 0.025 # Internal integration timestep to use +simConfig.seeds = {'conn': 1, 'stim': 1, 'loc': 1} # Seeds for randomizers (connectivity, input stimulation and cell locations) +simConfig.verbose = False # show detailed messages +simConfig.hParams = {'v_init': PYRcell['secs']['soma']['vinit']} + +# Recording +simConfig.recordCells = [] # which cells to record from +simConfig.recordTraces = {'Vsoma': {'sec': 'soma','loc': 0.5,'var': 'v'}} +simConfig.recordStim = True # record spikes of cell stims +simConfig.recordStep = 0.1 # Step size in ms to save data (eg. V traces, LFP, etc) + +# Saving +simConfig.filename = 'HHTut' # Set file output name +simConfig.saveFileStep = 1000 # step size in ms to save data to disk +simConfig.savePickle = False # Whether or not to write spikes etc. to a .mat file +simConfig.saveJson = True + +# Analysis and plotting +simConfig.analysis['plotRaster'] = {'saveData': 'raster_data.json', 'saveFig': True, 'showFig': True} # Plot raster +simConfig.analysis['plotTraces'] = {'include': [2], 'saveFig': True, 'showFig': True} # Plot cell traces +simConfig.analysis['plot2Dnet'] = {'saveFig': True, 'showFig': True} # Plot 2D cells and connections + + +simConfig.validateNetParams=True \ No newline at end of file diff --git a/examples/HHTut/src/cfg.py b/examples/HHTut/src/cfg.py deleted file mode 100644 index fff149cd4..000000000 --- a/examples/HHTut/src/cfg.py +++ /dev/null @@ -1,29 +0,0 @@ -from netpyne import specs - -cfg = specs.SimConfig() # object of class SimConfig to store the simulation configuration - -# Simulation parameters -cfg.duration = 1*1e3 # Duration of the simulation, in ms -cfg.dt = 0.025 # Internal integration timestep to use -cfg.seeds = {'conn': 1, 'stim': 1, 'loc': 1} # Seeds for randomizers (connectivity, input stimulation and cell locations) -cfg.verbose = False # show detailed messages -cfg.hParams = {'v_init': -71} - -# Recording -cfg.recordCells = [] # which cells to record from -cfg.recordTraces = {'Vsoma': {'sec': 'soma','loc': 0.5,'var': 'v'}} -cfg.recordStim = True # record spikes of cell stims -cfg.recordStep = 0.1 # Step size in ms to save data (eg. V traces, LFP, etc) - -# Saving -cfg.filename = 'HHTut' # Set file output name -cfg.saveFileStep = 1000 # step size in ms to save data to disk -cfg.savePickle = False # Whether or not to write spikes etc. to a .mat file -cfg.saveJson = True - -# Analysis and plotting -cfg.analysis['plotRaster'] = {'saveData': 'raster_data.json', 'saveFig': True, 'showFig': True} # Plot raster -cfg.analysis['plotTraces'] = {'include': [2], 'saveFig': True, 'showFig': True} # Plot cell traces -cfg.analysis['plot2Dnet'] = {'saveFig': True, 'showFig': True} # Plot 2D cells and connections - -cfg.validateNetParams=True \ No newline at end of file diff --git a/examples/HHTut/src/netParams.py b/examples/HHTut/src/netParams.py deleted file mode 100644 index ed4add11e..000000000 --- a/examples/HHTut/src/netParams.py +++ /dev/null @@ -1,38 +0,0 @@ -from netpyne import specs - -try: - from __main__ import cfg # import SimConfig object with params from parent module -except: - from cfg import cfg - -netParams = specs.NetParams() # object of class NetParams to store the network parameters - -# Population parameters -netParams.popParams['PYR'] = {'cellType': 'PYR', 'numCells': 200} # add dict with params for this pop - -# Cell parameters -## PYR cell properties -PYRcell = {'secs': {}} # cell rule dict -PYRcell['secs']['soma'] = {'geom': {}, 'mechs': {}} # soma params dict -PYRcell['secs']['soma']['geom'] = {'diam': 18.8, 'L': 18.8, 'Ra': 123.0} # soma geometry -PYRcell['secs']['soma']['mechs']['hh'] = {'gnabar': 0.12, 'gkbar': 0.036, 'gl': 0.003, 'el': -70} # soma hh mechanism -PYRcell['secs']['soma']['vinit'] = cfg.hParams['v_init'] # set initial membrane potential -netParams.cellParams['PYR'] = PYRcell # add dict to list of cell params - - -# Synaptic mechanism parameters -netParams.synMechParams['AMPA'] = {'mod': 'Exp2Syn', 'tau1': 0.1, 'tau2': 1.0, 'e': 0} - - -# Stimulation parameters -netParams.stimSourceParams['bkg'] = {'type': 'NetStim', 'rate': 10, 'noise': 0.5, 'start': 1} -netParams.stimTargetParams['bkg->PYR'] = {'source': 'bkg', 'conds': {'pop': 'PYR'}, 'weight': 0.1, 'delay': 'uniform(1,5)'} - - -# Connectivity parameters -netParams.connParams['PYR->PYR'] = { - 'preConds': {'pop': 'PYR'}, 'postConds': {'pop': 'PYR'}, - 'weight': 0.002, # weight of each connection - 'delay': '0.2+normal(13.0,1.4)', # delay min=0.2, mean=13.0, var = 1.4 - 'threshold': 10, # threshold - 'convergence': 'uniform(1,15)'} # convergence (num presyn targeting postsyn) is uniformly distributed between 1 and 15 \ No newline at end of file diff --git a/examples/LFPrecording/README.md b/examples/LFPrecording/README.md index 6fbde8dcc..60550c0f8 100644 --- a/examples/LFPrecording/README.md +++ b/examples/LFPrecording/README.md @@ -4,35 +4,35 @@ Two examples of using LFP recording in NetPyNE: 1) LFP recording for a single cell -- M1 corticospinal neuron -- with 700+ compartments (segments) and multiple somatic and dendritic ionic channels. The cell parameters are loaded from a .json file. The cell receives NetStim input to its soma via an excitatory synapse. Ten LFP electrodes are placed at both sides of the neuron at 5 different cortical depths. The soma voltage, LFP time-resolved signal and the 3D locations of the electrodes are plotted: -![lfp_cell](https://github.com/Neurosim-lab/netpyne/raw/lfp/examples/LFPrecording/lfp_cell.png) +![lfp_cell](https://raw.githubusercontent.com/suny-downstate-medical-center/netpyne/development/doc/source/figs/lfp_cell.png) 2) LFP recording for a network very similar to that shown in Tutorial 5. However, in this case, the cells have been replaced with a more realistic model: a 6-compartment corticostriatal neuron with multiple ionic channels. The cell parameters are loaded from a .json file. Cell receive NetStim inputs and include excitatory and inhibitory connections. Four LFP electrodes are placed at different cortical depths. The raster plot and LFP time-resolved signal, PSD, spectrogram and 3D locations of the electrodes are plotted: -![lfp_net](https://github.com/Neurosim-lab/netpyne/raw/lfp/examples/LFPrecording/lfp_net.png) +![lfp_net](https://raw.githubusercontent.com/suny-downstate-medical-center/netpyne/development/doc/source/figs/lfp_net.png) To record LFP just set the list of 3D locations of the LFP electrodes in the `simConfig` attribute `recordLFP` e.g. ``simConfig.recordLFP = e.g. [[50, 100, 50], [50, 200, 50]]`` (note the y coordinate represents depth, so will be represented as a negative value whehn plotted). The LFP signal in each electrode is obtained by summing the extracellular potential contributed by each segment of each neuron. Extracellular potentials are calculated using the "line source approximation" and assuming an Ohmic medium with conductivity sigma = 0.3 mS/mm. For more information on modeling LFPs see http://www.scholarpedia.org/article/Local_field_potential or https://doi.org/10.3389/fncom.2016.00065 . To plot the LFP use the ``sim.analysis.plotLFP()`` method. This allows to plot for each electrode: 1) the time-resolved LFP signal ('timeSeries'), 2) the power spectral density ('PSD'), 3) the spectrogram / time-frequency profile ('spectrogram'), 4) and the 3D locations of the electrodes overlayed over the model neurons ('locations'). See :ref:`analysis_functions` for the full list of ``plotLFP()`` arguments. -For further details see Tutorial 9: http://neurosimlab.org/netpyne/tutorial.html#recording-and-plotting-lfps-tutorial-9 +For further details see Tutorial 9: http://netpyne.org/tutorial.html#tutorial-9-recording-and-plotting-lfps ## Setup and execution Requires NEURON with Python and MPI support. -1. Type `nrnivmodl mod` to create a directory called either i686 or x86_64, depending on your computer's architecture. -2. To run the single cell example type: `python lfp_cell.py` or `mpiexec -np [num_proc] nrniv -mpi lfp_cell.py` -2. To run the network cell example type: `python lfp_net.py` or `mpiexec -np [num_proc] nrniv -mpi lfp_net.py` +1. Type `nrnivmodl mod` to compile .mod files (it will create a directory called either i686 or x86_64, depending on your computer's architecture). +2. To run the single cell example type: `python src/cell/init.py` or `mpiexec -np [num_proc] nrniv -mpi src/cell/init.py` +2. To run the network cell example type: `python src/net/init.py` or `mpiexec -np [num_proc] nrniv -mpi src/net/init.py` ## Overview of file structure: -* /cell_lfp.py: Code to run example of LFP recording in single cell. +* src/cell/init.py: Code to run example of LFP recording in single cell. -* /net_lfp.py: Code to run example of LFP recording in network. +* src/net/init.py: Code to run example of LFP recording in network. -* /PT5B_reduced_cellParams.json: Parameters of cell used for single cell example (netParams.cellParams). +* cells/PT5B_reduced_cellParams.json: Parameters of cell used for single cell example (netParams.cellParams). -* /IT2_reduced_cellParams.json: Parameters of cell used for network example (netParams.cellParams). +* cells/IT2_reduced_cellParams.json: Parameters of cell used for network example (netParams.cellParams). * /mod/: Mod files required for single cell and network models. diff --git a/examples/LFPrecording/src/cell/cfg.py b/examples/LFPrecording/src/cell/cfg.py index 9cca70b22..4a4b8325c 100644 --- a/examples/LFPrecording/src/cell/cfg.py +++ b/examples/LFPrecording/src/cell/cfg.py @@ -14,7 +14,5 @@ cfg.analysis['plotTraces'] = {'include': [('E',0)], 'oneFigPer':'cell', 'overlay': False, 'figSize': (5,6),'saveFig': True} # Plot recorded traces for this list of cells cfg.analysis['plotLFP'] = {'includeAxon': False, 'plots': ['timeSeries', 'locations'], 'figSize': (7.5,13.5), 'saveFig': True} -#simConfig.analysis['getCSD'] = {'timeRange': [10,45],'spacing_um': 150, 'vaknin': True} -#simConfig.analysis['plotCSD'] = {'timeRange': [10,45]} -#sim.analysis.getCSD(...args...) -#simConfig.analysis['plotCSD'] = {} +# cfg.analysis['plotCSD'] = {'timeRange': [10,45],'spacing_um': 150, 'vaknin': True, 'saveFig': True} +# cfg.analysis['plotCSDTimeSeries'] = {'saveFig': True} diff --git a/examples/LFPrecording/src/net/cfg.py b/examples/LFPrecording/src/net/cfg.py index 1ce8feab7..cf45fb3dc 100644 --- a/examples/LFPrecording/src/net/cfg.py +++ b/examples/LFPrecording/src/net/cfg.py @@ -22,5 +22,6 @@ cfg.analysis['plotLFP'] = {'includeAxon': False, 'figSize': (6,10), 'timeRange': [100,1000], 'saveFig': True} # optional: 'pop': 'E4' #simConfig.analysis['getCSD'] = {'spacing_um': 200, 'timeRange': [100,3000], 'vaknin': True} #simConfig.analysis['plotLFP'] = {'includeAxon': False, 'figSize': (6,10), 'timeRange':[100,900], 'minFreq': 10, 'maxFreq':60, 'norm':1, 'plots': ['spectrogram'], 'showFig': True} -#simConfig.analysis['plotLFP'] = {'includeAxon': False, 'figSize': (6,10), 'timeRange':[100,900], 'plots': ['spectrogram'], 'showFig': True} -#simConfig.analysis['plotCSD'] = True #{'timeRange':[100,200]} +# cfg.analysis['plotCSD'] = {'saveFig': True, 'spacing_um': 150, 'overlay': 'CSD'} +# cfg.analysis['plotCSDTimeSeries'] = {'saveFig': True} +# cfg.analysis['plotCSDPSD'] = {'saveFig': True} diff --git a/examples/M1detailed/cells/PT5B_full_cellParams.pkl b/examples/M1detailed/cells/PT5B_full_cellParams.pkl index 5006415fb..819681a7f 100644 Binary files a/examples/M1detailed/cells/PT5B_full_cellParams.pkl and b/examples/M1detailed/cells/PT5B_full_cellParams.pkl differ diff --git a/examples/NeuroMLImport/RS.mod b/examples/NeuroMLImport/RS.mod index 8e63abef6..7b061fa01 100644 --- a/examples/NeuroMLImport/RS.mod +++ b/examples/NeuroMLImport/RS.mod @@ -3,9 +3,9 @@ TITLE Mod file for component: Component(id=RS type=izhikevich2007Cell) COMMENT This NEURON file has been generated by org.neuroml.export (see https://github.com/NeuroML/org.neuroml.export) - org.neuroml.export v1.9.0 - org.neuroml.model v1.9.0 - jLEMS v0.10.7 + org.neuroml.export v1.10.0 + org.neuroml.model v1.10.0 + jLEMS v0.11.0 ENDCOMMENT @@ -13,7 +13,7 @@ NEURON { POINT_PROCESS RS - NONSPECIFIC_CURRENT i : To ensure v of section follows v_I + NONSPECIFIC_CURRENT i : To ensure v of section follows v_I RANGE v0 : parameter RANGE k : parameter RANGE vr : parameter @@ -24,13 +24,9 @@ NEURON { RANGE c : parameter RANGE d : parameter RANGE C : parameter - RANGE iSyn : exposure - RANGE iMemb : exposure - RANGE copy_v : copy of v on section - } UNITS { @@ -42,47 +38,45 @@ UNITS { (mV) = (millivolt) (mS) = (millisiemens) (uS) = (microsiemens) + (nF) = (nanofarad) (molar) = (1/liter) (kHz) = (kilohertz) (mM) = (millimolar) (um) = (micrometer) (umol) = (micromole) + (pC) = (picocoulomb) (S) = (siemens) } PARAMETER { - v0 = -60 (mV) - k = 7.0E-4 (uS / mV) - vr = -60 (mV) - vt = -40 (mV) - vpeak = 35 (mV) - a = 0.030000001 (kHz) - b = -0.002 (uS) - c = -50 (mV) - d = 0.1 (nA) - C = 1.0E-4 (microfarads) + v0 = -60 (mV) : was: -0.06 (voltage) + k = 7.0E-4 (uS / mV) : was: 7.0E-7 (conductance_per_voltage) + vr = -60 (mV) : was: -0.06 (voltage) + vt = -40 (mV) : was: -0.04 (voltage) + vpeak = 35 (mV) : was: 0.035 (voltage) + a = 0.030000001 (kHz) : was: 30.0 (per_time) + b = -0.002 (uS) : was: -2.0E-9 (conductance) + c = -50 (mV) : was: -0.05 (voltage) + d = 0.1 (nA) : was: 1.0E-10 (current) + C = 0.1 (nF) : was: 1.0E-10 (capacitance) } ASSIGNED { v (mV) - i (mA/cm2) - - copy_v (mV) - - v_I (nA) - - iSyn (nA) : derived variable + i (nA) : the point process current - iMemb (nA) : derived variable + v_I (mV/ms) : for rate of change of voltage + iSyn (nA) : derived variable + iMemb (nA) : derived variable rate_v (mV/ms) rate_u (nA/ms) } STATE { - u (nA) + u (nA) : dimension: current } @@ -101,7 +95,6 @@ BREAKPOINT { SOLVE states METHOD cnexp - copy_v = v i = v_I * C } diff --git a/examples/NeuroMLImport/poissonFiringSyn.mod b/examples/NeuroMLImport/poissonFiringSyn.mod index ed2686a6e..46f6146dd 100644 --- a/examples/NeuroMLImport/poissonFiringSyn.mod +++ b/examples/NeuroMLImport/poissonFiringSyn.mod @@ -3,9 +3,9 @@ TITLE Mod file for component: Component(id=poissonFiringSyn type=poissonFiringSy COMMENT This NEURON file has been generated by org.neuroml.export (see https://github.com/NeuroML/org.neuroml.export) - org.neuroml.export v1.9.0 - org.neuroml.model v1.9.0 - jLEMS v0.10.7 + org.neuroml.export v1.10.0 + org.neuroml.model v1.10.0 + jLEMS v0.11.0 ENDCOMMENT @@ -16,7 +16,6 @@ NEURON { RANGE averageRate : parameter RANGE averageIsi : parameter RANGE SMALL_TIME : parameter - RANGE i : exposure RANGE syn0_tauRise : parameter RANGE syn0_tauDecay : parameter @@ -24,9 +23,7 @@ NEURON { RANGE syn0_waveformFactor : parameter RANGE syn0_gbase : parameter RANGE syn0_erev : parameter - RANGE syn0_g : exposure - RANGE syn0_i : exposure RANGE iSyn : derived variable : Based on netstim.mod @@ -43,11 +40,13 @@ UNITS { (mV) = (millivolt) (mS) = (millisiemens) (uS) = (microsiemens) + (nF) = (nanofarad) (molar) = (1/liter) (kHz) = (kilohertz) (mM) = (millimolar) (um) = (micrometer) (umol) = (micromole) + (pC) = (picocoulomb) (S) = (siemens) } @@ -55,28 +54,24 @@ UNITS { PARAMETER { weight = 1 - averageRate = 0.05 (kHz) - averageIsi = 20 (ms) - SMALL_TIME = 1.0E-9 (ms) - syn0_tauRise = 0.5 (ms) - syn0_tauDecay = 10 (ms) - syn0_peakTime = 1.5767012 (ms) - syn0_waveformFactor = 1.2324 - syn0_gbase = 0.002 (uS) - syn0_erev = 0 (mV) + averageRate = 0.05 (kHz) : was: 50.0 (per_time) + averageIsi = 20 (ms) : was: 0.02 (time) + SMALL_TIME = 1.0E-9 (ms) : was: 1.0000000000000002E-12 (time) + syn0_tauRise = 0.5 (ms) : was: 5.0E-4 (time) + syn0_tauDecay = 10 (ms) : was: 0.01 (time) + syn0_peakTime = 1.5767012 (ms) : was: 0.0015767011966073639 (time) + syn0_waveformFactor = 1.2324 : was: 1.232399909181873 (none) + syn0_gbase = 0.002 (uS) : was: 2.0E-9 (conductance) + syn0_erev = 0 (mV) : was: 0.0 (voltage) } ASSIGNED { v (mV) - isi (ms) : Not a state variable as far as Neuron's concerned... - - syn0_g (uS) : derived variable - - syn0_i (nA) : derived variable - - iSyn (nA) : derived variable - - i (nA) : derived variable + isi (ms) : Not a state variable as far as Neuron's concerned... + syn0_g (uS) : derived variable + syn0_i (nA) : derived variable + iSyn (nA) : derived variable + i (nA) : derived variable rate_tsince (ms/ms) rate_tnextUsed (ms/ms) rate_tnextIdeal (ms/ms) @@ -86,11 +81,11 @@ ASSIGNED { } STATE { - tsince (ms) - tnextIdeal (ms) - tnextUsed (ms) - syn0_A - syn0_B + tsince (ms) : dimension: time + tnextIdeal (ms) : dimension: time + tnextUsed (ms) : dimension: time + syn0_A : dimension: none + syn0_B : dimension: none } diff --git a/examples/NeuroMLImport/syn0.mod b/examples/NeuroMLImport/syn0.mod index 264936489..35796ee8c 100644 --- a/examples/NeuroMLImport/syn0.mod +++ b/examples/NeuroMLImport/syn0.mod @@ -3,9 +3,9 @@ TITLE Mod file for component: Component(id=syn0 type=expTwoSynapse) COMMENT This NEURON file has been generated by org.neuroml.export (see https://github.com/NeuroML/org.neuroml.export) - org.neuroml.export v1.9.0 - org.neuroml.model v1.9.0 - jLEMS v0.10.7 + org.neuroml.export v1.10.0 + org.neuroml.model v1.10.0 + jLEMS v0.11.0 ENDCOMMENT @@ -17,9 +17,7 @@ NEURON { RANGE waveformFactor : parameter RANGE gbase : parameter RANGE erev : parameter - RANGE g : exposure - RANGE i : exposure @@ -36,23 +34,25 @@ UNITS { (mV) = (millivolt) (mS) = (millisiemens) (uS) = (microsiemens) + (nF) = (nanofarad) (molar) = (1/liter) (kHz) = (kilohertz) (mM) = (millimolar) (um) = (micrometer) (umol) = (micromole) + (pC) = (picocoulomb) (S) = (siemens) } PARAMETER { - tauRise = 0.5 (ms) - tauDecay = 10 (ms) - peakTime = 1.5767012 (ms) - waveformFactor = 1.2324 - gbase = 0.002 (uS) - erev = 0 (mV) + tauRise = 0.5 (ms) : was: 5.0E-4 (time) + tauDecay = 10 (ms) : was: 0.01 (time) + peakTime = 1.5767012 (ms) : was: 0.0015767011966073639 (time) + waveformFactor = 1.2324 : was: 1.232399909181873 (none) + gbase = 0.002 (uS) : was: 2.0E-9 (conductance) + erev = 0 (mV) : was: 0.0 (voltage) } ASSIGNED { @@ -60,18 +60,16 @@ ASSIGNED { v (mV) celsius (degC) temperature (K) - - g (uS) : derived variable - - i (nA) : derived variable + g (uS) : derived variable + i (nA) : derived variable rate_A (/ms) rate_B (/ms) } STATE { - A - B + A : dimension: none + B : dimension: none } diff --git a/examples/batchSGE/README.md b/examples/batchSGE/README.md new file mode 100644 index 000000000..0203b4c3a --- /dev/null +++ b/examples/batchSGE/README.md @@ -0,0 +1,55 @@ +# Running batch simulations using NetPyNE and SGE. + +An example of running NetPyNE batch with SGE +This model is a modification of tutorial 8: +http://www.netpyne.org/tutorial.html#tutorial-8-running-batch-simulations + +to run on an HPC implementing Sun Grid Engine +https://docs.oracle.com/cd/E19279-01/820-3257-12/n1ge.html +https://gridscheduler.sourceforge.net/htmlman/htmlman1/qsub.html + + +It implements each 4-core simulation on a compute node. + +## Requirements +- NEURON with MPI support +- NetPyNE + +## Quick steps to run batch sims + +1) python src/batch.py + +results are generated in tut8_data +simulation logs are generated in ~/qsub/ + +## Code structure: + +* src/netParams.py: Defines the network model. Includes "fixed" parameter values of cells, synapses, connections, stimulation, etc. Changes to this file should only be made for relatively stable improvements of the model. Parameter values that will be varied systematically to explore or tune the model should be included here by referencing the appropiate variable in the simulation configuration (cfg) module. Only a single netParams file is required for each batch of simulations. + +* src/cfg.py: Simulation configuration module. Includes parameter values for each simulation run such as duration, dt, recording parameters etc. Also includes the model parameters that are being varied to explore or tune the network. When running a batch, NetPyNE will create one cfg file for each parameter configuration. + +* src/init.py: Sequence of commands to run a single simulation. Can be executed via 'ipython init.py'. When running a batch, NetPyNE will call init.py multiple times, pass a different cfg file for each of the parameter configurations explored. + +* src/batch.py: Defines the parameters and parameter values explored in the batch, as well as the run configuration. +``` + b.runCfg = {'type': 'hpc_sge', + 'jobName': 'my_batch', + 'cores': 4, + 'script': 'init.py', + 'skip': True} +``` + +* arguments accepted by runCfg are described in netpyne/batch/grid.py +``` + sge_args = { + 'jobName': jobName, + 'walltime': walltime, + 'vmem': '32G', + 'queueName': 'cpu.q', + 'cores': 1, + 'custom': '', + 'mpiCommand': 'mpiexec', + 'log': "~/qsub/{}.run".format(jobName) + } +``` + diff --git a/examples/batchSGE/src/batch.py b/examples/batchSGE/src/batch.py new file mode 100644 index 000000000..8b6e19b87 --- /dev/null +++ b/examples/batchSGE/src/batch.py @@ -0,0 +1,30 @@ +from netpyne import specs +from netpyne.batch import Batch + +def batchTauWeight(): + # Create variable of type ordered dictionary (NetPyNE's customized version) + params = specs.ODict() + + # fill in with parameters to explore and range of values (key has to coincide with a variable in simConfig) + params['synMechTau2'] = [3.0, 5.0, 7.0] + params['connWeight'] = [0.005, 0.01, 0.15] + + # create Batch object with parameters to modify, and specifying files to use + b = Batch(params=params, cfgFile='src/cfg.py', netParamsFile='src/params.py',) + + # Set output folder, grid method (all param combinations), and run configuration + b.batchLabel = 'tauWeight' + b.saveFolder = 'tut8_data' + b.method = 'grid' + b.runCfg = {'type': 'hpc_sge', + 'jobName': 'my_batch', + 'cores': 4, + 'script': 'init.py', + 'skip': True} + + # Run batch simulations + b.run() + +# Main code +if __name__ == '__main__': + batchTauWeight() \ No newline at end of file diff --git a/examples/batchSGE/src/cfg.py b/examples/batchSGE/src/cfg.py new file mode 100644 index 000000000..10ce94aaa --- /dev/null +++ b/examples/batchSGE/src/cfg.py @@ -0,0 +1,21 @@ +from netpyne import specs + +# Simulation options +cfg = specs.SimConfig() # object of class SimConfig to store simulation configuration + +cfg.duration = 1*1e3 # Duration of the simulation, in ms +cfg.dt = 0.025 # Internal integration timestep to use +cfg.verbose = False # Show detailed messages +cfg.recordTraces = {'V_soma':{'sec':'soma','loc':0.5,'var':'v'}} # Dict with traces to record +cfg.recordStep = 0.1 # Step size in ms to save data (eg. V traces, LFP, etc) +cfg.filename = 'tut8' # Set file output name +cfg.saveJson = True +cfg.printPopAvgRates = True +cfg.analysis['plotRaster'] = {'saveFig': True} # Plot a raster +cfg.analysis['plotTraces'] = {'include': [0], 'saveFig': True} # Plot recorded traces for this list of cells + +cfg.saveDataInclude = ['simData', 'simConfig', 'netParams', 'net'] + +# Variable parameters (used in netParams) +cfg.synMechTau2 = 5 +cfg.connWeight = 0.01 diff --git a/examples/batchSGE/src/init.py b/examples/batchSGE/src/init.py new file mode 100644 index 000000000..45b8c3c6b --- /dev/null +++ b/examples/batchSGE/src/init.py @@ -0,0 +1,7 @@ +from netpyne import sim + +# read cfg and netParams from command line arguments if available; otherwise use default +simConfig, netParams = sim.readCmdLineArgs(simConfigDefault='src/cfg.py', netParamsDefault='src/params.py') + +# Create network and run simulation +sim.createSimulateAnalyze(netParams=netParams, simConfig=simConfig) diff --git a/examples/batchSGE/src/params.py b/examples/batchSGE/src/params.py new file mode 100644 index 000000000..d62822498 --- /dev/null +++ b/examples/batchSGE/src/params.py @@ -0,0 +1,36 @@ +from netpyne import specs, sim + +try: + from __main__ import cfg # import SimConfig object with params from parent module +except: + from tut8_cfg import cfg # if no simConfig in parent module, import directly from tut8_cfg module + +# Network parameters +netParams = specs.NetParams() # object of class NetParams to store the network parameters + +## Cell parameters/rules +PYRcell = {'secs': {}} +PYRcell['secs']['soma'] = {'geom': {}, 'mechs': {}} # soma params dict +PYRcell['secs']['soma']['geom'] = {'diam': 18.8, 'L': 18.8, 'Ra': 123.0} # soma geometry +PYRcell['secs']['soma']['mechs']['hh'] = {'gnabar': 0.12, 'gkbar': 0.036, 'gl': 0.003, 'el': -70} # soma hh mechanism +netParams.cellParams['PYR'] = PYRcell + +## Population parameters +netParams.popParams['S'] = {'cellType': 'PYR', 'numCells': 20} +netParams.popParams['M'] = {'cellType': 'PYR', 'numCells': 20} + +## Synaptic mechanism parameters +netParams.synMechParams['exc'] = {'mod': 'Exp2Syn', 'tau1': 0.1, 'tau2': cfg.synMechTau2, 'e': 0} # excitatory synaptic mechanism + +# Stimulation parameters +netParams.stimSourceParams['bkg'] = {'type': 'NetStim', 'rate': 10, 'noise': 0.5} +netParams.stimTargetParams['bkg->PYR'] = {'source': 'bkg', 'conds': {'cellType': 'PYR'}, 'weight': 0.01, 'delay': 5, 'synMech': 'exc'} + +## Cell connectivity rules +netParams.connParams['S->M'] = { # S -> M label + 'preConds': {'pop': 'S'}, # conditions of presyn cells + 'postConds': {'pop': 'M'}, # conditions of postsyn cells + 'probability': 0.5, # probability of connection + 'weight': cfg.connWeight, # synaptic weight + 'delay': 5, # transmission delay (ms) + 'synMech': 'exc'} # synaptic mechanism diff --git a/examples/dipoleRecording/src/cfg_cell.py b/examples/dipoleRecording/src/cfg_cell.py index 6b9ee4ccc..f0e6f7526 100644 --- a/examples/dipoleRecording/src/cfg_cell.py +++ b/examples/dipoleRecording/src/cfg_cell.py @@ -16,7 +16,4 @@ cfg.analysis['plotTraces'] = {'include': [('E',0)], 'oneFigPer':'cell', 'overlay': True, 'figSize': (5,3),'saveFig': True} # Plot recorded traces for this list of cells #cfg.analysis['plotLFP'] = {'includeAxon': False, 'plots': ['timeSeries', 'locations'], 'figSize': (5,9), 'saveFig': True} -#cfg.analysis['getCSD'] = {'timeRange': [10,45],'spacing_um': 150, 'vaknin': True} -#cfg.analysis['plotCSD'] = {'timeRange': [10,45]} -#sim.analysis.getCSD(...args...) -#cfg.analysis['plotCSD'] = {} \ No newline at end of file +#cfg.analysis['plotCSD'] = {'timeRange': [10,45], 'saveFig': True} diff --git a/examples/dipoleRecording/src/init_cell.py b/examples/dipoleRecording/src/init_cell.py index e7d48cf61..a0fce582c 100644 --- a/examples/dipoleRecording/src/init_cell.py +++ b/examples/dipoleRecording/src/init_cell.py @@ -1,5 +1,5 @@ from netpyne import sim -cfg, netParams = sim.loadFromIndexFile('index_cell.npjson') +cfg, netParams = sim.loadFromIndexFile('index.npjson') sim.createSimulateAnalyze(netParams = netParams, simConfig = cfg) -#sim.analysis.plotCSD() \ No newline at end of file +#sim.analysis.plotCSD() diff --git a/netpyne/__init__.py b/netpyne/__init__.py index 906d3fb50..61069316f 100644 --- a/netpyne/__init__.py +++ b/netpyne/__init__.py @@ -4,12 +4,13 @@ NetPyNE consists of a number of sub-packages and modules. """ + ############################################################################### # Note: this branch should NOT be merged to master/development # Updates for NeuroML support will be added on https://github.com/Neurosim-lab/netpyne/tree/neuroml_updates # and those changes pushed to development. This branch is for stable releases for # use on osbv2 -__version__ = '1.0.4.1+osbv2' +__version__ = '1.0.5+osbv2' ############################################################################### import os, sys diff --git a/netpyne/analysis/__init__.py b/netpyne/analysis/__init__.py index 95351848f..a719768f0 100644 --- a/netpyne/analysis/__init__.py +++ b/netpyne/analysis/__init__.py @@ -81,7 +81,6 @@ print('Warning: could not import interactive plotting functions; make sure the "bokeh" package is installed.') # Import CSD-related functions -# from .csd import getCSD, plotCSD from .csd import prepareCSD # Import dipole-related functions diff --git a/netpyne/analysis/csd.py b/netpyne/analysis/csd.py index 46fd43638..ba61f7100 100644 --- a/netpyne/analysis/csd.py +++ b/netpyne/analysis/csd.py @@ -19,10 +19,7 @@ import numpy as np import scipy - -import matplotlib -from matplotlib import pyplot as plt -from matplotlib import ticker as ticker +from numbers import Number import json import sys import os @@ -34,9 +31,109 @@ from .utils import exception, _saveFigData +def getBandpass( + lfps, + sampr, + minf=0.05, + maxf=300): + """ + Function to bandpass filter data + + Parameters + ---------- + lfps : list or array + LFP signal data arranged spatially in a column. + **Default:** *required* + + sampr : float + The data sampling rate. + **Default:** *required* + + minf : float + The high-pass filter frequency (Hz). + **Default:** ``0.05`` + + maxf : float + The low-pass filter frequency (Hz). + **Default:** ``300`` + + + Returns + ------- + data : array + The bandpass-filtered data. + + """ + + datband = [] + for i in range(len(lfps[0])):datband.append(bandpass(lfps[:,i], minf, maxf, df=sampr, zerophase=True)) + datband = np.array(datband) + return datband + +def vakninCorrection(x): + """ + Function to perform the Vaknin correction for CSD analysis + + Allows CSD to be performed on all N contacts instead of N-2 contacts (see Vaknin et al (1988) for more details). + + Parameters + ---------- + x : array + Data to be corrected. + **Default:** *required* + + + Returns + ------- + data : array + The corrected data. + + """ + + # Preallocate array with 2 more rows than input array + x_new = np.zeros((x.shape[0]+2, x.shape[1])) + + # Duplicate first and last row of x into first and last row of x_new + x_new[0, :] = x[0, :] + x_new[-1, :] = x[-1, :] + + # Duplicate all of x into middle rows of x_neww + x_new[1:-1, :] = x + + return x_new + +def removeMean(x, ax=1): + """ + Function to subtract the mean from an array or list + + Parameters + ---------- + x : array + Data to be processed. + **Default:** *required* + + ax : int + The axis to remove the mean across. + **Default:** ``1`` + + + Returns + ------- + data : array + The processed data. + + """ + + mean = np.mean(x, axis=ax, keepdims=True) + x -= mean + + + @exception def prepareCSD( sim=None, + timeRange=None, + electrodes=['avg', 'all'], pop=None, dt=None, sampr=None, @@ -49,6 +146,7 @@ def prepareCSD( getAllData=True, **kwargs ): + """ Function to prepare data for plotting of current source density (CSD) data @@ -58,6 +156,13 @@ def prepareCSD( sim : NetPyNE object **Default:** ``None`` + timeRange: list + List of length two, with timeRange[0] as beginning of desired timeRange, and timeRange[1] as the end + **Default:** ``None`` retrieves timeRange = [0, sim.cfg.duration] + + electrodes : list of electrodes to look at CSD data + **Default:** ['avg', 'all'] + pop : str Retrieves CSD data from a specific cell population **Default:** ``None`` retrieves overall CSD data @@ -88,7 +193,7 @@ def prepareCSD( norm : bool Subtracts the mean from the CSD data - **Default:** ``False`` + **Default:** ``False`` --> ``True`` saveData : bool Saves CSD data to sim object @@ -107,9 +212,10 @@ def prepareCSD( except: raise Exception('Cannot access sim') - ## Get LFP data from sim and instantiate as a numpy array - simDataCategories = sim.allSimData.keys() + ## Get LFP data from sim and instantiate as a numpy array + simDataCategories = sim.allSimData.keys() + if pop is None: if 'LFP' in simDataCategories: LFPData = np.array(sim.allSimData['LFP']) @@ -129,47 +235,43 @@ def prepareCSD( if dt is None: dt = sim.cfg.recordStep + # slice data by timeRange, if relevant + if timeRange is None: + timeRange = [0, sim.cfg.duration] + else: + LFPData = LFPData[int(timeRange[0] / sim.cfg.recordStep) : int(timeRange[1] / sim.cfg.recordStep), :] + # Sampling rate of data recording during the simulation if sampr is None: # divide by 1000.0 to turn denominator from units of ms to s sampr = 1.0 / (dt / 1000.0) # dt == sim.cfg.recordStep, unless specified otherwise by user + # Spacing between electrodes (in microns) if spacing_um is None: - spacing_um = sim.cfg.recordLFP[1][1] - sim.cfg.recordLFP[0][1] + # if not specified, use average spacing along y coord (depth) + yCoords = np.array(sim.cfg.recordLFP)[:,1] + spacing_um = (yCoords.max() - yCoords.min()) / (len(yCoords) - 1) # Convert spacing from microns to mm spacing_mm = spacing_um / 1000 - print('dt, sampr, spacing_um, spacing_mm values determined') - - ## This retrieves: - # LFPData (as an array) - # dt --> recording time step (in ms) - # sampr --> sampling rate of data recording (in Hz) - # spacing_um --> spacing btwn electrodes (in um) + # print('dt, sampr, spacing_um, spacing_mm values determined') - #################################### - - # Bandpass filter the LFP data with getbandpass() fx defined above - datband = getbandpass(LFPData, sampr, minf, maxf) - - # Take CSD along smaller dimension - if datband.shape[0] > datband.shape[1]: - ax = 1 - else: - ax = 0 + # Bandpass filter the LFP data with getBandpass() fx defined above + datband = getBandpass(LFPData, sampr, minf, maxf) + # now each row is an electrode - `datband` shape is (N_electrodes, N_timesteps) # Vaknin correction if vaknin: - datband = Vaknin(datband) + datband = vakninCorrection(datband) # norm data if norm: - removemean(datband, ax=ax) + removeMean(datband, ax=0) - # now each column (or row) is an electrode -- take CSD along electrodes - CSDData = -np.diff(datband, n=2, axis=ax) / spacing_mm**2 + # take CSD along electrodes dimension + CSDData = -np.diff(datband, n=2, axis=0) / spacing_mm**2 ##### SAVE DATA ####### # Add CSDData to sim.allSimData for later access @@ -181,105 +283,9 @@ def prepareCSD( sim.allSimData['CSDPops'][pop] = CSDData # return CSD_data or all data - if getAllData is True: + if getAllData: return CSDData, LFPData, sampr, spacing_um, dt - elif getAllData is False: - return CSDData - - -def getbandpass(lfps, sampr, minf=0.05, maxf=300): - """ - Function to bandpass filter data - - Parameters - ---------- - lfps : list or array - LFP signal data arranged spatially in a column. - **Default:** *required* - - sampr : float - The data sampling rate. - **Default:** *required* - - minf : float - The high-pass filter frequency (Hz). - **Default:** ``0.05`` - - maxf : float - The low-pass filter frequency (Hz). - **Default:** ``300`` - - - Returns - ------- - data : array - The bandpass-filtered data. - - """ - - datband = [] - for i in range(len(lfps[0])): - datband.append(bandpass(lfps[:, i], minf, maxf, df=sampr, zerophase=True)) - - datband = np.array(datband) - - return datband - - -def Vaknin(x): - """ - Function to perform the Vaknin correction for CSD analysis - - Allows CSD to be performed on all N contacts instead of N-2 contacts (see Vaknin et al (1988) for more details). - - Parameters - ---------- - x : array - Data to be corrected. - **Default:** *required* - - - Returns - ------- - data : array - The corrected data. - - """ - - # Preallocate array with 2 more rows than input array - x_new = np.zeros((x.shape[0] + 2, x.shape[1])) - - # Duplicate first and last row of x into first and last row of x_new - x_new[0, :] = x[0, :] - x_new[-1, :] = x[-1, :] - - # Duplicate all of x into middle rows of x_neww - x_new[1:-1, :] = x - - return x_new - - -def removemean(x, ax=1): - """ - Function to subtract the mean from an array or list - - Parameters - ---------- - x : array - Data to be processed. - **Default:** *required* - - ax : int - The axis to remove the mean across. - **Default:** ``1`` - - - Returns - ------- - data : array - The processed data. - - """ - - mean = np.mean(x, axis=ax, keepdims=True) - x -= mean + else: + from .lfp import prepareDataPerElectrode + CSDData = CSDData.T # to match the shape expected by prepareDataPerElectrode + return prepareDataPerElectrode(CSDData, electrodes, timeRange, sim) diff --git a/netpyne/analysis/csd_legacy.py b/netpyne/analysis/csd_legacy.py deleted file mode 100644 index 6c8f58eac..000000000 --- a/netpyne/analysis/csd_legacy.py +++ /dev/null @@ -1,651 +0,0 @@ -""" -Module with functions to extract and plot CSD info from LFP data - -""" - -from __future__ import print_function -from __future__ import division -from __future__ import unicode_literals -from __future__ import absolute_import - -from future import standard_library - -standard_library.install_aliases() - -try: - basestring -except NameError: - basestring = str - -import numpy as np -import scipy - -import matplotlib -from matplotlib import pyplot as plt -from matplotlib import ticker as ticker -import json -import sys -import os -from collections import OrderedDict -import warnings -from scipy.fftpack import hilbert -from scipy.signal import cheb2ord, cheby2, convolve, get_window, iirfilter, remez, decimate -from .filter import lowpass, bandpass -from .utils import exception, _saveFigData - - -def getbandpass(lfps, sampr, minf=0.05, maxf=300): - """ - Function to bandpass filter data - - Parameters - ---------- - lfps : list or array - LFP signal data arranged spatially in a column. - **Default:** *required* - - sampr : float - The data sampling rate. - **Default:** *required* - - minf : float - The high-pass filter frequency (Hz). - **Default:** ``0.05`` - - maxf : float - The low-pass filter frequency (Hz). - **Default:** ``300`` - - - Returns - ------- - data : array - The bandpass-filtered data. - - """ - - datband = [] - for i in range(len(lfps[0])): - datband.append(bandpass(lfps[:, i], minf, maxf, df=sampr, zerophase=True)) - - datband = np.array(datband) - - return datband - - -def Vaknin(x): - """ - Function to perform the Vaknin correction for CSD analysis - - Allows CSD to be performed on all N contacts instead of N-2 contacts (see Vaknin et al (1988) for more details). - - Parameters - ---------- - x : array - Data to be corrected. - **Default:** *required* - - - Returns - ------- - data : array - The corrected data. - - """ - - # Preallocate array with 2 more rows than input array - x_new = np.zeros((x.shape[0] + 2, x.shape[1])) - - # Duplicate first and last row of x into first and last row of x_new - x_new[0, :] = x[0, :] - x_new[-1, :] = x[-1, :] - - # Duplicate all of x into middle rows of x_neww - x_new[1:-1, :] = x - - return x_new - - -def removemean(x, ax=1): - """ - Function to subtract the mean from an array or list - - Parameters - ---------- - x : array - Data to be processed. - **Default:** *required* - - ax : int - The axis to remove the mean across. - **Default:** ``1`` - - - Returns - ------- - data : array - The processed data. - - """ - - mean = np.mean(x, axis=ax, keepdims=True) - x -= mean - - -# returns CSD in units of mV/mm**2 (assuming lfps are in mV) -def getCSD( - LFP_input_data=None, - LFP_input_file=None, - sampr=None, - dt=None, - spacing_um=None, - minf=0.05, - maxf=300, - norm=True, - vaknin=True, - save_to_sim=True, - getAllData=False, -): - """ - Function to extract CSD values from simulated LFP data - - Parameters - ---------- - LFP_input_data : list or numpy array - LFP data provided by user (mV). Each element of the list/array must be a list/array containing LFP data for an electrode. - **Default:** ``None`` pulls the data from the current NetPyNE sim object. - - LFP_input_file : str - Location of a JSON data file with LFP data. - **Default:** ``None`` - - sampr : float - Sampling rate for data recording (Hz). - **Default:** ``None`` uses ``1.0/sim.cfg.recordStep`` from the current NetPyNE sim object. - - dt : float - Time between recording points (ms). - **Default:** ``None`` uses ``sim.cfg.recordStep`` from the current NetPyNE sim object. - - spacing_um : float - Electrode contact spacing in units of microns. - **Default:** ``None`` pulls the information from the current NetPyNE sim object. If the data is empirical, defaults to ``100`` (microns). - - minf : float - Minimum frequency for bandpass filter (Hz). - **Default:** ``0.05`` - - maxf : float - Maximum frequency cutoff for bandpass filter (Hz). - **Default:** ``300`` - - norm : bool - Whether to subtract the mean from the signal or not. - **Default:** ``True`` subtracts the mean. - - vaknin : bool - Whether or not to to perform the Vaknin correction for CSD analysis. - **Default** ``True`` performs the correction. - - save_to_sim : bool - Whether to attempt to store CSD values in sim.allSimData, if this exists. - **Default** ``True`` attempts to store the data. - - getAllData : bool - True will have this function return dt, tt, timeRange, sampr, spacing_um, lfp_data, and CSD_data. - **Default** ``False`` returns only CSD_data. - - """ - - ### DEFAULT -- CONDITION 1 : LFP DATA COMES FROM SIMULATION ### - - # Get LFP data from simulation - if LFP_input_data is None and LFP_input_file is None: - - try: - from .. import sim - except: - raise Exception('No LFP input data, input file, or existing simulation. Cannot calculate CSD.') - - # time step used in simulation recording (in ms) - if dt is None: - dt = sim.cfg.recordStep - - # Get LFP data from sim and instantiate as a numpy array - sim_data_categories = sim.allSimData.keys() - - if 'LFP' in sim_data_categories: - lfp_data = np.array(sim.allSimData['LFP']) - else: - raise Exception('NO LFP DATA!! Need to re-run simulation with cfg.recordLFP enabled.') - - # Sampling rate of data recording during the simulation - if sampr is None: - # divide by 1000.0 to turn denominator from units of ms to s - sampr = 1.0 / (dt / 1000.0) # sim.cfg.recordStep --> == dt - - # Spacing between electrodes (in microns) - if spacing_um is None: - spacing_um = sim.cfg.recordLFP[1][1] - sim.cfg.recordLFP[0][1] - - ## This retrieves: - # lfp_data (as an array) - # dt --> recording time step (in ms) - # sampr --> sampling rate of data recording (in Hz) - # spacing_um --> spacing btwn electrodes (in um) - - ### CONDITION 2 : LFP DATA FROM SIM .JSON FILE ### # Note: need to expand capability to include a list of multiple files - - # load sim data from JSON file - elif LFP_input_data is None and '.json' in LFP_input_file: - - data = {} - with open(LFP_input_file) as file: - data['json_input_data'] = json.load(file) - - ## FOR MULTIPLE FILES - # for x in LFP_input_file: - # with open(x) as file: - # data[x] = json.load(file) - - # extract LFP data (only works in the 1 input file scenario; expand capability for multiple files) - for key in data.keys: - lfp_data_list = data[key]['simData']['LFP'] - - # cast LFP data as Numpy array - lfp_data = np.array(lfp_data_list) - - # get CSD data and relevant plotting params - csd_data = {} - for i in data.keys(): - csd_data[i] = {} - - if sampr is None: - csd_data[i]['sampr'] = 1.0 / ((data[i]['simConfig']['recordStep']) / 1000.0) - sampr = csd_data[i]['sampr'] - else: - csd_data[i]['sampr'] = sampr - - if spacing_um is None: - csd_data[i]['spacing_um'] = ( - data[i]['simConfig']['recordLFP'][1][1] - data[i]['simConfig']['recordLFP'][0][1] - ) - spacing_um = csd_data[i]['spacing_um'] - else: - csd_data[i]['spacing_um'] = spacing_um - - if dt is None: - csd_data[i]['dt'] = data[i]['simConfig']['recordStep'] - dt = csd_data[i]['dt'] - else: - csd_data[i]['dt'] = dt - - ## This retrieves: - # lfp_data (as a list) - # dt --> recording time step (in ms) - # sampr --> sampling rate of data recording (in Hz) - # spacing_um --> spacing btwn electrodes (in um) - - ### CONDITION 3 : ARBITRARY LFP DATA ### # NOTE: for condition 3 --> need to also retrieve the dt, sampr, and spacing_um !! - - # get lfp_data and cast as numpy array - elif len(LFP_input_data) > 0 and LFP_input_file is None: - lfp_data = np.array(LFP_input_data) - - #################################### - - # Convert spacing from microns to mm - spacing_mm = spacing_um / 1000 - - # Bandpass filter the LFP data with getbandpass() fx defined above - datband = getbandpass(lfp_data, sampr, minf, maxf) - - # Take CSD along smaller dimension - if datband.shape[0] > datband.shape[1]: - ax = 1 - else: - ax = 0 - print('ax = ' + str(ax)) - - # Vaknin correction - if vaknin: - datband = Vaknin(datband) - - # norm data - if norm: - removemean(datband, ax=ax) - - # now each column (or row) is an electrode -- take CSD along electrodes - CSD_data = -np.diff(datband, n=2, axis=ax) / spacing_mm**2 - - # noBandpass trial - datband_noBandpass = lfp_data.T - - if datband_noBandpass.shape[0] > datband_noBandpass.shape[1]: - ax = 1 - else: - ax = 0 - - if vaknin: - datband_noBandpass = Vaknin(datband_noBandpass) - - if norm: - removemean(datband_noBandpass, ax=ax) - - CSD_data_noBandpass = -np.diff(datband_noBandpass, n=2, axis=ax) / spacing_mm**2 - - # Add CSD and other param values to sim.allSimData for later access - if save_to_sim is True: - try: - from .. import sim - - sim.allSimData['CSD'] = {} - sim.allSimData['CSD']['sampr'] = sampr - sim.allSimData['CSD']['spacing_um'] = spacing_um - sim.allSimData['CSD']['CSD_data'] = CSD_data - sim.allSimData['CSD']['CSD_data_noBandpass'] = CSD_data_noBandpass - except: - print('NOTE: No sim.allSimData construct available to store CSD data.') - - # return CSD_data or all data - if getAllData is True: - return lfp_data, CSD_data, sampr, spacing_um, dt - if getAllData is False: - return CSD_data - - -# PLOTTING CSD - - -@exception -def plotCSD( - CSD_data=None, - LFP_input_data=None, - overlay=None, - timeRange=None, - sampr=None, - stim_start_time=None, - spacing_um=None, - ymax=None, - dt=None, - hlines=False, - layerLines=False, - layerBounds=None, - smooth=True, - fontSize=12, - figSize=(8, 8), - dpi=200, - saveFig=True, - showFig=True, -): - """ - Function to plot CSD values extracted from simulated LFP data - - Parameters - ---------- - CSD_data : list or array - CSD_data for plotting. - **Default:** ``None`` - - LFP_input_data : list or numpy array - LFP data provided by user (mV). Each element of the list/array must be a list/array containing LFP data for an electrode. - **Default:** ``None`` pulls the data from the current NetPyNE sim object. - - - overlay : str - Option to include LFP data overlaid on CSD color map plot. - **Default:** ``None`` provides no overlay - OPTIONS are 'LFP' or 'CSD' - - timeRange : list - Time range to plot [start, stop]. - **Default:** ``None`` plots entire time range - - sampr : float - Sampling rate for data recording (Hz). - **Default:** ``None`` uses ``1.0/sim.cfg.recordStep`` from the current NetPyNE sim object. - - stim_start_time : float - Time when stimulus is applied (ms). - **Default:** ``None`` does not add anything to plot. - **Options:** a float adds a vertical dashed line to the plot at stimulus onset. - - spacing_um : float - Electrode contact spacing in units of microns. - **Default:** ``None`` pulls the information from the current NetPyNE sim object. If the data is empirical, defaults to ``100`` (microns). - - ymax : float - The upper y-limit. - **Default:** ``None`` - - dt : float - Time between recording points (ms). - **Default:** ``None`` uses ``sim.cfg.recordStep`` from the current NetPyNE sim object. - - hlines : bool - Option to include horizontal lines on plot to indicate electrode positions. - **Default:** ``False`` - - layerLines : bool - Whether to plot horizontal lines over CSD plot at layer boundaries. - **Default:** ``False`` - - layerBounds : dict - Dictionary containing layer labels as keys, and layer boundaries as values, e.g. {'L1':100, 'L2': 160, 'L3': 950, 'L4': 1250, 'L5A': 1334, 'L5B': 1550, 'L6': 2000} - **Default:** ``None`` - - saveFig : bool or str - Whether and where to save the figure. - **Default:** ``True`` autosaves the figure. - **Options:** ``'/path/filename.ext'`` saves to a custom path and filename, valid file extensions are ``'.png'``, ``'.jpg'``, ``'.eps'``, and ``'.tiff'``. - - showFig : bool - Whether to show the figure. - **Default:** ``True`` - - """ - - print('Plotting CSD... ') - - # DEFAULT -- CONDITION 1 : GET CSD DATA FROM SIM - if CSD_data is None: - - from .. import sim - - LFP_data, CSD_data, sampr, spacing_um, dt = getCSD( - getAllData=True - ) # getCSD(sampr=sampr, spacing_um=spacing_um, dt=dt, getAllData=True) - - if timeRange is None: - timeRange = [0, sim.cfg.duration] - - tt = np.arange(timeRange[0], timeRange[1], dt) - - ymax = sim.cfg.recordLFP[-1][1] + spacing_um - - LFP_data = np.array(LFP_data)[int(timeRange[0] / dt) : int(timeRange[1] / dt), :] - CSD_data = CSD_data[:, int(timeRange[0] / dt) : int(timeRange[1] / dt)] - - CSD_data_noBandpass = sim.allSimData['CSD']['CSD_data_noBandpass'][ - :, int(timeRange[0] / sim.cfg.recordStep) : int(timeRange[1] / sim.cfg.recordStep) - ] - # CSD_data_noBandpass = CSD_data_noBandpass[:,int(timeRange[0]/dt):int(timeRange[1]/dt)] - - ### The problem with this setup is that it defaults to whatever was saved in .pkl !! - # sim_data_categories = sim.allSimData.keys() - - # if 'CSD' in sim_data_categories: - - # if timeRange is None: - # timeRange = [0,sim.cfg.duration] - - # dt = sim.cfg.recordStep - # tt = np.arange(timeRange[0],timeRange[1],dt) - - # spacing_um = sim.allSimData['CSD']['spacing_um'] - # #spacing_mm = spacing_um/1000 - - # ymax = sim.cfg.recordLFP[-1][1] + spacing_um - - # # get LFP data - # LFP_data = np.array(sim.allSimData['LFP'])[int(timeRange[0]/sim.cfg.recordStep):int(timeRange[1]/sim.cfg.recordStep),:] - - # # get CSD data - # CSD_data = sim.allSimData['CSD']['CSD_data'][:,int(timeRange[0]/sim.cfg.recordStep):int(timeRange[1]/sim.cfg.recordStep)] - - # # noBandpass trial - # CSD_data_noBandpass = sim.allSimData['CSD']['CSD_data_noBandpass'][:,int(timeRange[0]/sim.cfg.recordStep):int(timeRange[1]/sim.cfg.recordStep)] - - # else: - # raise Exception('No CSD data in sim.') - - # CONDITION 2 : ARBITRARY CSD DATA - elif CSD_data is not None: - if timeRange is None: - print('MUST PROVIDE TIME RANGE in ms') - else: - print('timeRange = ' + str(timeRange)) - - if dt is None: - print('MUST PROVIDE dt in ms') - else: - print('dt = ' + str(dt)) # batch0['simConfig']['recordStep'] - - if spacing_um is None: - print('MUST PROVIDE SPACING BETWEEN ELECTRODES in MICRONS') - else: - print('spacing_um = ' + str(spacing_um)) - - if ymax is None: - print('MUST PROVIDE YMAX (MAX DEPTH) in MICRONS') - else: - print('ymax = ' + str(ymax)) - - tt = np.arange(timeRange[0], timeRange[1], dt) - LFP_data = np.array(LFP_input_data)[int(timeRange[0] / dt) : int(timeRange[1] / dt), :] - - # PLOTTING - X = np.arange(timeRange[0], timeRange[1], dt) - Y = np.arange(CSD_data.shape[0]) - - # interpolation - CSD_spline = scipy.interpolate.RectBivariateSpline(Y, X, CSD_data) - Y_plot = np.linspace(0, CSD_data.shape[0], num=1000) - Z = CSD_spline(Y_plot, X) - - # plotting options - plt.rcParams.update({'font.size': fontSize}) - xmin = int(X[0]) - xmax = int(X[-1]) + 1 # int(sim.allSimData['t'][-1]) - ymin = 0 - extent_xy = [xmin, xmax, ymax, ymin] - - # set up figure - fig = plt.figure(figsize=figSize) - - # create plots w/ common axis labels and tick marks - axs = [] - numplots = 1 - gs_outer = matplotlib.gridspec.GridSpec( - 1, 1 - ) # (2, 2, figure=fig)#, wspace=0.4, hspace=0.2, height_ratios=[20, 1]) - - for i in range(numplots): - axs.append(plt.Subplot(fig, gs_outer[i * 2 : i * 2 + 2])) - fig.add_subplot(axs[i]) - axs[i].set_xlabel('Time (ms)', fontsize=fontSize) - axs[i].tick_params(axis='y', which='major', labelsize=fontSize) - axs[i].tick_params(axis='x', which='major', labelsize=fontSize) - - ## plot interpolated CSD color map - # if smooth: - # Z = scipy.ndimage.filters.gaussian_filter(Z, sigma = 5, mode='nearest')#smooth, mode='nearest') - - spline = axs[0].imshow( - Z, extent=extent_xy, interpolation='none', aspect='auto', origin='upper', cmap='jet_r', alpha=0.9 - ) - axs[0].set_ylabel('Contact depth (um)', fontsize=fontSize) - - # OVERLAY DATA ('LFP', 'CSD', or None) & Set title of plot - if overlay is None: - print('No data being overlaid') - axs[0].set_title('Current Source Density (CSD)', fontsize=fontSize) - - elif overlay is 'CSD' or overlay is 'LFP': - nrow = LFP_data.shape[1] - gs_inner = matplotlib.gridspec.GridSpecFromSubplotSpec( - nrow, 1, subplot_spec=gs_outer[0:2], wspace=0.0, hspace=0.0 - ) - subaxs = [] - - if overlay == 'CSD': - axs[0].set_title('CSD with time series overlay', fontsize=fontSize) - for chan in range(nrow): - subaxs.append(plt.Subplot(fig, gs_inner[chan], frameon=False)) - fig.add_subplot(subaxs[chan]) - subaxs[chan].margins(0.0, 0.01) - subaxs[chan].get_xaxis().set_visible(False) - subaxs[chan].get_yaxis().set_visible(False) - subaxs[chan].plot(X, CSD_data[chan, :], color='green', linewidth=0.3) # 'blue' - - elif overlay == 'LFP': - axs[0].set_title('CSD with LFP overlay', fontsize=fontSize) - for chan in range(nrow): - subaxs.append(plt.Subplot(fig, gs_inner[chan], frameon=False)) - fig.add_subplot(subaxs[chan]) - subaxs[chan].margins(0.0, 0.01) - subaxs[chan].get_xaxis().set_visible(False) - subaxs[chan].get_yaxis().set_visible(False) - subaxs[chan].plot(X, LFP_data[:, chan], color='gray', linewidth=0.3) - - else: - print('Invalid option specified for overlay argument -- no data overlaid') - axs[0].set_title('Current Source Density (CSD)', fontsize=fontSize) - - # add horizontal lines at electrode locations - if hlines: - for i in range(len(sim.cfg.recordLFP)): - axs[0].hlines(sim.cfg.recordLFP[i][1], xmin, xmax, colors='pink', linewidth=1, linestyles='dashed') - - if layerLines: - if layerBounds is None: - print('No layer boundaries given -- will not overlay layer boundaries on CSD plot') - else: - layerKeys = [] - for i in layerBounds.keys(): - axs[0].hlines(layerBounds[i], xmin, xmax, colors='black', linewidth=1, linestyles='dotted') - layerKeys.append(i) # makes a list with names of each layer, as specified in layerBounds dict argument - - for n in range(len(layerKeys)): # label the horizontal layer lines with the proper layer label - if n == 0: - axs[0].text( - xmax + 5, layerBounds[layerKeys[n]] / 2, layerKeys[n], color='black', fontsize=fontSize - ) - else: - axs[0].text( - xmax + 5, - (layerBounds[layerKeys[n]] + layerBounds[layerKeys[n - 1]]) / 2, - layerKeys[n], - color='black', - fontsize=fontSize, - verticalalignment='center', - ) - - # set vertical line at stimulus onset - if type(stim_start_time) is int or type(stim_start_time) is float: - axs[0].vlines(stim_start_time, ymin, ymax, colors='red', linewidth=1, linestyles='dashed') - - # save figure - if saveFig: - if isinstance(saveFig, basestring): - filename = saveFig - else: - filename = sim.cfg.filename + '_CSD.png' - try: - plt.savefig(filename, dpi=dpi) - except: - plt.savefig('CSD_fig.png', dpi=dpi) - - # display figure - if showFig: - plt.show() diff --git a/netpyne/analysis/dipole.py b/netpyne/analysis/dipole.py index 88f52d355..165467d4b 100644 --- a/netpyne/analysis/dipole.py +++ b/netpyne/analysis/dipole.py @@ -40,6 +40,10 @@ def plotDipole(showCell=None, showPop=None, timeRange=None, dpi=300, figSize=(6, p = sim.allSimData['dipolePops'][showPop] else: p = sim.allSimData['dipoleSum'] + + # if list (as a result of side-effect of some of save-load operations), make sure to convert to np.array + if isinstance(p, list): + p = np.array(p) p = p / 1000.0 # convert from nA to uA except: diff --git a/netpyne/analysis/lfp.py b/netpyne/analysis/lfp.py index da514d447..56dd2b7e6 100644 --- a/netpyne/analysis/lfp.py +++ b/netpyne/analysis/lfp.py @@ -38,7 +38,6 @@ def prepareLFP( electrodes=['avg', 'all'], pop=None, LFPData=None, - logy=False, normSignal=False, filtFreq=False, filtOrder=3, @@ -55,34 +54,20 @@ def prepareLFP( if not sim: from .. import sim - # create the output data dictionary - data = {} - data['electrodes'] = {} - data['electrodes']['names'] = [] - data['electrodes']['locs'] = [] - data['electrodes']['lfps'] = [] - # set time range if timeRange is None: timeRange = [0, sim.cfg.duration] - # create an array of the time steps - t = np.arange(timeRange[0], timeRange[1], sim.cfg.recordStep) - data['t'] = t - # accept input lfp data if LFPData is not None: # loading LFPData is not yet functional - lfp = LFPData[int(timeRange[0] / sim.cfg.recordStep) : int(timeRange[1] / sim.cfg.recordStep), :] + lfp = LFPData else: if pop and pop in sim.allSimData['LFPPops']: - lfp = np.array(sim.allSimData['LFPPops'][pop])[ - int(timeRange[0] / sim.cfg.recordStep) : int(timeRange[1] / sim.cfg.recordStep), : - ] + lfp = np.array(sim.allSimData['LFPPops'][pop]) else: - lfp = np.array(sim.allSimData['LFP'])[ - int(timeRange[0] / sim.cfg.recordStep) : int(timeRange[1] / sim.cfg.recordStep), : - ] + lfp = np.array(sim.allSimData['LFP']) + lfp = lfp[int(timeRange[0] / sim.cfg.recordStep) : int(timeRange[1] / sim.cfg.recordStep), :] # filter the signals if filtFreq: @@ -115,52 +100,56 @@ def prepareLFP( lfp[:, i] /= max(lfp[:, i]) # electrode preparation + data = prepareDataPerElectrode(lfp, electrodes, timeRange, sim) + + # data['electrodes']['lfps'] = np.transpose(np.array(data['electrodes']['lfps'])) + + return data + +def prepareDataPerElectrode(signalPerElectrode, electrodes, timeRange, sim): + + # create the output data dictionary + data = {} + data['electrodes'] = {} + data['electrodes']['names'] = [] + data['electrodes']['locs'] = [] + data['electrodes']['data'] = [] + + # create an array of the time steps + t = np.arange(timeRange[0], timeRange[1], sim.cfg.recordStep) + data['t'] = t + if 'all' in electrodes: electrodes.remove('all') electrodes.extend(list(range(int(sim.net.recXElectrode.nsites)))) - # if 'avg' in electrodes: - # electrodes.remove('avg') - # data['electrodes']['names'].append('avg') - # data['electrodes']['locs'].append(None) - # data['electrodes']['lfps'].append(np.mean(lfp, axis=1)) - for i, elec in enumerate(electrodes): - if isinstance(elec, Number) and (LFPData is not None or elec <= sim.net.recXElectrode.nsites): - lfpSignal = lfp[:, elec] + loc = None + if isinstance(elec, Number) and (elec <= sim.net.recXElectrode.nsites): + signal = signalPerElectrode[:, elec] loc = sim.cfg.recordLFP[elec] elif elec == 'avg': - lfpSignal = np.mean(lfp, axis=1) - loc = None - elif isinstance(elec, list) and ( - LFPData is not None or all([x <= sim.net.recXElectrode.nsites for x in elec]) - ): - lfpSignal = np.mean(lfp[:, elec], axis=1) - loc = None - - if len(t) < len(lfpSignal): - lfpSignal = lfpSignal[: len(t)] + signal = np.mean(signalPerElectrode, axis=1) + elif isinstance(elec, list) and (all([x <= sim.net.recXElectrode.nsites for x in elec])): + signal = np.mean(signalPerElectrode[:, elec], axis=1) data['electrodes']['names'].append(str(elec)) data['electrodes']['locs'].append(loc) - data['electrodes']['lfps'].append(lfpSignal) - - # data['electrodes']['lfps'] = np.transpose(np.array(data['electrodes']['lfps'])) - + data['electrodes']['data'].append(signal) return data @exception def preparePSD( LFPData=None, + CSD=False, sim=None, timeRange=None, electrodes=['avg', 'all'], pop=None, NFFT=256, noverlap=128, - nperseg=256, minFreq=1, maxFreq=100, stepFreq=1, @@ -181,24 +170,35 @@ def preparePSD( if not sim: from .. import sim - data = prepareLFP( - sim=sim, - timeRange=timeRange, - electrodes=electrodes, - pop=pop, - LFPData=LFPData, - logy=logy, - normSignal=normSignal, - filtFreq=filtFreq, - filtOrder=filtOrder, - detrend=detrend, - **kwargs - ) + if not CSD: + data = prepareLFP( + sim=sim, + timeRange=timeRange, + electrodes=electrodes, + pop=pop, + LFPData=LFPData, + logy=logy, + normSignal=normSignal, + filtFreq=filtFreq, + filtOrder=filtOrder, + detrend=detrend, + **kwargs + ) + else: + data = sim.analysis.prepareCSD( + sim=sim, + timeRange=timeRange, + electrodes=electrodes, + pop=pop, + getAllData=False, + **kwargs + ) + print('Preparing PSD data...') names = data['electrodes']['names'] - lfps = data['electrodes']['lfps'] + lfps = data['electrodes']['data'] allFreqs = [] allSignal = [] @@ -331,7 +331,7 @@ def prepareSpectrogram( if not timeRange: timeRange = [0, sim.cfg.duration] - lfps = np.array(data['electrodes']['lfps']) + lfps = np.array(data['electrodes']['data']) names = data['electrodes']['names'] electrodes = data['electrodes'] diff --git a/netpyne/analysis/mapping.py b/netpyne/analysis/mapping.py index d714dc008..2e9655c0a 100644 --- a/netpyne/analysis/mapping.py +++ b/netpyne/analysis/mapping.py @@ -245,4 +245,5 @@ def plotRaster( saveData=saveData, saveFig=saveFig, showFig=showFig, + syncLines=syncLines ) diff --git a/netpyne/analysis/network.py b/netpyne/analysis/network.py index ab4607893..068451b7f 100644 --- a/netpyne/analysis/network.py +++ b/netpyne/analysis/network.py @@ -51,6 +51,14 @@ def _plotConnCalculateFromSim( removeWeightNorm, logPlot=False, ): + # params validation + if groupBy == 'cell': + supportedFeatures = ['weight', 'delay', 'numConns'] + else: + supportedFeatures = ['weight', 'delay', 'numConns', 'probability', 'strength', 'convergence', 'divergence'] + if feature not in supportedFeatures: + print(f' Unsupported feauture "{feature}". Conn matrix with groupBy="{groupBy}" only supports features in {supportedFeatures}') + return None, None, None from .. import sim @@ -100,12 +108,9 @@ def list_of_dict_unique_by_key(seq, key): # Calculate matrix if grouped by cell if groupBy == 'cell': - if feature in ['weight', 'delay', 'numConns']: - connMatrix = np.zeros((len(cellGidsPre), len(cellGidsPost))) - countMatrix = np.zeros((len(cellGidsPre), len(cellGidsPost))) - else: - print(' Conn matrix with groupBy="cell" only supports features= "weight", "delay" or "numConns"') - return None, None, None + connMatrix = np.zeros((len(cellGidsPre), len(cellGidsPost))) + countMatrix = np.zeros((len(cellGidsPre), len(cellGidsPost))) + cellIndsPre = {cell['gid']: ind for ind, cell in enumerate(cellsPre)} cellIndsPost = {cell['gid']: ind for ind, cell in enumerate(cellsPost)} @@ -239,6 +244,8 @@ def list_of_dict_unique_by_key(seq, key): if feature == 'divergence': maxPreConnMatrix[popIndsPre[prePop], popIndsPost[postPop]] = numCellsPopPre[prePop] + preCellPops = {cell.gid: cell['tags']['pop'] for cell in cellsPre} + # Calculate conn matrix for cell in cellsPost: # for each postsyn cell @@ -254,8 +261,7 @@ def list_of_dict_unique_by_key(seq, key): if conn[preGidIndex] == 'NetStim': prePopLabel = conn[preLabelIndex] if preLabelIndex in conn else 'NetStim' else: - preCell = next((cell for cell in cellsPre if cell['gid'] == conn[preGidIndex]), None) - prePopLabel = preCell['tags']['pop'] if preCell else None + prePopLabel = preCellPops.get(conn[preGidIndex]) if prePopLabel in popIndsPre: if feature in ['weight', 'strength']: diff --git a/netpyne/analysis/spikes.py b/netpyne/analysis/spikes.py index bfc510199..39bb9c403 100644 --- a/netpyne/analysis/spikes.py +++ b/netpyne/analysis/spikes.py @@ -49,7 +49,7 @@ def prepareSpikeData( fileDesc=None, fileType=None, fileDir=None, - calculatePhase=False, + colorbyPhase=None, **kwargs ): """ @@ -157,6 +157,105 @@ def prepareSpikeData( sel = pd.concat([sel, ns]) numNetStims += len(spkindsNew) + # Calculating Phase of spikes + if isinstance(colorbyPhase,dict): + + try: + Signal = colorbyPhase['signal'] + except: + print("Importing signal for coloring spikes - No information about the signal") + Signal = 'LFP' + + if isinstance(Signal, basestring) or isinstance(Signal,np.ndarray): + from scipy import signal, stats, interpolate + from scipy.signal import hilbert + from math import fmod, pi + + rawSignal = [] + if isinstance(Signal, basestring): + # signal provided as a list packed in a pkl + if Signal.endswith('.pkl'): + import pickle + + with open(Signal, 'rb') as input_file: + rawSignal = pickle.load(input_file) + + try: + fs = colorbyPhase['fs'] + except: + print("Importing signal for coloring spikes - No frequency sampling provided") + fs = 1000.0 # assumes data sampled in ms + + time = np.linspace(0,len(rawSignal)*1000/fs,len(rawSignal)+1) # in milliseconds + + elif Signal == 'LFP': + try: + electrode = colorbyPhase['electrode'] + if electrode > sim.net.recXElectrode.nsites: + print('Wrong electrode number for coloring spikes according to LFP phase - Assigning first element in LFP recording setup') + electrode = 1 + except: + electrode = 1 + + rawSignal = [sim.allSimData['LFP'][n][electrode-1] for n in range(len(sim.allSimData['LFP']))] + fs = 1000.0/sim.cfg.recordStep + time = np.linspace(0,len(rawSignal)*1000/fs,len(rawSignal)+1) # in milliseconds + + else: + print('No signal recovered to color spikes according to its phase') + + # it is an array + else: + rawSignal = Signal.tolist() + try: + fs = colorbyPhase['fs'] + except: + print("Importing signal for coloring spikes - No frequency sampling provided") + fs = 1000.0 # assumes data sampled in ms + + time = np.linspace(0,len(rawSignal)*1000/fs,len(rawSignal)+1) # in milliseconds + + # Processing rawSignal + rawSignal.append(2*rawSignal[-1]-rawSignal[-2]) # extrapolation for the last element + rawSignal = stats.zscore(np.float32(rawSignal)) + rawSignal_ = np.r_[rawSignal[-1::-1],rawSignal,rawSignal[-1::-1]] # Reflect signal to minimize edge artifacts + + nyquist = fs/2.0 + + # parameters provided to filter the signal - setting defaults otherwise + try: + filtOrder = colorbyPhase['filtOrder'] + except: + filtOrder = 3 + + try: + filtFreq = colorbyPhase['filtFreq'] + except: + filtFreq = [1,500] + + if isinstance(filtFreq, list): # bandpass + Wn = [filtFreq[0]/nyquist, filtFreq[1]/nyquist] + b, a = signal.butter(filtOrder, Wn, btype='bandpass') + + elif isinstance(filtFreq, Number): # lowpass + Wn = filtFreq/nyquist + b, a = signal.butter(filtOrder, Wn) + + rawSignalFiltered_ = signal.filtfilt(b, a, rawSignal_) + + analytic_signal = hilbert(rawSignalFiltered_)[len(rawSignal):-len(rawSignal)] + amplitude_envelope = np.abs(analytic_signal) + instantaneous_phase = np.unwrap(np.angle(analytic_signal)) + instantaneous_phase_mod = [(fmod(instantaneous_phase[nn]+pi,2*pi)-pi)*(180/pi) for nn in range(len(instantaneous_phase))] + instantaneous_phase = np.r_[instantaneous_phase_mod] + + f_Rhythm = interpolate.interp1d(time,instantaneous_phase) + + sel['spkPhase'] = sel['spkt'].apply(f_Rhythm) + + else: + print('No signal recovered to color spikes according to its phase') + if len(cellGids) > 0 and numNetStims: ylabelText = ylabelText + ' and NetStims (at the end)' elif numNetStims: @@ -277,6 +376,12 @@ def prepareSpikeData( 'axisArgs': axisArgs, 'legendLabels': legendLabels, } + + if colorbyPhase: + spikeData.update({'spkPhases': sel['spkPhase'].tolist()}) + if 'include_signal' in colorbyPhase and colorbyPhase['include_signal']==True: + spikeData.update({'signal': analytic_signal}) + spikeData.update({'time': time}) if saveData: saveFigData(spikeData, fileName=fileName, fileDesc='spike_data', fileType=fileType, fileDir=fileDir, sim=sim) @@ -292,6 +397,7 @@ def prepareRaster( maxSpikes=1e8, orderBy='gid', popRates=True, + colorbyPhase=None, saveData=False, fileName=None, fileDesc=None, @@ -311,6 +417,7 @@ def prepareRaster( maxSpikes=maxSpikes, orderBy=orderBy, popRates=popRates, + colorbyPhase=colorbyPhase, saveData=saveData, fileName=fileName, fileDesc=fileDesc if fileDesc else 'raster_data', diff --git a/netpyne/analysis/spikes_legacy.py b/netpyne/analysis/spikes_legacy.py index 335ebd921..66e92bc2f 100755 --- a/netpyne/analysis/spikes_legacy.py +++ b/netpyne/analysis/spikes_legacy.py @@ -447,7 +447,7 @@ def plotRaster( lw=2, marker='|', markerSize=5, - popColors=None, + popColors={}, figSize=(10, 8), fontSize=12, dpi=100, @@ -1781,6 +1781,7 @@ def lognorm(meaninput, stdinput, binedges, n, popLabel, color): # ------------------------------------------------------------------------------------------------------------------- ## Plot spiking power spectral density # ------------------------------------------------------------------------------------------------------------------- +# @exception def plotRatePSD( include=['eachPop', 'allCells'], @@ -1795,7 +1796,8 @@ def plotRatePSD( smooth=0, norm=False, overlay=True, - popColors=None, + popColors={}, + yLogScale = True, ylim=None, figSize=(10, 8), fontSize=12, @@ -1884,6 +1886,11 @@ def plotRatePSD( **Default:** ``None`` uses standard colors **Options:** ``