Methods for frequency and correlation analyses of neurograms

Many physiological functions are based on motor rhythmic activities, among them breathing is a vital issue. The method presented here, or ‘temporal grid extraction’, aims at characterizing the temporal organization of such an activity. Beyond the measurement of the fundamental frequency, deﬁning the successive cycles, some signal processing tools are helpful in order to look for the presence of higher frequency components that potentially structure these cycles. The method is applied to neurograms recorded from frog brainstem preparations, where two cycle types, buccal and lung cycles, may alternate. It relies on: Using this method, the maps and proﬁles have revealed that a high frequency clock buccal lung


Introduction
Most of rhythmic motor behaviors are generated by central pattern-generators (CPG) whose mechanisms must be investigated in details in order to support the hypotheses implemented in computational models ( [3,6] ). Such CPG produce periodic excitatory inputs to motoneurons resulting in rhythmic muscle contractions. The rhythmic movements correspond to fundamental functions as breathing, chewing, suckling, walking, running, flying, swimming, etc... They can occur for short, long or permanent periods. Some rhythmic behaviors exhibit several modalities depending on feedback sensory inputs. For the locomotor rhythms, short term cycle-by-cycle feedback signals allows a fast adaptation to the environmental conditions and for the respiratory rhythms, the CPG is sensitive to hypercapnia and hypoxia. Moreover, in the case of breathing, longer term modulations occur in frequency and amplitude to meet the metabolism needs. Such rhythmic activities can be analyzed by signal processing tools which are adapted to extract specific features from the signals they generate. For the ventilatory system, the periodic breathing activity was recorded as neurograms from the phrenic nerve in mammals (for example in mice [7] or in piglet [1] ). In amphibians, the ventilatory activity is recorded in brainstem preparations, as neurograms from cranial motor roots. This fictive respiratory in vitro activity exhibit patterns corresponding to the respiratory movements recorded in vivo ( [15] ). In this paper, we present computational methods used to define several descriptors characterizing the rhythmic patterns displayed in neurograms recorded from brainstem preparations of amphibians ( [8,9] ). The motor activity recorded from cranial motor roots V or VII, shows quasiregular low-amplitude bursts associated with gill ventilation in frog tadpoles and identified as buccal cycles. In post-metamorphic tadpoles and adult frogs, interspersed with these buccal cycles, there are episodic, high amplitude bursts ensuring lung ventilation for gas exchange during air breathing ( [11,14] ). The aims of the signal processing presented here are: (1) to uncover an internal organization of the rhythmic activities and (2) to characterize rhythm stability and/or variability. These tools allow comparisons of the effect of different physiological conditions on these rhythm characteristics. We analyze these periodic signals in order to characterize them in both the frequency and the temporal domains in the aim to define key features of the motor cycles these signals are made of. We also define tools that allow comparison of these key features in order to evaluate the differences between preparations and/or changes in the recording conditions.

First step signal processing: filtering, segmenting, labelling
The data illustrating the methods presented in the next sections correspond to recordings from brainstem preparations. The neurograms are recorded at the level of V and VII motor roots, in two amphibian species, Lithobates catesbaianus and Pelophylax ridibundus , at different developmental stages, pre-and post-metamorphic, and in different conditions, as control condition, hypoxia or hypercapnia. The methods presented here focus on the signal processing tools only, since all the details concerning the experimental protocols, the resulting data and the physiological interpretation of these data are available in [9] for Pelophylax ridibundus , and in [8] for Lithobates catesbaianus . All these tools are implemented in Matlab, using standard functions that are mentioned at each corresponding step of the signal processing.
In the neurograms, two burst (or cycle) types are present accounting for two different functional modalities involved in ventilation, the lung bursts and the buccal bursts. The first step is to segment the signal into cycles, each cycle containing one burst.

