Jacob De Nobel

Diederick Vermetten

Hao Wang

Carola Doerr

Thomas Bäck

Tuning as a Means of Assessing the Benefits of New Ideas in Interplay with Existing Algorithmic Modules

Introducing new algorithmic ideas is a key part of the continuous improvement of existing optimization algorithms. However, when introducing a new component into an existing algorithm, assessing its potential benefits is a challenging task. Often, the component is added to a default implementation of the underlying algorithm and compared against a limited set of other variants. This assessment ignores any potential interplay with other algorithmic ideas that share the same base algorithm, which is critical in understanding the exact contributions being made. We explore a more extensive procedure, which uses hyperparameter tuning as a means of assessing the benefits of new algorithmic components. This allows for a more robust analysis by not only focusing on the impact on performance, but also by investigating how this performance is achieved. We implement our suggestion in the context of the Modular CMA-ES framework, which was redesigned and extended to include some new modules and several new options for existing modules, mostly focused on the step-size adaptation method. Our analysis highlights the differences between these new modules, and identifies the situations in which they have the largest contribution.

INTRODUCTION

With the continuous increase in interest for the field of optimization, many new algorithms get introduced every year. A large number of these algorithms are not completely novel, but instead add new algorithmic ideas to existing methods. Originally referring to one particular algorithm, CMA-ES has developed into a whole family of algorithms that are built around the core design of the original CMA-ES algorithm from [START_REF] Hansen | Completely Derandomized Self-Adaptation in Evolution Strategies[END_REF]. While this growth of the algorithm set helps to keep improving the state-of-the-art performance, it also raises a simple question: "How to assess the benefits of new algorithmic ideas?"

The naive way of performing such an assessment is to implement the algorithmic idea into a bare-bones version of the base algorithm, and to benchmark it against the default (and maybe some other variants). While this technique does manage to give an indication of the usefulness of the newly introduced component, the results are not always practical and hide important information, since they only consider the idea in isolation. Often, there tends to be an important interplay between algorithmic components, which is completely missed when doing the type of assessment described above.

We aim to provide in this work a roadmap for assessing these algorithmic ideas in a way which takes component interactions into account. This is achieved by considering the different algorithmic ideas as modules in a modular framework. Several of these types of frameworks have been developed over the years [START_REF] Birattari | Classification of Metaheuristics and Design of Experiments for the Analysis of Components[END_REF][START_REF] Lourenço | Evolving Evolutionary Algorithms[END_REF][START_REF] Tavares | On the Evolution of Evolutionary Algorithms[END_REF][START_REF] Sander Van Rijn | Evolving the structure of Evolution Strategies[END_REF][START_REF] Xu | Evaluating Component Solver Contributions to Portfolio-Based Algorithm Selectors[END_REF]. In this work, we present a use case for this roadmap by using the Modular CMA-ES (ModCMA), which is extended from the existing ModEA framework [START_REF] Sander Van Rijn | Evolving the structure of Evolution Strategies[END_REF], by both adding new modules and new options for existing modules (see Section 2 for details). With this modular framework, we show in this work how hyperparameter tuning can be used to assess the contributions of the newly implemented components.

We illustrate how this approach gives a detailed perspective on the benefits of new algorithmic ideas, by not only looking at pure performance metrics, but also considering the interplay with existing modules. We show, among other things, that the introduction of new Step-Size Adaptation (SSA) methods can be beneficial, but that it requires careful consideration of the interactions with other modules, such as the choice of recombination weights. We also discuss the limitations of this approach, and how to best use it to gain the most understanding about these new algorithmic ideas. redesign the framework to our specifications. The modifications we made rendered the name of the framework no longer befitting, as only CMA-ES variants can now be created using the framework, whereas the original framework also supported the design of other evolutionary algorithms. The new framework was dubbed the Modular CMA-ES (ModCMA) and is available as an open source Python package within the IOHprofiler [14] environment 1 . It is integrated with the IOHexperimenter, giving access to a broad set of benchmark problems, including a C++ implementation of the BBOB functions [START_REF] Hansen | Real-Parameter Black-Box Optimization Benchmarking 2009: Noiseless Functions Definitions[END_REF] from the COCO environment [START_REF] Hansen | COCO: A platform for comparing continuous optimizers in a black-box setting[END_REF]. In addition, this allows for easy data logging, which can be used directly with the interactive performance analysis and visualization from IOHanalyzer [START_REF] Wang | IOHanalyzer: Performance Analysis for Iterative Optimization Heuristic[END_REF].

Motivation. The primary goal behind redesigning the framework was to reduce its complexity, and to only include functionality compatible to the CMA-ES and its variants. The reasoning behind this is the fact that the framework mostly revolves around the CMA-ES. Other EAs are available in the framework, but are quite underdeveloped w.r.t. the CMA-ES. Moreover, introducing working interactions between the CMA-ES and operators from other EAs overly complicates the framework's structure. For example, ModEA contains a range of different methods for performing recombination. However, the canonical CMA-ES does not explicitly perform recombination. Instead, it updates its mean 𝑚 by taking a weighted average of the individuals in its current population, which it then uses to sample new individuals from a normal distribution. In other EAs, recombination occurs in a much more pronounced sense, for example by crossover. In order to make the modular algorithm of the CMA-ES function with these other forms of recombination, its original method for "recombination" had to be adapted. The CMA-ES however, is still only able to properly function with one of these recombination methods, the canonical one. As this pattern could be observed in other parts of the framework as well (i.e., mutation, selection), it was decided to remove these other methods all together and to focus solely on the CMA-ES.

The Modular CMA-ES

To design the Modular CMA-ES, we use the implementation from the popular CMA-ES tutorial [START_REF] Hansen | The CMA Evolution Strategy: A Tutorial[END_REF] as a starting point. This work provides a detailed description of the CMA-ES algorithm, including a practical guide to its implementation. From this basic design, we separate the CMA-ES in a number of functionally related blocks, in order to allow a customization of a specific part of the algorithm. This allows us to implement algorithmic variants of the CMA-ES as functional modules. From a user perspective, any of these modules could then be combined in order create a custom instantiation of the CMA-ES, by selecting an option for each available module.

In ModEA, eleven of such modules were already implemented. These were all reimplemented in the Modular CMA-ES, with a few changes to the structure of the options. Specifically, we removed the Pairwise Selection as a module. Instead, we incorporated this option in the Mirrored Sampling module as the option Mirrored Sampling with Pairwise Selection, converting this module from binary 1 https://github.com/IOHprofiler/ModularCMAES to ternary. This is done because the pairwise selection method is not suited for use without mirrored sampling [START_REF] Auger | Mirrored Sampling in Evolution Strategies with Weighted Recombination[END_REF].

We implemented a new module for performing boundary correction (see Section 2.2), and added five alternative options for performing step size adaptation (see Section 2.3). These two extensions to the framework will be the focus of our analysis through out this work. This set of changes give us the following list of modules for the redesigned Modular CMA-ES:

(1) Active Update: Bad candidate solutions are penalized in the covariance matrix update using negative weights [START_REF] Jastrebski | Improving Evolution Strategies through Active Covariance Matrix Adaptation[END_REF]. Note that in [START_REF] Hansen | The CMA Evolution Strategy: A Tutorial[END_REF], this is given as the default version, here we consider it to be optional. (2) Elitism: (𝜇 + 𝜆) -selection instead of (𝜇, 𝜆) -selection.

(3) Orthogonal Sampling: All the newly sampled points in the population are orthonormalized using a Gram-Schmidt procedure [START_REF] Wang | Mirrored Orthogonal Sampling with Pairwise Selection in Evolution Strategies[END_REF]. (4) Sequential Selection: Candidate solution are immediately ranked and compared with the current best solution. If improvement is found, no additional objective function evaluations are performed [START_REF] Brockhoff | Mirrored Sampling and Sequential Selection for Evolution Strategies[END_REF]. (5) Threshold Convergence: A method for balancing exploration with exploitation, scaling the mutation vectors to a required length threshold, which decays over time [START_REF] Piad-Morffis | Evolution strategies with thresheld convergence[END_REF]. [START_REF] Auger | Algorithms (X, sigma, eta): Quasi-random Mutations for Evolution Strategies[END_REF] Step-Size Adaptation: Supplementary to the default Cumulative Step size Adaptation (CSA), Two Point step size Adaption (TPA) [START_REF] Hansen | CMA-ES with Two-Point Step-Size Adaptation[END_REF] is implemented. TPA requires two additional objective function evaluations, used for evaluating both a shorter and a longer version of the population's center of mass. The version which shows the higher objective function value determines whether the step size should be increased or decreased. Five newly added mechanism for performing step size adaptation are implemented. They are described in detail in Section 2.3. (7) Mirrored Sampling: For every newly sampled point, its mirror image is added the population, by reversing its sign [START_REF] Auger | Mirrored Sampling in Evolution Strategies with Weighted Recombination[END_REF]. With Pairwise Selection, only the best point of each mirrored pair is used in recombination. (8) Quasi-Gaussian Sampling: Instead of performing the simple random sampling from the multivariate Gaussian, new solutions can alternatively be drawn from quasi-random sequences (a.k.a. low-discrepancy sequences) [START_REF] Auger | Algorithms (X, sigma, eta): Quasi-random Mutations for Evolution Strategies[END_REF]. We implemented two options for this module, the Halton and Sobol sequences. (9) Recombination Weights: Three options are implemented; 1) default weights (see [START_REF] Hansen | The CMA Evolution Strategy: A Tutorial[END_REF]), 2) equal weights: 𝑤 𝑖 = 1/𝜇, and 3) 𝑤 𝑖 = 1/2 𝑖 + 1/(𝜆2 𝜆 ) for 𝑖 = 1, 2, . . . , 𝜆. (10) Restart Strategy: When the optimization process stagnates, the CMA-ES can be restarted using a restart strategy. Two strategies are implemented. IPOP [START_REF] Auger | A restart CMA evolution strategy with increasing population size[END_REF] increases the size population after every restart by a constant factor. BIPOP [START_REF] Hansen | Benchmarking a BI-Population CMA-ES on the BBOB-2009 Function Testbed[END_REF] also changes the size of the population, but alternates between larger and smaller population sizes. [START_REF] Brockhoff | Mirrored Sampling and Sequential Selection for Evolution Strategies[END_REF] the search domain by applying a boundary correction operation. In Section 2.2, we describe six options for performing boundary correction which have been implemented.

