PyNeuroTrace - Python code for neural activity time series

Modern techniques for measuring neuronal activity using fluorescent biosensors and ultra-fast microscopy have allowed neuroscientists unprecedented access to neural information processing in vivo . The time series datasets generated from experiments sampling somatic action potentials from populations of neurons, or full-dendritic arbor sampling of populations of synapses, are becoming increasingly larger as new technologies allow for faster acquisition rates and higher temporal resolution of neural signals. Neuronal activities are sourced from an ever-expanding library of fluorescent indicators of distinct measures, including detectors of calcium, membrane voltage, and a range of neurotransmitters and neuromodulators. These biosensors are impacted by their unique molecular kinetics and inherent signal-to-noise properties. The quality of neural signal data sets are also impacted by acquisition instruments, which differ in sensitivity and sampling rate. All of these features, including underlying neural signals, biosensor properties, and microscope capabilities, must be considered during post-imaging signal processing with techniques that can scale to the size of modern neural datasets. To address this problem, here, we describe pyNeuroTrace , an open-source Python library developed to aid in processing neuronal signals from large fluorescent biosensor data sets, which allows dynamic control of filtering and signal processing with these unique aspects in mind before analyses of the underlying neuronal activity can be conducted.


Statement of need
Many neuroscience labs using optophysiological methods, such as sampling neural activity using two-photon microscopy or fiber photometry, typically must create and constantly adjust functions and filters to analyze raw recordings.Currently, there is limited standardization for approaches for signal processing, with techniques and algorithms scattered throughout the literature such as (Friedrich et al., 2017), requiring substantial framework overhead (Giovannucci et al., 2019) or implemented in various programming languages other than Python (Pérez-Ortega et al., 2024).Additionally, other popular tools used to extract calcium signals from imaging experiments often have limited signal processing utilities, lacking functions for calculating Δ /, forcing users to rely on independent methods (Pachitariu et al., 2017).Our library, pyNeuroTrace, seeks to address these problems by providing an analytic package written purely in Python specifically for neuronal activity data sets.Our package includes a collection of filters and algorithms relevant to optophysiological analysis that are implemented in a generalizable manner for time series data sets in either 1D arrays or a collection of recordings in 2D arrays.Additionally, with the increase in acquisition rates of new imaging techniques, we have implemented a subset of these algorithms using GPU-compatible code to significantly increase processing speeds to accommodate large datasets collected, such as those from large population sampling or ultra-fast kilohertz sampling rates.

DeltaF/F
There are several methods for calculating the change of intensity of a fluorescent indicator over time (Grienberger et al., 2022).We implemented the method described by Jia et al. for the calculation of Δ /, which normalizes the signal to a baseline, helping with bleaching or other changes that occur over time, influencing the detection or magnitude of events in the raw signal (Jia et al., 2010).This implementation includes several smoothing steps to mitigate shot noise (Jia et al., 2010).In short,  0 is calculated by finding the minimum signal in a window of the rolling average of the raw signal.Next, Δ is calculated as the difference in the raw signal and  0 , which is then divided by  0 to attain the trace for Δ / 0 .This Δ / 0 signal can be optionally smoothed using an exponentially weighted moving average (EWMA) to remove shot noise.Jia et al. defined their rolling average with the following equation:

𝐹 (𝜏 ) 𝑑𝜏
The variable  0 is defined using a second time constant,  2 , this parameter is the length of the rolling window to search for the minimum smoothed signal value to be used as a baseline: Thus Δ / is calculated using  0 and  where  is the original raw signal: Δ / =  () −  0  0 The two time constants,  1 and  2 , can be selected by users.An additional third parameter  0 is used to optimize the EWMA functions that smooths the final Δ / trace.Modifying these parameters will have a dramatic influence on the output signal.pyNeuroTrace uses defaults proposed by Jia et al. ( 0 = 0.2 s,  1 = 0.75 s, and  2 = 3 s), which work well for imaging at 30 Hz, these values should be adjusted for different imaging parameters.

Okada Filter
We also provide a Python implementation of the Okada Filter (Okada et al., 2016).This filter is designed to filter shot noise from traces in low-signal to noise paradigms, which is common for calcium imaging with two-photon microscopy where the collected photon count is low, and noise from PMT detectors can be nontrivial.This filter is defined by Okada et al. as: In this equation,   is the value in the neural activity trace at time .The value for , which is a coefficient, should be selected so that the product of   −  −1 and   −  +1 causes a sufficiently steep sigmoid curve which functions a binary filter in the equation.This function is equivalent to the following conditional states from Okada et al.: Essentially, in each trace the Okada filter replaces the point   with the average of adjacent values when the product of the differences in adjacent values is greater than zero.One useful characteristic of this smoothing algorithm is that it does not move the start position of events like other algorithms do (Okada et al., 2016)

