as systems neuroscience moves increasingly towards highly parallel physiological recording techniques , generation , management , and analysis of large complex data sets is becoming the norm .
the goal is to understand how neurons in visual cortex respond to visual stimuli , to the extent that the responses to arbitrary stimuli can be predicted .
neurons in close proximity are more likely to show functionally interesting interactions , and insights into how such localized populations work may help guide understanding of other parts of cortex , or even the brain as a whole . to this end
we need to record and analyse the simultaneous spiking behaviour of many neurons in response to a wide variety of visual stimuli .
we use 54-channel silicon polytrodes , in both rat and cat primary visual cortex , to extracellularly sample spiking activity constrained to roughly a cortical column ( figure 1a ) ( blanche et al . , 2005 ) .
time - locked visual stimuli are presented to the animal while simultaneously recording from dozens of neurons ( figure 1b ) .
waveforms are recorded continuously at a rate of 2.7  mb / s for up to 90  min ( 15  gb ) at a time .
a single animal experiment can last up to 3 days and generate hundreds of gb of data .
setting up our electrophysiology rig , with custom acquisition software written in delphi ( blanche , 2005 ) , was the first step .
although we had existing solutions in place for visual stimulation , waveform visualization and spike sorting , and spike train analysis , all three had limitations which were addressed by rewriting our software in python .
recording sites are closely spaced , such that a spike will typically appear on several sites at the same time ( see figure 3 ) .
stimuli are presented to the animal while stimulus information and extracellular voltage waveforms are acquired and saved to disk .
we chose the vision egg partly because of the language it was written in and written for : python .
we were thus introduced to python via one of its many packages , and the experience was so positive that it encouraged us to standardize on python for the spike sorting and spike train analysis projects to follow .
for one of us ( m. spacek ) , the switch to python has made programming a much more enjoyable and productive experience , and has resulted in greatly improved programming skills .
the benefits of python have been extolled at length elsewhere ( hetland , 2005 ; langtangen , 2008 ; lutz , 2006 ) .
briefly , python is a powerful , dynamically typed , interpreted language that  fits your brain  , with syntax akin to  executable pseudocode  .
some of its clarity stems from a philosophy to provide  one  and preferably only one  obvious way  to do a given task ( peters , 2004 ) , making features easy to remember .
its clarity is also due to a strong adherence to object - oriented programming principles [ chapter 7 of hetland ( 2005 ) is an excellent introduction ] . in python ,
this means that everything has attributes and methods ( methods are functions that are bound to and act on objects ) , and can thus be treated in a similar way .
a class can inherit attributes and methods from other classes hierarchically , allowing for substantial code reuse , and therefore less code to maintain .
python code is succinct compared to most other languages : a lot can be accomplished in only a few lines .
finally , python is free and open source , and encourages open source software development .
this is partly due to its interpreted nature : the source code and executable are typically one and the same .
python has a stable and feature - rich numeric library called numpy1 which provides an n - dimensional array object .
numpy arrays can be subjected to vectorized operations , most of which call static c functions , allowing them to run almost as fast as pure c code . yet
numpy turns python into an effective replacement for matlab ( the mathworks , natick , ma , usa ) , and is used extensively by dimstim , spyke , and neuropy .
while all three projects presented here were written in python , their use and implementation are very different .
spyke has a graphical user interface ( gui ) and looks like a native application , while neuropy is typically accessed from the python command line as a library . here
, we explore some of the features and benefits of python and its many add - on packages for the electrophysiologist , by introducing our own three packages as detailed working examples .
in our experiments , we needed a way to display and control a wide variety of stimuli with many different parameters , often shuffled with respect to each other in various ways .
since spike times are acquired at sub - millisecond temporal resolution , and since precise spike timing may play a role in neural coding ( mainen and sejnowski , 1995 ; vanrullen and thorpe , 2002 ) , we also wanted high temporal precision in the stimulus .
our prior stimulus software was written in fortran and ran under dos with a 32-bit extender .
the last graphics cards to support it were limited in the size and speed of movie frames they could draw to screen . moreover , these cards were limited to a screen refresh rate of 100  hz at our desired resolution .
we found significant artefactual phase - locking of responses in visual cortex at this frequency ( blanche , 2005 ) , which has been a concern reported elsewhere ( williams et al .
dimstim displays full - screen stimuli at a refresh rate of 200  hz , providing precise control of the display at 5  ms intervals without frame drops .
stimuli include manually controlled , drifting , and flashed bars and gratings , sparse noise , and m - sequence noise ( golomb , 1967 ) and natural scene movies .
stimulus parameters can be shuffled with or without replacement , independently or in covariation with each other .
parameters include spatial location and phase , orientation , speed , duration , size , mask , contrast , brightness , and spatial and temporal frequencies .
a copy of the script , and an index of the contents of the screen on each screen refresh , are sent to the acquisition computer , for simultaneous recording of stimulus and neuronal responses .
dimstim relies heavily on the vision egg2 library ( straw , 2008 ) to generate stimuli .
the vision egg uses the well - established opengl3 graphics language , which thanks to the demands of video games , is now supported by all modern video cards on all major platforms .
stimuli are displayed on a 19 '' iiyama hm903dtb and a 22 '' hm204dta crt monitor , two of only a handful of consumer monitors that are capable of 800    600 resolution at 200  hz .
unfortunately , like most other crts , these particular models have now been discontinued , but used ones may still be available .
hopefully the timing of lcd monitors will improve such that they can replace crts for temporally
multitasking operating systems ( oses ) present a challenge for real - time control of the screen . often , the os will decide to delay an operation to maintain responsiveness in other areas .
this can lead to frame drops , but can be mitigated by increasing the priority of the python process . setting the process and thread priorities to their maximum levels in the vision egg completely eliminated frame drops in windows xp , but with the unfortunate loss of mouse and keyboard polling . in dimstim , this meant that the user had no way of interrupting the stimulus script , other than by resetting the computer . moving to a computer with
a dual core cpu alleviated this problem , as the maximum priority python process was delegated to one core without interruption , while other os tasks such as keyboard polling ran normally on the second core .
dimstim communicates stimulus parameters on a frame - by - frame basis to the acquisition computer via a pci digital output board ( dt340 , data translations , marlboro , ma , usa ) , for simultaneous recording of stimulus timing alongside neuronal responses .
parameters are described by sending the row index of a large lookup table (  sweep table  ) on every screen refresh .
the sweep table contains all the combinations of the dynamic parameters , i.e. those stimulus parameters that can vary from one screen refresh to the next .
the digital output board is controlled by its driver 's c library . because python is written in c ( other implementations also exist ) , it has a c application programming interface ( api ) , and extensions to python can be written in c. we wrote such an extension to interact with the board 's c library , but today this is no longer necessary .
a new built - in python module called  ctypes  now allows interaction with a c library on any platform directly from within python code .
this is much simpler , as it removes the need to both write and compile c extension code using python 's somewhat tedious c api .
dimstim includes a demo ( olda_demo.py ) of how to use ctypes to directly interact with data translations ' open layers data acquisition library .
libraries for cards from other vendors ( such as national instruments ' ni - daqmx ) can be similarly accessed .
the photodiode signal , along with the raster signal from the video card and the digital outputs from the stimulus computer , were all recorded simultaneously .
we discovered that the contents of the screen always lagged by one screen refresh , due to opengl 's buffer swapping behaviour ( straw , 2008 ) .
this was corrected for by adding one frame time ( 5  ms ) to the timestamp of the digitized raster signal in the acquisition system .
several levels of uncorrected luminance were measured with a light meter ( minolta ls-100 ) and fit to a power law expression to determine the exponent corresponding to the gamma value of the screen ( blanche , 2005 ; straw , 2008 ) .
gamma correction can be set independently for each script , or globally across all scripts in dimstim 's config file .
natural scene movies used by dimstim were filmed outdoors with an ordinary compact digital camera ( canon powershot sd200 ) with 320    240 resolution at 60 frames per second ( fps ) .
unfortunately , this camera could record no more than 1  min of video at a time . to generate longer movies ,
multiple clips were filmed in succession , while keeping the camera as motionless as possible between the end of one clip and the start of the next .
concatenation of and conversion from multiple colour .avi files to a single uncompressed greyscale movie file was done using david mcnab 's y4m4 package .
processed movies were displayed in dimstim with the same visual angle subtended by the camera , at 67 fps ( three 5  ms screen refreshes per movie frame ) .
dimstim 's config file stores default values for a variety of generic parameters that apply to most stimuli .
these parameters include spatial location , size , orientation offset , and temporal and spatial frequencies . for simplicity ,
the config file can be edited by hand , but the typical procedure when optimizing parameters for the current neural population is to run a manually controlled bar or grating stimulus . for user convenience , the stimulus is shown simultaneously on two displays driven by two video outputs from the graphics card : one for the animal , and one for the user .
the parameters of the manual stimulus are controlled in real - time with the mouse and keyboard .
first , objects for storage of static and dynamic parameters are instantiated (  s  and  d  respectively , lines 56 ) . to these
are bound various different parameters as attributes ( denoted by a  .  ) . in this example
, most values are declared directly by the script , but two static parameters , grating orientation offset and gamma correction , are retrieved from their defaults in the config file , using the dimstim config parser object named  dc 
dynamic parameters , if assigned a list of multiple values , will iterate over those values over the course of the experiment . in this case ,
grating orientation , spatial frequency , and temporal frequency are all assigned multiple values ( lines 28 , 36 , 38 ) .
the rest remain constant for the duration of the experiment . in order to describe their interdependence and shuffling ,
each multiple - value dynamic parameter must be declared as a  variable  ( lines 5355 ) .
variables with the same dimension value (  dim  keyword argument ) covary with each other , and must therefore all have the same number of values and the same shuffle flag .
variables with different dimension values vary independently in a combinatorial fashion , with the lowest numbered dimension varying slowest , and the highest varying fastest .
this is implemented by dynamically generating a string object containing python code with the correct number of nested for loops ( equalling the number of independent variables specified in the script ) , and then executing the contents of the string with python 's exec ( ) function ( see the dimstim.core.sweeptable class ) .
next , the number of times to cycle through all combinations , and the frequency at which to insert a blank screen sweep ( for determining baseline firing rates ) are specified in their own objects ( lines 5758 ) .
finally , all these objects are passed together to the grating class ( which like all other dimstim stimuli , inherits from the experiment class ) to instantiate a grating experiment object , and the experiment is run ( lines 6265 ) . with 12 orientations , 6 spatial frequencies , and 4 temporal frequencies ,
each is presented four times for a total of 1152 stimulus sweeps , lasting 4  s each , for a total experiment time of about 77  min ( not including blank sweeps ) . a dimstim script describing a drifting sinusoidal grating .
such scripts may be edited at will , and are the primary way the user interacts with dimstim .
after some error checking , the script executes from the system 's command line , to which status messages are printed .
comments , denoted by # and " " " in python , are highlighted in red .
various checks are done to alert for any obvious errors in the user edited script .
this makes it possible to later reconstruct the sweep table for analysis , and even replay the entire experiment exactly , without the need for access to the original script on the stimulus computer . to ensure accurate timing , stimuli run only on the animal display , while the user display shows the system command line . in between experiments when no stimuli are running , a blank grey desktop is shown on the animal display .
dimstim 's config file stores default values for a variety of generic parameters that apply to most stimuli .
these parameters include spatial location , size , orientation offset , and temporal and spatial frequencies . for simplicity ,
the config file can be edited by hand , but the typical procedure when optimizing parameters for the current neural population is to run a manually controlled bar or grating stimulus . for user convenience ,
the stimulus is shown simultaneously on two displays driven by two video outputs from the graphics card : one for the animal , and one for the user .
the parameters of the manual stimulus are controlled in real - time with the mouse and keyboard . once the user is satisfied , the parameters are saved to the config file .
first , objects for storage of static and dynamic parameters are instantiated (  s  and  d  respectively , lines 56 ) . to these
are bound various different parameters as attributes ( denoted by a  .  ) . in this example
, most values are declared directly by the script , but two static parameters , grating orientation offset and gamma correction , are retrieved from their defaults in the config file , using the dimstim config parser object named  dc  ( lines 15 and 23 ) .
dynamic parameters , if assigned a list of multiple values , will iterate over those values over the course of the experiment . in this case , grating orientation , spatial frequency , and temporal frequency are all assigned multiple values ( lines 28 , 36 , 38 ) .
the rest remain constant for the duration of the experiment . in order to describe their interdependence and shuffling ,
each multiple - value dynamic parameter must be declared as a  variable  ( lines 5355 ) .
variables with the same dimension value (  dim  keyword argument ) covary with each other , and must therefore all have the same number of values and the same shuffle flag .
variables with different dimension values vary independently in a combinatorial fashion , with the lowest numbered dimension varying slowest , and the highest varying fastest .
this is implemented by dynamically generating a string object containing python code with the correct number of nested for loops ( equalling the number of independent variables specified in the script ) , and then executing the contents of the string with python 's exec ( ) function ( see the dimstim.core.sweeptable class ) .
next , the number of times to cycle through all combinations , and the frequency at which to insert a blank screen sweep ( for determining baseline firing rates ) are specified in their own objects ( lines 5758 ) .
finally , all these objects are passed together to the grating class ( which like all other dimstim stimuli , inherits from the experiment class ) to instantiate a grating experiment object , and the experiment is run ( lines 6265 ) . with 12 orientations , 6 spatial frequencies , and 4 temporal frequencies ,
each is presented four times for a total of 1152 stimulus sweeps , lasting 4  s each , for a total experiment time of about 77  min ( not including blank sweeps ) . a dimstim script describing a drifting sinusoidal grating .
such scripts may be edited at will , and are the primary way the user interacts with dimstim .
after some error checking , the script executes from the system 's command line , to which status messages are printed .
comments , denoted by # and " " " in python , are highlighted in red .
various checks are done to alert for any obvious errors in the user edited script .
this makes it possible to later reconstruct the sweep table for analysis , and even replay the entire experiment exactly , without the need for access to the original script on the stimulus computer . to ensure accurate timing , stimuli run only on the animal display , while the user display shows the system command line . in between experiments when no stimuli are running , a blank grey desktop is shown on the animal display .
once neural waveform and stimulus data were saved to disk by our acquisition system ( written in delphi ) , we needed a way to retrieve the data for visualization and spike sorting .
our existing program for this , also written in delphi , had some bugs and missing features . however
, the delphi environment required a license , the program would only run in windows , and the code was more procedural than object - oriented .
in particular , some of the code had blocks ( if statements , for / while statements ) that were nested many layers deep , making it difficult to follow . 
flat is better than nested  ( peters , 2004 ) is another python philosophy .
several short , shallow blocks of code are easier to understand and manage than one long deep block .
spyke has a cross - platform gui with native widgets for data visualization and navigation , and spike sorting ( figure 3 ) . spike waveforms are displayed in two ways : spatially according to the polytrode channel layout ( spike window ) , and vertically in chart form ( chart window ) .
local field potential ( lfp ) waveforms are also displayed vertically in chart form ( lfp window ) .
polytrode channels are closely spaced ( 4375 m ) over two or three columns ( figure 1a ) . a single spike
can generate a signal on multiple channels , hence the need to visualize waveforms according to their polytrode channel layout .
channels are colour - coded to make them easy to distinguish and align across windows .
spyke looks and behaves like a native gui application , with menus , buttons , and resizable windows .
navigation is mouse and keyboard based . a horizontal slider and combo box at the top of the main spyke window control file position in time .
left and right arrow keys , and page up and page down keys step through the data with single timepoint or 1 ms resolution respectively .
clicking on any data window ( spike , chart , or lfp ) centres all three windows on that timepoint .
holding ctrl and scrolling the mouse wheel over a data window zooms it in or out in time .
holding ctrl and clicking on a channel enables or disables it . hovering the mouse over a data window
displays a tooltip with the timestamp , channel , and voltage currently under the mouse cursor .
main spyke window ( top ) , with data windows ( bottom ) showing high - pass waveforms in polytrode layout ( left ) and chart layout ( middle ) .
a third data window shows the low - pass lfp waveforms ( right ) concurrently recorded from a subset of channels ( colour coded ) .
the shaded region in the middle of both the chart and lfp windows represents the time range spanned by the window to its left .
this is a python interface to the wxwidgets c++ gui library which generates widgets on windows , linux , and osx .
now well over a decade old ( rappin and dunn , 2006 ) , wxpython is a stable library that has adapted to changing oses .
widgets include everything from windows , menus , and buttons , to more complex list and tree controls .
wxpython has a big advantage over other gui libraries in its use of widgets that are native to the os the program is running on , such that they look and behave identically to normally created widgets in that os .
itself a wxpython based gui application , wxglade takes the programmer 's visual layout and automatically generates the corresponding layout code in python .
this code can then be included in the programmer 's own code base , typically by defining a class that inherits from the automatically generated code .
although wxglade is not necessary for writing a gui with wxpython , we found it much faster and easier than writing all of the layout code by hand .
unfortunately , some widgets are inherently different on different oses . writing and testing a wxpython gui on only one os will therefore not guarantee perfect functionality on another . to do so would require checking for the current os , and implementing certain things differently depending on the os .
spyke does not currently do this , and has so far only been thoroughly tested in windows .
although wxpython is one of the best ( rappin and dunn , 2006 ) , it has bugs7  some of them longstanding  that had to be worked around in spyke .
it has two interfaces : one that mimics the familiar plotting commands of matlab , and another that is much more object - oriented .
scaling of plots is handled automatically by matplotlib , such that when the wxpython window is resized by dragging its corner or edge , the plotted traces inside resize accordingly .
another benefit of matplotlib is its antialiasing abilities , providing beautiful output with subpixel resolution .
there is some performance penalty for using such a high level drawing library , but performance is fast enough on fairly ordinary hardware ( pentium m 1.6 ghz notebook ) , even when scrolling through 54 channels of data with thousands of data points on screen at a time . more importantly
the data acquisition files are complex , with different types of data multiplexed throughout the file . on opening
, the file must be parsed to determine the number and offset values of hundreds of thousands of records in the file .
for multi gb files , this can take up to a few minutes . to deal with this
this is done using python 's pickle module , which can take a snapshot of almost any python object in memory , serialize it , and save it to disk as a  pickle  .
a pickle can then later be restored ( unpickled ) to memory as a live python object , even on a different platform . in this case , a custom written file object containing all of the parse information is saved to disk as a .parse file of only a few mb in size . restoring from the .parse file is about an order of magnitude faster than reparsing the entire acquisition file .
segments of waveform data are loaded from the acquisition file , nyquist interpolated , and sample - and - hold delay ( shd ) corrected on the fly as needed ( blanche and swindale , 2006 ) .
interpolation is performed to improve spike detection , and nyquist interpolation is the optimum method of reconstructing a bandwidth - limited signal at arbitrary resolution .
to do so , a set of sinc function kernels is generated ( one kernel per interpolated data point , each kernel with a different phase offset ) and convolved with the data . for shd correction , a different set of kernels is generated for each channel . correcting for each channel 's shd requires appropriate modification of the phase offset of each kernel for that channel .
for example , interpolating from 25 to 50  khz with shd correction requires two appropriately phase corrected kernels per channel . each kernel is separately convolved with the data ( using numpy.convolve ( ) ) , and the resulting data points are interleaved to return the final interpolated waveform .
the  dynamic multiphasic  method searches for two consecutive threshold crossings of opposite polarity within a defined period of time . the second crossing 's threshold is dynamically set according to the amplitude of the first phase of the spike . for both methods , primary thresholds are calculated separately for each channel , based on the standard deviation or median noise level of either the entire recording or of a narrow sliding window thereof .
spatiotemporal detection lockouts prevent double triggering off of the same spike , while minimizing the chance of missed spikes .
some algorithms , such as these event detection methods , can not be easily vectorized and require a custom loop .
due to its dynamic typing and interpreted nature , long loops are slow to execute in python . for the majority of software development ,
developer time is usually much more valuable than cpu time ( hetland , 2005 ) , but numerically intensive software is the exception . writing fast python extensions in c has always been possible , but the c interface code required by python 's api is tedious to write , and writing in c eliminates the convenience of working in python syntax . to get around this ,
the cython9 package ( a fork of the pyrex package ) specifies a sublanguage almost identical to python , with some extra keywords to declare loop variables as static c types . after issuing the standard python setup.py build command ,
such code is automatically translated into an intermediary c file including all of the tedious interface code .
this is subsequently compiled into object code and is accessible as a standard c extension module from within python , just as a handwritten c extension would be .
this yields the computational speed of c loops when needed , with the developmental speed , convenience and familiarity of python syntax to implement them .
cython was used to write the custom loop that iterates over timepoints and channels for each of the event detection methods . for 25  khz sampled waveform data on 54 channels ,
this amounts to 1.35 million iterations per second of data . on an average single - core notebook computer ( pentium m 1.6  ghz )
the data is partitioned into blocks ( typically 1  s long ) , and each is searched independently , allowing multiple core cpus to be exploited .
( gil ) in the c implementation of python , multiple processes must typically be used instead of multiple threads to take advantage of multiple cores .
unfortunately , a process can require significantly more memory and more time to create than a thread .
there are ways around the gil , but the best solution for spyke is not yet clear .
search options are controlled in the  detect  tab in spyke 's main window ( figure 3 ) .
searches can be limited to specific time ranges in the file , in the number of events detected , and whether to search linearly or randomly .
random sampling is important to build up a temporally unbiased collection of detected events with which to build templates .
searching for the next or previous spike relative to the current timestamp can be done quickly using the keyboard .
searches are restricted to enabled channels , allowing for a targeted increase in the number of events belonging to a spatially localized template .
when a search completes , the sort window ( figure 4a ) opens and is populated with any newly detected events .
the user then visually sorts the detected events ( typically only a fraction of all spikes in the recording ) into templates corresponding to isolated neurons .
any number of event or template mean waveforms can be overplotted with each other . although the mouse may be used , keyboard commands are more efficient for toggling the display of events and templates , and moving events and keyboard focus around between the sorted template tree ( left column ) and unsorted event list ( right column ) .
the event list has sortable columns for event i d , maximum channel , timestamp , and match error .
all the events in the list can be matched against the currently selected template , and those match errors populate the error column .
sorting the event list by maximum channel or match error makes manual template generation much easier , because it clusters similar events close to each other in the unsorted event list .
templates and their member spikes are represented in the tree ( left ) , and unsorted detected events in the list ( middle ) . selecting a template or event in either the tree or the list plots its waveform ( right ) .
the tree currently has keyboard focus , making its selections more distinctly coloured than those of the list .
unsorted events have colour coded channels , while each template ( and its member spikes ) has a single identifying colour . here
, template 0 ( red ) , a putative neuron near the top of the polytrode , has 6 member spikes , and its mean waveform is being overplotted with an unsorted event ( # 1260 , multicoloured ) , which fits quite well .
template 1 ( orange ) and all of its member spikes are plotted near the middle of the polytrode . also plotted further down is another unsorted event ( # 1150 , multicoloured ) , which obviously does not fit either template .
( b ) the integrated pyshell window exposes all of spyke 's objects and functionality at the python command line .
template 0 's dictionary ( a mapping from names to values ) of its 6 member events is referenced and returned on lines 12 .
once templates have been generated , a full event detection is run across the whole recording , and the templates are matched against each detected event . or , each template can be slid across the recording and matched against every timepoint in the recording ( blanche et al . , 2005 ) .
either way , matching to target and non - target spikes or noise generally yields a non - overlapping bimodal error distribution . for each template
, a threshold is manually set at the trough between the two peaks in the distribution , and events whose
match errors fall below this threshold are classified as spikes of that template . at any point in the sorting process , the entire  sort  session object , which among other information includes detected events , generated templates , and sorted spikes , can be saved to disk as a .sort file , again using python 's pickle module .
sort sessions can then be restored from disk and sorting can resume in spyke , or their sorted spike times can be used for spike train analysis ( see neuropy section ) .
this increases the file size , but allows for review of detected and sorted spikes without the need to access the original multi gb continuous data acquisition file .
integrated into spyke is patrick o'brien 's pyshell ( figure 4b ) , an enhanced python command line that is part of the wxpython package .
this was , and continues to be , a very useful tool for testing existing features and for developing new ones .
neuropy ( or almost any other python package ) can be imported and used directly from this command line .
for example , spike sorting validation is not yet implemented in spyke 's gui , but all of neuropy 's functionality including autocorrelograms ( to check refractory periods ) can be accessed by typing import neuropy in spyke 's pyshell .
the  dynamic multiphasic  method searches for two consecutive threshold crossings of opposite polarity within a defined period of time . the second crossing 's threshold is dynamically set according to the amplitude of the first phase of the spike . for both methods , primary thresholds are calculated separately for each channel , based on the standard deviation or median noise level of either the entire recording or of a narrow sliding window thereof .
spatiotemporal detection lockouts prevent double triggering off of the same spike , while minimizing the chance of missed spikes .
some algorithms , such as these event detection methods , can not be easily vectorized and require a custom loop .
due to its dynamic typing and interpreted nature , long loops are slow to execute in python . for the majority of software development ,
developer time is usually much more valuable than cpu time ( hetland , 2005 ) , but numerically intensive software is the exception . writing fast python extensions in c has always been possible , but the c interface code required by python 's api is tedious to write , and writing in c eliminates the convenience of working in python syntax . to get around this ,
the cython9 package ( a fork of the pyrex package ) specifies a sublanguage almost identical to python , with some extra keywords to declare loop variables as static c types . after issuing the standard python setup.py build command , such code
is automatically translated into an intermediary c file including all of the tedious interface code .
this is subsequently compiled into object code and is accessible as a standard c extension module from within python , just as a handwritten c extension would be .
this yields the computational speed of c loops when needed , with the developmental speed , convenience and familiarity of python syntax to implement them .
cython was used to write the custom loop that iterates over timepoints and channels for each of the event detection methods . for 25  khz sampled waveform data on 54 channels ,
this amounts to 1.35 million iterations per second of data . on an average single - core notebook computer ( pentium m 1.6  ghz )
the data is partitioned into blocks ( typically 1  s long ) , and each is searched independently , allowing multiple core cpus to be exploited .
( gil ) in the c implementation of python , multiple processes must typically be used instead of multiple threads to take advantage of multiple cores .
unfortunately , a process can require significantly more memory and more time to create than a thread .
there are ways around the gil , but the best solution for spyke is not yet clear .
search options are controlled in the  detect  tab in spyke 's main window ( figure 3 ) .
searches can be limited to specific time ranges in the file , in the number of events detected , and whether to search linearly or randomly .
random sampling is important to build up a temporally unbiased collection of detected events with which to build templates .
searching for the next or previous spike relative to the current timestamp can be done quickly using the keyboard .
searches are restricted to enabled channels , allowing for a targeted increase in the number of events belonging to a spatially localized template .
when a search completes , the sort window ( figure 4a ) opens and is populated with any newly detected events .
the user then visually sorts the detected events ( typically only a fraction of all spikes in the recording ) into templates corresponding to isolated neurons .
although the mouse may be used , keyboard commands are more efficient for toggling the display of events and templates , and moving events and keyboard focus around between the sorted template tree ( left column ) and unsorted event list ( right column ) .
the event list has sortable columns for event i d , maximum channel , timestamp , and match error .
all the events in the list can be matched against the currently selected template , and those match errors populate the error column . sorting the event list by maximum channel or match error makes manual template generation much easier , because it clusters similar events close to each other in the unsorted event list .
templates and their member spikes are represented in the tree ( left ) , and unsorted detected events in the list ( middle ) . selecting a template or event in either the tree or the list plots its waveform ( right ) .
the tree currently has keyboard focus , making its selections more distinctly coloured than those of the list .
unsorted events have colour coded channels , while each template ( and its member spikes ) has a single identifying colour . here
, template 0 ( red ) , a putative neuron near the top of the polytrode , has 6 member spikes , and its mean waveform is being overplotted with an unsorted event ( # 1260 , multicoloured ) , which fits quite well .
template 1 ( orange ) and all of its member spikes are plotted near the middle of the polytrode .
also plotted further down is another unsorted event ( # 1150 , multicoloured ) , which obviously does not fit either template .
( b ) the integrated pyshell window exposes all of spyke 's objects and functionality at the python command line .
template 0 's dictionary ( a mapping from names to values ) of its 6 member events is referenced and returned on lines 12 .
once templates have been generated , a full event detection is run across the whole recording , and the templates are matched against each detected event . or , each template can be slid across the recording and matched against every timepoint in the recording ( blanche et al . , 2005 ) .
either way , matching to target and non - target spikes or noise generally yields a non - overlapping bimodal error distribution . for each template
, a threshold is manually set at the trough between the two peaks in the distribution , and events whose match errors fall below this threshold are classified as spikes of that template . at any point in the sorting process , the entire  sort  session object , which among other information includes detected events , generated templates , and sorted spikes , can be saved to disk as a .sort file , again using python 's pickle module .
sort sessions can then be restored from disk and sorting can resume in spyke , or their sorted spike times can be used for spike train analysis ( see neuropy section ) .
this increases the file size , but allows for review of detected and sorted spikes without the need to access the original multi gb continuous data acquisition file .
integrated into spyke is patrick o'brien 's pyshell ( figure 4b ) , an enhanced python command line that is part of the wxpython package .
this was , and continues to be , a very useful tool for testing existing features and for developing new ones .
neuropy ( or almost any other python package ) can be imported and used directly from this command line .
for example , spike sorting validation is not yet implemented in spyke 's gui , but all of neuropy 's functionality including autocorrelograms ( to check refractory periods ) can be accessed by typing import neuropy in spyke 's pyshell .
after spike sorting , we needed a way to analyse spike trains and their relation to stimuli .
our initial decision was to use matlab for spike train analysis , and we soon developed a collection of matlab scripts for the job , with one function per .m file .
for example , one .m file would load each neuron 's data from disk and return all of them in a cell array of structures .
furthermore , the code became difficult to manage as each additional function required an additional .m file .
we were also faced with out of memory errors , limited gui capabilities , and a high licensing cost .
although matlab 's toolboxes are a major benefit , scipy10 ( jones et al . , 2001 ) , an extensive python library of scientific routines , provides most of the equivalent functionality .
much of scipy is a wrapper for decades - old , highly tested and optimized fortran code .
another package , mlabwrap11 , allows a licensed matlab user to access all of matlab 's functionality , including all of its toolboxes , directly from within python .
although in the end we did not need to use mlabwrap , its existence erased any remaining hesitations about switching to python for analysis .
a data - centric object hierarchy ( figure 5a ) quickly emerged as a natural way to organize neuropy .
each object in the hierarchy has an attribute that references its parent object , as well as all of its child objects .
specifically ,  data  is an abstract object from which all  animals  are accessible .
each animal has polytrode  tracks  , each track has  recordings  , and each recording has both  sorts  ( spike sorting sessions ) and  experiments  ( which describe stimuli ) .
finally , each sort contains a number of  neurons  , one of whose attributes is a numpy array of spike times .
( b ) example code using neuropy to plot the spike - triggered average ( sta ) of two neurons in response to an m - sequence noise movie ( see text for details ) .
each row corresponds to a neuron , and each column corresponds to the sta within a fixed time range following the m - sequence white noise stimulus . on responses
neuropy relies on a hierarchy of data folders on the disk with a fairly rigid naming scheme , such that animal , track , recording , experiment , and sort ids can be extracted from file and folder names .
all objects have a unique i d under the scope of their parent , but not necessarily under the scope of their grandparent .
all data can be loaded in at once by creating an instance of the data class and then calling its .load ( ) method . however , most often only a subset of data is needed , such as only the data from a given animal , track , or recording . for example
, an object representing recording 92 from the default track of the default animal can be instantiated by typing recording(92 ) at the command line .
this recording 's data can then be loaded from disk into the object by calling its .load ( ) method .
a recording loads the neurons from its default sort , which can also be modified .
some analyses are written as simple methods of one of the data objects , but most have their own separate class which is instantiated by a data object 's method call .
many analyses generate plots , some of them interactive ( such as the population spike raster plot ) , again using matplotlib and wxpython .
currently implemented analyses include interspike interval histograms , instantaneous firing rates and their distributions , cross - correlograms and autocorrelograms , and spike - triggered averages ( stas ) ( dayan and abbott , 2001 ) .
more specialized analyses include binary codes of population spike trains , their correlation coefficient distributions , maximum entropy ising modelling of such codes ( using scipy.maxent ) , and several other related analyses ( schneidman et al . , 2006 ; shlens et al . , 2006 ; spacek et al . , 2007 ) . because of the data - centric organization , new analyses are easy to add .
neuropy is used interactively as a library from the python command prompt , usually in an enhanced shell such as pyshell ( figure 4b ) or the more widely used ipython12 .
an example of neuropy use is shown in figure 5b , which calculates and plots the sta of neurons 2 and 5 of the default animal and track .
the sta estimates a neuron 's spatiotemporal receptive field by averaging the stimulus ( in this case , an m - sequence noise movie ) at fixed time intervals preceding each spike .
recording 92  was recorded during m - sequence noise movie playback , and is used in this example . line  1 imports all of neuropy 's functionality into the local namespace .
next , an object representing recording 92 is instantiated and bound to the name r92 for convenience , and its data is loaded from disk ( lines 23 ) .
its dictionary of available experiments is requested and printed out ( lines 45 ) ; only one experiment is available , with i d 0 .
stas are calculated with respect to this experiment by calling its .sta ( ) method and passing the ids of the desired neurons ( line 6 ) .
 object , which upon further inspection contains two  sta  objects , one per requested neuron ( lines 710 ) .
finally , the stas object 's .plot ( ) method is called with default options , displaying the result for both neurons ( figure 5c ) .
it allows the user to quickly discover what methods and attributes are available for any given object , eliminating the need to recall them from memory ( figure 4b ) .
instead of immediately returning the raw result or plotting it , most analyses in neuropy return an analysis object , which usually has .calc
the .calc ( ) method is run automatically on instantiation , and the results are stored as attributes of the analysis object .
. these can be modified , and .calc ( ) can be called again to update the result attributes .
once satisfied with the calculation , the user can call the .plot ( ) method .
each time a new plot is generated , it does so from the existing results , saving on unnecessary recalculation time .
we have described python packages for three tasks pertinent to systems neuroscience : visual stimulus generation , waveform visualization and spike sorting , and spike train analysis .
python allowed us to meet these software challenges with a level of performance not normally associated with a dynamically typed interpreted language .
performance challenges included time - critical display and communication of visual stimuli , parsing and streaming of multiplexed data from gb sized files , on the fly nyquist interpolation and shd correction , fast execution of non - vectorizable algorithms , and parallelization .
other challenges , whose solutions were simpler than in a statically compiled language , included a cross - platform native gui , the storage and retrieval of relatively complex data structures to and from file ( .parse and .sort files ) , and a command line environment for interactive data analysis .
spyke is the most recent and remains under heavy development , while new analyses are added to neuropy as needed . as with most other python packages , all
neuropy relies on the stimulus description and timing signals generated by dimstim , and on the spike sorting results from spyke .
these three packages depend on many other open source packages , which themselves rely on yet other packages ( e.g. the vision egg currently depends on pyopengl and pygame ) .
because it greatly encourages object - oriented programming , python code is easier to organize and reuse than matlab code .
this is important for scientific code which tends to continually evolve as new avenues are explored .
often , scientific code is quickly written and bug - tested , used once or twice , and then forgotten about , with little chance of re - use outside of copying and pasting .
dimstim , spyke , and neuropy have roughly 3000 , 5000 , and 4000 lines of code respectively ( excluding comments and blank lines ) . fewer lines
make code maintenance easier , not just because there is less code to maintain , but also because each line is closer to all other lines , making it easier to navigate .
we encourage others in neuroscience to consider python for their programming needs , and hope that our three examples ( available at http://swindale.ecc.ubc.ca/code ) may be of use to others , whether directly or otherwise . rallying around a common open - source language may help foster efforts to increase sharing of data and code , efforts deemed necessary ( teeters et al . ,
the authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest .
 electrophysiology is increasingly moving towards highly parallel recording techniques which generate large data sets . 
 we record extracellularly in vivo in cat and rat visual cortex with 54-channel silicon polytrodes , under time - locked visual stimulation , from localized neuronal populations within a cortical column . to help deal with the complexity of generating and analysing these data 
 , we used the python programming language to develop three software projects : one for temporally precise visual stimulus generation (  dimstim  ) ; one for electrophysiological waveform visualization and spike sorting (  spyke  ) ; and one for spike train and stimulus analysis (  neuropy  ) . 
 all three are open source and available for download ( http://swindale.ecc.ubc.ca/code ) . 
 the requirements and solutions for these projects differed greatly , yet we found python to be well suited for all three . here 
 we present our software as a showcase of the extensive capabilities of python in neuroscience . 
Introduction
Dimstim: Visual Stimulus Generation
Usage
Spyke: Waveform Visualization and Spike Sorting
Spike sorting
Neuropy: Spike Train Analysis
Conclusion
Conflict of Interest Statement