In Table 1, an overview is given of all currently implemented modules and their options in the Modular CMA-ES framework.

Boundary Correction

In the original framework, a boundary correction function taken from [START_REF] Li | Mixed-Integer Evolution Strategies for Parameter Optimization and Their Applications to Medical Image Analysis[END_REF] was implemented, and always applied after each mutation. In some cases, however, this operator can degrade the performance of the algorithm quite drastically. We therefore decided to make the boundary correction optional, and to implement it as a module, for it to only be used when beneficial. A number of different boundary correction strategies were implemented, taken from [START_REF] Caraffini | Infeasibility and structural bias in differential evolution[END_REF]:

Step-Size Adaptation

In this work, we consider a number of alternative step size adaptation mechanisms for new options for the Modular CMA-ES. We take inspiration from [START_REF] Krause | Qualitative and Quantitative Assessment of Step Size Adaptation Rules[END_REF], which provides a qualitative evaluation of multiple step size adaptation mechanisms used in ES. In addition the CSA and TPA step size adaptation methods, which were already implemented, we implemented the following procedures:

(1) Median success rule (MSR) [START_REF] Ouassim | A Median Success Rule for Non-Elitist Evolution Strategies: Study of Feasibility[END_REF]: The MSR mechanism adapts the step-size 𝜎 as follows: it firstly computes a success rate by checking the number of current individuals that are better than some user-defined quantile of the function values in the previous population, then accumulates such success rates in every iteration, and finally decides to increase the step-size if the cumulated values is bigger than 1/2 and decrease it otherwise. (2) Population success rule (PSR) [START_REF] Loshchilov | A Computationally Efficient Limited Memory CMA-ES for Large Scale Optimization[END_REF]: determines the success rate of the current population using a rank-based approach. It firstly sorts all individuals in the current and previous population together, then retrieves the set of ranks of individuals belonging to the current iteration and the one for the previous iteration, and finally calculates the average rank difference between those two sets as the population success rate, which controls the step-size updates.

(3) xNES step-size adaptation (xNES) [START_REF] Glasmachers | Exponential Natural Evolution Strategies[END_REF][START_REF] Krause | Qualitative and Quantitative Assessment of Step Size Adaptation Rules[END_REF][START_REF] Wierstra | Natural Evolution Strategies[END_REF]: calculates the length of each standardized mutation vector and subtracts from it the expected length of the standard Gaussian vector. The resulting difference is then scalarized using the same weights used in the recombination, which is finally fed into an exponential function to generate a multiplicative coefficient to modify the step-size. (4) mean-xNES step-size adaptation (m-XNES) [START_REF] Krause | Qualitative and Quantitative Assessment of Step Size Adaptation Rules[END_REF]: functions similarly to xNES, with the exception that it takes the standardized differential vector between current center of mass and the one in the previous iteration and compares it to the expected length of the standard Gaussian vector. ( 5) xNES with log normal prior Step size adaptation (p-xNES) [START_REF] Krause | Qualitative and Quantitative Assessment of Step Size Adaptation Rules[END_REF]: resembles the principle of self-adaptation for step-sizes, where 𝜆 trial step-sizes are generated from a lognormal distribution which takes the current step-size as its mean and each trial step-size is used to sample a candidate point. To determine the new step-size, this method calculates the weighted sum of the log-transformed trial step-sizes, where those assigned to their corresponding candidate points in the recombination.

INCREMENTAL ASSESSMENT OF MODULE PERFORMANCE

With the introduction of these new module settings, we have a clear use-case for the assessment of algorithmic ideas within the CMA-ES algorithm. Since these options are implemented into a framework with many existing modules, it will not suffice to look at them in isolation. Instead, we should carefully consider the potential interactions with the existing modules and investigate their impact on the empirical performance of ModCMA. Previous work [START_REF] Sander Van Rijn | Algorithm Configuration Data Mining for CMA Evolution Strategies[END_REF] used data from a complete enumeration of all module settings to analyze the contribution of each individual module. However, such an approach becomes intractable when we are confronted with a huge set of algorithmic variants, or more importantly if we aim to obtain the contribution of some new modules implemented incrementally to an existing portfolio of algorithms, which we have investigated previously. In addition, this complete enumeration approach ignores entirely the configuration of continuous strategy parameters, e.g., 𝑐 1 , 𝑐 𝜇 , and 𝑐 𝑐 , which have been shown to significantly impact the per-instance performance of the resulting configurations [START_REF] Belkhir | Per instance algorithm configuration of CMA-ES with limited budget[END_REF].

To properly address the problem of determining the contribution of a single module setting to an existing portfolio of modules, we make use of hyperparameter optimization, which has previously been shown to achieve results comparable to the complete enumeration method, while being much more easily extendable to other hyperparameters [START_REF] Vermetten | Integrated vs. sequential approaches for selecting and tuning CMA-ES variants[END_REF]. We propose the following roadmap to formalize this procedure, which is designed to be generic, so that it can function with any modular algorithm, hyperparameter tuner, and performance metric:

(1) Select a modular implementation of the base algorithm to which the new module has been added, a hyperparameter optimizer and a performance metric. (2) Collect a list of the existing modules and relevant hyperparameters (without the new module to assess). This will be the search space for the hyperparameter optimization. (3) Run the selected hyperparameter optimizer on this search space, ideally for a wide set of relevant benchmark functions. This data will then serve as the baseline performance. (4) Extend the original search space by including the new module to assess, and run the hyperparameter optimization on this extended search space (using the exact same setup as the baseline). ( 5) Compare the data from the baseline to the experiment with the extended search space. This should not only be done from a performance perspective, but also from the resulting configurations themselves. This allows for the analysis of potential interactions between modules.

Performance Measure

Assuming a set of optimization algorithms A = {𝐴 1 , 𝐴 ], and adopt a performance measure which aggregates the hitting times of each of these targets; the Area Under the ECDF Curve (AUC):

AUC(𝐴, 𝑓 , V) = ∫ 𝐵 1 𝐹 (𝑡; 𝐴, 𝑓 , V)d𝑡, 𝐹 (𝑡; 𝐴, 𝑓 , V) = 1 𝑁 |V | 𝑣 ∈V 𝑁 𝑖=1 1(𝑇 (𝐴, 𝑓 , 𝑣, 𝑖) ≤ 𝑡),
where 1 is the characteristic function. We note that most hyperparameter tuning methods are built with minimization in mind. As such, we use the Area Over the Curve (AOC) instead of AUC, since we know AOC(𝐴, 𝑓 , V) = 𝐵 -AUC(𝐴, 𝑓 , V).

In order to collect the AOC measure from the runs of the Mod-CMA, we integrated it into the IOHprofiler [14], which provides ease-of-use logging functionality required to calculate the AOC of each run.

Experimental Overview

In this paper, we use the irace [START_REF] López-Ibáñez | The irace package: Iterated racing for automatic algorithm configuration[END_REF][START_REF] López-Ibáñez | The irace package, Iterated Race for Automatic Algorithm Configuration[END_REF] hyperparameter optimizer. Irace 2 is based on the principle of iterated racing, in where each 2 Implemented in R, freely available at [START_REF] López | The irace Package: Iterated Race for Automatic Algorithm Configuration[END_REF].

race 3 repeatedly executes several configurations until there is a statistically significant reason to discard enough of them to move to the next race (thus inherently allocating more runs to more promising configurations).

Four runs of irace are performed for each of the 24 objective functions in the BBOB single objective noiseless problem suite [START_REF] Hansen | COCO: A platform for comparing continuous optimizers in a black-box setting[END_REF][START_REF] Hansen | Real-Parameter Black-Box Optimization Benchmarking 2009: Noiseless Functions Definitions[END_REF], of which the first function instance is used in 5D. Each run of Irace is given a budget of 1 000 algorithm evaluations, which themselves have a budget of 10 000 • 𝐷 function evaluations. We use the AOC attained by a run of a given configuration as the objective function value for irace. Irace will designate one or more configurations as elites, which are the best configurations found by irace. We validate the performance of these elite configurations by performing 25 validation runs, with the same random seeds for all configurations. We use the results of these runs to assess the final performance.

Conform our roadmap, we define a baseline by tuning the existing modules from ModCMA, which are shown (plain text) in Table 1. In addition, we tune four continuous hyperparameters 𝑐 1 , 𝑐 𝜇 , 𝑐 𝑐 , and 𝑐 𝜎 , which control the dynamics of the adaption of the covariance matrix (𝑐 1 , 𝑐 𝜇 , and 𝑐 𝑐 ) and of the step-size (𝑐 𝜎 ).

We compare two experiments to our baseline where in addition to the existing modules, 1) a number of new SSA methods (see Section 2.3) are included, and 2) a new boundary correction module (see Section 2.2) is added to the tuned parameters 4 . Both of these experiments are using the same experimental setup as the baseline experiment (excluding the tuned parameters). Note that in the boundary correction experiment, the new SSA methods cannot be selected and vice versa.

Single Module Performance

Before considering our proposed method, we run a basic benchmarking experiment on each of the individual module options (including the new options). This is similar to the common approach of benchmarking a new module against a set of other algorithm variants. We show the resulting best single-module configurations (a.k.a. the virtual best solver, VBS for short) relative to the default CMA-ES in Table 2. In this table, we see that among the new modules, only two have been selected: MSR for F23 and m-XNES for F5. We can further look at the over-all contributions of the newly introduced step-size settings by plotting the ECDF-curves over all functions, as done in Figure 1. In this figure, we can clearly see that most methods are quite competitive, with the only exception being xNES, which has a overall worse performance than the others. Overall, the MSR method seems to be quite effective, but there is no strict domination over the other settings.