Nonnegative Deconvolution
When the kinetics of biosensors are known, nonnegative deconvolution (NND) can be useful to deconvolve raw or Δ / traces back into the most likely discrete triggers of increases in the sensors.pyNeuroTrace includes an implementation of an algorithm that can efficiently compute NND in linear time (Podgorski & Haas, 2012).This can also be useful for denoising, as well particularly good at finding smaller magnitude events in fluorescent imaging that are often obfuscated by machine noise (Podgorski & Haas, 2012).

Event Detection
The event detection module uses several strategies to identify neuronal activity events in time series datasets.These methodologies have been previously discussed and compared by Sakaki et al. (Sakaki et al., 2018).These include two generalizable methods and one that requires prior knowledge of recorded event shape.The generalizable methods include filtering the signal through an exponentially weighted moving average (ewma) or cumulative sum of movement above the mean (cusum).The final filter is a matched filter that finds the probability of the trace matching a previously defined shape, such as one described by an exponential rise and decay of calcium signal generated by a genetically encoded calcium indicator (GECI).
Computationally efficient algorithms were selected, allowing for the possibility of applying them in real-time during experiments.

Visualization
pyNeuroTrace has several built-in visualization tools depending on the format of the data.2D arrays of neuronal timeseries can be displayed as heat maps Figure 1 or as individual traces Figure 2. The heatmap is a useful visualization tool for inspecting many traces at once; additionally, at the bottom of the plot, the stimuli timing is displayed if provided Figure 1.This functionality allows for quick visual inspection of from a population of neurons or signals sampled across a neuronal structure, such as a dendritic arbor.
One of the in-built visualizations is specific to the data structure generated by a custom acousto-optic random access multi-photon (RAMP) microscope (Sakaki et al., 2020) Figure 1.This microscope uses acousto-optic deflectors (AODs) to perform inertia-free scanning between preselected points of interest, allowing for extremely fast acquisition rates for sampling neuronal activity throughout complex 3D neural structures.The scan engine of the microscope allows for random-access sampling for imaging activity across the entire dendritic arbor morphology of a single neuron (Sakaki et al., 2020).This type of imaging does not generate a traditional image.The microscope instead links acquired neuronal traces to points of interest organized into a hierarchical tree structure representing the neuronal morphology in a complex data file.pyNeuroTrace will use additional information from this microscope to link structure and function by color-coding where neuronal activity comes from in a sampled neuron.For individual or small numbers of activity traces, pyNeuroTrace has a line plot feature Figure 2.This is an ideal option for inspecting the shape of events, which may be difficult to appreciate from the colormaps in the heatmap visualization.Dotted lines are plotted vertically across the traces of neural activity to indicate when stimulus presentation occurred during an experiment.Additionally, if a record of stimulus or trigger times is provided pyNeuroTrace can plot the average evoked response in a recording.Figure 3 shows the average evoked response to a full-field OFF and ON stimuli from the data in Figure 1.This graph shows the evoked synaptic weights differ for the different visual stimuli presented to the animal.

GPU Acceleration
Several of the filters in pyNeuroTrace have been rewritten to be almost entirely vectorized in their calculations.The benefit is more noticeable when comparing the difference in performance while using large data sets, such as those generated using a longer time series or faster acquisition rates.These vectorized implementations gain further speed by being executed on a GPU using the Cupy Python library (Okuta et al., 2017).The GPU-accelerated filters can be imported from the pyneurotrace.gpu.filtersmodule, and a CUDA-compatible graphics card is required for their execution.This functionality is becoming increasingly crucial as acquisition rates increase for kilohertz imaging of activity (Zhang et al., 2019), which can generate arrays of hundreds of thousands of data points in just a few minutes of recording.
Figure 4 shows the difference in calculating arrays of various sizes using either the CPU or vectorized GPU-based approach of the Δ / function.The CPU used in these calculations was an Intel i5-9600K with six 4.600GHz cores; the GPU was an NVIDIA GeForce RTX 4070 with CUDA Version 12.3.When compared, the GPU calculation was on average in the order of 100 times faster on time series up to 100,000 points in size.To vectorize the functions several were modified.For example the EMWA used to smooth the Δ / signal as described by Jia et al. was changed to an approximation using convolution with an exponential function.The kernel used to perform this is defined as: Where  is defined as: is a user-selected time constant in seconds, which is translated into the number of samples using the acquisition rate used to acquire the data. is a window parameter for the kernel calculated using : This filters for smaller values that have a minuscule influence on the weighted average.The kernel needs to be normalized such that the sum of its elements is 1, to produce smoothing with the same approximate value as the non-vectorized implementation: The normalized kernel is then convolved with the Δ / signal, d: This convolved signal,  is then normalized to the cumulative sum of the exponential kernel: To demonstrate the differences between the CPU and GPU implementations of the EWMA calculations were performed on an array of random values Figure 5.These were generated from the same array using the respective decays for either implementation using the time constant of 50 milliseconds and a sampling rate of 2kHz.Depending on user parameters, the difference between the two outputs typically ranges in magnitude from 1e-16 to 1e-12.These discrepancies can also be attributed to differences in floating-point number accuracy between CPU and GPU calculations.

