{'fraction_non_alphanumeric': 0.02651668541432321, 'fraction_numerical': 0.010048743907011624, 'mean_word_length': 4.374929463925836, 'pattern_counts': {'":': 0, '<': 23, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 2, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Characterizing mixed mode oscillations shaped by noise and bifurcation structure

Many neuronal systems and models display a certain class of mixed mode oscillations (MMOs) consisting of periods of small amplitude oscillations interspersed with spikes. Various models with different underlying mechanisms have been proposed to generate this type of behavior. Stochastic versions of these models can produce similarly looking time series, often with noise-driven mechanisms different from those of the deterministic models. We present a suite of measures which, when applied to the time series, serves to distinguish models and classify routes to producing MMOs, such as noise-induced oscillations or delay bifurcation. By focusing on the subthreshold oscillations, we analyze the interspike interval density, trends in the amplitude and a coherence measure. We develop these measures on a biophysical model for stellate cells and a phenomenological FitzHugh-Nagumo-type model and apply them on related models. The analysis highlights the influence of model parameters and reset and return mechanisms in the context of a novel approach using noise level to distinguish model types and MMO mechanisms. Ultimately, we indicate how the suite of measures can be applied to experimental time series to reveal the underlying dynamical structure, while exploiting either the intrinsic noise of the system or tunable extrinsic noise.


I. INTRODUCTION
Mixed mode oscillations (MMOs), composed of welldefined periods of subthreshold oscillations (STOs) separating spikes, appear in biological rhythms, neural dynamics, chemical oscillations, and network oscillators, with many models proposed to capture this phenomenon in these and other applications 5 . One question in applications is the significance of the STOs -that is, do the interspike intervals (ISIs) and STOs encode important information for the system (e.g., Refs. 8 or 23). An element in the answer is identifying the mechanisms that drive the STOs. A challenge in this identification is that noisy time series generated by structurally different models can appear hardly distinguishable from each other, even quantitatively. Experimentally observed behavior can be captured with different types of models with different noise levels, so that model identification and calibration requires consideration of several classes of models over wide parameter ranges.
Previous studies indicate that a minimum of three degrees of freedom is necessary to produce MMOs in deterministic models (see examples in Ref. 5), while stochastic van der Pol and FitzHugh-Nagumo-type models illustrate how noise can drive MMOs in simple 2D models 26,27 . Similarities in MMO signals often stem from an underlying Hopf bifurcation and canard structure in the system 7,26,30 . In periods where the system is in some sense a) Electronic mail: peterphysik@gmail.com b) Electronic mail: rachel@math.ubc.ca close to the Hopf bifurcation, STOs can be both supported and disrupted by noise. In these models it is possible to choose parameters so that the MMOs are driven by the deterministic dynamics, while for other parameter choices the appearance of MMOs relies on stochastic fluctuations.
While on the one hand noise can make it difficult to distinguish between different models with various routes for producing stochastic MMOs, it is nevertheless possible to use noise as a way to identify underlying mechanisms. In this paper we compare a suite of measures at different noise levels to extract the structure of the model, and thus the mechanism for the STOs and MMOs.
The choice of measures is partly motivated by the presence of multiple time scales, an important feature that allows both well-defined periods of STOs and more than one mechanism for MMO generation. As described in more detail in the next section, the common substructure of the models we consider corresponds to a subsystem with a slowly varying control parameter. Depending on the combination of the other model parameters, that control parameter can slowly vary through a Hopf bifurcation, in which case the STO-spike transition exhibits the dynamics of bifurcation delay (Ref. 6 and references therein). For other parameter combinations corresponding to quiescence in the deterministic system, noise excites subthreshold coherent oscillations with a frequency near that of the Hopf bifurcation. The weak damping of these oscillations is due to the presence of multiple time scales, observed for parameters near a Hopf bifurcation. These noise-induced STOs are just one type of oscillation appearing in the context of coherence resonance (CR). In its broadest sense CR is observed when a range of noise levels excites coherent oscillations in a system that is quiescent without noise, with a maximum coherence at an optimal noise level. This phenomenon is most commonly cited for transitions between steady equilibria with large excursions, such as relaxation oscillations 22,28 , but in this paper it refers to the noise-induced STOs near a HB, also exhibiting CR-type behavior.
Whether the STOs appear via slow passage through a Hopf bifurcation (SPHB) or via a CR-type mechanism, the impact of the noise on the MMOs is concentrated in the interspike interval (ISI) where the STOs arise and eventually transition to a spike. It is no accident that this is also the interval where the slow dynamics are most prominent, and it is well known that noise can have a significant impact in such intervals 4,12,24 . The suite of measures used in this paper focuses on the STOs and the ISI behavior, capitalizing on previous analytical and computational studies of noise sensitivity of SPHB 4 and CR-driven oscillations 17,36 . These previous analyses show how noise can enhance or suppress the STOs through interactions with the fast-slow dynamics inherent in the MMOs, and these characteristics are captured by the behavior of the suite of measures. Focused on the noise sensitive STOs in the ISI, these measures can isolate differences between models and identify mechanisms for the MMOs as noise levels and bifurcation and control parameters vary.
An additional advantage of these measures is that they are based on features that are easy to analyze in time series data, making them amenable for use on experimental and simulated data. Furthermore, recognition of the multiple time scales of the MMOs suggests ways to experimentally introduce fluctuations into the system as a tunable extrinsic noise. When scaled appropriately by considering the ISI dynamics, this extrinsic noise can be used to mimic the effects of other intrinsic noise sources. Then this introduction of noise, combined with the suite of measures, can be applied both to identify the likely mechanism for MMOs and limit the variety of models one must consider for calibration.