ANALYSIS AND RESULTS

In this section, we present the results of our hyper-parameter tuning experiment. We consider two paths to analyze the contributions of the newly introduced modules: the performance-perspective and the perspective of the selected modules. We start by examining our baseline. This is followed by an analysis for the performanceperspective and a deeper analysis into the selected modules.

Baseline

As mentioned in Section 3.2, we conduct a baseline tuning experiment. Since we run 4 runs of irace for each function, this results in 4 sets of elites (each set has up to 5 configurations), for which we then perform the verification runs. We plot the distribution of the AOC for each of these configurations in Figure 2, in addition to this, the AOC of the default CMA-ES and the VBS is shown. From this figure, it is clear that the tuning of all parameters at once is much better than simply selecting a single-module variant, as is to be expected. This plot also highlights the variance in performance of the final found configurations. There are two main reasons for Figure 3: Module counts of all elites found in the baselineexperiment, over all 24 BBOB-functions. The option numbers correspond to those in Table 1 this fact: the inherent stochasticity of the CMA-ES itself, and the large impact of the initially generated configurations of irace. We discuss these challenges in detail in Section 5.

From this baseline data, we can also study the resulting configurations themselves. This can be done by aggregating the modules which have been selected in the final elite configurations in the separate irace runs, as is visualized in Figure 3. In this figure, we can see that there is a large variability in the selected module options, which seems to indicate that they are all usable for at least some functions. One notable exception is the weights option "equal", which is chosen in less than 1% of configurations.

Performance analysis

First, we visualize the distributions of the AOC of the single best configuration found in each run of irace (based on the verification runs) in Figure 4. In this plot, we can see that the effect of introducing the new modules is quite mixed. For some functions, performance decreases (e.g., on F8) after introducing new modules, while for others we see the desired improvement (e.g. on F23).

In order to better show these differences, we show in Figure 6 the AOC of the single best configurations found in both the SSA and bound-experiments relative to the best configuration from the baseline. Here we observe a generally negative trend, with outliers in both directions. This seems to indicate that these new modules are not always beneficial to the final performance. For example, we can consider F12, where the configuration found by the baseline has an average AOC of 1 159, while the best configuration found when including the new SSA-methods in the search space reaches an average AOC of 1 480. We show the expected running time of these two configurations in Figure 5, where we can clearly observe this difference. However, we can observe a large variance between runs, which can partly explain poor performance. Indeed, if we look at the average AOC as found during the irace run (instead of the later verification runs), the difference between these two configurations is only 7%, even though the distance between them in the verification runs is much larger. This leads to an important observation about the assessment of the new algorithmic modules: when judging results purely from the average performance measures, it is necessary to also consider the overall variability of the experiment, as well as the inherent stochasticity of the base algorithm. We perform the same procedure for the boundary correction methods. The impact of this module is expected to be smaller, since for most of the "easier" functions, the boundary condition rarely violated. For some of the more challenging functions however, the penalty value given by BBOB function itself might not be sufficient to "guide" the algorithm back in bounds, but an explicit boundary correction could be beneficial in these cases. We can see that this seems to indeed be the case in Figure 6, where on the more complex functions, e.g., F21, the performance is improved when the boundary correction module is tuned.

In Figure 6, we also see that the inclusion of the new SSA methods manages to improve the overall performance for some functions. As an example, on F23 we saw an improvement of 17.1% over the best baseline configuration. If we consider all four elite configurations and compare the average performance differences, the average improvement is even higher, at 22.3%. The stability of this improvement is promising, but in order to fully grasp how the inclusion of the new SSA mechanisms leads to this improvement, we need to analyze the selected modules across these different experiments.

Module Analysis

We have seen that the performance of the elite configuration found on F23 improves when we include the new SSA modules in the search space. In order to identify what this performance can tell us about the new modules themselves, we should study the configurations in more detail. The obvious way to see the difference is by looking at how often the new module options have been selected in the final elite configurations. Over 20 elites, the PSR update was selected 14 times, MSR once, and CSA five times. This shows that these new modules are indeed used in the successful configurations. To see how the inclusion of these module options changes the interactions with the other modules, we look at the combined module activation plot, which is shown in Figure 7. From this figure, we can see that there are some interesting differences between the two sets of configurations: the options for the restart and mirrored module are not as uniform when using the new SSA methods, and the weights option is changed completely. These observations show that there is a clear interplay between these modules.

In addition to module activation, we can look into the distributions of configured continuous hyperparameters. To illustrate this, we study F3, and show the pairwise relations between the four continuous hyperparameters and the final AOC value in Figure 11. From the marginal distribution (shown on the diagonal), we can see that the optimized setting of 𝑐 𝜎 differs the most across the SSA, boundary correction, and the baseline experiments. This is a direct result of the introduction of the new SSA methods, each of which prefers slightly different setting for this parameter. This indicates that even though the final performance of the elite configurations is The lower the line, the better its performance, scaled within each band according to the AOC. The option numbers correspond to those in Table 1. similar between the baseline and the SSA-experiment, the inclusion of new SSA methods clearly changes the selected configurations. We can extend this module analysis to all functions by aggregating the most important differences found between the baseline and SSA-experiments. First, we can plot how often each new module option is selected in the elites for each function, as is done in Figure 8. We can use the same principle to study the interaction with the other modules. For the binary modules, we can directly capture the module difference by looking at which modules occur more or less often in the final set of elites, as is visualized in Figure 9. While this does not directly generalize to modules with more settings, we can create a similar plot for the other modules by considering the overlap in selected module distributions. This is visualized in Figure 10. From these figures, it becomes clear that the elites on some functions are barely affected by the inclusion of the new modules, while others require completely different module settings to properly exploit the changed search space.

The cost of tuning. We should note that only considering the final elite configurations does not tell the full story of a module's contribution. As noted previously, introducing a new module increases the size and complexity of the search space, which has a large impact on the hyperparameter tuning task. If a module is very dependent on the settings of other hyperparameters, this can lead to deterioration of the final results, since the initially sampled configurations are likely to have worse performance than those in the baseline. This is visualized in Figure 12, where this is clearly seen on function F5. This is a linear slope function, but the BBOBspecification does not include a sufficient penalty for leaving the search space. As a result, an algorithm which quickly leaves the search space will reach the required objective value very quickly. Thus, when adding boundary correction methods, 5 6 random configurations are not able to abuse this loophole, leading to a worse initial performance. While for F5, the function is simple enough that the good configurations can still be found (and the inclusion of the default CMA-ES settings in the initial population means that there is always at least one good configuration present), the same issue exists to a lesser extent in other functions. Figure 12 also shows that the "tunability" of modules on different functions varies widely. For instance, on functions F16 -F18, the spread of AOC values is significantly larger than those on functions F19 -F21, suggesting that it is relatively more difficult to tune the modules in the latter since the tuner will very likely take a considerably larger budget to identify optimal configurations. Also, while on some functions it is trivial to get improvement (e.g., F7) over the default CMA-ES, it is a lot more challenging on others, for example on functions F16 -F18.

CHALLENGES

We discuss three key challenges for the module assessment procedure based on hyperparameter optimization as identified in this work.

Influence and stochasticity of the hyperparameter tuning: While we showed that assessing the impact of an algorithmic component by using a hyperparameter tuning approach provides useful insights, there are several factors which can complicate this approach. Since hyperparameter tuning is a very challenging problem, with many different approaches to solving it, the kind of tuner used will have a large impact on the resulting assessment [START_REF] Vermetten | Integrated vs. sequential approaches for selecting and tuning CMA-ES variants[END_REF]. In this paper, we used irace, which tends to focus on converging to a single configuration, instead of covering a large set of different solutions. This necessitates running multiple repetitions of the irace procedure itself, as the initialization might otherwise have too much impact on the final configurations. This can quickly become computationally expensive.

Algorithm-inherent stochasticity: As we discussed in the results, we need to take care when drawing conclusions from the performance of the different CMA-ES configurations. Since CMA-ES is inherently stochastic, the amount of variance of the configurations on a certain function has a large impact on the search procedure of irace. Since we end up selecting elites based on the average performance, we are inherently underestimating the AOC of the final configuration. Even though irace largely mitigates this by using statistical testing in the races to decide when to discard configurations, there will always be some degree of underestimation of the performance (the median performance in the verification runs is 3.4% worse than predicted from the irace runs).

Limits of the per-instance analysis: In the current setup, the performance assessment is done on an per-instance basis. While this can be preferred over tuning for large sets of functions/instances [START_REF] Belkhir | Per instance algorithm configuration of CMA-ES with limited budget[END_REF], it does have some drawbacks. Specifically, if a module is designed to have a good performance over a wide set of functions, but other settings exist for each individual function which outperform it, this new module would not be seen as beneficial. Because of this, we argue that module assessment by hyperparameter tuning should not replace the traditional assessments, but rather complement it for more in-depth, per-instance analysis.

DISCUSSION AND FUTURE WORK

We introduced a roadmap for assessing the performance of individual algorithmic ideas, which takes into account the interplay with other existing settings by comparing the results of hyperparameter tuning. Since this approach requires a modular design to function as intended, we use the Modular CMA-ES framework, which we extended with new modules. Our analysis showed that the newly added step size adaptation mechanisms are not always useful, but do provide clear benefits in several functions. The results also showed that SSA is most useful when combined with a different weights option.

