Microelectrode recordings from human peripheral and cranial nerves provide a means to study both afferent and efferent axonal signals at different levels of detail, from multi- to single-unit activity. Their analysis can lead to advancements both in diagnostic and in the understanding of the genesis of neural disorders. However, most of the existing computational toolboxes for the analysis of microneurographic recordings are limited in scope or not open-source. Additionally, conventional burst-based metrics are not suited to analyze pathological conditions and are highly sensitive to distance of the microelectrode tip from the active axons. To address these challenges, we developed an open-source toolbox that offers advanced analysis capabilities for studying neuronal reflexes and physiological responses to peripheral nerve activity. Our toolbox leverages the observation of temporal sequences of action potentials within inherently cyclic signals, introducing innovative methods and indices to enhance analysis accuracy. Importantly, we have designed our computational toolbox to be accessible to novices in biomedical signal processing. This may include researchers and professionals in healthcare domains, such as clinical medicine, life sciences, and related fields. By prioritizing user-friendliness, our software application serves as a valuable resource for the scientific community, allowing to extract advanced metrics of neural activity in short time and evaluate their impact on other physiological variables in a consistent and standardized manner, with the final aim to widen the use of microneurography among researchers and clinicians.

Among the available methods of accessing the nervous system, the technique of microneurography has distinct benefits, such as being minimally invasive and having high spatial and temporal selectivity.1 This electrophysiological method, in which a tungsten microelectrode is inserted percutaneously into an accessible peripheral or cranial nerve, is utilized in over 50 laboratories across the world to directly assess neural activity within the peripheral nervous system, allowing one to record the activity of postganglionic sympathetic axons supplying muscle or skin2–4 as well as somatosensory afferents originating from muscles, joints, or skin.1,5,6 In the following, we will first review the main approaches and analytical tools used to quantify neural activity and study its effects on other response variables and then we will introduce our novel toolbox. Crucially, unlike most available toolboxes or custom software, our toolbox is accessible to novice users in biomedical signal processing. Furthermore, it increases flexibility for specialized analyses with no need of coding skills, decreasing the computational time also for experienced researchers. In this context, the toolbox holds relevance especially in the domain of sympathetic microneurography. Its adaptability across diverse user backgrounds and different recordings, to both multi-unit and single-unit data, extends its utility beyond basic research, offering valuable insights for future clinical applications, where understanding sympathetic neural activity is of paramount importance for elucidating underlying pathophysiological mechanisms and developing targeted therapeutic interventions.

A brief summary of the existing computational toolboxes for analyzing complex microneurographic datasets in human physiology or pathophysiology is reported here. LabChart (ADInstruments, Australia)7 offers a user-friendly interface and a comprehensive set of tools for data processing and analysis, enabling multi-unit burst identification with extraction of key parameters, such as amplitude, duration, and inter-burst intervals, and generating informative graphs and reports, as well as allowing the analysis of single-unit data following extraction of individual action potentials by template-matching. LabChart is accessible to researchers with few coding skills but may have limited options for customizing analysis algorithms or implementing advanced signal processing techniques, restricting flexibility for specialized analyses. Spike28 provides advanced features for processing microneurographic signals, including spike detection and clustering. Spike2 offers a wide range of tools for spectral analysis, correlation analysis, and other advanced analysis techniques. NeuroExplorer9 is specifically designed for neural signal analysis. It offers a range of features for analyzing microneurographic signals, including burst identification, spectral processing, and correlation analysis. While Spike2 and NeuroExplorer are excellent for spike analysis, they may have limited functionality for other types of analyses commonly used for assessing neural activity and its effects at the systemic level. In addition, they require coding skills and some important investment in learning the software. Other software options for analyzing microneurographic recordings include the ZOOM/SC software developed at the University of Umea and the 662C-3 Nerve Traffic Analysis System developed at the University of Iowa. However, these programs are limited in their capacity to perform advanced postprocessing analyses, primarily focusing on data acquisition, visualization, and preliminary examinations, and lack open-access availability. Among the open-source solutions, three complementary options are APTrack,10 openMNGlab,11 and PyDapsys.12 APTrack is very specific to data visualization and preprocessing with identification of action potentials. openMNGlab presents a software customized for microneurography data analysis. It is designed to load data from major data acquisition solutions such as Spike213 and DAPSYS,14 facilitating data visualization, fiber tracking, and feature extraction for subsequent analysis. However, using openMNGlab for data analysis requires specific hardware instrumentation during the acquisition phase for data format compatibility as well as a certain level of proficiency in programming languages for effective utilization. PyDapsys addresses the challenge of accessing electrophysiology data recorded with the proprietary DAPSYS system. By providing open access to DAPSYS-recorded files and enabling conversion to commonly used open research formats, PyDapsys ensures data interoperability, long-term accessibility, and compatibility with other analysis tools. Nevertheless, researchers may still face challenges associated with Python programming prerequisites and some effort for initial setup in this scenario as well.