Summary of results
We use a combination of amplitude trend, ISI density, and coherence measure based on power spectral density (PSD) to classify and differentiate mechanisms for MMOs. We compare phenomenological and physiological models to explore the types of STO and MMO behaviors that can be captured by simplified or reduced stochastic models over a range of noise levels. We then analyze these results within intra-model comparisons to get a deeper understanding of those MMO characteristics related to model structure, those dictated by deterministic dynamics, and those driven by noise. We focus particularly on cases where CR and/or SPHB are key mechanisms. In the conclusion we give a detailed list of identifying characteristics of MMO mechanisms obtained from this suite of measures. Understanding the dynamical features reveals the sensitivity in model calibration on bifurcation structure and on model choice, and illustrates ways that noise can be used to help to differentiate between models.
For STOs of the CR-type we observe longer tails in the ISI density together with stronger peaks in the coherence measure vs. noise, and no trend in the average STO amplitude leading up to the spike. In contrast, SPHB-dominated STOs have different behaviors in all three measures.
Simplified low dimensional integrate and fire (IF) models can be used to capture some STO and ISI behavior over a range of noise levels through adjustment of the speed of the control variable and reset. However, they can not mimic differences observed from underlying bifurcation structure, particularly in the case where a combination of CR and SPHB are at play.
Models with voltage-dependent control variables can have a rich variety of underlying deterministic behavior, allowing possibilities for both CR-and SPHB-driven MMOs. These differences can drive significantly different stochastic behavior in the ISI and PSD behavior as compared to simple IF models.
Related to this last observation, our analysis also illustrates how refractory dynamics, dynamics following the spike at the initial stage of the ISI, can influence ISI density for increasing noise, and disrupt the coherence of the STOs in MMOs. The influence of the reset or return mechanism, as well as the underlying deterministic behavior, observed in the measures considered here is consistent with the influence of these modeling aspects observed for noise-driven clusters of spikes, that is, repeated spikes without STOs in the ISI 19 .
The article is organized as follows: In Section II we first give a brief overview of two mathematical models in the literature -one physiological, the other phenomenological -that have been used to explain MMOs in neural systems, and discuss different mechanisms for MMOs in these stochastic systems. In the same section we present a number of different computational measures for analyzing important characteristics of the time series in order to identify the underlying model structure. In Section III we generate and classify MMOs with the physiological model and find matching MMOs generated by the phenomenological model. We then apply the measures to these time series and provide a set of features that -in combination -make the time series distinguishable. Secs. IV and V focus on the comparison of time series generated by the same model with varying parameters, to further highlight differences in the computational measures that appear with different model or bifurcation structure. We also highlight the dramatic effect of weak noise on different families of MMOs that are close in parameter space. In Sec. VI we outline how tunable extrinsic noise can be introduced to mimic the effect of intrinsic noise, based on the inherent multiple time scales and the fact that the noise impact is concentrated in the ISI. We also illustrate the applicability of our analysis to another MMOgenerating model.

II. MODELS AND MEASURES
We review the setup and dynamics of two models whose STO dynamics will be analyzed in this paper: A detailed biophysical model for MMOs in stellate cells ('SC' model) as well as a modified FitzHugh-Nagumo model ('MFN'). The two main MMO-generating mechanisms in these noisy models are slow passage through a Hopf bifurcation (SPHB) and coherence resonance (CR).

A. Biophysical model -SC
We consider a biophysical model for the MMOs observed experimentally in the layer II stellate cells of the medial entorhinal cortex introduced in Ref. 1. It consists of a voltage equation that includes the three usual Hodgkin-Huxley currents for sodium, potassium and a leak current together with a persistent sodium current and a hyperpolarisation-activated current. Together with the dynamical equations for six gating variables, this model is a system of seven coupled ordinary differential equations ('7DSC'), given in App. B.
Since we focus on the STO dynamics, for most of this paper we analyze a reduced, three-dimensional version (3DSC) of the full model, introduced in Ref. 30, consisting of only the equations for the voltage and two gating variables r f and r s (In the equation for r s , we use the term from the original work 1 (see footnote in Ref. 30).). This 3DSC model provides a good approximation of the ISI dynamics of the full 7D model 30 . Throughout this article, we refrain from giving units in the equations or for parameters (see App. B for the original units and parameter values). The three equations considered for the reduced model are: Eq. 1 contains a noise term η(t) representing white noise, corresponding to fluctuations in the persistent sodium current. According to Ref. 35, this current gives the main stochastic contribution in stellate cells, due to the low number of channels. The model (Eqs. 1-3) is treated as an integrate and fire (IF) model, where spiking and refractory dynamics of the full (7D) model are replaced by a reset value once the trajectory of the system crosses a threshold defined in terms of V . The value of reset is chosen such that the trajectory is placed close to a post-spike value of the full deterministic system (V = −80, r f = r s = 0 as in Ref. 30). The deterministic (D = 0) version of this reduced model has been analyzed in detail 30,31,34 , and there have been also studies of the noisy (D = 0) case 19 As discussed in Ref. 30, the 3DSC system can be viewed as a 2D subsystem in V -r f with slowly varying control parameter r s . In the nondimensionalized version of the 3DSC model analyzed in Ref. 31, the ratio of time scales for V and r f is approximately 0.02, with the time scale for r s even slower than that of r f by a factor of 0.3.
From that viewpoint, Fig. 1 shows the relevant part of the bifurcation diagram (with fixed I app ) of the underlying 2D subsystem with r s treated as the bifurcation parameter. For the parameters considered, there is a subcritical Hopf bifurcation at r s,H and a canard transition at r s,c (r s,H ≈ 0.08437, r s,c ≈ 0.08271 for I app = −2.45 and r s,H ≈ 0.09241, r s,c ≈ 0.09078 for I app = −2.58obtained with XPPAUT 10 ).
In the stochastic 3DSC model (D > 0), MMOs are observed for a broader range of I app values than listed above, generated by two mechanisms: a slow passage through a Hopf bifurcation and coherence resonance. For intermediate values of I app and noise, mixtures of both of these mechanisms are observed. One of the purposes of this paper is to provide measures that characterize and identify these different mechanisms in various settings. The top panel of Fig. 1 shows the stochastic dynamics of the 3DSC model when the deterministic dynamics corresponds to a SPHB, or delay bifurcation. The dependent variable V is attracted to a steady state for r s < r s,H . As the variable r s slowly ramps through r s,H , the full system does not immediately make the transition to spiking. Instead, there is a delayed transition, as STOs increase gradually about the unstable steady state. It is well known from previous studies that for delay bifurcations, or more generally for certain systems with alternating slow/fast dynamics, certain characteristics of the slow dynamics are exponentially sensitive to noise. That is, for noise levels above those exponentially small in terms of the parameter of the slow time scale, the behavior of the slow dynamics is qualitatively changed by the noise, as has been discussed for examples of delay bifurcations and resonance [2][3][4]15,27 and in other slow systems 12,21,24 . In the example in Fig. 1 the stochastic model typically has shorter intervals of STOs as compared with the deterministic model, as even small noise typically drives faster transitions to spiking.
The example in the bottom panel of Fig. 1 corresponds to STOs of the CR-type. The underlying deterministic system has a steady state corresponding to r s = r s,0 < r s,H , so STOs are damped. In the stochastic system coherence resonance (CR) occurs when noise amplifies the STOs in the presence of weak damping (see Refs. 11,17,20,36 and references therein; note that this type of CR is different to what is described, e.g., in Refs. 22 or 28). These STOs have a frequency corresponding to that of the Hopf bifurcation at r s,H of the 2D subsystem. Results in Ref. 36 indicate that CR drives STOs in the range of noise levels where coherence is optimal, that is, where the PSD is relatively narrow with significant power. The analysis in Ref. 36 shows that the amplification factor of the STOs is inversely proportional to the proximity to the Hopf bifurcation, so that the STOs are a prominent and prolonged feature in the ISI for r s near r s,H . Variability in the amplitude of these STOs yields a nontrivial probability of spiking, leading to frequent appearance of MMOs, even for I app < I app,H if noise is strong enough 30 .