The current version of the Modular CMA-ES framework is a good step in the direction of complete modularization of the CMA-ES algorithm, but some further enhancements can still be made. This would allow for even more precise control over each of the individual components, leading to an ideal testbed for new algorithmic ideas, which can then be evaluated using the approach outlined in this paper. However, since this can be computationally intensive, we should aim to share and reuse data as much as possible, by developing and maintaining a well-organized repository for this type of benchmark data. This does not only reduce the amount of computation needed to test new modules, but it also gives rise to the possibility of testing methods to re-use data from other experiments. Ideally, this would allow for the usage of methods from transfer learning to significantly shorten the time needed to assess a module's performance, even within a large modular search space.

Additionally, we note that while the proposed module assessment is inherently dependent on the used hyperparameter tuning method, the overall procedure remains the same no matter which tuner is used. As a result, the analysis of the results should take into account the particularities of the tuner, such as the way configurations are generated. Further research should still be done into different hyperparameter optimization methods (e.g., SMAC [START_REF] Hutter | Sequential modelbased optimization for general algorithm configuration[END_REF], MIP-EGO [START_REF] Wang | Cooling Strategies for the Moment-Generating Function in Bayesian Global Optimization[END_REF], SPOT [START_REF] Bartz-Beielstein | SPOT: An R Package For Automatic and Interactive Tuning of Optimization Algorithms by Sequential Parameter Optimization[END_REF], GGA [START_REF] Ansótegui | Model-based Genetic Algorithms for Algorithm Configuration[END_REF], hyperband [START_REF] Li | Hyperband: A novel bandit-based approach to hyperparameter optimization[END_REF], etc.). Moreover, one could further investigate a module's contribution to the portfolio by developing a credit assignment scheme, e.g. using Shapley values [START_REF] Fréchette | Using the Shapley Value to Analyze Algorithm Portfolios[END_REF]. Additionally, an analysis pipeline for this type of benchmarking could be designed within existing tools like the IO-Hanalyzer [START_REF] Wang | IOHanalyzer: Performance Analysis for Iterative Optimization Heuristic[END_REF], which would greatly reduce the amount of effort needed to assess new algorithmic ideas.

( 1 )

 1 None: No correction is applied to infeasible coordinates of solutions. (2) Uniform Resample (UR): Replaces all infeasible coordinates of a solution with new coordinates sampled uniformly at random within the search space. (3) Mirror Correction Strategy (MCS): Mirrors all infeasible coordinates of a solution with respect to its closest boundary. (4) Complete One-tailed Normal Correction Strategy (COTN): All infeasible coordinates are replaced with new coordinates inside the search space according to a rescaled one-sided normal distribution centered on the boundary. (5) Saturation Correction Strategy (SCS): All infeasible coordinates is set to the closest corresponding bound. (6) Toroidal Correction Strategy (TCS): All infeasible coordinates get reflected off the opposite boundary.

  of (run, target, ...) pairs

Figure 1 :

 1 Figure 1: ECDF-curve of all single-module stepsize options. Figure generated using IOHanalyzer [41].our baseline. This is followed by an analysis for the performanceperspective and a deeper analysis into the selected modules.

Figure 2 :

 2 Figure 2: Distribution of the area over the ECDF curve for the final elite configuration of the baseline irace runs. All AOC's are averages of 25 verification runs. The VBS singlemodule configurations can be seen in Table2.

Figure 4 :

 4 Figure 4: Distribution of the single best elites from the baseline and the tuning with the additional modules. AOC values are the result of averaging over 25 verification runs.

Figure 5 :

 5 Figure 5: Comparison of the Expected Running Time of the best configurations found on F12 by both the baseline and the SSA experiments. Shaded areas indicate the outer quantiles (20-80).

Figure 6 :

 6 Figure 6: The relative improvements per function of the best configuration found by irace relative to the baseline experiment's best configuration.

Figure 7 :

 7 Figure 7: Combined module activation plot for the elites found in the baseline and SSA experiments, for function 23.The lower the line, the better its performance, scaled within each band according to the AOC. The option numbers correspond to those in Table1.

Figure 8 :

 8 Figure 8: Heatmap showing the fraction of the elite configuration in which each of the options for either SSA (top) or boundary correction (bottom) are active.

Figure 9 :

 9 Figure 9: Heatmap showing difference in the fraction of the elite configuration in which each the of the binary modules are active, between the baseline and the SSA experiment. Positive values indicate a module is turned on more often in the SSA experiments.

Figure 10 :

 10 Figure 10: Heatmap showing difference in the distribution of the ternary modules selected in the final elites, between baseline and the SSA experiment. 0 indicates that the distribution is identical, while 1 indicates that there is no overlap at all.

Figure 11 :

 11 Figure 11: Distribution of the continuous hyperparameters from the elite configurations found in all three experiments.

Figure 12 :

 12 Figure 12: Distribution of the relative AOC values found in the initial race of irace (relative to the default CMA-ES configuration; positive values equate to lower AOC.)

Table 1 :

 1 The modules available for the Modular CMA-ES.

	Boundary Correction: If candidate solutions are sampled
	outside the search domain, they can be transformed back into

Table 2 :

 2 Table

showing the AOC of the best singlemodule configuration for each function (VBS), compared to that of the default CMA-ES. The name of the solver corresponds to the module which is active, e.g. <mod-ule_name>_<option_value>. Note that these values does not include benefits from tuning the continuous hyperparameters, which are set to the default values for all configurations in this table.

Table 2 .

 2

										Option 0
										Option 1
		400								Option 2
	Module Count	200 300
		100
		0	Ac tiv e	El iti st	O rt ho go na l Se qu en tia l Th re sh ol d	SS A	M irr or ed	Sa m pl er	W ei gh ts	Re st ar t

The initial iteration of irace consists of random configurations and the default CMA-ES setting.

All of the code used in these experiments, and the resulting data, is available in[START_REF] Jacob De Nobel | Data and Code from: Tuning as a means of assessing the benefits of new ideas in interplay with existing algorithmic modules[END_REF]


Milan Djilas
email: milan.djilas@inserm.fr

Christine Azevedo-Coste

David Guiraud

Ken Yoshida

Spike Sorting of Muscle Spindle Afferent Nerve Activity Recorded with Thin-Film Intrafascicular Electrodes

Afferent muscle spindle activity in response to passive muscle stretch was recorded in vivo using thin-film longitudinal intrafascicular electrodes. A neural spike detection and classification scheme was developed for the purpose of separating activity of primary and secondary muscle spindle afferents. The algorithm is based on the multiscale continuous wavelet transform using complex wavelets. The detection scheme outperforms the commonly used threshold detection, especially with recordings having low signal-to-noise ratio. Results of classification of units indicate that the developed classifier is able to isolate activity having linear relationship with muscle length, which is a step towards online model-based estimation of muscle length that can be used in a closed-loop functional electrical stimulation system with natural sensory feedback.

Introduction

Functional electrical stimulation (FES) is one solution for restoring movement in paralyzed limbs. In order to correct for disturbances and unavoidable musculoskeletal modeling errors, closed-loop FES is desirable. With the advent of advanced implanted prosthetic interfaces, natural sensors are being explored as an alternative source for feedback information. By directly interfacing the peripheral nerves, it is possible to record signals from natural sensors distributed throughout the body [START_REF] Haugland | Slip information provided by nerve cuff signals: application in closed-loop control of functional electrical stimulation[END_REF][START_REF] Hoffer | Neural signals for command control and feedback in functional neuromuscular stimulation: a review[END_REF][START_REF] Sinkjaer | Integrating sensory nerve signal into neural prosthesis devices[END_REF][START_REF] Strange | Restoration of use of paralyzed limb muscles using sensory nerve signals for state control of FES-assisted walking[END_REF]. Muscle spindles are one type of natural sensor. Their main function is to signal changes in the length of the muscle within which they reside [START_REF] Kandel | Principles of Neural Science[END_REF]. Changes in the muscle length are associated with changes in the angles of the joints that the muscles cross. Therefore, their activity can be used to provide feedback information about muscle state in a closed-loop FES system [START_REF] Yoshida | Closed-loop control of ankle position using muscle afferent feedback with functional neuromuscular stimulation[END_REF], as these sensors remain intact and active below the level of lesion in spinal cord injured patients [START_REF] Mortimer | Motor prostheses[END_REF]. In a number of studies muscle spindle afferent activity was recorded using cuff electrodes to provide natural sensory feedback [START_REF] Haugland | Slip information provided by nerve cuff signals: application in closed-loop control of functional electrical stimulation[END_REF][START_REF] Strange | Restoration of use of paralyzed limb muscles using sensory nerve signals for state control of FES-assisted walking[END_REF][START_REF] Jensen | Improving signal reliability for on-line joint angle estimation from nerve cuff recordings of muscle afferents[END_REF][START_REF] Micera | Neuro-fuzzy extraction of angular information from muscle afferents for ankle control during standing in paraplegic subjects: an animal model[END_REF][START_REF] Riso | Nerve cuff recordings of muscle afferent activity from tibial and peroneal nerves in rabbit during passive ankle motion[END_REF]. Despite their chronic stability and noise immunity, the use of cuff electrodes is ultimately limited by their low selectivity, which requires them to be distributed throughout the body and in close proximity to the sensory end organs targeted for use in feedback [START_REF] Upshaw | Digital signal processing algorithms for the detection of afferent nerve activity recorded from cuff electrodes[END_REF]. An alternative to the cuff is the longitudinal intrafascicular electrode (LIFE). It is designed to be implanted longitudinally within the peripheral nerve where it can record activity from a relatively small population of nerve fibers [START_REF] Chemineau | A modeling study of the recording selectivity of longitudinal intrafascicular electrodes[END_REF][START_REF] Goodall | Separation of action potentials in multiunit intrafascicular recordings[END_REF][START_REF] Goodall | Analysis of single-unit firing patterns in multi-unit intrafascicular recordings[END_REF][START_REF] Malmstrom | Recording properties and biocompatibility of chronically implanted polymer-based intrafascicular electrodes[END_REF][START_REF] Yoshida | Characterization of signals and noise rejection with bipolar longitudinal intrafascicular electrodes[END_REF].