All these tools are focused on the analysis of axonal signals in the peripheral nervous system, with single-unit recordings consisting of sequences of action potentials of uniform amplitude and morphology and varying temporal firing distributions. Here, we concentrate on the analysis of multi-unit recordings of Muscle Sympathetic Nerve Activity (MSNA) or Skin Sympathetic Nerve Activity (SSNA) as these are some of the most common signals investigated.15 A conventional method used for quantifying MSNA or SSNA involves counting the number of multi-unit bursts generated per minute (burst frequency) or the number of bursts generated in 100 heartbeats, from the integrated or RMS (root mean square)-processed neurogram.16 However, this approach has its limitations, particularly in pathological conditions characterized by elevated sympathetic activation, such as heart failure or hypertension, where these metrics tend to approach saturation levels (100 bursts per 100 heartbeats). Moreover, it loses much of the information contained in a multi-unit recording. Other conventional methods employ amplitude-based measurements from processed microneurographic signals,17,18 but their reliability is compromised by their susceptibility and dependence on proximity of the electrode tip to the active axon, posing challenges for both intra-subject and inter-subject comparisons. In addition, these measurements lack the ability to distinguish between the recruitment of new units and increased firing of active fibers. As detailed below, our approach instead makes use of the amplitude of the raw neurogram, containing spatiotemporal information, to execute unit sorting or clustering.19,20 To interpret the functional role of neural signals, it is necessary to study physiological responses from the observation of action potential temporal sequences.21,22 For a long time, peripheral neurophysiology has been studied with analysis of compound action potentials, time-locked to a stimulus or event, or by assessing changes in firing; but, time-averaging methods neglect the information contained in the exact timing of spikes. Thus, a viable alternative is to consider the precise action potential timing of interest composing the train of spikes, rather than losing important information by signal averaging/integration. Along these lines, several strategies (Fig. 1) can be used to extract useful indexes for decoding the neuronal code, as briefly discussed hereafter. (i) Spiking-latency: the position in time of a spike in relation to previous specific events (i.e., ventricular systoles or touch contact). With this approach, only one action potential is considered. Therefore, this index could be sensitive to individual spikes that may originate spontaneously asynchronously with other fibers in multi-unit analysis. (ii) Spatiotemporal pattern and synchronization: firing of the neuron could be temporally related to the activity of other units. Distinct patterns of spatiotemporal activity could represent certain features of neuromodulation. Synchrony may be important because it promotes temporal summation of the released transmitter, increasing the efficiency of transmission to the neuroeffector by enhancing coordinated activation of post-synaptic receptors. (iii) Spiking-phase: the latency of a single action potential could be assessed relative to a periodic oscillation. For example, sympathetic postganglionic fibers could convey information in the discharge phase relative to the cardiac cycle.23 An increase in baroreceptor (type I) discharge during systole inhibits the excitatory pattern of preganglionic nerves.24–26 This evidence supports the hypothesis that the baroreceptor acts as a phasic forcing input that entrains a central oscillator within the brainstem. Consequently, as a consequence of lack of entrainment, a continuous shift is observed in the phase between the cardiac cycle and MSNA recordings when baroceptor activity reduces.

FIG. 1.

Different metrics for decoding the neuronal code respect to the ECG signal. In the first plot, the raw MSNA signal is displayed, revealing three distinct groups of identified spikes. The second plot features the ECG signal with labeled R-wave peaks. Various latency and phase metrics were derived based on these two signals. Latencies are computed as temporal differences between neural events and the R-wave peaks of the ECG signal (e.g., Unit 1 to R peaks) or among different unit groups (e.g., Unit 1 to Unit 2). Phase indices are calculated relative to the RR interval of the cardiac cycle. These indices offer valuable insights into reflex behaviors (e.g., baroreflex), providing a nuanced understanding of neural activity in relation to cardiac biomechanics. Importantly, these metrics can be easily extended and computed with respect to other oscillatory signals.

FIG. 1.

Different metrics for decoding the neuronal code respect to the ECG signal. In the first plot, the raw MSNA signal is displayed, revealing three distinct groups of identified spikes. The second plot features the ECG signal with labeled R-wave peaks. Various latency and phase metrics were derived based on these two signals. Latencies are computed as temporal differences between neural events and the R-wave peaks of the ECG signal (e.g., Unit 1 to R peaks) or among different unit groups (e.g., Unit 1 to Unit 2). Phase indices are calculated relative to the RR interval of the cardiac cycle. These indices offer valuable insights into reflex behaviors (e.g., baroreflex), providing a nuanced understanding of neural activity in relation to cardiac biomechanics. Importantly, these metrics can be easily extended and computed with respect to other oscillatory signals.

Close modal

Physiological oscillations play a crucial role in various neural processes, including peripheral nerve transduction and information transmission.27 Power spectral analyses, including methods such as fast Fourier transform (FFT),28,29 autoregressive (AR) modeling,30 and maximum entropy methods,31 decompose time series data into frequency components. While these methods may not account for dynamic frequency changes, wavelet analyses32,33 are employed to capture transient fluctuations. Understanding the dynamics of physiological rhythms and their response to perturbations is essential for unraveling the underlying mechanisms and their functional implications. By using phase response curves,34,35 it is possible to quantify the dynamic changes within the phase domain over time. This approach allows for the characterization of phase transitions, phase locking, or phase synchronization, providing a more comprehensive understanding of complex systems. In the context of cardiovascular analysis, for example, phase indexes take into account changes in the duration of the heart cycle (cardiac interval), accounting for variations in heart rate (HRV).36 Therefore, in addition to latencies, which are calculated with respect to a single time event, phase response indexes consider the dynamic characteristics of the reference inherent cycle. It is worth noting that, despite their conceptual clarity, these methods are currently absent in existing computational toolboxes for analyzing the activity recorded in peripheral nerves.

In light of this gap, we developed a practical framework, integrating solutions to address these requirements in our program. Moreover, we integrated cross-wavelet analysis23 searching for variations in amplitude and phase between signals (with different lag times). Although the interpretation of these results may be complex and require an initial investment, standardizing these metrics can pave the way for new advanced levels of analysis, leading to neurophysiological modeling and clinical translatability enabled by bioengineering methods.