B. Modified FitzHugh-Nagumo model -MFN
In Ref. 26, a FitzHugh-Nagumo-type (FN) model is given by the following pair of coupled nonlinear stochastic differential equations: We follow Ref. 26 and use ǫ = 0.005, a = 0.9 and g(x) = 7x 2 + 0.08(1 − e −x/0.08 ), chosen to give a ratio of the time scales of the STOs and the spikes (relaxation oscillations) that is close to that typically observed experimentally and in the SC model. The parameter b is a control or bifurcation parameter. We refer to this model as the modified FN model (MFN). Taking a = −1, D = 0, adding a constant term in Eq. 4 and a v-dependent term in Eq. 5 yields the general FN model 16 . With a = −1, g(x) = x, b = D = 0 and a simple rescaling Eqs. 4 and 5 are the original van der Pol (vdP) model 14 . Fig. 2 shows the corresponding bifurcation diagram: a supercritical Hopf bifurcation with a canard transition to The influence of noise on the MMO dynamics of the vdP model has been studied recently in detail 27 , providing scaling relationships between the proximity to the Hopf value and the noise levels that lead to different types of behavior. Here we use the model given by Eqs. 4 and 5 as a basis for comparison with the biophysical model for two reasons. First, the use of g(x) gives control over the time scale and shape of the spike (see above). The other important effect that the nonlinearity has on the system is a shift of the canard point, enlarging the parameter region of stable STOs (i.e., b c −b H ) by a factor of roughly 1.5. This leads to MMOs with well-defined STOs over significant parameter ranges even in the presence of stronger noise values, in contrast to the vdP model studied in Ref. 27. One other difference, compared to the vdP model, is that the choice of g(x) yields a less abrupt transition from STOs to relaxation oscillations, allowing for a greater probability of observing spikes with reduced amplitude, an effect rarely observed in the SC model.
As in Ref. 19 we augment the MFN model (Eqs. 4 and 5) with a slow variation in the control parameter given by an additional dynamical equation for b. The presence of a slowly varying control parameter is a feature common to a number of MMO models including the SC model, and is necessary to reproduce some of the features observed in IF models. Additional comments about the model of Eqs. 4 and 5 relative to the augmented models are included in Subsec. V A.
The first augmented model we consider is an IF model with a linear ramp for b plus reset at a fixed threshold value for u:ḃ Once the trajectory surpasses u th , b is reset to b rs and we choose to reset u and v such that the system is on the nullcline (u rs = b rs and v rs = b rs (b rs − a)(1 − b rs )) or very close to it (see Subsec. V A for a discussion of the dynamics when reset values are off of the nullcline). This system, which we reference as the linearly augmented modified FN model (LMFN), generates only MMOs within the parameter values considered. The stochastic LMFN model exhibits features of both a SPHB and CR, depending on the value of ǫ 2 as well as where the trajectory is reset. This allows for a comparison with simple IF models, highlighting the effect of noise on the dynamics for different ramp speeds of the control parameter b and reset values.
The second augmentation of the MFN model considered is a nonlinear u-dependent equation for b(u) with a functional form similar to that in the SC model:

C. Measures
For the analysis of the MMOs and specifically the STOs in the ISIs, we use a number of measures: ISI density, amplitude of the STOs preceding a spike, power spectral distribution (PSD) of the STOs, and a coherence measure. The time series were generated by numerically integrating the coupled differential equations given in the preceeding subsections using a simple Euler-forward routine with constant time step. For the output, it was ensured that the data was sampled at a high enough rate 1/∆t, such that the STOs could clearly be identified (40∆t T STO ).
The density of the interspike intervals (ISI) was obtained in the following way: In IF-type models (3DSC and LMFN), ISI duration is measured from the reset to the crossing of the threshold corresponding to initiation of the spike. The typically long refractory period without STOs in the 3DSC model was removed. In spiking models (7DSC, NLMFN), it is measured as the time between two consecutive crossings of a threshold (i.e., the time between two spikes including a spike).
To compare the different models (with different time scales), we measure the ISI duration in T STO , the approximate mean STO period. For simplicity, we use only one T STO for each of the model versions and obtain these approximately from one time series. We use the following approximate values for T STO : 3DSC: 107.5, NLMFN: 0.47, LMFN: 0.45, but note that the actual average T STO vary by about ±5% among the classes introduced in the next section. Average ISI duration is used as a coarse measurement to calibrate FN-type models to the biophysical (SC) model, and the ISI density is used to compare models and parameter choices, including influence of the noise.
A second characteristic of the STOs used to calibrate FN-type models to the biophysical model is average STO amplitude trend before a spike. For this, the maxima of the low-pass filtered time series were extracted and averaged over many ISIs (for details see App. C).
To obtain the power spectral distribution (PSD) of the STOs, the spikes were removed from the time series and it was ensured that the normalization of the PSD is consistent within each model. From the PSD, we compute a measure of coherence, namely 11 where f p is the frequency of the peak in the PSD corresponding to the STOs, h p is its associated maximal power and ∆f is its full width at half maximum. Since the different models use different units, a comparison of the absolute value of β across models was not performed. For details regarding the computation of the PSD and β, see App. C.