In a previous study a simple corrector controller was implemented to follow a desired joint angle trajectory in the presence of externally applied disturbances [START_REF] Yoshida | Closed-loop control of ankle position using muscle afferent feedback with functional neuromuscular stimulation[END_REF]. The approach proved to be applicable as feedback in online closed-loop control in restrictive conditions (limited motion speed and range), probably due to the unaccounted variation of dynamic sensitivity of the two different types of muscle spindle sensory endings. Type Ia sensory fibers predominately encode information about the rate of change of muscle length, and type II sensory fibers predominately encode information about the muscle length. The former introduce a component in nerve response that makes the relationship between nerve activity and muscle length velocity dependant. The new generation of LIFE, the thin-film LIFE (tfLIFE) [START_REF] Stieglitz | Implantable biomedical microsystems for neural prostheses[END_REF], was used in a recent study in which effort was made towards modeling the neural response of muscle afferents to passive muscle stretch [START_REF] Djilas | Interpretation of muscle spindle afferent nerve response to passive muscle stretch recorded with thin-film longitudinal intrafascicular electrodes[END_REF]. Even though the multichannel tfLIFE provides improved selectivity compared to the older single-channel LIFE, the same effect of velocity dependence in the neural response was observed. Solving this issue would require the decomposition of the multiunit recorded signals into components originating from different fiber types. Isolation of the activity of type II sensory fibers could allow for a linear model approximation to be used to track muscle length variations [START_REF] Azevedo | Towards a model-based estimator of muscle length and force using muscle afferent signals for real time FES control[END_REF][START_REF] Djilas | Interpretation of muscle spindle afferent nerve response to passive muscle stretch recorded with thin-film longitudinal intrafascicular electrodes[END_REF]. In the following we give a brief overview of the state of the art in real-time neural spike sorting of multiunit ENG recorded with LIFE.

Voltage threshold triggering is a common way for detecting neural spikes [START_REF] Lewicki | A review of methods for spike sorting: the detection and classification of neural action potentials[END_REF]. This method requires minimal computing power for signal processing and it is also easy to implement. On the other hand it does not always provide acceptable isolation between nerve and noise spikes. A method using the discrete wavelet transform (DWT) for signal denoising was developed by Donoho [START_REF] Donoho | De-noising by soft-thresholding[END_REF] to detect action potentials buried in noise. It involves thresholding of the detail coefficients in the wavelet decomposition. After denoising, spikes are detected by the simple voltage thresholding method. Eventually, the signals are reconstructed by applying the inverse DWT. The procedure was applied on ENG signals by Diedrich et al. [START_REF] Diedrich | Analysis of raw microneurographic recordings based on wavelet de-noising technique and classification algorithm: wavelet analysis in microneurography[END_REF] and recently on intraneural signals recorded using thin-film LIFE [START_REF] Citi | On the use of wavelet denoising and spike sorting techniques to process electroneurographic signals recorded using intraneural electrodes[END_REF]. In both studies the authors used the Symlet 7 wavelet because of its similarity to typical action potential waveforms in the recordings. No quantitative analysis was given to justify the choice of wavelet and the denoising threshold was chosen empirically. In the latter study cycle spinning was used, a method developed by Coifman and Donoho [START_REF] Coifman | Translation-Invariant De-Noising[END_REF], to reduce but not completely remove the effects of DWT translation variance. Signal detection can also be performed in wavelet space, without the need for reconstructing the signal. That way the required computing power is reduced which is an important issue for real-time implementation. Only being interested in analysis, and not synthesis, also relaxes the constraints on wavelet choice, since not all wavelets are suitable for signal reconstruction. K. H. Kim and S. J. Kim [START_REF] Kim | A wavelet-based method for action potential detection from extracellular neural signal recording with low signal-to-noise ratio[END_REF] have suggested using multiscale wavelet analysis as an equivalent to multiple approximations of matched filters. Their method utilizes the point-wise product of wavelet transform coefficients over several selected scales. The detection method proposed does not require "quantitative" a priori information on either the target signal or background noise, and only involves qualitative information that is common to the neural signal recordings, that is, spike waveform shapes that are most common. Unlike the discrete wavelet transform, the continuous wavelet transform (CWT) can operate at every scale, from that of the original signal up to some maximum scale that one determines by trading off the need for detailed analysis with available computational horsepower. During computation the analyzing wavelet is shifted smoothly over the full domain of the analyzed function. This eliminates the problem of translation variance when using the DWT. The stationary DWT and its equivalent methods reduce but not completely eliminate this problem [START_REF] Misiti | Wavelets and Their Applications[END_REF]. The complex wavelet transform (CoWT) is a complexvalued two-dimensional extension of the standard wavelet transform. It provides a suitable framework to incorporate two wavelets into one transformation, and one wavelet being the real part, and the other the imaginary part of the transform. A recent study demonstrated that spike detection in the auditory nerve based on this approach outperforms a matched filter approach [START_REF] Bourien | Comparison of three spike detectors dedicated to single unit action potentials of the auditory nerve[END_REF].

Detected neural action potentials can be classified using different sets of features. Action potentials waveform shapes depend on neuron type, electrode construction, electrode placement relative to the neuron, and the local properties of the tissue surrounding the electrode. Action potentials from different neurons can be distinctively different, but can also be quite similar. Moreover, high noise levels on weak nerve signals make rapid and accurate classification of spikes challenging. Classification using Fourier transform coefficients as features gives higher classification error rates when compared to methods using time domain features [START_REF] Goodall | Identification of single units in multi-unit recordings from peripheral nerves[END_REF][START_REF] Schmidt | Computer separation of multi-unit neuroelectric data: a review[END_REF]. One method for choosing features automatically is with principal component analysis [START_REF] Glaser | Separation of neuronal activity by waveform analysis[END_REF][START_REF] Glaser | On-line separation of interleaved neuronal pulse sequences[END_REF]. The idea behind principal component analysis (PCA) is to find an ordered set of orthogonal basis vectors that capture the directions in the data of largest variation. Real-time implementation of PCA is also possible using transversal filter structures [START_REF] Schmidt | Computer separation of multi-unit neuroelectric data: a review[END_REF]. To the best of our knowledge, there are no results published about the use of PCA on intrafascicular ENG recordings. McNaughton and Horch reported that classification using artificial neural networks (ANNs) outperformed both methods using time-domain features and methods using Fourier transform coefficients [START_REF] Mcnaughton | Action potential classification with dual channel intrafascicular electrodes[END_REF]. In another study the ANN approach managed to classify 6 out of a total of 10 units present in the signals. With more than 10 units present, the number of separable units dropped [START_REF] Mirfakhraei | Classification of action potentials in multi-unit intrafascicular recordings using neural network pattern-recognition techniques[END_REF]. In these studies a three-layer feedforward ANN was used. The main drawback of using ANNs is that the results of classification are very sensitive to the data set used to train these networks, and obtaining a "good" training set is difficult, especially in the case with low SNR [START_REF] Kim | A wavelet-based method for action potential detection from extracellular neural signal recording with low signal-to-noise ratio[END_REF]. Moreover, there are no criteria for defining the appropriate structure and size of an ANN. Many trials are necessary to obtain good results. At the end, it is uncertain whether acceptable results are obtained or not. Moreover, amplitude and shape of neural spikes change over time, due to electrode drift and fibrous encapsulation of the electrode recording sites (tissue reaction to implantation). Classifiers based on artificial neural networks would require tedious and supervised relearning procedures which are conducted in laboratory conditions [START_REF] Mirfakhraei | Recognition of temporally changing action potentials in multiunit neural recordings[END_REF]. This would not be very practical in real-life working conditions of neuroprosthetic devices.

Independent component analysis is a method that can improve source separation by taking advantage of multiple channels available. However, one restrictive assumption of this approach is that the minimum number of channels must equal the number of sources [START_REF] Lewicki | A review of methods for spike sorting: the detection and classification of neural action potentials[END_REF]. It will be shown later that the number of units picked up by one channel of the LIFE is greater than the number of available recording sites on the electrode. Methods for source separation based on the relative difference between spike features recorded by different channels [START_REF] Mcnaughton | The stereotrode: a new technique for simultaneous isolation of several single units in the central nervous system for multiple unit records[END_REF] are not applicable in the case of the LIFE either. Analysis of recorded data shows that the different channels need not pick up activity from the same units.

In this paper we present a novel approach to neural spike detection and classification based on the continuous wavelet transform using complex wavelets.

Methods

Acute Animal Experiments.

Acute rabbit experiments were conducted on 10 New Zealand white rabbits (median weight 4 kg and standard deviation 0.24 kg) to acquire experimental data for validating the spike sorting algorithm.

2.1.1. Animal Preparation. Anesthesia was induced and maintained throughout the experiments with periodic intramuscular doses of a cocktail of 0.15 mg/kg Midazolam (Dormicum, Alpharma A/S), 0.03 mg/kg Fetanyl and 1 mg/kg Fluranison (combined in Hypnorm, Janssen Pharmaceutica). In order to immobilize the left leg of the rabbit, it was anchored at knee and ankle joints to a fixed mechanical frame using bone pins placed through the distal epiphyses of the femur and tibia. The common calcaneal tendon was attached to the arm of a motorized lever system (Dualmode system 310B Aurora Scientific Inc.) using a yarn of polyaramid fibers (Kevlar 49, Goodfellow Cambridge Ltd). The motorized lever system provided both the actuation and measuring. Pulling the Kevlar fibers produced ankle extension and releasing tension on the Kevlar fibers resulted in ankle flexion (stretched muscle returning to its resting state due to its intrinsic elasticity). Fixation to the mechanical frame insured the elimination of mechanical vibration that might have resulted from a free swinging foot. A tripolar cuff electrode was implanted around the sciatic nerve. It was used to find the length-tension curve for the medial gastrocnemius (MG) muscle. Electrical charge was delivered using a stimulation unit (Grass Technologies SD9), coupled with a photoelectric isolation unit (Grass Technologies PSIU6). The nerve was stimulated with 300 microseconds pulses and a pulse repetition frequency of 2 Hz. The stimulation intensity was set to the level that produced maximal nerve twitch response (maximal amplitude of compound action potential). Keeping the stimulation level constant, muscle length was varied in small incremental steps. Isometric force produced by the stimulated muscle was simultaneously monitored.