To address the limitations of the available tools, reported above, we developed a user-friendly Graphical User Interface (GUI) with open-source computational core to assist researchers in advanced analyses of microneurographic recordings. The application has been made publicly accessible through a dedicated GitHub repository (https://github.com/Neuro-Robotic-Touch-Laboratory/Microneurography-Toolbox). By providing this software tool, we aim to empower researchers with a reliable and user-friendly solution for analyzing microneurographic recordings in complex datasets. The toolbox is implemented with a GUI, released in a compiled version, under an open Creative Commons license. Simultaneously, the source code of the computational core is made available for the sake of transparency, flexibility, and future upgrades and fostering collaboration and adaptation for MATLAB users. Our GUI provides an intuitive, interactive, and accessible platform for conducting advanced analysis tasks related to neural activity within the peripheral nervous system. Additionally, the underlying MATLAB-based code can be upgraded and customized by users with more advanced coding skills, allowing for further enhancements and modifications. The GUI is designed for broad compatibility with diverse acquisition systems, ensuring adaptability to various experimental setups. In addition to flexibility and compatibility, another seemingly contradictory strength is specificity. While existing software, such as Spike2 and NeuroExplorer, are designed for general neural signals, we have fine-tuned and tested advanced capabilities of our toolbox specifically for microneurography. Consequently, the software does not require extensive parameter adjustments for spike detection, regardless of the experimental setup and electrode impedance, ensuring consistent analyses across different microneurographic recordings. Furthermore, the use of this semi-automated program can significantly reduce the time required to analyze each acquired dataset. Moreover, microneurography-derived time series are extracted, including the burst integral function, which can be further used for subsequent processing operations, such as spectral, time–frequency analysis or modeling. The software is suitable for both multi-unit and single-unit data analysis and, by employing cross-wavelet analysis and phase/latency response curves, the program quantifies dynamic changes within time, frequency, and phase domain, characterizing phase locking, transitions, and synchronizations and evaluating prominent effects on other physiological variables over subsequent and previous cyclic events in multichannel recordings.20,37 The software testing encompassed a range of experimental conditions, including studies involving breathing stressors, such as voluntary apneas and periodic breathing, and chemoreflex challenges, such as hypoxic-normocapnic ventilatory responses (HVRs), and hypercapnic-normoxic ventilatory responses (HCVRs), thereby elucidating the dynamic aspects of sympathetic transduction to breathing disorders and chemoreflex challenges. Furthermore, the program can extend to study physiological perturbations akin to those encountered in established protocols such as the Ewing battery or other reflex assessments, i.e., baroreflex with drug challenges (phenylephrine) or ergoreflex.38,39 Finally, the program allows the automatic generation of reports with outputs, graphs, and statistical analyses, eliminating the need to use additional statistical software, which can be costly and require specialized expertise. This feature facilitates access and interpretation of results by various users and the standardization of reports contributes to improve transparency, reproducibility, and collaboration within the field of peripheral neurophysiology. Therefore, another potential advantage of the program lies in facilitating data comparison and archiving. Enabling traceability from experiments to analysis enhances opportunities for data sharing and fosters better understanding among researchers.

Alongside the article, in the GitHub repository of the toolbox, we offer an informative video tutorial serving as a practical guide, effectively showcasing the utilization of the software and its capabilities.

We developed a user-friendly GUI using MATLAB to facilitate interaction with advanced signal processing libraries and algorithms. Subsequently, we created an executable file and an installer to enable the use of the GUI application on computers without a MATLAB license, requiring only the free license MATLAB runtime environment. The software was tested using 35 recordings of MSNA data acquired through a standardized technique described in the literature. These recordings were collected at Fondazione Toscana Gabriele Monasterio, in Pisa, Italy (Ethics Protocol Code: SNAP-HF - Study No. 16314). Some anonymized test data, along with the code, are available on the GitHub repository. Briefly, during the recording the patient maintained a semi-recumbent position on a bed, with the backrest inclined at 45° and the legs horizontal. A cushion was placed under the knee to achieve an angle of 30°–45°. After identifying the possible course of the common peroneal nerve at the level of the head of the fibula, this was confirmed by delivering external electrical impulses (2–10 mA; Stimulus Isolator, ADInstruments). A sterile, insulated tungsten microelectrode (20 mm length, resistance 1 MΩ, Frederick Haer and Co., Bowdoin, ME) was inserted into the common peroneal nerve at the fibular head and an uninsulated microelectrode about 2 cm away served as the reference electrode. Low-intensity electrical impulses (0.2 ms, 1 Hz, 0.01–1.0 mA) were delivered to optimize the electrode placement in the nerve until motor axons were activated (evidenced by muscle contractions in the foot without paresthesia) for values ≤0.02 mA. The microelectrode was further manipulated until observing bursts of MSNA, defined as pulse-synchronous and increasing with a maximal inspiratory apnea. The electrical signal, adequately amplified and filtered (20 000X, bandwidth 0.3–5.0 kHz; NeuroAmpEX, ADInstruments), was digitized (10 kHz), recorded, processed, and analyzed using LabChart's specific functions.

To begin using the toolbox, users should install and launch the application “MNG_toolbox.exe.” The interface comprises several panels, including Start, Burst and Spikes, Time–frequency analysis (spectral and wavelet analysis), Modelling (correlation and stepwise regression), Transduction analysis (annotate and transduction), Report, and Acknowledgements (Fig. 2).

FIG. 2.

Software macro-workflow with the six panels of the GUI: Start, Burst and Spikes, Time-frequency analysis, Modeling, Transduction, and Report.

FIG. 2.

Software macro-workflow with the six panels of the GUI: Start, Burst and Spikes, Time-frequency analysis, Modeling, Transduction, and Report.

Close modal

In the Start panel, users can select the recording to analyze, which can be in “*.ADICHT” (LabChart), “*.MAT” (MATLAB), or “*.h5” (Hierarchical Data Format version 5, HDF5) format. The first two supported file formats are both proprietary, and without the software/license, data cannot be accessed or created. However, utilizing open-source libraries, such as scipy.io.savemat in Python,40 allows for the creation of MATLAB files. The latter is an open-source file format able to support complex heterogeneous data. Instructions on the MATLAB structure required for data storage are provided on the GitHub repository, ensuring that data remain accessible across different platforms. The toolbox allows for channel selection and the generation of a figure displaying the signals from the chosen channels. The interface supports zooming and time axis modifications. Users can specify or create the folder for saving output figures and files. After selecting the file, users can choose the desired channels for analysis and the display order of the selected signals. The start panel allows optical inspection of all the selected signals, including markers, and comments inserted during the acquisition phase (Fig. 3). Within the start panel, the interval of interest can be chosen by selection of comments made during the recording, entering timestamps, manual selection of specific signal attributes, or any combination of the above. Furthermore, the app allows the selection of subintervals for statistical comparisons within the interval of interest. Importantly, all the analyses described below can be conducted separately for each subinterval, allowing for statistical comparisons using parametric or non-parametric methods such as t-tests (paired, unpaired), Mann–Whitney or Wilcoxon, and ANOVA (also for repeated measurements). By switching from the start panel to other functions of the application, several derived signals are calculated if the related fundamental signals are selected. The derived signals comprise detection of R-waves in the ECG including continuous heartrate and HRV, detection of peak inspiration and expiration from respiration belt or spirometry signals including calculation of continuous respiration rate, and extraction of systolic, diastolic, and notch values from blood pressure signals, end tidal CO2 values from CO2 signals, and minute ventilation from spirometry flow signals. For artifact or interference prone signals, the detected R-waves and peak in/expiration can be controlled and manually edited. The software enables versatility by processing the single-channel neural data even in the absence of additional signals such as ECG or respiration signals.

FIG. 3.

(Left) Channels selection to subdivide the dataset to be analyzed; (Right) Start panel for displaying the selected channels with comments and markers (red lines) inserted during the experimental session.

FIG. 3.

(Left) Channels selection to subdivide the dataset to be analyzed; (Right) Start panel for displaying the selected channels with comments and markers (red lines) inserted during the experimental session.

Close modal

The Burst and Spikes panels enable compound and single-unit analysis.

Burst analysis provides various standard MSNA metrics based on burst detection. Users can manually add or remove bursts as necessary. The output consists of visual plots (Fig. 4) and an Excel file (Fig. 5) with burst information. The raw MSNA neurogram is averaged with time constant of 0.2 s, getting the root mean square (RMS) signal. Two parameters predefined by the user are used to identify bursts: a threshold and a time window. A burst to be identified must exceed a certain amplitude threshold and remain above it for a predefined time. To mitigate undesirable baseline variations of the RMS signal, for instance due to the interference of external noise components, the signal can be adjusted through the detrending checkbox. Based on burst identification, various conventional indicators are extracted, such as burst frequency, amplitude, duration, area or integral, and other derived indices. In the burst panel (Fig. 4—top), the raw signal is displayed in the first subplot, the processed RMS signal overlaid with green lines indicating the locations of the ECG R-waves in the second subplot, and a temporal series of each burst area calculated in the third subplot. At the bottom, the distributions of the burst integral, burst amplitude, burst duration, and inter-burst intervals are reported. The application allows one to manually remove bursts that do not meet the criteria for neural activity, such as short-duration artifacts or electromyographic (EMG) interference from contracting muscle. In addition to generating the plotted results, the software also generates an Excel file (.xls) with different sheets in the corresponding analysis folder. A sheet in the Excel file contains the following for each selected interval: number of bursts, burst rate, burst incidence, the median, mean, and standard deviation of the burst integral, burst amplitude, burst duration, RR interval, heart rate, and HRV parameters (namely, SDNN, RMSSD, pNN50, pLF, pHF, LF/HF ratio, VLF, LF, and HF). The other Excel sheets correspond to inference outcomes of statistical hypothesis tests chosen by the user.

FIG. 4.

Burst and spike analysis: (Top) Figure from burst analysis displaying the raw signal, processed RMS signal, and distributions of burst-based metrics; the raw signal is displayed in the first subplot, the processed RMS signal overlaid with green lines indicating the locations of the ECG R-waves in the second subplot, and a temporal series of each burst area calculated in the third subplot. At the bottom the distributions of the burst integral, burst amplitude, burst duration and inter-burst intervals. (Bottom) An output figure from spikes analysis including the raw signal, identified spike events, spikes sorted shapes, and first spike phases for each unit relative to the cardiac cycle.

FIG. 4.

Burst and spike analysis: (Top) Figure from burst analysis displaying the raw signal, processed RMS signal, and distributions of burst-based metrics; the raw signal is displayed in the first subplot, the processed RMS signal overlaid with green lines indicating the locations of the ECG R-waves in the second subplot, and a temporal series of each burst area calculated in the third subplot. At the bottom the distributions of the burst integral, burst amplitude, burst duration and inter-burst intervals. (Bottom) An output figure from spikes analysis including the raw signal, identified spike events, spikes sorted shapes, and first spike phases for each unit relative to the cardiac cycle.

Close modal
FIG. 5.

Burst and spike outputs: (Top left) The window of the app allows for selecting the statistical tests to compare different subintervals of the recording; (Bottom left) An example of an Excel sheet with the statistical results and generated boxplots of the extracted metrics; (Right) An example of an Excel sheet with means and standard deviations of the extracted metrics for different subintervals of the recording.

FIG. 5.

Burst and spike outputs: (Top left) The window of the app allows for selecting the statistical tests to compare different subintervals of the recording; (Bottom left) An example of an Excel sheet with the statistical results and generated boxplots of the extracted metrics; (Right) An example of an Excel sheet with means and standard deviations of the extracted metrics for different subintervals of the recording.

Close modal

Spike analysis focuses on the precise timing of action potentials within the cycle of another oscillating signal. Users can choose to perform spike sorting or clustering, or consider all units in a single cluster. The toolbox detects spikes, calculates firing rates, and visualizes spike phases and latencies relative to specific events. By activating the spikes sorting checkbox, it is possible to perform spike sorting using Wave_clus algorithms,41 alternatively units can be clustered using the spikes clustering option, or considered in a single cluster. Wave_clus has been identified as a top-performing spike sorting algorithm in comparative studies,42 and parameters were extensively fine-tuned and tested with microneurographic data, contributing to result accuracy and reliability. Spikes clustering consists in binning the amplitudes of the three peaks in the detected triphasic action potential templates (the main peak, corresponding to the hyperpolarization of the neuron, is at the center), followed by a combinatorial approach to classify the morphologies.

In each cluster, the firing rate is calculated both beat-to-beat and within a 500 ms temporal window. Additionally, the analysis involves determining the first spike phase, average firing phase, and last spike phase, all relative to the cardiac cycle. The spike analysis panel [Fig. 4(b)] displays identified spike time events highlighted for each unit in the first subplot, the raw signal with spike events, and firing rate for the selected unit in the second subplot. The phases of the first spikes relative to the cardiac cycle are visually represented by a circle ranging from zero to two pi, along with the corresponding phase-locking values.

The toolbox provides methods for investigating the dynamic and spectral features of the recorded signals (Fig. 6). Users can perform time–frequency analysis, such as wavelet transform or spectrogram analysis, to explore the dynamic changes of signals and cross-correlations in phase and amplitude domains. The spectral analysis panel displays the results of a variety of methods for spectral estimation, encompassing both parametric methods, such as autoregressive (AR) models, and non-parametric methods such as the fast Fourier transform (FFT), periodogram, and Welch function. Using the Yule–Walker, Burg, Covariance, and Modified Covariance methods for AR modeling offers a robust approach, considering different algorithmic strengths and adaptability to various signal characteristics and noise types.43,44 This enables the user to choose the most suitable approach for specific analysis needs and signal characteristics. On the wavelet panel, wavelet and cross-wavelet analyses can be performed, using the wavelet coherence toolbox.45 Wavelet and cross-wavelet analyses can be performed on signals acquired during the experimental session or time series derived from previous processing steps. Cross-wavelet analysis can be performed with three different time lags chosen by the user. Wavelet-based analyses provide the ability to compare the dynamic behavior between signals in amplitude and phase in a frequency-specific way.

FIG. 6.

Time-frequency analysis: (Top) The spectral analysis panel displays the results of a variety of methods for spectral estimation, encompassing both non-parametric methods, such as the Fast Fourier Transform (FFT), periodogram, and Welch functions, and parametric methods, using Yule Walker, Burg, Covariance and Modified Covariance methods for AR modeling; (Bottom) Figure from wavelet analysis displays the selected signals, wavelet estimation of each signal, and cross-wavelet correlation with three different time lags chosen by the user.

FIG. 6.

Time-frequency analysis: (Top) The spectral analysis panel displays the results of a variety of methods for spectral estimation, encompassing both non-parametric methods, such as the Fast Fourier Transform (FFT), periodogram, and Welch functions, and parametric methods, using Yule Walker, Burg, Covariance and Modified Covariance methods for AR modeling; (Bottom) Figure from wavelet analysis displays the selected signals, wavelet estimation of each signal, and cross-wavelet correlation with three different time lags chosen by the user.

Close modal

The toolbox offers capabilities for linear cross-correlations and multiple regression between acquired or derived signals (Fig. 7). Additionally, stepwise regression algorithms aid in the systematic selection of predictor variables, facilitating the creation of robust regression models. Users are provided with several preprocessing options to enhance data quality and optimize the output models. These options include selecting the time lag for each signal to be used in the linear model, downsampling and applying a moving average to each signal based on user-defined parameters and identifying and removing outliers from the resulting time series. An Excel file (.xls) is also generated in the corresponding analysis folder with the model and regressors of the predicted variable.

FIG. 7.

Modeling: users can optimize data quality and output models through preprocessing options, such as choosing signal time lags, downsampling, applying moving average, and removing outliers. (Top) Linear cross correlation analysis; (Bottom) Linear regression and Bland-Altman plots with predictor variables selected by stepwise algorithms.

FIG. 7.

Modeling: users can optimize data quality and output models through preprocessing options, such as choosing signal time lags, downsampling, applying moving average, and removing outliers. (Top) Linear cross correlation analysis; (Bottom) Linear regression and Bland-Altman plots with predictor variables selected by stepwise algorithms.

Close modal

The Transduction panels offer annotation capabilities for different physiological signals and evaluate the prominent effects of neural activity on other physiological signals by plotting the response variable associated with different neural metrics and tracking changes over specific time intervals (Fig. 8).

FIG. 8.

(Top) The annotate panel displays MSNA raw neurogram signal, beat-by-beat arterial pressure, electrocardiogram, and respiratory belt signals with corresponding annotated markers (R-wave peaks of the ECG channel; blood pressure end-diastolic/systolic peaks; inspiration and expiration peaks from respiration belt). (Bottom) Plots with Pearson correlations between extracted neural indices (latencies and phases of the first spikes) and response variables (blood pressure and RR interval) evaluated in previous and subsequent cardiac cycles.

FIG. 8.

(Top) The annotate panel displays MSNA raw neurogram signal, beat-by-beat arterial pressure, electrocardiogram, and respiratory belt signals with corresponding annotated markers (R-wave peaks of the ECG channel; blood pressure end-diastolic/systolic peaks; inspiration and expiration peaks from respiration belt). (Bottom) Plots with Pearson correlations between extracted neural indices (latencies and phases of the first spikes) and response variables (blood pressure and RR interval) evaluated in previous and subsequent cardiac cycles.

Close modal

The annotate panel displays MSNA raw neurogram signal, beat-by-beat arterial pressure, electrocardiogram, and respiratory belt signals with corresponding annotated markers (i.e., R-wave peaks of the ECG channel and blood pressure end-diastolic/systolic peaks). The Transduction panel displays pressure and heart rate responses to MSNA spike-based indexes extracted in the phase and time domains. The timing of the neural events is extracted from the MSNA signal using the same method as the previous spike analysis and therefore offers the ability to perform spike sorting or clustering or to consider all spikes as one cluster. Spikes are then evaluated for their latencies and phases with respect to the previous and subsequent R-waves. In particular, the firing rate, the first spike phase, the average firing phase, and the last spike phase are calculated across each cardiac cycle. Important information on sympathetic regulation of the circulation is provided by the pressor and heart rate responses to these spikes-based indexes. The program plots the response variable (i.e., arterial pressure or RR interval) for each heartbeat associated with a spike index and tracks changes over the next and previous 10 cardiac cycles.

The Report panel (Fig. 9) allows users to generate a comprehensive report summarizing the analysis results. The report includes figures, tables, and statistical measures obtained during the analysis process. It provides a concise overview of the recorded data, burst and spike analysis, time–frequency, modeling, and transduction analysis. A significant feature of the Report panel is the ability to perform statistical analysis on the extracted metrics. This functionality not only supports in-depth comparisons within different subintervals of a single dataset (intra-subject studies) but also extends to the comparative analysis across diverse recordings (inter-subjects or prospective studies). The report is saved in PDF format to ease sharing and documentation.

FIG. 9.

The Report panel allows users to generate a comprehensive report summarizing the analysis results, and includes controls for executing statistical tests, comparing different subintervals of a single dataset and/or across diverse recordings.

FIG. 9.

The Report panel allows users to generate a comprehensive report summarizing the analysis results, and includes controls for executing statistical tests, comparing different subintervals of a single dataset and/or across diverse recordings.

Close modal

The open-source toolbox presented in this work provides a valuable computational resource for analyzing microneurographic recordings of neural activity within the peripheral nervous system. The software encompasses two seemingly contradictory strengths: specificity for microneurography and flexibility for further enhancements. Moreover, it offers an intuitive user-friendly GUI, innovative metrics for in-depth microneurographic data analysis while being a cost-effective solution, in both financial and temporal terms. By offering advanced analysis methods and user-friendly interfaces, the toolbox facilitates the exploration of neuronal discharge patterns and their relationship with physiological responses. Its accessibility makes it suitable for researchers and healthcare professionals with diverse backgrounds, promoting collaboration and knowledge sharing in the field of peripheral neurophysiology, bioengineering, and clinical application. Furthermore, future versions of the toolbox will include expanded functionality for evaluating other systems, such as tactile sensory afferents in other multichannel datasets.

The toolbox can be customized and expanded by the user community, ensuring its continuous development and adaptability to future research needs. Depending on the availability of supporting grants, the toolbox will be maintained by the developers and additional groups that may join this initiative.

The development and curation of tools for open science in health sciences and engineering require long-lasting support from institutional grants. The authors would thus like to thank the agencies that contributed to this aim. (i) This project primarily stems from the TUNE-BEAM (TUscany NEtwork for BioElectronic Approaches in Medicine: AI-based predictive algorithms for fine-tuning of electroceutical treatments in neurological, cardiovascular and endocrinological diseases, Cod. H14I20000300002) project funded by the Tuscany Region within the Bando Ricerca Salute 2018 call. Moreover, the development of additional functions and the curation of the toolbox were enabled by complementary grants: (ii) the project is funded by European Union – NextGenerationEU. PNRR – “MNESYS – A Multiscale integrated approach to the study of the nervous system in health and disease” project, Cod. PE_00000006; (iii) the project is funded by European Union – NextGenerationEU. PNRR – “THE - Tuscany Health Ecosystem”, Spoke 8 project, Cod. ECS_00000017; (iv) the project is funded by European Union - NextGenerationEU. PNRR – BRIEF “Biorobotics Research and Innovation Engineering Facilities” PNRR M4C2, Cod. IR0000036, CUP J13C22000400007; (v) the project is funded by the Italian Ministry of Universities and Research via the PRO3 project CRONONC-Lab, CUP J53C22000160005; and (vi) the project is funded by Fundação para a Ciência e a Tecnologia via the FCT project UIDB/00645/2020 (https://doi.org/10.54499/UIDB/00645/2020). Moreover, the authors would like to thank Professor Silvestro Micera for scientific support and stimulating discussions and Dr. Domenico Camboni for technical assistance.

The authors have no conflicts to disclose.

Giacomo D'Alesio: Conceptualization (equal); Data curation (lead); Formal analysis (lead); Investigation (supporting); Methodology (equal); Resources (equal); Software (lead); Validation (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (lead). Lars Ingmar Stumpp: Conceptualization (supporting); Data curation (equal); Formal analysis (equal); Investigation (supporting); Methodology (equal); Software (lead); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (lead). Paolo Sciarrone: Data curation (equal); Investigation (supporting); Methodology (equal); Resources (supporting); Software (supporting); Validation (supporting); Writing – review & editing (equal). Alessandro Navari: Data curation (equal); Investigation (supporting); Methodology (equal); Resources (supporting); Software (supporting); Writing – review & editing (equal). Francesco Gentile: Data curation (equal); Investigation (supporting); Methodology (equal); Resources (supporting); Software (supporting); Validation (supporting); Writing – review & editing (equal). Chiara Borrelli: Investigation (supporting); Methodology (equal); Resources (supporting); Writing – review & editing (equal). Sara Ballanti: Investigation (supporting); Methodology (equal); Resources (supporting); Software (supporting); Writing – review & editing (equal). Eleonora Degl'Innocenti: Investigation (supporting); Methodology (equal); Resources (supporting); Writing – review & editing (equal). Adrian Carrasco: Methodology (supporting); Resources (equal); Software (supporting); Writing – review & editing (equal). Ana Catarina Costa: Methodology (supporting); Resources (equal); Software (supporting); Writing – review & editing (equal). Alexandre Andrade: Methodology (supporting); Resources (equal); Software (supporting); Writing – review & editing (equal). Andrea Mannini: Funding acquisition (supporting); Methodology (supporting); Project administration (lead); Resources (equal); Software (supporting); Writing – review & editing (equal). Vaughan Gary Macefield: Funding acquisition (supporting); Investigation (supporting); Methodology (equal); Resources (supporting); Software (supporting); Supervision (supporting); Validation (supporting); Writing – review & editing (equal). Michele Emdin: Conceptualization (supporting); Investigation (supporting); Methodology (supporting); Resources (lead); Software (supporting); Validation (supporting); Writing – review & editing (equal). Claudio Passino: Conceptualization (supporting); Investigation (supporting); Methodology (supporting); Resources (lead); Software (equal); Validation (supporting); Writing – review & editing (equal). Alberto Mazzoni: Conceptualization (supporting); Data curation (equal); Formal analysis (supporting); Funding acquisition (supporting); Investigation (supporting); Methodology (equal); Project administration (supporting); Resources (equal); Software (equal); Supervision (supporting); Validation (supporting); Visualization (supporting); Writing – original draft (equal); Writing – review & editing (supporting). Alberto Giannoni: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (lead); Project administration (equal); Resources (lead); Software (supporting); Supervision (supporting); Validation (equal); Visualization (supporting); Writing – original draft (supporting); Writing – review & editing (equal). Calogero Maria Oddo: Conceptualization (lead); Data curation (lead); Formal analysis (lead); Funding acquisition (lead); Investigation (lead); Methodology (lead); Project administration (lead); Resources (lead); Software (equal); Supervision (lead); Validation (lead); Visualization (equal); Writing – original draft (lead); Writing – review & editing (lead).

The data that support the findings of this study are openly available in the Microneurographic Toolbox via GitHub at https://github.com/Neuro-Robotic-Touch-Laboratory/Microneurography-Toolbox (Ref. 46).

1.
Å. B.
Vallbo
, “
Microneurography: How it started and how it works
,”
J. Neurophysiol.
120
(
3
),
1415
1427
(
2018
).
2.
Å. B.
Vallbo
,
K. E.
Hagbarth
, and
B. G.
Wallin
, “
Microneurography: How the technique developed and its role in the investigation of the sympathetic nervous system
,”
J. Appl. Physiol.
96
(
4
),
1262
1269
(
2004
).
3.
R. F.
Rea
and
B. G.
Wallin
, “
Sympathetic nerve activity in arm and leg muscles during lower body negative pressure in humans
,”
J. Appl. Physiol.
66
(
6
),
2778
2781
(
1989
).
4.
A.
Shamsuzzaman
, “
Microneurography
,” in
Encyclopedia of the Neurological Sciences
(
Elsevier
,
2014
), pp.
6
10
.
5.
R.
Ackerley
and
R. H.
Watkins
, “
Microneurography as a tool to study the function of individual C-fiber afferents in humans: Responses from nociceptors, thermoreceptors, and mechanoreceptors
,”
J. Neurophysiol.
120
(
6
),
2834
2846
(
2018
).
6.
T.
Mano
,
S.
Iwase
, and
S.
Toma
, “
Microneurography as a tool in clinical neurophysiology to investigate peripheral neural traffic in humans
,”
Clin. Neurophysiol.
117
(
11
),
2357
2384
(
2006
).
7.
See https://www.adinstruments.com/products/labchart for “
LabChart - Life Science Data Acquisition Software | Data Analysis
” (accessed Aug. 7, 2023).
8.
See https://ced.co.uk/products/spkovan for “
CED Spike2: Analysis
” (accessed Aug. 7, 2023).
9.
See https://www.neuroexplorer.com/ for “
NeuroExplorer
” (accessed Aug. 7, 2023).
10.
A. P.
Nickerson
et al, “
Open-source real-time closed-loop electrical threshold tracking for translational pain research
,”
J. Vis. Exp.
194
,
e64898
(
2023
).
11.
F.
Schlebusch
,
F.
Kehrein
,
R.
Röhrig
,
B.
Namer
, and
E.
Kutafina
, “
openMNGlab: Data analysis framework for microneurography – A technical report
,”
Stud. Health Technol. Inform.
283
,
165
171
(
2021
).
12.
P.
Konradi
et al, “
PyDapsys: An open-source library for accessing electrophysiology data recorded with DAPSYS
,”
Front. Neuroinform.
17
,
1250260
(
2023
).
13.
See https://ced.co.uk/products/spkovda for “
CED Spike2: Data acquisition
” (accessed Mar. 6, 2024).
14.
See http://www.dapsys.net/ for “
DAPSYS Data Acquisition Processor System
” (accessed Mar. 6, 2024).
15.
V. G.
Macefield
, “
Sympathetic microneurography
,”
Handb. Clin. Neurol.
117
,
353
364
(
2013
).
16.
G.
Sundlöf
and
B. G.
Wallin
, “
The variability of muscle nerve sympathetic activity in resting recumbent man
,”
J. Physiol.
272
(
2
),
383
397
(
1977
).
17.
V. G.
Macefield
, “
Recording and quantifying sympathetic outflow to muscle and skin in humans: Methods, caveats and challenges
,”
Clin. Auton. Res.
31
(
1
),
59
75
(
2021
).
18.
K.
Notay
et al, “
Validity and reliability of measuring resting muscle sympathetic nerve activity using short sampling durations in healthy humans
,”
J. Appl. Physiol.
121
(
5
),
1065
(
2016
).
19.
S. A.
Klassen
et al, “
Asynchronous action potential discharge in human muscle sympathetic nerve activity
,”
Am. J. Physiol. Heart Circ. Physiol.
317
(
4
),
H754
H764
(
2019
).
20.
B. M.
Shafer
et al, “
Action potential amplitude and baroreflex resetting of action potential clusters mediate hypoxia-induced sympathetic long-term facilitation
,”
J. Physiol.
600
(
13
),
3127
3147
(
2022
).
21.
V. G.
Macefield
and
B. G.
Wallin
, “
Physiological and pathophysiological firing properties of single postganglionic sympathetic neurons in humans
,”
J. Neurophysiol.
119
(
3
),
944
956
(
2018
).
22.
V. G.
Macefield
,
B. G.
Wallin
, and
A. B.
Vallbo
, “
The discharge behaviour of single vasoconstrictor motoneurones in human muscle nerves
,”
J. Physiol.
481
(
Pt 3
),
799
809
(
1994
).
23.
R. W.
De Boer
and
J. M.
Karemaker
, “
Cross-wavelet time-frequency analysis reveals sympathetic contribution to baroreflex sensitivity as cause of variable phase delay between blood pressure and heart rate
,”
Front. Neurosci.
13
,
447984
(
2019
).
24.
S. M.
Barman
and
G. L.
Gebber
, “
“Rapid” rhythmic discharges of sympathetic nerves: Sources, mechanisms of generation, and physiological relevance
,”
J. Biol. Rhythms
15
(
5
),
365
379
(
2000
).
25.
S. C.
Malpas
, “
The rhythmicity of sympathetic nerve activity
,”
Prog. Neurobiol.
56
(
1
),
65
96
(
1998
).
26.
S.
Zhong
,
S. Y. I.
Zhou
,
G. L.
Gebber
, and
S. M.
Barman
, “
Coupled oscillators account for the slow rhythms in sympathetic nerve discharge and phrenic nerve activity
,”
Am. J. Physiol.
272
(
4 Pt 2
),
R1314
R1324
(
1997
).
27.
C. C.
Canavier
, “
Phase-resetting as a tool of information transmission
,”
Curr. Opin. Neurobiol.
31
,
206
213
(
2015
).
28.
R. D.
Berger
,
J. P.
Saul
, and
R. J.
Cohen
, “
Transfer function analysis of autonomic regulation. I. Canine atrial rate response
,”
Am. J. Physiol. Heart Circ. Physiol.
256
(
1
),
H142
H152
(
1989
).
29.
P. B.
Persson
, “
Spectrum analysis of cardiovascular time series
,”
Am. J. Physiol.
273
(
4
),
R1201
R1210
(
1997
).
30.
H.
Akaike
, “
Power spectrum estimation through autoregressive model fitting
,”
Ann. Inst. Stat. Math.
21
(
1
),
407
419
(
1969
).
31.
T. J.
Ulrych
and
T. N.
Bishop
, “
Maximum entropy spectral analysis and autoregressive decomposition
,”
Rev. Geophys.
13
(
1
),
183
200
, https://doi.org/10.1029/RG013i001p00183 (
1975
).
32.
A.
Stefanovska
,
M.
Bracic
, and
H. D.
Kvernmo
, “
Wavelet analysis of oscillations in the peripheral blood circulation measured by laser Doppler technique
,”
IEEE Trans. Biomed. Eng.
46
(
10
),
1230
1239
(
1999
).
33.
C.
Torrence
and
G. P.
Compo
, “
A practical guide to wavelet analysis
,”
Bull. Am. Meteorol. Soc.
79
(
1
),
61
78
(
1998
).
34.
B.
Kralemann
et al, “
In vivo cardiac phase response curve elucidates human respiratory heart rate variability
,”
Nat. Commun.
4
,
2418
(
2013
).
35.
R. M.
Smeal
,
G. B.
Ermentrout
, and
J. A.
White
, “
Phase-response curves and synchronized neural networks
,”
Philos. Trans. R. Soc. B
365
(
1551
),
2407
(
2010
).
36.
R.
Tiwari
,
R.
Kumar
,
S.
Malik
,
T.
Raj
, and
P.
Kumar
, “
Analysis of heart rate variability and implication of different factors on heart rate variability
,”
Curr. Cardiol. Rev.
17
(
5
),
10
(
2021
).
37.
B. M.
Shafer
et al, “
Acute hypoxia elicits lasting reductions in the sympathetic action potential transduction of arterial blood pressure in males
,”
J. Physiol.
601
(
3
),
669
687
(
2023
).
38.
D. J.
Ewing
and
B. F.
Clarke
, “
Diagnosis and management of diabetic autonomic neuropathy
,”
Br. Med. J.
285
(
6346
),
916
(
1982
).
39.
M.
Alam
and
F. H.
Smirk
, “
Observations in man upon a blood pressure raising reflex arising from the voluntary muscles
,”
J. Physiol.
89
(
4
),
372
(
1937
).
40.
See https://docs.scipy.org/doc/scipy/reference/generated/scipy.io.savemat.html for “
scipy.io.savemat—SciPy v1.12.0 Manual
” (accessed Mar. 8, 2024).
41.
R. Q.
Quiroga
,
Z.
Nadasdy
, and
Y.
Ben-Shaul
, “
Unsupervised spike detection and sorting with wavelets and superparamagnetic clustering
,”
Neural Comput.
16
(
8
),
1661
1687
(
2004
).
42.
J.
Wild
,
Z.
Prekopcsak
,
T.
Sieger
,
D.
Novak
, and
R.
Jech
, “
Performance comparison of extracellular spike sorting algorithms for single-channel recordings
,”
J. Neurosci. Methods
203
,
369
376
(
2011
).
43.
J.
Pardey
,
S.
Roberts
, and
L.
Tarassenko
, “
A review of parametric modelling techniques for EEG analysis
,”
Med. Eng. Phys.
18
(
1
),
2
11
(
1996
).
44.
A. S.
Yilmaz
,
A.
Alkan
, and
M. H.
Asyali
, “
Applications of parametric spectral estimation methods on detection of power system harmonics
,”
Electric Power Syst. Res.
78
(
4
),
683
693
(
2008
).
45.
A.
Grinsted
,
J. C.
Moore
, and
S.
Jevrejeva
, “
Application of the cross wavelet transform and wavelet coherence to geophysical time series
,”
Nonlinear Process Geophys.
11
(
5/6
),
561
566
(
2004
).