III. ANALYSIS OF MMO CHARACTERISTICS: INTER-MODEL COMPARISON
In this section we present an analysis of the characteristics of MMOs with computational measures, focusing on the STOs, produced by the different models as introduced in Section II. STO amplitude trend before a spike and average ISI provide a coarse classification of the MMOs in the subsections below. For each class of behavior the parameters of the 3DSC model are chosen to produce a type of STO observed in the 7DSC model, and the parameters for the FN-type models are chosen such that they reproduce the mean ISI and amplitude trend similar to the 3DSC. The results illustrate features of the MMOs that can be reproduced by both the biophysical and phenomenological models. The measures introduced in Subsec. II C provide a finer comparison of the three broad classes. We introduce different noise levels to highlight the mathematical elements of both the underlying deterministic model and stochasticity that influence the MMO characteristics.
A. Class 1: Short ISI, increasing STO amplitude In Class 1, we analyze MMOs with short ISIs during which the amplitude of the STOs increases. Examples for this class of MMOs are shown in Fig. 1 (top panel) for the 3DSC model and Fig. 3 for the LMFN model. The comparison of the amplitude dynamics of the STO and ISI density is shown in Fig. 4, for 3DSC and LMFN, where the ISI density and mean agree up to a shift due to different reset conditions. This Class 1 behavior is not readily reproducible within the NLMFN model where the choices of c b needed for short ISIs lead to STOs with a larger amplitude throughout the ISI.
As shown in Sec. II, the underlying dynamics for both the 3DSC and the LMFN model corresponds to a slow passage through a Hopf bifurcation. In the SC model, the increasing amplitude of the STOs are generated when the system spirals away from an unstable fixed point for r s values beyond the subcritical HB. In the LMFN model, the trajectory passes through the underlying supercritical HB at b H followed by an increasing STO amplitude, with escape from STOs to a spike appearing for relatively large b (around 0.34) beyond the canard point b c . The significance of the speed of variation of b in the LMFN model is discussed further in Subsec. V A.
The noise levels shown in Figs. 1 (top panel) and 3 are among the lowest at which stochastic effects are observed, with the STO dynamics dominated by the deterministic trend of increasing amplitude. Increasing the noise strengths has very similar effects in both models.
• Typical for a SPHB, increased noise drives an earlier escape, indicated by shorter ISIs (Fig. 6) and lower average amplitude before escape (Fig. 5). Differences due to the nature of the underlying Hopf bifurcation (sub-/supercritical) are evident only in the amount of reduction of the amplitude of the STOs. This reduction is more significant in the 3DSC, where the STOs are unstable for the underlying 2D system with r s > r s,H , as compared to the LMFN model, where the supercritical Hopf bifurcation yields stable STOs for b > b H .
• Increased noise shifts the ISI density to include peaks at shorter durations, consistent with earlier escape described above. Contributions at longer ISIs are due to a small fraction of trajectories for which the noise disrupts the transition and lengthens the ISI. Class 1 behavior survives only for low to moderate noise levels, as larger noise regularly eliminates the STOs. This is illustrated by the considerable spread of both the ISI density and the PSD (Fig. 7) for large noise levels.
• The coherence measure β is not clearly defined for MMOs of Class 1, since the strong amplitude trend of the deterministic dynamics within the short ISI periods yields PSDs with many well-defined peaks (one of them at the STO frequency -see Fig. 7). The SC model shows more peaks due to larger amplitudes of the STOs across the ISI. The strong peak at low frequencies comes from the spiking (i.e., thresholding and concatenation).  (bottom panel) and 8 for the 3DSC, NLMFN, and LMFN, respectively. These time series were calibrated for similar average amplitude trend and average ISI durations. Fig. 9 shows that for low noise levels, it is possible to find parameter values for LMFN that compare well in amplitude behavior with 3DSC for both Class 2A and Class 2B. For NLMFN, the amplitude behavior is similar to 3DSC for Class 2B but differs in Class 2A for STOs well before the spike. In the refractory period following the spike in the NLMFN the system relaxes near but not on the left null cline with b < b H , yielding slowly damped STOs in the first part of the ISI. For the LMFN, the reset value for u and v is chosen very close to the fixed point, avoiding this part of STOs with decaying amplitude. Comparable average ISIs are easily obtained in the LMFN with ǫ 2 chosen an order of magnitude smaller than in Class 1, and for NLMFN with c b near unity for weak noise. Differences in the ISI density are apparent from Fig. 10 and are related to differences in the deterministic dynamics as discussed in Secs. IV and V.
As we see below, by tuning parameters and noise levels, the FN-type model can capture some but not all of the amplitude, ISI and coherence behavior for Class 2. Results for different noise levels point to differences in the underlying structure of the models.
• For noise levels increased by an order of magnitude, both subclasses for 3DSC show roughly a 20-30% decrease in the average amplitude of the largest STO immediately preceding a spike, while the FN-type models show a slightly greater decrease (Fig. 11). For Class 2B, the entire average amplitude curve of 3DSC shifts to lower   levels with increased noise, while both FN-type models show both shifts and changes in the average amplitude curve shape. For Class 2A, both 3DSC and LMFN show an increase of the average amplitude of STOs well before the spike for increased noise, indicating STOs driven both by CR and by SPHB. For the NLMFN model, the return sets b < b H following an earlier spike, so that the combination of damped STOs for b < b H and CR effects yields robust STOs with less variation with noise level than for for the LMFN model.
• In Class 2B, the 3DSC model has longer tails in the ISI density for both small and intermediate noise levels. These longer tails are not seen in the FN-type models due to differences in the underlying deterministic dynamics (see Sec. IV). For weak noise, the multi-modal ISI density for the NLMFN is due to the stronger attraction to the underlying deterministic MMO or STO dynamics, with larger STOs initially in the ISI. In the LMFN model the reset onto the nullcline avoids damped STOs initially in the ISI, allowing a unimodal ISI density. For both Class 2A and 2B, 3DSC and LMFN show a shift to shorter ISIs and a similar shape of the ISI density for larger noise, with a larger variance for LMFN in Class 2A. For NLMFN there is no shift in the ISI density with increased noise, just an increased variance both for Class 2A and 2B, reflecting the dynamic return mechanism following the spike. • Fig. 13 shows the differences between the coherence measure β (Subsec. II C) for the three different models for various noise levels. For 3DSC in Class 2B the STOs are driven by CR and there is a peak in β (indicating CR), while for Class 2A noise reduces β in general, since the STOs are driven by a SPHB. For low noise levels in NLMFN, the robust STOs are represented by a multipeaked PSD (equivalent to the ISI density -see above), for which β is not defined. For larger noise levels in Class 2, β is defined, but noise only disrupts the underlying STOs, so β decreases with noise. For LMFN, a moderate noise level can enhance the STOs, particularly at the beginning of an ISI, yielding an optimal noise level for coherence in both types of Class 2. The influence of the underlying bifurcation structure on the coherence measure is discussed further in Sec. IV. Here we analyze STOs with an average ISI similar in length to Class 2, but with a constant average amplitude preceding a spike. Fig. 14 shows the trajectories from the 3DSC, NLMFN, and LMFN calibrated for similar average amplitude behavior and average ISI duration. Class 3 can be observed in 3DSC only at higher noise levels compared with Classes 1 and 2. Fig. 15 shows that the average STO amplitude for Class 3 is constant for at least the last ten periods of STOs with a slight amplitude increase in transition to the spike. This nearly constant average amplitude results from strong variation in amplitude with a constant ensemble average. A clear single peak in the PSD (data not shown) confirms well-defined STOs with variability in phase and amplitude.
The amplitude and average ISI characteristics of these MMOs are readily reproduced with the FN-type models at stronger noise levels with the average amplitude lower by roughly a factor of two compared to Classes 1 and 2. Class 3 behavior can also be obtained for the MFN model (Eqs. 4 and 5) with constant control parameter b as in Ref. 26. The ISI densities shown in Fig. 15 are similar to each other, with a slightly more significant tail in the ISIs for 3DSC. Figs. 16 and 17 show the trend of the STO amplitude and the density of ISI when noise levels increase.
• For stronger noise the increase in amplitude directly before the spike in the FN-type models is eliminated, with stronger noise levels rather than deterministic dynamics dominating the transition to spikes. At higher noise levels, distinguishing stochastic fluctuations from regular STOs is less dependable, so that the amplitude dynamics for the largest noise levels are not shown in Fig. 16. In the 3DSC model, there is a slight increase in average amplitude during the ISI which is due to our algorithm (Subsec. II C) not fully eliminating the trend of increasing V in steady state with increasing r s for r s < r s,H . • The ISI densities for 3DSC and LMFN are more concentrated at shorter ISI durations, while the ISI density for NLMFN remains spread over a large range of ISIs. This is due to the underlying stable STOs for NLMFN, together with the nonlinear return mechanism that typically returns b to values well below b H following a spike, thus allowing for a range of ISI durations. With the ISI density concentrated at low values, clustered spikes are more frequent in 3DSC and LMFN with nearly tonic spiking at the higher noise level (cf. Ref. 19).
• For Class 3, β for the 3DSC model shows a similar behavior as for Class 2B (Fig. 13), with a clear maximum. The values of β in Class 3 are lower than in Class 2 for 3DSC, mainly because of smaller amplitudes of the STOs (being inversely proportional to the distance from r s,H for CR 36 -also see Subsec. II C). Similarly, the trend for β in Class 2B and 3 is the same in LMFN. For the NLMFN model, the value of c b is the same for Class 2B and Class 3 and the behavior of β is shown in Fig. 13.