A tfLIFE structure was implanted in the tibial branch of the sciatic nerve innervating the MG muscle of the rabbit's left hind limb. It was located 3 cm distal to the cuff electrode implantation site. The tfLIFE enabled the monitoring of multichannel ENG from the fascicle in which the structure was implanted. By having the electrode implanted very close to the muscle, chances of having anything except muscle spindle activity recorded are minimized. Moreover, in order to have purely muscle afferent activity in the recordings, the sciatic nerve was crushed proximally of the cuff and tfLIFE implantation sites using a pair of forceps. It should also be mentioned that with increasing levels of anesthesia, the effects of reflex mechanisms diminish, and decreased tonic stiffness of muscles is observed [START_REF] Pages | Effect of halothane anesthesia on responses from primary endings of muscle spindles in the cat[END_REF][START_REF] Sherrington | On reciprocal innervation of antagonist muscles. Fourth note[END_REF]. As for Golgi tendon organ activity that encodes muscle force information, it accounts for only a small part of the muscle afferent signal under passive conditions. A histological study reports that they account for less than 5% of the total number of receptors in tibialis anterior muscle [START_REF] Wei | Joint angle signaling by muscle spindle receptors[END_REF].

Animals were euthanized at the end of the experiments. All procedures used in experiments were approved by the Danish Committee for the Ethical use of Animals in Research.

Data Acquisition

System. The amplification system consisted of a low-noise preamplifier (AI402, Axon Instruments), followed by a gain-filter amplifier (Cyberamp 380, Axon Instruments). Signals were recorded using a custom modified multichannel digital tape recorder (ADAT-XT, Alesis). ENG data were band pass filtered (4th order Bessel, corners at 0.1 Hz and 10 kHz), amplified (gain 5000), and acquired with a sampling rate of 48 kHz per channel. Out of the eight available electrode sites on the tfLIFE, only the four having the lowest background noise level were selected to be recorded from, due to the limited number of available channels on our recording system.

Signals for driving the motorized lever system (Aurora Scientific) were generated on a portable computer using LabVIEW (National Instruments). Before feeding the analog signal into the lever system input, the signal passed through a low-pass filtering stage (corner at 100 Hz) in order to remove any quantization noise resulting from the D/A conversion. Quantization noise would translate into vibration of the lever arm which could have induced activity of muscle spindles [START_REF] Goodwin | Proprioceptive illusions induced by muscle vibration: contribution by muscle spindles to perception?[END_REF][START_REF] Jansen | The central control of the dynamic response of muscle spindle receptors[END_REF]. Length and force signals were recorded simultaneously with the ENG.

Muscle Stretch Protocol.

The recorded nerve activity is a mixture of activity from two sensory fiber types. A convenient and common way for studying muscle spindle afferent response is applying sinusoidal extensions to a muscle and simultaneously recording the muscle receptor afferent ENG [START_REF] Goodwin | Proprioceptive illusions induced by muscle vibration: contribution by muscle spindles to perception?[END_REF][START_REF] Jansen | The central control of the dynamic response of muscle spindle receptors[END_REF]. With recordings made with such a protocol, it is later possible to analyze the contributions of the two components to the aggregate recorded activity. In our experiments the MG muscle was passively stretched by rotating the ankle in the extension/flexion plane using the lever arm. The initial muscle length was set to the muscle length at which the produced isometric force was maximal. Ankle position was set so it was flexed 90 • and then finely adjusted by experimentally finding the maximum of the length/tension curve.

The muscle was presented with sinusoidal stretches of 2 frequencies: 0.01 Hz and 0.25 Hz. Both stretch profiles had peak-to-peak amplitudes of 4 mm which covered a large portion of the normal range of motion of the ankle. The two frequencies were selected to model slow movements, such as slow walking (0.25 Hz) and correction for postural control (0.01 Hz). They were selected to be below a critical frequency corresponding to a system pole related to the activation and conduction delays of the neuromuscular system at 1-2 Hz. The maximum frequency in the pass band of the system is about 1/5 the pole frequency [START_REF] Yoshida | Closed-loop control of ankle position using muscle afferent feedback with functional neuromuscular stimulation[END_REF]. The faster sine frequency was chosen accordingly. According to Kralj and Bajd, FES assisted walking with minimal energy cost is between 0.35 and 0.56 m/s [START_REF] Kralj | Functional Electrical Stimulation: Standing and Walking after Spinal Cord Injury[END_REF]. If we use FES gait with stride length of 1 m, then the frequency is between 0.35 and 0.56 Hz. Regarding body sway, it has been reported that most of the energy in able-bodied body sway is in the band between 0.02 and 0.2 Hz [START_REF] Demura | Comparison of power spectrum characteristics of body sway during a static upright standing posture in healthy elderly people and young adults[END_REF][START_REF] Yoneda | Frequency analysis of body sway in the upright posture. Statistical study in cases of peripheral vestibular disease[END_REF]. It is more difficult to estimate low frequencies because of the lack of velocity sensitivity and lower levels of neural activity. We chose lower frequencies to demonstrate a worst case scenario. The durations of the recordings were 2 minutes for the slower stretch (to allow for one full cycle of the sinusoid to complete), and just over 1 minute for the faster stretch (4 cycles of the sinusoid). Simultaneous recordings were made from the four intrafascicular electrode sites, together with force and position recordings from the muscle lever system. The muscle was not stimulated while being stretched.

Neural Spike Detection.

The methodology for spike detection we developed is an expansion of the idea of using complex wavelets, so it covers a range of temporal scales. Contrary to the matched filtering approach, where a priori information is necessary about action potential shapes, using the multiscale complex wavelet approach provides multiple approximations of matched filters, making it a more generic approach and probably more robust when addressing the issue of changes in the spike waveforms. A training set of action potentials waveforms is therefore not necessary, which is an advantage of the wavelets approach over using artificial neural networks.

Choice of Complex Wavelet Family and Scale Factors.

The following equation defines the wavelet transform W:

W(α, τ) = +∞ -∞ x(t) 1 √ α Ψ t -τ α dt, ( 1
)
where the real numbers α and τ denote scale and translation, respectively. The wavelet transform essentially performs a correlation analysis between the input signal x and the translated and dilated version of a reference signal called the mother wavelet Ψ. Hence, it would be expected that the output would have local maxima where the input signal most closely resembles the analysis template, that is, the wavelet function. Some wavelet basis functions are similar in shape to neural action potentials. In addition, the basis function is dilated over a range of scales. If the scales are well chosen, the wavelet transform can act as a number of effective approximations of the matched filter, even though the exact action potential waveforms are not known. In the case of complex wavelets the mother wavelet function Ψ is complex and the wavelet transform is also complex. In order to find the optimal complex wavelet, around 30 action potentials with different waveforms were visually identified by inspecting the recorded ENG data and extracted. Five action potentials with distinctly different shapes are shown on Figure 1. The CoWT was computed for all extracted waveforms using a series of complex wavelet families available in the MATLAB wavelets toolbox (The Mathworks). The computations were done using a range of scales to find the optimal scales, which produced the wavelet coefficient with the maximal magnitude. Optimal scales were selected to be the ones that produced at least one coefficient with a magnitude larger than 95% of the maximal CoWT response among all scales. An example of timescale representations of two action potentials with different shapes is shown on Figure 2. In this example, the CoWT was computed using the cgau1 wavelet using scales from 1 to 16. Panel (a) shows coefficient magnitudes for the first action potential waveform. Panel (b) shows only the coefficients above the 95% threshold. On panels (c) and (d) are the corresponding plots for the other action potential waveform. The complex Gaussian family and the complex Morlet family of wavelets both produced well localized peaks in CoWT space, that is, these peaks appeared with a small or zero delay relative to the position of the waveform peaks in the time domain. In terms of the magnitude of the output coefficients, the complex Morlet family produced wavelet coefficients with a 50% lower magnitude compared to the complex Gaussian family. This was consistent for all extracted action potentials. With the complex Shannon family both the peak localization and the magnitude of the output coefficients were poor. Therefore, only the complex Gaussian family of wavelets was considered in the process of selecting the optimal wavelet to be used for spike detection.

Within the complex Gaussian family, as a general rule, higher-order wavelets produced wavelet coefficients with lower magnitudes, for example, cgau1 produced a larger response than cgau2, and cgau2 produced a larger response than cgau3 and so on. Consequently, the cgau1 wavelet was chosen as optimal for neural spike detection. The real and imaginary parts of the cgau1 wavelet are shown on Figure 3. The optimal scales range was from 1 to 6.

For comparison, the above analysis was also performed using non-complex wavelets that support CWT, from which the db2 wavelet produced wavelet coefficients with maximal magnitudes. Compared to the cgau1 wavelet, the range of scales for the db2 wavelet was larger more than two times (scales from 2 to 16). Comparing the inter-quartile ranges of the scale factors, the cgau1 wavelet requires 3 times less scale factors in order to cover all action potential waveforms. This is important when later looking into the implementation of the algorithm. If the number of scale factors were the same, the cgau1 wavelet would require double processing time compared to using the db2 wavelet. As the number of scale factors for the db2 wavelet is more than double, implementing the algorithm using the cgau1 wavelet would require less processing time.

Another benefit of having a smaller range for scale factors is the fact that it should result in better detector specificity. Using the analogy that low scale factors correspond to low frequencies in the signal spectrum and higher scale factors to higher frequency components in the signal spectrum, then a wider range of scale factors would correspond to a wider frequency bandwidth of the transforms and thus there would be more noise influence on detection performance. We will show this to be correct later on.