Animals
Free-swimming albino Xenopus laevis tadpoles were reared in 10% Steinberg's solution and kept in a 12-hour light/dark cycle at 22°C.All experiments were performed on Stage 48 tadpoles.

In vivo Single-Cell Electroporation of DNA and Dye
Single-cell electroporation (SCE) was used to label individual cells with fluorescent calcium indicators (Haas et al., 2001).A single tectal cell was transfected in vivo using SCE with a GCaMP6m plasmid Figure 1.In brief, the tadpole was anesthetized using 0.02% MS-222 in Steinberg's solution.A pulled glass pipette filled with 1 / of plasmid DNA was stimulated using an Axoporator 800A electroporator (Molecular Devices).The electroporator parameters were set to 1 ms, -40 V pulses at 300 Hz for 300 ms.In another tadpole, an individual tectal cell was loaded with 2 mM of CAL-590 3000 Dextran (AAT Bioquest) Figure 4.The parameters for the dye electroporation were set to 1 ms, -40 V pulses at 900 Hz for 50 ms.

In vivo Imaging of Neuronal Activity
Neuronal activity recordings were performed using two-photon microscopy.All animals were reversibly paralyzed by bath application of 2 mM pancuronium dibromide for 5 minutes (Tocris) before being transferred to a custom imaging chamber (Sakaki et al., 2020).The GCaMP6m expressing neuron was imaged using a custom AOD-based RAMP microscope previously described (Sakaki et al., 2020).A Ti:Sapphire laser (Chameleon Vision II, Coherent) tuned to 910 nm was used for two-photon excitation, and emitted light was captured using a 60x, 1.1 NA water-immersion objective (LUMPL, Olympus).Tectal cells filled with CAL-590 AM were imaged using a custom microscope built from a commercial kit (Podgorski, 2021).The excitation source was a 1030 nm custom Ytterbium-doped fiber laser with an average power of 68 W, 150 fs pulse width, and 40 MHz repetition rate (Tangerine-SP, Amplitude).Full-field raster imaging was collected at 13.5 Hz through a 25x 1.0 NA water-immersion objective (HC IRAPO, Leica) Figure 2. The microscope's integrated scan mode was used to capture calcium signals at 2 kHz from a branch of a neuron labeled with CAL-590 3k Dextran Figure 4.All signal processing was performed using pyNeuroTrace.

Ethics Statement
All animal data were collected from experiments reviewed and approved by the University of British Columbia Animal Care Committee under protocols A19-0297 and A24-0049.These experiments were conducted in accordance with the guidelines established by the Canadian Council on Animal Care.

Figure 1 :
Figure 1: An example of a heatmap generated by pyNeuroTrace.A) An example of a lab-specific visualization of points of interest across a dendritic arbor imaged in vivo with an AOD random-access two-photon microscope.Plotted data was recorded from a single dendritic branch in vivo using a random access two-photon microscope.The indicator is a GCaMP6m, a plasmid encoding the GECI was single-cell electroporated (Haas et al., 2001) B) Heatmap of visually evoked responses imaged at 148 Hz, showing intensity (yellow being highest), in addition to stim indicators along the X axis, and branch indicator along the Y.

Figure 2 :
Figure 2: Six traces from neuron cell bodies imaged in vivo with a SLAP2 microscope.Tectal cells were bulk loaded with CAL-590-AM, a calcium sensitive dye.Plotting individual traces better highlights the event shape

Figure 3 :
Figure 3: Average stimulus-locked responses from the same in vivo imaging experiment of the dendritic branch in Figure 1 .The response for each branch node is plotted for two visual stimuli, full field OFF and ON

Figure 4 :
Figure 4: Comparison between Δ / with EWMA calculations for different array sizes using either the CPU (blue) or GPU (orange).

Figure 5 :
Figure 5: Overlay of the EWMA calculations using the CPU implementation and GPU approximation in red and blue.The difference in values from the output is also plotted.The EWMA traces diverge after the initial window length of the GPU version's kernel, but never by more than a factor of 1e-11.