IV. INTRA-MODEL COMPARISONS FOR STELLATE CELLS
The comparisons in the previous section point to the role of the underlying bifurcation structure on the characteristics of the STOs and MMOs. We summarize the results for the 3DSC model from Sec. III, comparing within that model. We also discuss alternative values of I app for the 7DSC model in order to compare with previous studies of the full model.

A. Comparison for different ranges of Iapp
Within the 3DSC model (and similarly the 7DSC model) differences in the stochastic behavior can be related to the differences in the deterministic behavior for different ranges of I app . For I app < I app,H ≈ −2.575 in 3DSC (see Subsec. II A), the stable steady state has a value of r s = r s,0 well below r s,H . For I app > I app,H there is no steady state for the deterministic system, and the STOs are driven by SPHB with r s as the slowly varying control parameter. In the case I app ≈ I app,H , characteristics of both CR and SPHB are observed.
• Amplitude: For I app well below I app,H , r s,0 is well below r s,H , so only high noise levels can drive MMOs, and STOs are purely of the CR-type without any trend in the average amplitude (as in Class 3). In contrast in Class 1, with I app clearly above I app,H the STOs due to the SPHB have an increasing trend in amplitude before the spike. Well-defined periods of STOs survive only for lower noise levels, as the slow passage is sensitive to very weak noise. For values of I app closer to I app,H , as in Class 2, both CR and SPHB can affect the dynamics, so that periods of STOs with increasing amplitude trend can survive over a large range of noise levels.
• ISI density: Significant differences are seen at lower noise levels, where longer tails in the ISI density correspond to STOs of the CR-type, not observed in STOs driven by SPHB.
• Coherence measure β: CR-driven STOs in Class 2B and Class 3 are responsible for the more pronounced peak in the coherence measure β in Fig. 13.
Intermediate values of I app = −2.5, −2.65 were also analyzed (not shown), and as expected intermediate behavior between the three classes is observed. For I app = −2.5, STO and ISI behavior was closer to that of Class 2B than Class 1, without contributions to longer ISIs, and for intermediate to larger noise levels very weak coherence is observed with a mild peak in β around D ≈ 10 −5 .

B. Families of MMOs for Iapp > Iapp,H at weak noise
Here we consider the effect of noise on MMO families in the full 7DSC model in the range I app > I 7D app,H (i.e., within the regime of deterministic MMOs), as illustrated through the ISI density. The 3DSC model shows a similar behavior, with a shift in the values of I app . Recall that for any fixed I app > I 7D app,H ≈ −2.702, the underlying deterministic dynamics is a slow passage of r s through r s,H . We restrict our study to very low noise levels and consider trajectories that fit into Classes 1 and 2A (Sec. III) when stronger noise is added. For the sake of comparison with the analysis in Ref. 34, we use an alternative equation for the dynamics of r s . Instead of Eq. 3, we use Eq. B9 given in App. B, which alters the numerical values of the model but not the qualitative features. Fig. 18 shows examples of the ISI density for different families of MMOs. The 1 s families (for notation, see Ref. 34) have MMOs with s STOs in each ISI, yielding ISI densities with a single point mass. The MMOs for a second, more complex family, J i (J j ) l , have an ISI with i STOs followed by l ISIs each with a different type of j STOs that increase in amplitude and period for each of the subsequent l ISIs. The J i (J j ) l families then have 1 + l point masses in the ISI density. We consider only those MMO families that are stable in the deterministic model. In the range −2.55 < I app < −2.4 stable 1 s and J i (J j ) l families are found for alternating ranges of I app . For I app −2.55, the stable solutions are only 1 s families densely filling the space in I app 34 . The time series analyzed earlier as Class 1 in this notation would be a mixture of 1 2 and 1 3 .
For values of I app near the center of the stability range for a particular type of MMO, well separated from other families, weaker noise broadens the ISI density and stronger noise shifts the density towards shorter ISIs when noise drives an early escape typical in a slow passage through a HB, as shown in Fig. 18 for J 0 (J 1 ) 1 . The bottom panel of Fig. 18 shows that for values of I app near the edge of a stability range, very weak noise can lead to a mixture of close families of MMOs. As noise drives a faster escape to spiking the peak corresponding to the longest ISI and largest STO of J 0 (J 1 ) 4 is reduced as the nearby MMOs of J 0 (J 1 ) 3 occur, appearing as stronger nearby peaks in the ISI density. There is also considerable sensitivity in the ISI density when stability ranges are small and densely packed. Noise can drive a sampling of a mixture of families close in terms of the bifurcation parameter, yielding clearly defined peaks that are not present in the deterministic case, as shown in Fig. 19. In the top panel, the deterministic MMO is 1 1 , and noise excites a second peak for a family 1 2 or a mixture of families between 1 1 and 1 2 . At I app = −2.66 the stable solution is 1 12 and weak noise drives a sampling of the six nearby MMO solutions 1 s with 10 ≤ s ≤ 15. At stronger noise the longer ISIs are lost and the ISI density spreads out and shifts towards shorter ISIs.
In conclusion, even very weak noise can alter the dynamics of a MMO-generating system significantly, partic-ularly if the deterministic solutions are close in parameter space. Weak noise drives a sampling between a few or many of these solutions. Then it can be difficult to distinguish between complex deterministic solutions like the J i (J j ) l families and a noisy trajectory sampling a few different MMO families.

V. INTRA-MODEL COMPARISONS FOR THE MODIFIED FN MODEL
As in the SC model, the underlying bifurcation structure and deterministic behavior influence the source and characteristics of the STOs. For the FN-type models, the underlying bifurcation is a supercritical Hopf, with stable STOs in the 2D reduced system for u-v with constant b H < b < b c .

A. Linearly augmented modified FN model
In the LMFN model with appropriate reset the slowly varying control variable sweeps through an underlying HB and canard transition, which always yields MMOs in the deterministic setting. We consider two aspects of the control parameter that can be varied to change the characteristics of the STOs, the rate of change in b,ḃ = ǫ 2 , and its reset b rs , following a spike/threshold crossing.
In the different classes studied in Sec. III, we analyzed four different combinations of ǫ 2 and b rs : slower variation in b given by smaller ǫ 2 = 0.00055, with reset b rs < b H (Class 3) or b rs ≈ b H (Class 2B), and larger values of ǫ 2 = 0.001 (Class 2A) and ǫ 2 = 0.0147 (Class 1) both with b rs ≈ b H .