Algorithm Implementation.

Prior to detection, ENG signals were band pass filtered to remove noise and artifacts. Implementing a matched filter to remove powerline noise in real time is also a possibility. The algorithm needs to be robust enough to adapt to the changing parameters of the powerline noise. Not only does the amplitude of the harmonics change, due to multipath propagation, but so does the frequency [START_REF] Olguín | Adaptive notch filter for EEG signals based on the LMS algorithm with variable step-size parameter[END_REF]. Fitting a sine wave of 50 Hz and a number of its harmonics on the raw data and subtracting the fit give good results. Fitting and subtraction of noise harmonics up until the 5th or 6th harmonic of the noise was sufficient to remove slow baseline oscillations. Performing the fit on 20 millisecond windows (corresponding to 50 Hz) insureed that the fitting algorithm locks onto the phase of the noise. This value is also suitable for real-time implementation since it is below the delay that can be tolerated in closed-loop FES control.

The detection algorithm consisted of finding peaks in the signal transformed into wavelet space that crossed a preset threshold. Wavelet coefficients were computed only at optimal scales. In order to have the scales independent of length of the signal being processed, a windowed continuous wavelet transform was implemented. The transform was computed using a 20 millisecond moving window, matching the window duration from the noise removal step.

Depending on the shape of an action potential, there were cases where computation of the wavelet transform resulted in maxima from different scales to appear at different time instants-if more than one wavelet coefficient had a magnitude larger that a preset detection threshold, multiple detections of the same neural spike occurred. The effect of multiple peaks in wavelet space is illustrated on Figure 2, where for one action potential the maxima are localized at a single time instant (one vertical bar on the panel (b)), while for the other waveform the maxima appear in 2 time instants (2 vertical bars on panel (d)). In the latter case there would be two events detected above the threshold for only one action potential. In order to avoid this, a refractory period was introduced into the algorithm: when a spike is detected, another event can be registered only after the expiration of the preset time interval. (Not to be confused with the nerve fiber refractory period. This refractory period, implemented in the algorithm, is a period during which any other peaks crossing the detector threshold are ignored). Exploratory data analysis on all extracted action potentials showed that the scattering of wavelet coefficient maxima was never larger than 146 microseconds. The refractory period was therefore set to this value.

Classification of Action Potentials.

It is important to note that our objective was not to classify all the action potentials to eventually have accurate information about single-unit activity, but rather to isolate activity from subsets of fibers that would be usable for closed-loop FES. In other words, we were interested in classes that would provide a linear relationship between neural firing and muscle length and therefore independent of muscle stretch velocity.

In principle, there are two approaches when wavelets are used in pattern recognition problems. The first is to use one wavelet to represent all spike shape variations, and the other is to use different wavelets for each of the spike templates. In the former case one wavelet may not be sufficient if there is a large variation between the action potentials waveforms. This is very probable if a larger number of units are present in the recording. In the latter case, the best representation of different spike waveforms would be achieved by designing new wavelets for continuous wavelet transform. The procedure consists of approximating a given pattern using least squares optimization under constraints leading to an admissible wavelet [START_REF] Misiti | Wavelets and Their Applications[END_REF]. We have tried this on spike waveforms extracted from experimentally recorded ENG, and good fits could not be found for all waveforms because of the imposed constraints. Only certain biphasic action potentials produced good fits. Even if a set of wavelets could be found that represents all spike waveforms, the computational power required for real-time parallel computation of CoWT coefficients and processing would exceed the computing power available today. A compromise between the two approaches could be the representation of the full set of action potential waveforms by a reduced set of wavelets, for example, using complex wavelets. Two different action potential waveforms would be represented by the real and imaginary parts of a complex wavelet. Estimates of the number of units from which the electrode picks up activity at the point when the muscle is maximally stretched. Results are shown for each channel for both sinusoidal frequencies, for all 10 experiments.

The multiscale CoWT has an advantage that it also offers a framework for classifying the detected neural spikes. Action potentials differ in their shape and amplitude and it was necessary to choose a feature set and a distance metric with which they would be distinguishable. Exploratory data analysis on the extracted neural spike waveforms indicated that the CoWT coefficients computed using the same range of scale factors as in the detection algorithm could be suitable as classification features. Visual inspection of time-scale plots, like those shown on Figure 2, indicate that they are different for different spike waveforms. Classification could therefore be performed by creating feature vectors using the computed wavelet coefficients and then clustering the data using the Euclidean distance metric. Feature vectors were created by concatenating rows of the time-scale plots for the real and imaginary parts of the CoWT transform, where each row consisted of the CoWT coefficients from a particular scale.

In our context, there are two steps in the classification of detected action potentials. The first is the calibration phase, where the classification is performed offline. The purpose of the calibration is to identify action potentials that encode relevant information for closed-loop control. In other words, we look for classes of action potentials having a linear relationship between firing rate and muscle length. Once they are identified, the second step of the classification is to recognize and track these particular action potentials online in order to estimate muscle state.

Action potentials were classified using k-means. In order to avoid local minima in the optimization, clustering was repeated 50 times (replicates), each time using different starting points. Clustering using less replicates sometimes produced different clustering results using the same data set. After the classification was complete, firing rates of each class were computed. Eventually, the linearity of the relationship between firing rates and muscle length was checked for each class.

One requirement for using k-means algorithm is to know in advance the total number of classes. In order to estimate the number of units each site of the tfLIFE picks up activity from, aggregate afferent firing rates of the postprocessed ENG signals were computed by counting the number of peaks above threshold in a 1-ms moving time window. This window duration matches approximately the absolute refractory period of mammalian sensory nerve fibers [START_REF] Von Brücke | Recovery of responsiveness in motor and sensory fibers during the relative refractory period[END_REF]. Since there cannot be two spikes originating from the same axon within this period, an estimation of the minimal number of axons can be made by counting the number of spikes in the window. Miscounts can occur in cases where two or more spikes overlap. By extending the above analysis onto a number of periods of the faster sine wave, the probability of spike overlap becomes smaller, as it is unlikely that the same subset of spikes will overlap in each period at the same phase. The number of units being picked up by the recording electrode changes depending on how much the muscle is stretched. The more the muscle is stretched, the more units are firing. Values in the Table 1 are estimates of the number of units at the time the muscle is maximally stretched. No significant differences in the numbers are found between the slower and faster sine wave data. On the other hand, there is a large variation when comparing between rabbits and in some cases between different channels of one electrode. Statistically, the median number of units picked up at maximal stretch is 8, with a standard deviation of 2.7. To account for the changing number of classes, the total number of classes for off-line k-means clustering was set to 10, allowing the algorithm to create empty classes.

Evaluation.

In order to be able to evaluate the performance of the spike sorting algorithm, knowledge of the exact timing and class of each action potential in the ENG signal is needed. Because of uncertainty about this information in experimental data, artificial signals based upon recorded action potentials were synthesized. Five action potentials with distinctly different shapes were chosen to represent 5 neural spikes originating from different axons, that is, different spike classes. These are shown on Figure 1. The waveforms were normalized and used to synthesize spike trains. Spike train firing rates were randomly chosen from within the range found in the literature [START_REF] Matthews | The sensitivity of muscle spindle afferents to small sinusoidal changes of length[END_REF]. With the exception of burst firing, muscle spindle afferents can fire with a rate up to about 75 Hz. Burst firing was not considered because the amplitude of an action potential firing in burst mode can vary often as much as 50% [START_REF] Mcnaughton | The stereotrode: a new technique for simultaneous isolation of several single units in the central nervous system for multiple unit records[END_REF]. Spike amplitudes were scaled by integer values ranging from 3 to 6 standard deviations of the background noise level, which corresponds to the range of values found by inspecting the recorded data. These 4 scaling factors represented different SNR levels for which the analysis was performed. Therefore, the SNR is defined here as the ratio of the peak amplitude of the noise-free action potential and the standard deviation of background noise. Signals were synthesized by adding the spike trains onto experimentally recorded background noise. Signals with up to 10 units firing simultaneously were synthesized. Spikes having the same waveform and amplitude were considered to be from the same axon (belonging to the same class). A total of 900 signals were generated: 100 signals with 2 units active, another 100 with 3 units, and so on until 100 signals with 10 units. Steps in the

Results

Performances of the wavelet-based detector and the detector using simple amplitude thresholding are shown on Figure 5 in the form of receiver operating characteristics, or ROC curves. These curves are graphical representations of detector sensitivity versus specificity using a range of detection thresholds. On the whole range of SNR levels the waveletbased detector outperforms the detector based on amplitude thresholding; that is, for any given specificity, the corresponding sensitivity is greater for the wavelet-based detector. The performance gap becomes especially prominent with low SNR.

Compared to detection using the non-complex db2 wavelet, detection using the cgau1 wavelet shows better specificity in ROC space (Figure 6(a)). This is most probably due to the wider range of wavelet transform scale factors required for the db2 wavelet (Figure 6(b)). A wider scales range translates into a wider frequency bandwidth, as explained earlier. For the cgau1 wavelet the scales range was from 1 to 6, and for the db2 wavelet it was from 2 to 14.

Classification results are shown in the form of classification error rates which are ratios of the number of misclassified spikes to the total number of spikes being classified. The wavelet-based classification is compared to two other methods of classification: principal components analysis (PCA) and template matching. Results are shown on Figure 7 starting from the case when only two different spike classes are present in the signal up to the case where 10 units are simultaneously firing. Classification based on template matching produced the highest classification error rates, while wavelet-based and PCA-based approaches showed similar results.

The spike sorting technique was eventually applied on experimentally recorded muscle spindle afferent nerve activity. Only flexion periods of ankle joint motion (stretch periods of the MG muscle) were analyzed. The detection threshold was chosen to be seven times the standard deviations of the background noise level (in wavelet space). Throughout all the trials, this threshold value corresponded to the point on the ROC curves, where the specificity starts to rapidly deteriorate while at the same time there is little improvement in sensitivity.