Filtering with RMS and Segmenting the signals into cycles
Each neurogram is sampled at 20 0 0 Hz, which means, according to Nyquist-Shannon, that it contains information in the frequency range [0-10 0 0] Hz. As the respiratory frequencies are around 1 Hz, the neurogram is low-pass filtered using a zero-phase moving root mean square (RMS) ( [2] ) with a rectangular window of width W 0 adapted to the fundamental signal frequency, i.e. the buccal . In C2 and C3, local minima of S 200 are marked by blue vertical small bars indicating the segmentation of the signal into individual cycles. This segmentation is also applied to S 10 (B2 and B3). In B3, red stars indicate the maximal amplitude recorded for each cycle. In this example, the small burst corresponds to a buccal burst, and the high one corresponds to a lung burst.
frequency. In the case of the neurograms recorded in L. catesbaianus , W 0 is set at 200 ms producing a S 200 signal ( Fig. 1 C). The cut-off frequency of this low pass filter is 1.6 Hz which allows low frequency analysis around the fundamental frequency of 1 Hz. In the case of P. ridibundus neurograms, whose typical fundamental frequency is around 2 Hz, i.e. about twice the fundamental frequency in L. catesbaianus , the value of W 0 is set at 100 ms leading to a S 100 signal (the matlab function used for this filter process is "filtfilt").
This resulting smoothed filtered signal, S 200 or S 100 , where all the high frequency noises haves been eliminated, allows the segmentation of the signal into individual successive cycles (see Fig. 1 C2 and C3) by automatic detection of the local minima of this signal (the matlab function used for this detection process is "findpeaks").
This segmentation process defines all the cycles of the neurogram. In order to analyze the intra-cycle temporal organization, the previous cycle segmentation is applied to a new signal computed by RMS filtering of the neurogram, S 10 signal, using a window width set at 10 ms ( Fig. 1 B2). This new filter has a cut-off frequency of 32 Hz.