Characteristics for ǫ2, brs in Section III
The main underlying feature that influences amplitude dynamics for these different parameter combinations, is time spent in the parameter range b > b H . • Amplitude increase is more pronounced for low noise levels and smaller values of ǫ 2 . For low noise levels, the system follows the underlying STO dynamics, less likely to be driven to spiking by noise. For smaller values of ǫ 2 , the control parameter behaves almost as a constant relative to the oscillation frequency. This allows attraction to larger STOs for slowly increasing b H < b < b c , with the possibility of reaching values of b > b c before spiking, so that larger amplitudes are typical before escape. For larger values of ǫ 2 , as in Class 1 or Class 2A, b increases through the STO region less slowly. Then the trend of increasing amplitude consistent for b H < b < b c survives but increased variation in b limits the attraction to the stable STOs as b approaches and exceed b c . The behavior is closer to a series of transients that is more susceptible to noise-induced escape to spiking. For larger noise levels, the escape to spiking is noise-driven in or near the stable STO region. An increase in STO amplitude farther from the spike, for b ≈ b H is observed due to the role of CR together with the stable STOs. An earlier escape limits the opportunity for amplitude increase closer to the spike.
• The relation between ǫ 2 and b rs and ISI density is as expected: smaller values of ǫ 2 or b rs lead to greater likelihood of longer ISIs.
• The coherence measure β shown in Fig. 13 increases for long stretches of large amplitude oscillations, and decreases with fluctuating amplitude. Trajectories with the reset close to the HB b rs ≈ b H show the largest value of β for moderate noise levels. Slower variation in b in these cases also reduces fluctuation in amplitude, thus increasing coherence. For reset b rs < b H there is fluctuation in amplitude due to oscillations driven by both CR and STOs, so β is reduced somewhat overall. For larger values of noise coherence is destroyed in all cases. Fig. 20 shows the analysis of an additional parameter combination, a very slow sweep through the HB and reset before the HB (ǫ 2 = 0.0001, b rs = 0.305), whose coherence measure β was included in Fig. 13.

Characteristics for slower increase in b with brs < bH
• The difference in amplitude dynamics for lower and higher noise levels is as described above. For lower noise levels and very slow dynamics of b, the system is attracted to larger amplitude STOs before spiking, so that larger amplitudes are typical for more STOs before escape.
• The ISIs are very long, and only for stronger noise are they comparable to Class 3 (cf. Fig. 15). In contrast to Class 2, only for stronger noise levels is the ISI density concentrated at shorter ISIs, since the trajectory takes longer to reach b H .
• The shape of the curve β vs. D is qualitatively similar to the classes studied in Sec. III. However, at very low noise values, β is larger than other cases, caused by CR-driven STOs and slower variation in b as described above.
The last panel in Fig. 20 shows the ISIs for an intermediate case, with larger ǫ 2 and b rs < b H . The ISI density then shows characteristics between Class 2 and Class 3, with a narrow ISI density for smaller noise levels, and longer tails for stronger noise, as oscillations of the CRtype appear for b < b H .
We add one additional note about resetting off of the nullcline. The choice of reset used in Sec. III (reset very close to the nullcline), can be relaxed to a reset near the nullcline without qualitative changes in the measures used above. A choice of reset away from the nullcline modifies the amplitude dynamics of the STO at the beginning of the ISI. For reset with b rs well below b H , STOs with considerable amplitude follow the reset and relax as b approaches b H , similar to NLMFN. For Classes 1 and 2 this increased amplitude supports an earlier escape to spiking and reduced ISIs, particularly in Class 1. For low noise levels in Class 2 the ISIs are then narrower, since the larger STOs are less susceptible to small noise. For Class 3 the stronger noise dominates, so that variation in the reset has limited effect.
As noted in Subsec. II B, the model of Eqs. 4 and 5 has a larger distance between b H and b c , as compared with vdP, yielding a larger region in parameter space where well-defined periods of STOs are observed in the presence of noise. This behavior is confirmed by comparisons of PSDs for MFN and for a rescaled version of vdP to effectively increase the canard value. This comparison illustrates that the distance between b H and b c plays a critical role in the robustness of STOs and MMOs, in addition to the reset and ramp speed. that support phenomenological behavior observed in the SC model. Here we also discuss cases that differ qualitatively from those studied in Sec. III, and compare the amplitude, ISI density, and coherence measure within this model. • Increasing amplitude before the spike over a range of noise levels.
• Multi-peaked ISI density for lower noise levels and concentrated ISI density for larger noise.
• Decreasing trend in β, defined only for larger values of the noise.
For values of c b < 0.95 the return value of b following a spike takes values between b H and b c , thus shortening the ISI significantly. As discussed above, we do not classify this behavior as Class 1, since the return value yields larger amplitude STOs and limited increase for low noise levels, similar to those observed for a related model in Subsec. VI B.
In contrast, for c b = 1.1, the deterministic system exhibits stable STOs where b oscillates near b c . The stochastic behavior is similar to that of the case with c b = 0.95, but differences due to the underlying deterministic dynamics are observed for lower noise levels and c b = 1.1: • Increased noise drives an earlier escape to spiking. This results in a reset of b < b H , yielding both CR-type and SPHB-driven oscillations, and a limited opportunity for amplitude increase just before the spike.
• The ISI density is broader for c b = 1.1 (compared to c b = 0.95).

Intermediate values
The attracting deterministic behavior for b consists of small oscillations midway between b H and b c , so characteristics include some elements from Class 2 for smaller noise, and similarities to Class 3 and c b = 1.5 (see below) for larger noise: • The amplitude dynamics is similar to Class 2 for lower noise, with a weaker increase. • The ISI behavior has characteristics closer to Class 2.
• The coherence measure β vs. D behaves similarly to c b = 1.5, with larger values of β for small noise. For larger values of c b , the attracting state of the deterministic system takes values of b near b H . For c b = 1.6 the long term deterministic behavior is quiescent with the steady state corresponding to b just below b H . For c b = 1.5 the deterministic system has stable STOs, with small oscillations in b just above b H . To drive MMOs in either of these cases, noise must be strong enough to cause excursions beyond the canard transition at b c , certainly stronger than in the cases for c b ≈ 1 above. ISI and amplitude measures can therefore only be derived for these stronger values of noise, leading only to Class 3 MMOs with the following characteristics (Fig. 21): • Primarily constant trend in average amplitude before a spike.
• Broad ISI density with long tails.
• A clear difference between c b = 1.5 and 1.6 due to the underlying deterministic dynamics is seen in the coherence measure β (Fig. 13). For c b = 1.5 all noise levels disrupt the coherence of the STOs produced deterministically, so that β strictly decreases with noise. For c b = 1.6, β has a maximum indicating an optimal noise level typical of CR-driven oscillations 28,36 .