The detected units were classified into 10 clusters. The analysis was performed on data from all rabbits. Two to three spike classes per rabbit showed a linear relationship between their computed neural firing rate and instantaneous muscle length. Since this relationship was not linear when using the aggregate afferent firing rate, the result is an indication that the algorithm is capable of isolating activity of units less sensitive to muscle stretch velocity. Results from one rabbit are shown on Figure 8 of muscle stretch. A linear regression analysis performed on both shows that the fit on (b) is better.

Discussion

Results show that the CoWT is the preferred method for neural spike detection. Even though this wavelet-based classification does not show improvement in error rates compared to the PCA-based algorithm, the advantage of using the wavelet-based approach is that it provides a unique framework for both spike detection and classification, that is, after computing the complex wavelet coefficients in the detection stage, no additional computation is required in the subsequent classification.

Cluster Centroid Comparison for Different Stretch Rates.

It was of interest to compare activities of the same units for the slow and fast muscle stretch rates in order to see if there is a change in their velocity sensitivity. However, it is difficult to identify same units in the 2 stretching conditions. Even if units have a similar shape, there is no guarantee that they are from the same sensory neuron. As we have seen from the extracted waveforms, there are fewer distinctive shapes than there are action potentials in a recording. Nevertheless, comparison of cluster centroids from the 2 data sets was done. In order to graphically present cluster inter-distance, PCA was used to reduce the dimensionality of the timescale signatures of cluster centroids. Results for one rabbit are shown on Figure 9. Full-circle and empty-triangle markers correspond to cluster centroids from the slower and faster sinusoidal muscle stretch, respectively. Pairs of centroids are easily identifiable, which is a strong indication that they correspond to the same sensory neuron. When relationships between firing rate and muscle length are plotted for the cluster pairs, firing rates from the slower muscle motion exhibit linear or close to linear relationships with length, while firing rates for the faster muscle motion show obvious nonlinear behavior. Under the assumption that the cluster pairs correspond to the same sensory neurons, one could conclude that classes showing different behavior at different muscle stretch rates originate from group Ia sensory fibers, while classes showing linear behavior regardless of the stretch rate originate from group II sensory fibers. There is however not enough evidence to conclude that the units that keep their linear behavior during the faster stretch rate are velocity insensitive. Even the 0.25 Hz muscle stretching is a very slow rate. If a definite conclusion is to be made, the behavior of these units needs to be studied using faster stretch rates than those used in this study.

Aggregate Compared to Single-Cluster Activity.

In cases where the relationship between the aggregate firings rate and muscle length from the faster sinusoidal stretch is nonlinear, activity of single classes exhibits different behavior. A few distinctive cases are shown on Figure 10. Class 1 shows the typical behavior of a sensory neuron with a low activation threshold and saturation point before maximal muscle extension is reached. Activity is registered even at minimal muscle stretch and saturation occurs even before the half of the muscle stretch range. The second class illustrates the difference between saturation and velocity sensitivity. Instead of maintaining a constant firing rate after saturation, with further increase in muscle length the rate starts dropping at maximal stretch, where the velocity of the sinusoidal movement decreases. The third plot is an example of a class with a linear firing rate throughout the whole range of motion. These are the action potentials we are interested in.

The fourth plot most probably represents 2 action potentials with similar waveforms classified in one single class. The first starts firing as soon as the muscle starts stretching, after which it quickly saturates. As the muscle is stretched further, the second starts firing and continues firing until maximal muscle extension is reached. Classes with this kind of behavior could also be useful for control purposes as long as their activity can be modeled as a piece-wise linear function.

Application in Closed-Loop FES.

Results show that the spike sorting algorithm may be useful in closed-loop FES using natural sensory feedback. The spike sorting scheme seems to be capable of isolating the activity of secondary sensory endings from the aggregate neural activity of muscle spindle afferents, making it possible to establish a linear relationship between muscle length and neural firing rate. This result is a step towards an online model-based estimator of muscle length. The more classes having linear behavior are found, the more robust the estimation of muscle state would eventually be. Information about cluster centroids from the calibration step would be used as initial values for the classifier and each detected spike would be assigned to one of these initial clusters. Cluster centroid for that class would be updated, taking into account the signature of the new class member. Updating insures the algorithm can adapt to any slow changes in the shape of action potentials, resulting from electrode drift or fibrous encapsulation of the electrode recording sites. Lastly, not all classes would be needed. Only classes relevant for feedback purposes would be tracked.

To give an idea of the computational requirements, off-line computation of complex wavelet coefficients at 16 scales using the cgau1 wavelet took approximately 4 seconds for 100 ms of data (update time in closed-loop control). It is expected that the processing time would be reduced approximately tenfold when the algorithm is implemented using a low-level programming language, and even more if implemented using an application-specific integrated circuit (ASIC). Therefore, it is expected that realtime implementation is possible.

Regardless of the spike algorithm, the overall performance of spike sorting is highly dependant on the quality of signal recording. Apart from signal processing methods to minimize the effects of noise, work on novel signal acquisition techniques [START_REF] Djilas | Improving the signal-to-noise ratio in recordings with thinfilm longitudinal intrafascicular electrodes using shielding cuffs[END_REF] would be beneficial.

Conclusions

The CoWT offers a convenient framework for neural spike detection and classification: (1) it is time invariant, unlike the undecimated DWT, (2) it needs 3 times less scale factors compared to the CWT to achieve the same performance in terms of spike detection and localization, (3) it requires less computing power (50%), and (4) it is covers a narrower frequency band, resulting in processing that is less sensitive to noise. Contrary to the previous work [START_REF] Citi | On the use of wavelet denoising and spike sorting techniques to process electroneurographic signals recorded using intraneural electrodes[END_REF][START_REF] Kim | A wavelet-based method for action potential detection from extracellular neural signal recording with low signal-to-noise ratio[END_REF], no inverse wavelet transform and no additional computation is necessary in the classification stage as we propose a classification scheme based on the wavelet coefficients used to detect spikes. By doing so, computation cost is reduced, and detection and classification are unified in a single data processing flow.

Evaluation on synthesized multiunit ENG shows the wavelet-based neural spike detection outperforms the threshold detection method, especially in cases of low SNR. Results are even better when comparing the CoWT to CWT.

We demonstrated that the complex wavelet-based classification is able to isolate afferent muscle spindle activity having a linear relationship with muscle length. Results are consistent in all experiments. This is a step towards an online model-based estimator of muscle length that can be used in a closed-loop FES system with natural sensory feedback. We show that by using a single signal processing method based on the complex wavelet transform, it is possible to estimate the length of a muscle based on ENG recordings.

The algorithm could be further optimized in terms of computational speed. Reducing the number of scaling factors used to compute the CoWT could perhaps produce negligible deterioration in performance. We were conservative when making the initial choice in order to demonstrate feasibility and to focus on the method itself. This potential improvement will be investigated in more detail before proceeding with the real-time hardware implementation of the spike sorting scheme.

Figure 1 :

 1 Figure 1: Five distinctive action potential waveforms extracted from experimentally recorded data. Length of each trace is 1 milliseconds.

Figure 2 :Figure 3 :

 23 Figure 2: Time-scale plots of two action potential waveforms with different waveform shapes (a, c) and the same plots showing only wavelet coefficients larger than 95% of the coefficient with the maximal magnitude (b, d).

  classes together, no noise All classes together with added background noise

Figure 4 :

 4 Figure 4: Steps in the creation of a synthetic signal. Top 3 traces are the generated action potential trains for 3 units. The onset of the first spike train is at t = 1 second, the second train at t = 2 second, and so on. The fourth trace is the superposition of traces 1-3. Background noise is eventually added resulting in a signal on the bottom trace.

Figure 5 :Figure 6 :

 56 Figure 5: ROC curves for four SNR levels. Performances of a simple threshold detector (empty-circle line) and the wavelet-based detector (full-triangle line) are compared.

Figure 7 :

 7 Figure 7: Classification error rates depending on the number of units simultaneously active, compared for three classification approaches: template matching (black), principal components analysis (gray), and wavelet-based (white).

Figure 8 :

 8 Figure 8: Afferent neural firing rate versus muscle length. Aggregate activity of all detected spikes (a) and activity from 2 clusters having a good linear fit to the data (b) were used to compute the firing rates. Linear regression analyses are the full lines. Muscle length was normalized by 4 mm.

Figure 9 :

 9 Figure 9: Comparison of cluster centroids for the slow (full circles) and fast (triangles) muscle stretch rates.

10 ComputationalFigure 10 :

 1010 Figure 10: Typical profiles of firing rates versus muscle length for single clusters. On each plot the abscissa is the normalized muscle length and the ordinate is the firing rate in spikes per second.

Table 1 :

 1 Estimates of the number of units.

		0.01 Hz sinusoidal stretch	0.25 Hz sinusoidal stretch
	Rabbit Ch 1 Ch 2 Ch 3	Ch 4	Ch 1 Ch 2 Ch 3 Ch 4
	1	1 0	9	1 0	8	9	9	9	6
	2	13	13	10	12	12	12	9	12
	3	7	8	5	6	7	7	6	7
	4	7	8	5	6	7	7	4	6
	5	1 5	1 5	5	3	1 4	1 4	4	3
	6	8	9	8	8	7	7	8	8
	7	6	4	7	10	6	3	6	8
	8	5	11	8	12	6	10	10	11
	9	1 1	9	8	8	1 0	9	8	8
	10	9	7	7	9	8	8	7	9

Acknowledgments

The authors would like to thank Henrik Barlebo, Ole Sørensen, Torben Madsen, Jens Sørensen and the staff at the Department of Pathology at Aarhus University Hospital Aalborg for providing technical support for the acute animal work. The financial support for this work was provided by the EADS Corporate Foundation.