Labelling the cycles and building separate buccal and lung signals
Since it is known that, in vivo and in vitro, the lung bursts have higher amplitude than the buccal ones ( [14,15] ), the maximal amplitude value is recorded ( Fig. 1 B3) for each cycle represented as a portion of S 10 signal in order to identify this cycle as buccal or lung. Taking a segmented S 10 signal computed from a given neurogram, all the maximal amplitude values are recorded (one maximal amplitude for each cycle) and are ranked from the smaller to the higher in a monotonic increasing function ( Fig. 2 A1). If an abrupt increase appears in this function, it means that there are lung cycles in this neurogram: the amplitude that corresponds to the inflection point of the abrupt increase defines the amplitude threshold between buccal and lung bursts. Correspondingly, the amplitude distribution is bimodal ( Fig. 2 A2). Each cycle can thus be labelled either 'buccal' or 'lung' depending on its maximal amplitude when compared to the threshold. In order to analyze separately the buccal cycle population and the lung cycle population, once segmented and labelled, all the buccal bursts are concatenated, ranked in their chronological order, in order to build the 'buccal S 10 signal' (blue in Fig. 2 B). Similarly, the 'lung S 10 signal' (green in Fig. 2

Analysis in frequency domain: continuous Wavelet Transform and frequency profile
Even in a stable condition, there is an intrinsic variability of this type of physiological recordings, to which it is pertinent to apply time-frequency analysis. Such an analysis is applied to a S W signal with Continuous Wavelet Transform (CWT) ( [1,7,10] ). For W 0 = 100 ms or 200 ms, S W 0 was first subsampled at 20 Hz (after applying an anti-aliasing low pass filter eliminating all the frequencies above 10 Hz), looking for frequencies between 0.1 and 10 Hz, with a complex Morlet wavelet basis function including five oscillations ( [5] ). The coefficients of each time-frequency amplitude map are computed for 100 values linearly distributed between 0.1 and 10 Hz, and at a time resolution corresponding to the sampling at 20 Hz. They are normalized between 0 and 1 by setting the maximal value at 1. Such a time-frequency map displays the temporal evolution of a signal S 200 dominant frequencies in color scale ( Fig. 3 B, for frequency domain [0.1 4] Hz). Since a neurogram is recorded in a given condition (baseline, hypercapnia or hypoxia), the CWT map reveals a stable band of the dominant fundamental frequency. Thus it is relevant to compute a profile in low frequency domain for a S 200 signal as the mean value over time, at each frequency, of the corresponding coefficients of the amplitude map ( Fig. 3 B, white curve). The prominent peak of this profile corresponds to the band of fundamental frequency, i.e the ventilation frequency. In order to allow a comparison between frequency profiles, independently of the signal durations and amplitudes, a normalization is applied to each frequency profile which consists in setting the value of the maximal amplitude of the profile at 1, i.e. at the dominant frequency. A similar CWT analysis can be performed on S 10 signals, in frequency domain [1,100] Hz. In this case, the subsampling of S 10 is performed at 200 Hz, and the time-frequency map is computed for 100 frequency values distributed between 1 and 100 Hz using a logarithmic grid, and at a time resolution corresponding to the sampling at 200 Hz. A profile in the low and high frequency domain is computed from the time-frequency map as previously described. Each profile is normalized by setting its integral at 1. Fig. 4 illustrates two examples of the frequency profiles: a profile in the low frequency domain (A) and a profile in the low and high frequency domain (B). The two profiles in (A), computed on S 200 buccal and lung L. catesbaianus signals, built from a same neurogram, indicate that a common fundamental frequency close to 1 Hz characterizes both signals. This common fundamental frequency suggests a coupling between the buccal and the lung oscillators, leading to a rhythmic regularity whatever the type, lung or buccal, of the successive bursts, as shown in figure Fig. 2 B. In Fig. 4 B, the frequency profiles in the frequency domain [1,100] Hz, computed on S 10 buccal signals of P. ridibundus pre-metamorphic tadpoles, exhibit, in both conditions normocapnia and hypercapnia, a narrow peak in the low frequency domain, corresponding to the fundamental ventilatory buccal rhythm, and a wide peak with a maximum around 25 Hz. This wide peak corresponds to intra-burst oscillations. There is no hypercapnia effect on the wide peak around 25 Hz, but the fundamental frequency is shifted from 1.6 Hz to 1.2 Hz ( [9,12,13] ).

Amplitude maps, amplitude profiles and oscillation profiles
The aim of this signal processing is to get a temporal alignment of the bursts by crosscorrelation, in order to detect a potential common time organization of the neuronal activity inside the bursts of a given neurogram.

Amplitude maps
In a given S 10 signal, a reference cycle is chosen, for instance the cycle that is the most correlated with all the other cycles (see Fig. 5 A a1, the red burst; the matlab function used for the computation of the correlation is "xcorr"). Then we position all the other cycles according to the maximal value of the crosscorrelation function (sliding dot product) between each cycle and the reference cycle ( [4] ). In a first step, the amplitude vectors representing the S 10 cycles are padded with zeros in order to get an identical length for all these vectors. In a second step, a symmetric matrix of crosscorrelation coefficients is computed, whose generic term of indices i and j is the maximal value of the crosscorrelation function between the amplitude vectors of cycle i and cycle j. The reference cycle is then the cycle whose mean correlation with all the other cycles is the highest: it can be considered as the cycle the most representative of all the cycles. At a third step, the time lag, which corresponds to the maximum of the crosscorrelation function between a cycle and the reference cycle, defines the time position of this cycle relatively to the reference cycle when building the amplitude matrix. In ( Fig. 5 A a2), a cycle amplitude is coded in gray scale. In ( Fig. 5 B b1) all the cycles are stacked building an amplitude map. When the cycles are horizontally positioned according to their time lag, the bursts are aligned ( Fig. 5 B b2). The cycles can then be ranked in the vertical stack according to their maximal amplitude ( Fig. 5 B b3). This representation of the cycles is the amplitude map. Another sorting rule can be chosen for the vertical ranking of the bursts in the amplitude map, for instance using the crosscorrelation between each cycle with the reference cycle.

Amplitude and oscillation profiles
The last step consists in representing all the buccal bursts (or lung bursts) of a S 10 signal, by a typical oscillatory pattern, i.e. an amplitude profile to which the low frequency component is eliminated. An amplitude profile is computed as the mean value of an amplitude map ( Fig. 5 C). Once a reference cycle is chosen, this amplitude profile is the same whatever the vertical sorting rule that defines the amplitude map ( Fig. 5 C c2 and c3). A filtered version of an amplitude profile (black dashed line in Fig. 5 C c3) is computed using a zero-phase, amplitude maintaining Savitzky-Golay filter with a polynomial function of degree 3 and a moving window of 201 points (i.e. 0.1 sec; the "sgolayfilt" function of matlab is used for this filtering process). The filtered profile is then subtracted from the amplitude profile, resulting in an oscillation profile (thick line in Fig. 5 C c3) that exhibits the high frequency components only. Each high amplitude peak of this profile corresponds to 'columns' of black points in the corresponding amplitude map (i.e. high amplitude peaks of the cycles). In Fig. 6 , 2 and 3, the amplitude maps are built using the crosscorrelation order, i.e. the cycles are stacked and ranked according to the decreasing crosscorrelation value with the reference cycle from bottom to top, as indicated by the plots of these values beside each map. The analysis illustrated in Fig. 6 allows to exhibit the shape modifications induced on the amplitude profile by metamorphosis and by hypercapnia in P. ridibundus tadpoles. Metamorphosis reduces both the buccal burst duration in each cycle and the number of oscillations inside each burst. In normocapnia, the burst shape in post-metamorphic tadpoles is more stable than in pre-metamorphic tadpoles, as shown by higher values of the crosscorrelation coefficients between the cycles and their respective reference cycles (see the blue curves in Fig. 6 , 2, A and B). At both stages, hypercapnia induces oscillations that are more marked than in normocapnia. The common feature of these profiles is the presence of oscillations that reveals a temporal grid. This grid structures the neuronal activity inside each burst at a frequency of about 25 Hz, i.e. in the range of the high frequency peak observed in Fig. 4 A.

Comparisons between oscillation profiles
A spectral analysis by Fast Fourier Transform (FFT) of the oscillation profiles illustrated in Fig. 7 A exhibits a same frequency range (18-28 Hz) for the dominant frequencies of each of them (using "fft" function of matlab), explaining the apparent similarity between these profiles computed on buccal S 10 signals of six different L. catesbaianus pre-metamorphic preparations. A careful comparison of these profiles using a crosscorrelation study shows that each oscillation profile has its own specific features. In order to evaluate how much a single oscillation profile is representative of the whole buccal burst population of a given neurogram, two oscillation profiles are computed from two distinct sets of the same number (half the total number) of buccal bursts randomly selected in the corresponding S 10 signal. Each cycle of the S 10 signal belongs either to one set or to the other set, but both profiles are based on the same reference cycle. The crosscorrelation coefficient, i.e. the maximal value of the crosscorrelation function obtained from these two oscillation profiles, is computed: it indicates the similarity degree between these two oscillation profiles. This crosscorrelation coefficient is called intra-individual crosscorrelation coefficient. Fig. 7 B illustrates the result of this computation for the six S 10 signals of pre-metamorphic L. catesbaianus tadpole preparations in normocapnia. The mean intra-individual crosscorrelation coefficient is 0.79. In order to evaluate the specificity of an oscillation profile for a given preparation, it is possible to compute an inter-individual crosscorrelation coefficient between two oscillation profiles from two respective S 10 signals. Each profile is now computed on the basis of all the buccal cycles of the corresponding signal. The matrix of such coefficients between all the pairs of the six oscillation profiles represented in Fig. 7 A, is illustrated in Fig. 7 C in color scale. The mean inter-individual crosscorrelation coefficient between two distinct pairs is 0.31, i.e. far lower than the mean intra-individual crosscorrelation coefficient. In the same way, in order to evaluate the similarity between the lung oscillation profile and the buccal oscillation profile for a given recording that exhibits both buccal and lung activities, we define the lung-buccal crosscorrelation coefficient. Each oscillation profile is computed using a common buccal reference cycle, from respective buccal or lung S 10 signal built as described in section 2.2 . Such pairs of superimposed buccal-lung profiles (buccal in blue, lung in green), centered and normalized in amplitude, are shown in Fig. 8 for L. catesbaianus post-metamorphic tadpoles and adult frogs, in control conditions and in hypoxia. For each of the 8 pairs of oscillation profiles, the buccal and lung profiles exhibit between 6 and 11 dominant peaks temporally coincident, suggesting a common time grid. The crosscorrelation coefficient value (in the range of 0.67-0.89) computed for each buccal-lung profile pair underlines the high degree of similarity between these two profiles, thus it confirms the presence of a common time grid on which both buccal and lung cycles are temporally organized.

Autocorrelogram maps and profiles
In order to evaluate the duty cycle shape stability/variability cycle by cycle, an S 10 signal is analyzed using the autocorrelation tool. Fig. 9 illustrates how this tool is applied to recordings from P. ridibundus tadpoles before and after metamorphosis, in control condition and in hypercapnia. The S 10 signal is segmented in successive pairs of two consecutive cycles, the last cycle of each pair is the first cycle of the next pair (red rectangle in Fig. 9 A). An autocorrelation function, or autocorrelogram, is computed for each pair. All these autocorrelograms are color coded ( Fig. 9 B) and stacked to build an autocorrelogram map as shown in Fig. 9 C and D. Such autocorrelogram maps exhibit the changes occurring with the stage (compare pre-metamorphic and post-metamorphic maps in Fig. 9 C and D, left and right columns respectively) and with the condition (normocapnia in Fig. 9 C and hypercapnia in D). Note that, for the pre-metamorphic tadpole, the color contrast is stronger in hypercapnia than in normocapnia, corresponding to a more stable duty cycle shape. It is the opposite for the post-metamorphic tadpole, where the color contrast is stronger, i.e. the shape stability is higher in normocapnia than in hypercapnia. From each autocorrelogram map, an autocorrelogram profile is computed as the mean value of all the lines of the map. The maximal possible value of the amplitude of the second peak of an autocorrelogram profile (black and red arrows in Fig. 9 E) is 0.5, which corresponds to a signal where all the cycles are identical. By comparing the value of the amplitude of this second peak with 0.5, we evaluate the cycle feature stability: the closer to 0.5 the higher the stability. In Fig. 9 E, the superimposition of the autocorrelogram profiles summarizes the condition effects on the cycle feature stability: for the pre-metamorphic tadpole, this stability is higher in hypercapnia than in normocapnia, while the situation is reverse for the post-metamorphic tadpole. In Fig. 9 F, the distributions of the second peak amplitude measured on each autocorrelogram show the differences between the two conditions: for the pre-metamorphic tadpole, the mean value of the distribution is 0 . 1 ± 0 . 05 in normocapnia and 0 . 26 ± 0 . 04 in hypercapnia, and for the post-

Conclusion
The frequency, amplitude, oscillation and autocorrelation profiles presented here, computed using time-frequency and crosscorrelation tools, are applied for neurograms recorded in amphibian brainstem preparations. These signal representations allow quantification of the short term and longterm stability/variability of rhythmic motor activities and of their similarity/dissimilarity. Numerous pseudo-periodic physiological signals can be represented using such profiles, for instance the mammalian breathing and locomotion.