VI. APPLICATION TO OTHER MODELS
A. The 3DSC model with noise in Iapp The stochastic nature of the SC system results from noise in ion channels and noisy external signals. The model we consider above follows Ref. 35, including the primary noise source in the gating variable of the persistent sodium channel. Here, we compare our results with those obtained with a different noise source in the 3DSC model, namely, in the applied current I app . The consideration of this additive noise is motivated by several factors. First, systems with multiple time scales often can be divided into regions of time or space where noise plays a more or less significant role. This raises questions about where or when noise plays a critical role, and whether the type of noise can make a significant difference. Second, since I app is the only variable in the system that can be easily experimentally controlled, this raises the question of whether varying I app can be used as a controllable input to probe the dynamics. Providing noisy input is a technique used in experimental neuroscience (e.g., Ref. 25) to investigate the structure of the complex dynamics. Finally, understanding whether different types of noise sources have similar or different effects contributes both to the understanding of the mechanism for the dynamics and to model identification.
To test whether fluctuations in I app can reproduce dynamics similar to those generated by noise in the persistent sodium current, we consider the voltage equation in the interval near the end of the ISI, when the dynamics of V is slow. This is the stage at which STOs are generated for parameters near the Hopf bifurcation of the reduced system. In that part of the cycle, V can be approximated by a constant V 0 in Eq. 1. This approximation suggests that even though the noise is parametric, at the end of the ISI when the STOs are generated, the noise behaves as additive noise to leading order with a coefficient D ′ ≈ (0.15G p (V 0 − E Na )) 2 D and Eq. 1 can be approximated bẏ Eq. 9 can be interpreted as having a noisy applied current with I app + √ 2D ′ η(t) and D ′ ≈ 68D (for V 0 ≈ −55) instead of the original noisy gating variable. Then one would expect that the effects of different noise levels on the STOs as described in Sec. III could be reproduced by varying the noise level in I app , as long as the results are scaled in a manner consistent with Eq. 9. Indeed for the ranges of noise levels considered in this paper, the results for stochastic I app are essentially the same as described in the classification of Sec. III. With the above scaling, the plots for the coherence measure β (as in the top panel of Fig. 13) overlap almost perfectly. We also reproduced similar behaviors of the ISI density and amplitude dynamics as analyzed for the 3DSC model in Sec. III (data not shown).

B. A self-coupled FN-type model
Here, we analyze another MMO-generating model, a self-coupled modified FN model derived and used for the study of coupled Hodgkin-Huxley neurons 7,33 . Similar to the MFN model presented earlier, v represents a voltage variable, h is gating variable, and s is a dynamic coupling between neurons. Compared to Refs. 7 and 33, we add a noise term in the first differential equation: H(v) is the Heaviside function. This system of differential equations has multiple time scales similar to the models analyzed so far: the dynamic variable s plays the role of a slowly varying control parameter, decreasing and passing through a HB at s H for the underlying v-h system. Also, the spikes are included in the model (not just a threshold and reset) but this is done using a sharp switch (H(v)) rather than a smooth nonlinear function for the control variable as in our NLMFN. in the NLMFN case with c b = 1.5. Another noticeable difference between Eqs. 10-12 and the MFN models is that the dynamics of s yields a return value well below the Hopf point s H following a spike. After the return there is a significant interval during which s approaches s H , which has implications for both the amplitude and ISI behavior. Class 1 time series are obtained with δ = 2.5 and D = 10 −7 , yielding both the amplitude dynamics and the ISI density similar to those shown in Fig. 4. One difference for the model of Eqs. 10-12 is that increasing noise strength typically increases the amplitude of the STOs before a spike. This is consistent both with the slowing of s near s H and with the opportunity for CRtype STOs in the interval where s > s H . As in Class 1, we see the appearance of new well-defined peaks in the ISI density with low noise levels. For this particular choice of parameters, the first of these new peaks appears at longer ISIs, similar to the behavior in the upper panel of Fig. 19 for the 7DSC model.
For time series with the characteristics of Class 2, we find δ = 0.3 and D = 3 × 10 −6 in Eqs. 10-12 yielding dynamics comparable to Sec. III in terms of average ISI and amplitude trend (Fig. 23). For the STOs well before the spike, an increase of the amplitude results from the algorithm that subtracts the average v as v has a decreasing trend. The amplitude before the spike increases with noise level, once again due to the combination of CR-driven STOs before s crosses s H , and a slowing of s. Larger values of D shift the ISI density towards smaller values and broadens it, yet a certain minimal ISI is maintained due to the return of s well beyond s H . The effect of such a return mechanism was also observed in the ISI density in the NLMFN.
Reduced speed of s (δ = 0.17) combined with stronger noise (D = 3 × 10 −4 ) leads to time series that fit into our Class 3. The amplitude dynamics is very similar, but the ISI density is different from those shown in Fig. 15. At this noise level, the ISI density is still centered around 15 T STO , again due to the return of s well above s H . Only for the highest noise levels considered (D = 3 × 10 −3 ) the noise drives an early escape, reflected in the ISI density concentrated at shorter ISI lengths.
In Fig. 24 we show the coherence measure β for the two parameter sets used for generating Class 2 and Class 3 time series. As with the earlier analyses, obtaining β for Class 1 time series is not possible due to a multipeaked PSD. For Classes 2 and 3, we neglect a strong peak in the PSD resulting from the concatenation and compute β as in previous examples. There is a weak peak in β at intermediate noise levels, consistent with previous examples that exhibit MMOs related to both CR and SPHB. In Eqs. 10-12 this combination of STO mechanisms results from the return of s and its slower variation near the HB.

VII. CONCLUSION
Attention towards oscillatory neural dynamics has expanded to include coherent subthreshold oscillatory activity in the voltage of some neuron classes as well as spikes 9 . Combined, the sequence of subthreshold oscil-lations (STOs) followed by a spike comprises a class of mixed-mode oscillations MMOs. As these MMOs have been recognized more frequently in both experimental and modeling settings there is an increasing challenge to distinguish between different MMO-generating mechanisms. For stochastic systems there are additional routes for MMOs that include noise-induced oscillations and transients that are sustained with noise.
In this article, we analyzed MMOs of this type both in a biophysical model for the dynamics of stellate cells and in augmented versions of phenomenological FN-type models. Our analysis focuses on the STO part of the signal where the slow dynamics are prominent and noise has the greatest impact on the STOs. In both model types we observe signatures of two distinct oscillation-generating mechanisms: slow passage through a Hopf bifurcation (SPHB) and noise-induced coherent oscillations with frequencies close to that of the Hopf bifurcation. The latter of these appears in certain types of coherence resonance (CR).
For noise levels that are in the range of what is observed experimentally, features of the underlying deterministic models can be hidden or transformed in the stochastic time series. A further complication in identifying the MMO mechanism in stochastic models is that substantially distinct models can easily be tuned to produce very similar time series. Then model calibration based on time series alone must search through a larger parameter space or range of stochastic models needed for real world data.
We use a suite of measures to be able to distinguish between different MMO-and STO-generating mechanisms in the presence of noise, differentiating between the types of underlying models as well as analyzing the influence of certain general model parameters. Identifying possible MMO mechanisms in this way limits the broad range of parameter or model space for finding appropriate MMO mechanisms or calibrating a biophysical model. Furthermore this type of comparison can also identify appropriate classes of reduced models, with appropriate stochastic behavior over a range of parameters, to be used to approximate a full biophysical model. Such reduced models are used both for simplicity and for computational speed within larger modeling frameworks, as, e.g., in Ref. 13.
The suite of these measures was chosen to focus on the STOs in the ISI, where noise has the greatest impact on the character of the MMOs. We show that the measures can also be used to exploit the noise to identify the route to the MMOs. Given the variety of routes to MMOs, more than one such measure is needed, and we use three distinct measures (interspike interval, trend of the STO amplitude and noise-dependent coherence). While the focus of these measures is on bifurcation parameters, noise levels, and slowly varying control parameters, we found that these measures also reveal information about the refractory behavior or reset in the models. Furthermore, we have focused on using an approach that can easily be applied to time series, so that it can be used in both experimental and simulation settings.
We summarize the main characteristics obtained from these measures for different MMO generating mechanisms.
STOs of the CR-type have distinct features in these measures in comparison with those dominated by SPHB: 1. MMOs driven primarily by the deterministic behavior of SPHB display a strong trend in increasing STO amplitude, while those driven by CR have a weaker trend or no trend at all, for average amplitude.
2. For small noise ISI densities are highly concentrated for families driven by SPHB, in contrast to the STOs of the CR-type that have ISI densities with long tails. For larger noise values, these ISI densities are more similar and may depend on the return mechanism (see items 7-9 below).
3. The coherence measure has a clear optimal noise level for CR-driven STOs. PSDs are typically multipeaked for SPHB for smaller noise, so that a coherence measure is not well defined for small noise. For larger noise values, the coherence measure for SPHB is strictly decreasing.
While each model has a HB for the underlying 2D subsystem, differences in the criticality of these bifurcations can be observed depending on the underlying attracting deterministic dynamics.
4. For STOs dominated by SPHB, an increase in noise level typically drives a greater reduction in STO amplitude before the spike for a subcritical HB than a supercritical HB, since the latter has attracting STOs.
5. The underlying deterministic dynamics influence different behaviors of β as a function of noise level. Systems with a supercritical HB may have underlying attracting behavior of small oscillations that can support increased coherence. This type of deterministic behavior is typically not seen for subcritical HB.
6. For STOs dominated by SPHB in the subcritical case, even very low noise levels can drive a sampling from a variety of MMO families, together with a shift and spread of ISI values, making it difficult to distinguish from MMOs driven by other mechanisms.
In addition to identifying features related to noise level and bifurcation or slowly varying control parameters, we identified a number of additional behaviors that are related to the reset or refractory dynamics following a spike.
7. Reset and rate of increase of the control parameter in the simple IF model can be varied to capture the amplitude and ISI behavior of the different classes of time series in the physiological model. However, a simple IF model may not have the flexibility to capture different behaviors of the coherence measure.
8. For the dynamic voltage-dependent control parameters, an early escape to spiking can translate into a lower return value. In that case the ISI density does not shift with larger noise but is just less concentrated. This behavior distinguishes it from IF-type models with a fixed reset or models where the return mechanisms are inde-pendent of the spike. 9. For return mechanisms with limited damping of STOs following the spike, the coherence measure shows only a decrease in the coherence measure or no welldefined coherence over the ISI.
One important element in these comparisons is that characteristics of the time series, captured by the suite of measures, can change in distinct ways when the noise level is varied. A direct way of varying noise levels in neuroscience experiments is by injecting a noisy current. Recognizing that the noise has its greatest impact in the ISI where the slow dynamics are prominent leads to the proposal that tunable extrinsic noise can be introduced through an applied current to mimic the effects of the intrinsic noise. We show that this is indeed the case with an appropriate scaling of the noise in the physiological model. The distinctions between model features highlighted by these measures have been observed in an additional measure related to spike clusters, repeated spikes without STOs in the ISI 19 . There it was shown that the dynamic return mechanism of NLMFN, where an earlier spike can result in a reset value farther from the HB, can result in more robust STOs in the ISI, consistent with the ISI density behavior of NLMFN. The spike cluster frequency increases much more dramatically with noise for systems with SPHB-driven STOs, as compared with the case where STOs are CR-driven, consistent with the amplitude and ISI density results observed here. Also, small perturbations in the reset value can translate into a large increase in spike cluster frequency in stochastic systems, particularly where the STOs are SPHB-driven for an underlying subcritical HB.
Our analysis in this paper is focused on a specific type of MMOs, namely a combination of small amplitude oscillations and spikes relevant for neural systems. We have proposed measures that are focused on characteristics of MMOs that are particularly sensitive to the noise, due to the presence of multiple scales. This suggests that understanding the presence of multiple time scales and noise-sensitive characteristics of the underlying bifurcation structures would provide a solid basis for identifying measures that characterize MMOs in broader or more generic settings.
Appendix C: Details for the measures used to characterize MMOs 1. Amplitude trend before a spike First, the data (time series) was low-pass filtered by convolving it with a triangle function (e.g., Ref. 32) of a length that is on the order of half a STO period. A simple algorithm then finds the spikes, chooses a window before each spike, removes the average and finds the local maxima in each window going backwards in time. The amplitude at each consecutive maximum before a spike is averaged over all ISIs in the time series. In this analysis the choice of the low pass filter ('triangle length') is crucial, to avoid reduction of amplitudes or miscounting of the order of the maxima resulting from the time scale of the filter being too short or long.

Power spectral distribution
For the spiking models, spikes were removed from the time series by deleting data points corresponding to a typical spike duration (including recovery) after the respective variable (V or u) crossed a certain value. In the IF-type models, a similar but much shorter stretch of data was removed. The remaining data points were concatenated. Typically, PSDs computed from the full time series (including spikes and recovery) show significant power in broad frequency ranges that differ strongly between different models due to the different shapes of the spikes. Removing the spikes makes the power contribution of the STOs more clearly visible. To remove low frequency content in the PSDs, the average of the time series was removed. The routine spctrm from numerical recipes 29 was used to compute an estimate of the PSD. The normalization is such that the total power in the PSD is equal to the mean squared amplitude of the time series. Our 'standard' PSD for the MFN models is obtained from time series sampled at ∆t = 0.02 with a Bartlett window of size 4096 data points with overlap 29 . The 'standard' PSD for the SC models is obtained from time series sampled at ∆t = 2.5 with the same window function. For time series sampled at a different ∆t, we scale the PSD accordingly. We scale the PSDs such that the frequencies are expressed in 1/T STO (see Subsec. II C).

Coherence measure β
To compute the coherence measure β according to Eq. 8, we obtain the relevant values of the PSD by fitting a non-normalized Cauchy-Lorentz distribution to the peak of the PSD corresponding to STOs.
The fit is generally good for intermediate and strong values of noise. For weaker noise values the PSDs are often multi-peaked making it difficult to obtain a meaningful coherence measure. The coherence measure β has units of the PSD (squared unit of the considered variable when using T STO as a time scale) and depends strongly on the sampling rate of the original time series and the windowsize of the PSD-estimator. It is therefore difficult to compare absolute values of β between time series with different units, scales and sampling rates.