Parameter Estimation in the Age of Degeneracy and Unidentifiability

Abstract

Parameter estimation from observable or experimental data is a crucial stage in any modeling study. Identifiability refers to one’s ability to uniquely estimate the model parameters from the available data. Structural unidentifiability in dynamic models, the opposite of identifiability, is associated with the notion of degeneracy where multiple parameter sets produce the same pattern. Therefore, the inverse function of determining the model parameters from the data is not well defined. Degeneracy is not only a mathematical property of models, but it has also been reported in biological experiments. Classical studies on structural unidentifiability focused on the notion that one can at most identify combinations of unidentifiable model parameters. We have identified a different type of structural degeneracy/unidentifiability present in a family of models, which we refer to as the Lambda-Omega ($\mathsf{\Lambda }$-$\mathsf{\Omega }$) models. These are an extension of the classical lambda-omega ($\lambda$-$\omega$) models that have been used to model biological systems, and display a richer dynamic behavior and waveforms that range from sinusoidal to square wave to spike like. We show that the $\mathsf{\Lambda }$-$\mathsf{\Omega }$ models feature infinitely many parameter sets that produce identical stable oscillations, except possible for a phase shift (reflecting the initial phase). These degenerate parameters are not identifiable combinations of unidentifiable parameters as is the case in structural degeneracy. In fact, reducing the number of model parameters in the $\mathsf{\Lambda }$-$\mathsf{\Omega }$ models is minimal in the sense that each one controls a different aspect of the model dynamics and the dynamic complexity of the system would be reduced by reducing the number of parameters. We argue that the family of $\mathsf{\Lambda }$-$\mathsf{\Omega }$ models serves as a framework for the systematic investigation of degeneracy and identifiability in dynamic models and for the investigation of the interplay between structural and other forms of unidentifiability resulting on the lack of information from the experimental/observational data.

1. Introduction

Mathematical models are useful tools that can be used to interpret experimental or observational data, make quantitative predictions that are amenable for experimental testing, and identify the mechanisms that underlie the generation of pattern of activity in terms of the interactions among the system’s components [1,2,3]. An important step in connecting models with experimental or observational data is to estimate the model parameters by fitting the model outputs to the available data (Figure 1).

The ability of models to make predictions, provide mechanistic explanations, and be useful for decision making all depends on the accuracy and reliability of the parameter estimation process. A large number of parameter estimation tools are available to scientists as well as methods to discover data-drive non-linear dynamic equations [4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20] (and references therein) and tools to link data with models continue to develop.

A key feature of these tools is the minimization of an error function, a metric of the difference between the model output (simulations using estimated parameters) and the available (target) experimental/observed data. These data can be continuous (data collected with a high-frequency sampling, as compared to the scale of the process of interest, so that one can make the continuous approximation for mathematical purposes) or discrete (point process capturing the occurrence of events of interest). Error functions can be constructed by using all the data available in a point-to-point fashion (e.g., Figure 2A) or by computing one or more attributes that characterize these data (e.g., oscillation frequency, oscillation amplitude, oscillation duty cycle, stationary states, characteristic raise and decay times; Figure 2B).

There are several difficulties associated with the implementation of parameter estimation tools. These difficulties are data related (lack of access to all state variables, inconsistent gaps across trials), computational (algorithmic nature), statistical (data is noisy and therefore one can at best expect to estimate distributions of parameter values around a “true” mean), and structural (degeneracy, mathematical nature).

In an ideal situation, there would be a unique parameter set that fits the experimental/observational data (e.g., Figure 1A). We refer to the underlying model as identifiable [24]. In mathematical terms, the function linking parameter sets and data is bijective. In practice, these systems are subject to fluctuations from uncertain sources and one obtains a distribution of parameter sets (in a high-dimensional spaces whose dimensions are equal to the number of parameters of interest) centered at the “true” parameter set. The former has been termed structural identifiability, while the second has been termed practical identifiability or estimability [25]. Roughly speaking, the narrower the estimated parameters distribution (the smaller the variance), the more precise the parameter estimation process. However, this assumes that the center of the distribution is a good approximation to the “true” value. When necessary, the necessary validation of this assumption and the circumstance under which it is true can be achieved by using ground truth data generated by mathematical models. When continuous data is too noisy to be amenable for parameter estimation purposes, one might be able to extract useful discrete data in the form of attributes that capture the most relevant aspects of these data or point processes that are embedded in the data (e.g., Figure 2C and Figure 3). However, this is performed at the expense of loosing information and therefore has consequences for the accuracy of the parameter estimation process. Additional issues emerge when not all state variables are accessible or if gaps in discrete datasets are inconsistent across trials, requiring the use of imputation tools.

Unidentifiability, the opposite of identifiability [25,26,27,28,29,30,31], is associated with the concept of degeneracy. Degeneracy refers to these situations where multiple sets of parameter values can produce the same observable output (e.g., Figure 1B and Figure 2D1 if the output is the period), therefore making the inverse problem (of finding parameters given the data) ill posed. In mathematical terms, the function linking parameter sets and data is not injective (at most surjective). (There is in fact an implicit assumption that these links are surjective.) Degeneracy is not a problem associated with the statistical uncertainty in the knowledge of a unique parameter set from which the data is generated referred to above, but a structural problem inherent to mathematical models where the same patterns (e.g., temporal, spatial) can be obtained from multiple parameter sets. Degeneracy gives rise to the concept of level sets in parameter space (e.g., Figure 2D1) [32,33]. These are geometric objects (curves, surfaces, hypersurfaces) along which a given attribute of activity remains constant (see additional examples and a detailed discussion in the Appendix A), and are the geometric instantiation of the idea that a number of well-defined combinations of unidentifiable parameters can form an identifiable set. In other words, one can at most identify combinations of unidentifiable parameters.

Identifiability and unidentifiability have been studied in a number of systems. Structural identifiability has been studied using different approaches, including differential algebra [27,34,35,36,37,38,39], Taylor series [40], similarity transformations [41,42] and the Fisher information matrix [26,43,44,45] (see also [46,47,48]). In recent work [32], we have used numerical simulations and dynamical systems tools to characterize the frequency and duty cycle level sets and explain what are the dynamic mechanisms responsible for their generation in the FitzHugh–Nagumo (FHN) [49] model and the Morris–Lecar model [50,51] in the oscillatory regimes. Within a given attribute level set (e.g., period), the oscillatory patterns are non-identical (similarly to Figure 2D), and therefore dynamic balances operate to maintain the attribute of interest. The Morris—Lecar model is a neuronal model describing the interplay of two ionic currents (calcium and potassium), a leak (linear) current and a capacitive current. The FHN model is a caricature model that has been used to capture the oscillatory behavior in a number of fields (see Appendix A.3). In the relaxation oscillations regime, the oscillatory behavior of the FHN model is qualitatively similar to the one exhibited by the Morris–Lecar model and other models of biological relevance [2,52,53]. When the time scales are not well separated, the FHN model exhibits oscillations qualitatively similar to the HIV model in Figure 2.

The use of non-trivial attributes makes the comparison between data and model output computationally less expensive; instead of having to use a large set of points (e.g., the histogram envelope captured by the red curve in Figure 3B or the full time courses in Figure 2) one can use a significantly smaller number of attributes. In general, attributes refer to the parameters that characterize the data according to the way one chooses the data to be organized. On one extreme, the attributes are the whole set of data points and one ideally has access to all the model’s state variables. On the other extreme, the attributes are the minimal number of “numbers” necessary to fully characterize the available data. For example, frequency and amplitude are enough to describe a sinusoidal oscillation. A third attribute would be necessary if this oscillation is not centered at zero. Additional attributes such as the time constant governing the transition from an active (up) to a passive (down) phases would be necessary to discriminate between sinusoidal and relaxation oscillations. In the example in Figure 2, the oscillations can be characterized by the frequency, amplitude and duty cycle. In the example in Figure 3, the data consists of the interspike (or interpeak) intervals (ISIs or IPIs) histogram, which releases us from making detailed spike-time comparisons that could be strongly affected by noise and therefore complicates the error minimization process. This histogram can be characterized by four attributes (Figure 3B): the histogram peak $F_{max}$, the characteristic rise and decay constants ${\tau }_{rise}$ and ${\tau }_{decay}$, respectively, and the value of $\tau$ at which the histogram “begins”, ${\tau }_0$. Additional attributes are the mean and variance of the ISI distribution. (Neuronal spiking behavior typically approximates a Poisson process [54,55] and therefore the mean ISI is equal to its standard deviation.) Similar type of histograms could be computed for the frequency, amplitude and duty cycle in the HIV-1 viral kinetics example in Figure 2C. Note that attributes are used in connection with the data, while parameters are used in connection to the model used to generate or describe these data. A detailed discussion about the relationship between attributes (of patterns of activity) and model parameters in simple models including the prototypical FitzHugh–Nagumo oscillator [49] and a Hepatitis C viral model [56] (see also [57]) is presented in the Appendix A.

Oscillations are ubiquitous in dynamical systems, particularly in biological systems [2,21,52,53,58,59,60,61,62,63] where degeneracy has been observed, not only in models, but also in experiments [64,65,66] (see also [32,33]). The mechanisms underlying the presence of degeneracies in biological oscillations are not well understood. It is also unclear how the presence of degeneracies affect the parameter estimation process for models describing biological oscillations and oscillations in general.

We have identified a family of canonical models for unidentifiability/degeneracy where the limit cycle trajectories are unidentifiable and the degeneracy can be analytically computed. There are infinitely many parameter sets that produce identical stable oscillations, except possible for a for a phase shift, which depends on the initial phase. In other words, all the attributes of the of the oscillatory patterns are degenerate and have the same level sets, which can be analytically computed. We refer to the generic formulation of this model as the $\mathsf{\Lambda }$-$\mathsf{\Omega }$ model. It is a generalization of the so-called $\lambda$-$\omega$ system [67,68], also referred to as the Poincaré oscillator [69,70]. These are real-valued special cases of the complex Ginzburg–Landau equation [71,72,73,74] and their role for modeling biological systems are discussed in detail in [58] (see references therein). The general formulation can exhibit bistability and a number of bifurcations. The classical $\mathsf{\Lambda }$-$\mathsf{\Omega }$ exhibits sinusoidal oscillations. Extended versions of this model produce oscillations with more complex waveforms, including relaxation and spike-like oscillations, and a rich repertoire of dynamic behavior, including bistability and even multistability. In contrast to other models showing structural unidentifiability due to redundancy in the number of parameters with physical meaning (see Appendix A for a detailed discussion), the number of parameters in the $\mathsf{\Lambda }$-$\mathsf{\Omega }$ models cannot be reduced without reducing the model complexity. To our knowledge, the issues of unidentifiability and degeneracy in the family of $\mathsf{\Lambda }$-$\mathsf{\Omega }$ models has not been raised before.

Because of the properties described above, the family of $\mathsf{\Lambda }$-$\mathsf{\Omega }$ models are an ideal system to ask fundamental questions about degeneracy of dynamical systems and their relationship with parameter estimation unidentifiability. The goal of this paper is to address these issues. Our results contribute to the development of a framework for the investigation of unidentifiability in dynamic models.

2. Methods

2.1. A Viral Infection Kinetic Model

For illustrative purposes, we will use the HIV viral infection model presented in [21] (see also [22,23]), which, as noted in [23] belongs to a more general class of viral infection models [75]. The model is described by the following equations for the concentration T of the uninfected T cells, the concentration $T^{*}$ of productively infected T cells, and the concentration V of free virus particles in the blood

$$\frac{dT}{dt}=\delta -\alpha T+pT(1-T/T_{max})-kVT,$$

(1)

$$\frac{d{T}^{*}}{dt}=-\beta T^{*}+kVT,$$

(2)

$$\frac{dV}{dt}=-\gamma V+N\beta T^{*}-ikVT,$$

(3)

where $\alpha$ (day${}^{-1}$), $\beta$ (day${}^{-1}$) and $\gamma$ (day${}^{-1}$) are the death rates of the uninfected T cells, the infected T cells and the virus particles, respectively, k (mm${}^3$ day${}^{-1}$) is the contact rate of the uninfected T cells and the viral particles, $\delta$ (mm${}^{-3}$ day${}^{-1}$) represents the constant production of the T cells, N is the average number of virus particles produced by an infected T-cell, and p (day${}^{-1}$) and $T_{max}$ (mm${}^{-3}$) are the growth rate and the carrying capacity, respectively, associated with the logistic growth of uninfected T cells in the absence of virus particles, infected T cells and other factors. As shown in Figure 2, this model displays oscillatory behavior. Oscillations are typically observed in a number of virus dynamical systems [21,76,77,78].

2.2. Lambda-Omega ($\mathsf{\Lambda }$-$\mathsf{\Omega }$) Models: General Formulation

The $\mathsf{\Lambda }$-$\mathsf{\Omega }$ systems have the general form

$$\frac{dx}{dt}=\mathsf{\Lambda }(r)x-\mathsf{\Omega }(r)y,$$

(4)

$$\frac{dy}{dt}=\mathsf{\Omega }(r)x+\mathsf{\Lambda }(r)y,$$

(5)

with

$$r^2=x^2+y^2,$$

(6)

where x and y are state variables and t is time. The functions $\mathsf{\Lambda }$ and $\mathsf{\Omega }$ are typically chosen to be polynomials of degree two or four [67,69,70], but they could in principle be well-behaved enough functions with the necessary properties to support oscillations as discussed in this paper.

2.3. Modified $\mathsf{\Lambda }$-$\mathsf{\Omega }$ Models: General Formulation

The modified $\mathsf{\Lambda }$-$\mathsf{\Omega }$ models have the general form

$$\frac{dx}{dt}=\mathsf{\Lambda }(r)x-\mathsf{\Omega }(r)y,$$

(7)

$$\frac{dy}{dt}=\widehat{\mathsf{\Omega }}(r)x+\widehat{\mathsf{\Lambda }}(r)y,$$

(8)

where the functions $\widehat{\mathsf{\Lambda }}$ and $\widehat{\mathsf{\Omega }}$ are different from the functions $\mathsf{\Lambda }$ and $\mathsf{\Omega }$, but belong to the same class. We use this formulation to investigate the consequences of the break of symmetry present when $\widehat{\mathsf{\Lambda }}=\mathsf{\Lambda }$ and $\widehat{\mathsf{\Omega }}=\mathsf{\Omega }$.

2.4. Parameter Estimation Algorithms

To study the impacts of the $\mathsf{\Lambda }$-$\mathsf{\Omega }$ model’s degenerate mathematical structure on the estimation of its parameters, the behavior of three parameter estimation algorithms, covering a range of methods, were characterized. These include a genetic algorithm (GA) [10,11,79,80], gradient descent algorithm [6,81], and the sequential neural posterior estimation (SNPE) method [15,82,83]. The purpose of each of them is to find parameter values that fully minimize the difference between estimated solutions produced by the algorithm’s output and a dataset representing the target (optimal) solution to the model. In practical applications of parameter estimation, experimental (or observational) data typically acts as this set of target data; however, for the purpose of this study, when working with the $\lambda$-$\omega$ model, a simulated output using a predetermined set of parameters ($\lambda =1$, $\omega =1$, $a=1$, $b=1$) acts as a ground truth dataset (GTD) to stand in as a proxy for experimental data. To categorize the impacts of noise on the identifiability of the model’s parameters, each parameter estimation algorithm was tested using datasets with variable levels of noise applied to the GTD. We used Gaussian white noise with mean zero and variance D ($D=0,100,250$), added to the both equations in the $\lambda$-$\omega$ model. This approach seeks to activate the transient dynamics of the system by perturbing its trajectory from the limit cycle. Under certain conditions, the activation of transients is believed to resolve some types of unidentifiability issues. Specifically in the case of the $\lambda$-$\omega$ model, if one is able to estimate $\lambda$ from the transient dynamics, then the estimation of b is straightforward, eliminating the unidentifiability of the two parameters. In our simulations, we used $\Delta t=0.01$.

2.4.1. Stochastic Optimization: Genetic Algorithms (GA)

The GA is an evolutionary algorithm based on the principles of natural selection, in which only “fit” individuals of a given generation are used to create subsequent generations [79,80] (see also [10,11,12,13]). The fitness of parameter sets within the population are those that produce an output most closely resembling the GTD. In our numerical experiments, we quantify this as the difference between the ground truth and the simulated differences and amplitudes. However, other implementations of the algorithm can consider the attributes independently or simply use the time series produced by each parameter set as a whole. Both of these approaches were tested and showed no improvement on the identifiability of parameters. We used a multiobjective GA, where the estimation process begins by creating an initial population of randomly distributed parameter sets [79]; for the purposes of this study, populations were created using random values within the range of $[0,2]$. In each generation, only the top half of the population that most optimally matches the GTD are used as “parents” to build the next generation. The next generation consists of both crosses between two independent parent sets and the original parent sets themselves with random mutations being added to parameter values at a rate of 1 over the number of parameters estimated. The number of generations produced and population size are specified by the user, with executions in our study using 500 generations and populations of 250 individuals. To ensure the adequacy of these conditions, trials with a larger population size (up to 2000 individuals) and/or higher number of generations (up to 1000 generations) were tested and shown to have no impact on the performance of the GA.

2.4.2. Sequential Neural Posterior Estimation (SNPE)

We implemented SNPE using the simulation-based inference (SBI) developer toolbox [15,82,83]. SNPE falls under a larger category of algorithms called Sequential Neural Processes, which use Bayesian inference to estimate parameter values [15,82,83]. This technique is a statistical method for calculating likelihood of events based on prior observations, and it is applied to determine which parameter values are most likely to produce the behavior observed in the GTD [15,82,83]. Similarly to the GA, Sequential Neural Processes start with an initial distribution of parameters, which is then refined, in this case using Bayesian inference. For this application of SNPE, an initial uniform distribution on $[0,2]$ is created as a prior distribution. Then, 10,000 simulations of the model are created by sampling from parameter sets within this distribution in order to infer a posterior distribution based on which sets produce an output most closely matching the GTD.

3. Results

3.1. Unidentifiability and Degeneracy: A Preliminary Discussion Using Simple Models

Here, we introduce some ideas related to unidentifiability and degeneracy in dynamical systems, including the concepts of activity attributes, level sets and parameter redundancy in models and their outputs. We do this in the context of relatively simple models. A detailed mathematical description of these models and the unidentifiability analysis are presented in the Appendix A.

In the 1D minimal model example discussed in the Appendix A.1, the minimal model has no level sets. The number of minimal model parameters is equal to the number of attributes that characterize the pattern. In fact, the minimal model parameters are themselves the attributes. In the 2D minimal model example discussed in the Appendix A.2, there are no level sets for the time constant $\tau$ and this attribute is determined uniquely by the minimal model parameter b in Equation (A6). In contrast, the damped oscillations frequency $\omega$ (A7) is a combination of the minimal model parameters b and c and has level sets in this parameter space. Clearly, if one can estimate the two attributes $\tau$ and $\omega$, one can identify the minimal model parameters b and c. If one can only estimate $\omega$, then the model has degeneracy. Typically, by design, the minimal model parameters have no physical (or biological) meaning, but only dynamic meaning. They are a minimal set of model parameters in that eliminating one of them reduces the dynamic complexity.

The situation changes when models are created to reproduce or explain experimental or observational data. The number of physically meaningful model parameters is larger than the number of minimal model parameters needed to maintain the model dynamic complexity. The level sets and degeneracy that we find in the realistic models in both cases (Appendix A.1 and Appendix A.2) emerge because the minimal model parameters (with dynamic meaning) are combinations of (realistic) parameters with physical meaning. From that perspective, the realistic parameters generate a dynamic redundancy, which is reflected by the fact that the minimal model parameters are combinations of physical parameters. This type of degeneracy cannot be resolved unless one has access to additional data.

The FHN model discussed in the Appendix A.3 is a non-linear minimal model in the oscillatory regime. In the example presented in the Appendix A.3, there are level sets, which emerge as two or more minimal model parameter sets dynamically balance each other to maintain an attribute constant (e.g., frequency, see Figure A3A,B) [32] similarly to the HIV model oscillations discussed in the context of Figure 2. Attribute level sets such as frequency, amplitude, and duty cycle do not coincide except perhaps in some specific cases [32]. Additional degeneracies emerge in realistic models (e.g., the Morris–Lecar model [32]) as the number of parameters with biological meaning increases with the consequent increase in parameter redundancies. However, it is not always possible to reduce non-linear models to minimal models.

Figure A3C in the Appendix A.1 illustrates the fact that additional attributes may be necessary to characterize oscillatory patterns. The non-identical sinusoidal and square-wave functions have identical frequency, amplitude and duty cycle. This raises the need of identifying the appropriate attributes to characterize a given pattern and their optimal number.

3.2. The $\mathsf{\Lambda }$-$\mathsf{\Omega }$  Systems Are Canonical Models for Unidentifiability/Degeneracy

Here, we discuss the role of the $\mathsf{\Lambda }$-$\mathsf{\Omega }$ models generically described by Equations (4)–(6) as canonical models for unidentifiability where a limit cycle with the same frequency can be generated by multiple combinations of model parameters.

A change of coordinates from Cartesian to polar

$$x=rcos(\theta )\mathrm{and}y=rsin(\theta ),$$

(9)

transforms system (4) and (5) into

$$\frac{dr}{dt}=r\mathsf{\Lambda }(r),$$

(10)

$$\frac{d\theta }{dt}=\mathsf{\Omega }(r).$$

(11)

The equation for the radius of the limit cycle is independent of the the angular velocity $\theta$. Therefore, the limit cycles are circles around zero and their radii $\overline{r}$ are the solutions of

$$\mathsf{\Lambda }(\overline{r})=0.$$

(12)

The number of limit cycles $\overline{r}$ depends on the form and complexity of the function $\mathsf{\Lambda }(r)$. Since $r>0$, a limit cycle of radius $\overline{r}$ is stable if and only if ${\Delta }^{\prime }(\overline{r})<0$. The frequency of the solution along this limit cycle is $\mathsf{\Omega }(\overline{r})$ and the period is $2\pi /\mathsf{\Omega }(\overline{r})$. Because, if a stable limit cycle exists, in the limit of $t\to \infty$ the dynamics are described completely by $\mathsf{\Omega }(\overline{r})$, $\mathsf{\Lambda }$-$\mathsf{\Omega }$ systems are often referred to as phase oscillators. Note that $\overline{r}=0$ corresponds to a fixed point and therefore the $\mathsf{\Lambda }$-$\mathsf{\Omega }$ systems always have a fixed point whose stability properties also depend on the sign of ${\Delta }^{\prime }(0)$.

Limit cycle and frequency degeneracy emerges when the function $\Delta$ is defined by two or more parameters. Additional degeneracy in the frequency results from the parameters defining the function $\mathsf{\Omega }$.

3.3. Lambda-Omega Systems of Order 2: Limit Cycle Degeneracy

$\mathsf{\Lambda }$-$\mathsf{\Omega }$ systems of order 2 ($\mathsf{\Lambda }$-${\mathsf{\Omega }}_2$) correspond to the following quadratic choices for $\mathsf{\Lambda }(r)$ and $\mathsf{\Omega }(r)$.

$$\mathsf{\Lambda }(r)=\lambda -b{r}^2\mathrm{and}\mathsf{\Omega }(r)=\omega +a{r}^2,$$

(13)

where $\omega$, $\lambda$, a and b are constant. Equations (4) and (5) become

$$\frac{dx}{dt}=(\lambda -b{r}^2)x-(\omega +a{r}^2)y,$$

(14)

$$\frac{dy}{dt}=(\omega +a{r}^2)x+(\lambda -b{r}^2)y.$$

(15)

Following [67], we refer to these models as $\lambda$-$\omega$ models. They have been also referred to as Poincaré oscillators [69] and are real-valued special cases of the complex Ginzburg–Landau equation as mentioned above [71,72,73,74].

The parameters $\lambda$ and $\omega$ control the linear dynamics and the parameters a and b control the contribution of the non-linear terms to the system’s dynamics. In this sense, the $\lambda$-$\omega$ model is minimal. As we will see, the parameters $\lambda$ and $\omega$ play a stronger role in controlling the transient dynamics than the parameters a and b.

3.3.1. Limit Cycles, Fixed Points, Stability and Oscillation Frequency

These $\lambda$-$\omega$ systems have a single limit cycle for

$$\overline{r}=\sqrt{\frac{\lambda }{b}},$$

(16)

if $\lambda /b>0$, which is stable if and only if $\lambda <0$. Otherwise, the fixed point $(\overline{x},\overline{y})=(0,0)$ is stable. Along the limit cycle, the frequency is constant, and given by

$$\overline{\omega }=\omega +a\frac{\lambda }{b}.$$

(17)

In Cartesian coordinates, along the limit cycle

$$\frac{dx}{dt}=-\overline{\omega }y,\mathrm{and}\frac{dy}{dt}=\overline{\omega }x,$$

(18)

whose solution satisfying $x(0)=0$ and $y(0)=1$ is given by

$$x(t)=\sqrt{\frac{\lambda }{b}}sin(\overline{\omega }t)\mathrm{and}y(t)=\sqrt{\frac{\lambda }{b}}cos(\overline{\omega }t).$$

(19)

Solutions for other initial conditions are phase shifts of this solution. This implies that limit cycles are identical and trajectories moving along the limit cycle have the same frequency. Figure 4a illustrates the x and y traces for representative parameter values. Figure 4b shows the corresponding phase-plane diagram, including the circular limit cycle characteristic of sinusoidal trajectories. Finally, Figure 4c shows the speed diagram, also characteristic of sinusoidal trajectories.

3.3.2. Limit Cycle Degeneracy: The Level Sets for All Possible Attributes Coincide

Along the limit cycle

$${\overline{r}}^2=\frac{\lambda }{b}\mathrm{and}\frac{d\theta }{dt}=\overline{\omega }.$$

(20)

The limit cycles have the same amplitude for the infinitely many combinations of values $\lambda$ and b for which their quotient is constant and they have the same frequency for the infinitely many combinations of the four parameters such that $\overline{\omega }$ is constant. The amplitude level sets are curves in the two-dimensional $\lambda$-b parameter space and the frequency level sets are hypersurfaces in the four-dimensional $\lambda$-b-$\omega$-a parameters space. Note that even along an amplitude level set there is frequency degeneracy along the limit cycle.

Figure 5 illustrates the identical limit cycles for different combinations of parameter values. The different phase-plane diagrams supporting this limit cycle degeneracy indicate that the transient dynamics are different for the different cases. This is also clear from Equation (10). The time constant associated with the evolution of the envelope oscillations approaching the limit cycle (or moving away from it, if it is unstable) is $\lambda$ dependent.

3.3.3. Limit Cycle Degeneracy: Symmetries in the Phase-Plane Diagram

The x- and y-nullclines for the $\lambda$-$\omega$ system are given by the solutions $y=F(x)$ and $x=G(y)$ of

$$a{y}^3+bx{y}^2+(\omega +a{x}^2)y-(\lambda -b{x}^2)x=0$$

(21)

and

$$a{x}^3-b{x}^{2}y+(\omega +a{y}^2)x+(\lambda -b{y}^2)y=0,$$

(22)

respectively.

The phase-plane diagram in Figure 4b shows the nullclines and trajectory for the solution to the $\lambda$-$\omega$ system for representative parameter values whose solution is shown in Figure 4a. We express the y-nullcline as a function of y, and not of x, because as a function of x this nullcline is not invertible (Figure 4b).

System (14) and (15) is preserved under the following changes of coordinates: (i) $(x,y)\to (-x,-y)$, (ii) $(x,y)\to (-y,x)$, and (iii) $(x,y)\to (y,-x)$, but not under $(x,y)\to (y,x)$ and $(x,y)\to (-y,-x)$. The preservation property (i) reflects the fact that the two nullclines are odd-functions: $F(x)=-F(-x)$ and $G(y)=-G(-y)$. The preservation properties (ii) and (iii) reflect the fact that the two nullclines are the negative of one another $F(x)=-G(x)$. Taken together, these properties reflect the fact that the two nullclines are orthogonally identical in the sense that they can be obtained from one another by a $\pi /2$ rotation (Figure 4b). The limit cycles are preserved if both nullclines are rotated by the same angle [74]. Figure 5 illustrates that these asymmetries are preserved for parameter values for different values of the amplitude and frequency of the limit cycle.

3.4. Compensation Mechanisms of Generation of Level Sets in the $\lambda$-$\omega$ Model

Here, we investigate the dynamic mechanisms of generation of level sets in terms of the model parameters by focusing on the geometric properties of the phase-plane diagrams in Cartesian coordinates. We aim to identify the geometric constraints in the phase-plane diagrams for the parameter values that belong to the same level sets and give rise to the limit cycle degeneracy.

In our analysis, we assume to have no information about the properties of the limit cycles and level sets from the polar coordinates analysis above. The reason for doing so is to develop tools to analyze the models when the symmetries of the $\lambda$-$\omega$ model are broken and the level sets for the different attributes are split. In these cases, Cartesian coordinates are more amenable for analysis than polar coordinates.

Because of the model symmetries reflected in the geometric properties of the x- and y-nullclines, it is natural to assume that if a limit cycle exists, then it has circular shape and it is centered at the origin (Figure 6).

3.4.1. Level Sets in the $\lambda$-b Parameter Space

Setting $y=0$ in the x-nullcline (21) shows that it crosses the x-axis at $x=\pm \sqrt{\lambda /b}$, independently of the values of $\omega$ and $\lambda$. Likewise, setting $x=0$ in the y-nullcline (22) shows that it crosses the y-axis at $y=\pm \sqrt{\lambda /b}$, independently of the values of $\omega$ and $\lambda$ (Figure 7a).

The assumption that the limit cycle is circular and centered at the origin implies that the limit cycle trajectory is vertical on the x-axis and horizontal on the y-axis. This requires limit cycle trajectory to intersect the x-nullcline when the latter intersects the x-axis (at $x=\pm \sqrt{\lambda /b}$) and to intersect the y-nullcline when the latter intersects the y-axis (at $y=\pm \sqrt{\lambda /b}$). In other words, the intersection of the x-nullcline with the x-axis (or the intersection of the y-nullcline with the y-axis) determines the radius of the limit cycle trajectory without the need of the transformation to polar coordinates.

In general, changes in $\lambda$ and b change the shape of the nullclines (compare Figure 5a–c). However, when $\lambda$ and b change in such a way that $\lambda /b$ remains constant, then the intersection of the x-nullcline with the x-axis is constant and therefore the radius of the limit cycle is constant. In other words, along a given level set, $\lambda$ and b compensate for each other so as to maintain the intersection of the x-nullcline with the x-axis constant (or, by symmetry, the intersection of the y-nullcline with the y-axis constant).

Along these level sets, the frequency of the oscillations on the limit cycle is constant since from (14) and (15) with $\lambda -b{r}^2=0$ it follows that

$$\frac{dx}{dt}=-(\omega +a\frac{\lambda }{b})y,$$

(23)

$$\frac{dy}{dt}=(\omega +a\frac{\lambda }{b})x,$$

(24)

whose solution is a sinusoidal function of frequency $\omega +a\lambda b^{-1}$ and phase ${tan}^{-1}(1)$. This implies that solutions belonging to the same level set are identical.

3.4.2. Frequency Level Sets in the $\omega$-a Parameter Space: Static Compensation Mechanisms When $\lambda$ and $\omega$ Belong to the Same Level Set

Computation of the derivative of the x-nullcline (21) at $(\pm \sqrt{\lambda b^{-1}},0)$ gives

$$\frac{dy}{dt}=\mp \frac{2\lambda }{\omega +a\lambda /b}.$$

(25)

Similarly, computation of the derivative of the y-nullcline (22) at $(0,\pm \sqrt{\lambda b^{-1}})$ gives

$$\frac{dy}{dt}=\pm \frac{\omega +a\lambda /b}{2\lambda }.$$

(26)

Therefore, the nullclines for parameter values on the same frequency level set ($\omega +ab/\lambda =cont$) for the same value of $\lambda$ are tangent at their point of intersection (Figure 7b). This constraint is lost when $\lambda$ changes, even though the parameter values may belong to the same amplitude level set (Figure 7c).

3.4.3. Frequency Level Sets in the $\omega$-a Parameter Space: Dynamic Compensation Mechanisms When $\lambda$ and b Do Not Belong to the Same Level Set

Here, we focus on frequency level sets ($\omega +ab/\lambda =const$) for parameter values such that $\lambda /b\ne const$. In these cases, the oscillation amplitudes are different (Figure 8A,B). The x- and y-nullclines still intersect the corresponding limit cycles at the x- and y-axis, respectively, but, in order to maintain a constant frequency, the speed of the limit cycle trajectory needs to be greater than the limit cycle radius (Figure 8C).

3.5. Unfolding of Attribute Level Sets for Different Attributes Due to Symmetry Breaking

The complete degeneracy of the limit cycle oscillations’ in the $\lambda$-$\omega$ model results from a number of symmetries present in the model, and a break in this symmetry is expected to disrupt the degeneracy. However, from previous work [32,33] and our discussion above in the context of the HIV model (Figure 2), we expect degeneracies to remain for the attributes characterizing the oscillation patterns. Here, we address this issue by using the modified $\mathsf{\Lambda }$-$\mathsf{\Omega }$ model (7) and (8) where all the functions are quadratic and r is given by (6). The resulting modified $\lambda$-$\omega$ model reads

$$\frac{dx}{dt}=(\lambda -b{r}^2)x-(\omega +a{r}^2)y,$$

(27)

$$\frac{dy}{dt}=(\widehat{\omega }+\widehat{a}{r}^2)x+(\widehat{\lambda }-\widehat{b}{r}^2)y,$$

(28)

and reduces to the $\lambda$-$\omega$ model discussed above for $\widehat{\lambda }=\lambda$, $\widehat{\omega }=\omega$, $\widehat{a}=a$ and $\widehat{b}=b$.

Geometrically, from the dynamical systems point of view, there are four primary ways of breaking the symmetry between the x- and y-nullclines: (i) generating of a time scale separation between the two variables, (ii) expanding or shrinking the linear component of one equation with respect to the other, (iii) rotating one nullcline while keeping the other fixed, and (iv) displacing one nullcline while keeping the other fixed.

Here, we systematically analyze the effects of these symmetry breaking actions by choosing the appropriate values of the parameters in Equation (28) while keeping the parameter values in Equation (27) fixed. For simplicity, we focus on the amplitude and frequency level sets that emerge for each one of the models. In all cases, the canonical level sets are not preserved and new level sets emerge, which are different for different attributes and the four different ways of breaking the model symmetry. As expected, degeneracy remains for these attributes. As the parameter values move away from symmetry, both the shapes of the limit cycles and the shapes of the speed diagrams change, indicating that new dynamic compensation mechanisms to maintain the attributes’ degeneracy arise.

3.5.1. Onset of Frequency and Amplitude Level Sets Due to the Presence of a Time Scale Separation

A time scale separation between the first and second equations is generated by multiplying each one of the parameters in the y-equation by a parameter $ϵ$

$$\widehat{\omega }=ϵ\omega ,\widehat{\lambda }=ϵ\lambda ,\widehat{a}=ϵa,\widehat{b}=ϵb.$$

(29)

Figure 8A shows the frequency and amplitude level sets for representative parameter values. The presence of a time scale separation does not break the rotational symmetry of the two nullclines, but only the relative dynamics of the two variables, which is reflected in shapes of the limit cycles (Figure 8B2,C2) and the non-uniform speed along the limit cycle (Figure 8B3,C3) as well as the shapes of the oscillatory patterns (Figure 8B1,C1). The speed diagrams show the presence of dynamic compensation mechanisms responsible for the generation of the level sets.

3.5.2. Onset of Frequency and Amplitude Level Sets Due to Changes in the Linear Properties of the y-Nullcline

The region of linear behavior around the origin of the y-nullcline can be expanded or shrunk by multiplying the linear parameters $\omega$ and $\lambda$ by a constant $c>1$, while leaving the non-linear parameters a and b unchanged

$$\widehat{\omega }=c\omega ,\widehat{\lambda }=c\lambda ,\widehat{a}=a,\widehat{b}=b.$$

(30)

Figure 9 shows our results for three values of c: $c=1$ (symmetry), $c=1.25$ and $c=1.5$. The amplitude level sets almost coincide, but the frequency level sets are more separated (Figure 9a). This is reflected in the shapes of the limit cycles (Figure 9b2,c2) and speed diagrams (Figure 9b3,c3). The speed diagrams show the presence of dynamic compensation mechanisms responsible for the generation of the level sets. These mechanisms are different from the ones discussed above.

3.5.3. Onset of Frequency and Amplitude Level Sets Due to Rotation of the y-Nullcline with Respect to the x-Nullcline

In order for the y-nullcline to be rotated by an angle $\varphi$, the parameters $\widehat{\omega }$, $\widehat{\lambda }$, $\widehat{a}$ and $\widehat{b}$ must be given by (see Appendix B.1)

$$\widehat{a}=acos(\varphi )-bsin(\varphi ),$$

(31)

$$\widehat{b}=asin(\varphi )+bcos(\varphi ),$$

(32)

$$\widehat{\omega }=\omega cos(\varphi )+\lambda sin(\varphi ),$$

(33)

$$\widehat{\lambda }=\lambda cos(\varphi )-\omega sin(\varphi ).$$

(34)

Figure 10 shows our results for three values of $\varphi$: $\varphi =0$ (symmetry), $\varphi =\pi /6$ and $\varphi =\pi /3$. The frequency and amplitude level sets do not coincide. The speed diagrams show the presence of dynamic compensation mechanisms responsible for the generation of the level sets. Comparison with the previous cases shows that these mechanisms are different from the ones discussed above. In particular, the speed diagrams in Figure 10b3,c3 have prominent intersections, which are almost entirely absent in the cases discussed above.

3.5.4. Onset of Frequency and Amplitude Level Sets Due to the Displacement of the y-Nullcline with Respect to the x-Nullcline

The modified $\lambda$-${\omega }_2$ system where the y-nullcline is displaced $x_c$ units in the horizontal direction from the position determined by system (14) and (15) is given by

$$\frac{dx}{dt}=(\lambda -b{r}^2)x-(\omega +a{r}^2)y,$$

(35)

$$\frac{dy}{dt}=(\omega +a{r}^2)x+(\lambda -b{r}^2)y+P(x,y),$$

(36)

where

$$P(x,y)=-3a{x}_c{x}^2-a{x}_c{y}^2+2b{x}_{c}xy+3a{x}_c^{2}x-b{x}_c^{2}y-a{x}_c^3-\omega x_c.$$

(37)

Figure 11 shows our results for three values of $x_c$: $x_0=0$ (symmetry), $x_c=0.2$ and $x_c=0.4$. These results are similar to the previous cases, in that the frequency and amplitude level sets do not coincide and the speed diagrams show the presence of dynamic compensation mechanisms responsible for the generation of the level sets with prominent intersections (Figure 11b3,b3) as is only observed in Section 3.5.3.

3.6. Lambda-Omega Systems of Order 4: Degeneracy, Hopf Bifurcations and the Emergence of Bistability

$\mathsf{\Lambda }$-$\mathsf{\Omega }$ systems of order 4 ($\mathsf{\Lambda }$-${\mathsf{\Omega }}_4$) correspond to the following quartic choices for $\mathsf{\Lambda }(r)$ and $\mathsf{\Omega }(r)$

$$\mathsf{\Lambda }(r)=\lambda -b{r}^2+c{r}^4\mathrm{and}\mathsf{\Omega }(r)=\omega +a{r}^2+d{r}^4.$$

(38)

They are the normal forms of the Hopf bifurcation [84], and they have been used to model the dynamics of oscillatory systems in response to external inputs [85,86].

The limit cycles, if they exist are given by

$$\overline{r}=\sqrt{\frac{b\pm \sqrt{b^2-4\lambda c}}{2\lambda }}.$$

(39)

The fixed points and limit cycles are stable (unstable) if $-b+2c{\overline{r}}^2<0$ ($-b+2c{\overline{r}}^2>0$).

Bistability emerges when there is a stable fixed point, a stable limit cycle and an unstable limit cycle in between. There are a number of alternative scenarios that include the presence of one or two limit cycles. Along the limit cycles, the frequency is given by

$$\overline{\omega }=\omega +a{\overline{r}}^2+d{\overline{r}}^4.$$

(40)

The degeneracies that emerge in this system, particularly in the bistable case, are more complex than in the $\lambda$-$\omega$ system discussed above and the level sets live in higher-dimensional parameter spaces. The analysis of these degeneracies and the compensation mechanisms that give rise to them goes along the same lines of the discussion above for the $\lambda$-$\omega$ system; however, a detailed analysis is beyond the scope of this paper.

3.7. Degeneracy in Extended $\lambda$-$\omega$ Systems with Generic Waveforms: From Sinusoidal to Square and Spike-like Waveforms

The limit cycles in the standard $\lambda$-$\omega$ systems have sinusoidal waveforms. However, non-sinusoidal waveforms are ubiquitous in biological systems. Variations of the $\lambda$-$\omega$ system showing square-wave and spike-like waveforms have been introduced in [87]. Here, we discuss the presence of degeneracy, in particular the complete limit cycle degeneracy in these systems.

The extended lambda-omega $\lambda$-$\omega$ systems (or extended Poincare oscillators) have the form

$$\frac{dx}{dt}=f(x,y;p),$$

(41)

$$\frac{dy}{dt}=g(x,y;p),$$

(42)

where

$$f(x,y;p)=\lambda x-\omega y-(bx+ay)(x^2+y^2)-cy\left(1+\frac{x}{\sqrt{x^2+y^2}}\right)-d\frac{y^3}{x^2+y^2},$$

(43)

and

$$g(x,y;p)=\omega x+\lambda y+(ax-by)(x^2+y^2)+cy\left(1+\frac{x}{\sqrt{x^2+y^2}}\right)+d\frac{x{y}^2}{x^2+y^2}.$$

(44)

The parameters p ($\omega$, $\lambda$, a, b, c and d) are non-negative constants. For $c=d=0$, system (4) and (5) reduces to the standard $\lambda$-$\omega$ system. The parameters $\lambda$ and $\omega$ control the linear dynamics, the parameters a and b control the cubic-like non-linear terms (the non-linearities in the standard $\lambda$-$\omega$ system) and the parameters c and d control the remaining non-linear terms.

A change of coordinates from Cartesian to polar

$$x=rcos(\theta )\mathrm{and}y=rsin(\theta ),$$

(45)

transforms system (4) and (5) into

$$\frac{dr}{dt}=(\lambda -b{r}^2)r,$$

(46)

$$\frac{d\theta }{dt}=\omega +a{r}^2+2c{cos}^2(\theta /2)+d{sin}^2(\theta ),$$

(47)

where r is given by (6)

3.7.1. Sinusoidal, Spike-like and Square-Wave Sustained Oscillatory Patterns for the Extended $\lambda$-$\omega$ System

From Equation (10), the extended $\lambda$-$\omega$ system (4) and (5) has a circular limit cycle of radius

$$\overline{r}=\sqrt{\frac{\lambda }{b}}$$

(48)

provided $\lambda /b>0$ (Figure 4). The limit cycle is stable if and only if $\lambda <0$. Otherwise, the fixed point $(\overline{x},\overline{y})=(0,0)$ is stable. The properties of the right-hand side of Equation (A16)

$$F(\theta )=\omega +a\frac{\lambda }{b}+2c{cos}^2(\theta /2)+d{sin}^2(\theta )$$

(49)

determine the various types of temporal patterns along these limit cycles.

For $c=d=0$, $F(\theta )$ is constant (Figure 12A2), and therefore $\lambda$-$\omega$ evolves linearly with time (Figure 12A3). The solutions are sinusoidal (Figure 12A1) with period

$$T=\frac{2\pi }{\mathsf{\Omega }}$$

(50)

where

$$\mathsf{\Omega }=\omega +a\frac{\lambda }{b}.$$

(51)

For $d=0$, $F(\theta )$ peaks at $\theta =0$ (Figure 12B2), causing $\theta$ to evolve in a fast-slow manner (Figure 12B3) giving rise to spike-like patterns (Figure 12B1). The period is given by (see Appendix B)

$$T=\frac{2\pi }{\sqrt{\mathsf{\Omega }}\sqrt{\mathsf{\Omega }+2c}}.$$

(52)

For $c=0$, $F(\theta )$ peaks twice at $\theta =\pm \pi /2$ (Figure 12C2), causing $\theta$ to evolve in a fast slow manner, changing from fast to slow and vice versa twice within a cycle (Figure 12C3) and give rise to square-wave patterns (Figure 12C1). The period is given by (see Appendix B)

$$T=\frac{2\pi }{\sqrt{\mathsf{\Omega }}\sqrt{\mathsf{\Omega }+d}}.$$

(53)

3.7.2. Amplitude and Period Level Sets

The amplitude level sets are as for the standard $\lambda$-$\omega$ system discussed above. The frequency level sets for the extended $\lambda$-$\omega$ system depend on the values of c and d. It is illustrative to consider two separate cases: (i) $c=0$ and (ii) $d=0$ for which the period T can be easily computed.

From (52), for $d=0$, the period level sets are given by the solution to

$$\sqrt{\mathsf{\Omega }}\sqrt{\mathsf{\Omega }+2c}=K_{freq}.$$

(54)

By solving one gets

$$\omega +a\frac{\lambda }{b}=-c+\sqrt{c^2+K_{freq}^2}.$$

(55)

Similarly, from (52), for $c=0$, the period level sets are given by the solution to

$$\sqrt{\mathsf{\Omega }}\sqrt{\mathsf{\Omega }+d}=K_{freq}.$$

(56)

By solving one gets

$$\omega +a\frac{\lambda }{b}=\frac{-d+\sqrt{d^2+4K_{freq}^2}}{2}$$

(57)

Similarly to the standard $\lambda$-$\omega$ system ($c=d=0$), for all combinations of parameter values satisfying ${\overline{r}}^2=\lambda /b$ and either (55) or (57) the system has the same amplitude and frequency level sets.

3.8. Parameter Estimation Algorithms Recover at Most Parameter Values on a Level Set

One potential way to disambiguate degenerate models is to obtain information about the transient behavior of trajectories. For example, for the $\lambda$-$\omega$ models, the parameter $\lambda$ controls the evolution of the envelope of the oscillations with increasing amplitude converging to the limit cycle (with amplitude $\sqrt{\lambda /b}$). However, in realistic situations one does not have clean access to the transient evolution of trajectories unless one is able to perform perturbations to the oscillatory patterns. It has been argued that the presence of noise is useful to improve parameter estimation results since noise causes trajectories to explore wider regions of the phase space and therefore more information about the dynamics, particularly the transient behavior of trajectories, is available to estimate parameters.

We used the two parameter estimation algorithms described in Section 2.4. In Figure 13, we fixed the values of $\omega =a=1$ and attempted to estimate the values of $\lambda$ and b using ground truth data (GTD) generated by using $\lambda =b=\omega =1$. In Figure 14, we used the same GTD and attempted to estimate the four parameters.

Our results show that relatively small levels of noise are not enough to provide useful information, while higher levels of noise affect the efficacy of the parameter estimation algorithm. While these results are not conclusive, they illustrate that the degeneracy present in the models poses an obstacle to the ability to estimate the model parameters.

4. Discussion

Parameter estimation from experimental or observational data is a critical stage in any modeling study. There are several issues that conspire against the accurate estimation of model parameters, particularly the parameters in dynamic models (e.g., described by differential equations). Data are noisy and therefore one can at best expect to be able to identify a distribution of model parameters centered around the “true” parameter set (if it exists and is unique). In addition, one may not have access to data for all the state variables assumed to describe the process. Moreover, one may have experimental/observational access to discrete events (e.g., sequence neuronal spike times, record of a viral infection), but not the data that underlie the generation of these events (e.g., membrane potential, concentrations of the virus, infected and uninfected cells). Besides, if the data are too noisy or exhibit strong irregularities, one may only have reliable access to a number of activity attributes (e.g., oscillation preferred frequency, neuronal spiking mean firing rate) that may not be enough to capture the complete dynamic structure of these data (e.g., sinusoidal and relaxation oscillations may have the same frequency, amplitude and duty cycle). Furthermore, there may be gaps in the datasets and these gaps may be inconsistent across trials. Finally, even if one has clean experimental/observational access to all the state variables and the data are noiseless, the models may exhibit parameter degeneracy (multiple combinations of parameter give rise to the same pattern), rendering them structurally unidentifiable.

Classical structural unidentifiability is associated with the notion that one can at most identify combinations of unidentifiable model parameters [25,26]. This is due to an excess in the number of parameters with physical meaning as compared to the number of parameters necessary to compute the attributes that describe these data. We discussed these ideas in detail using a number of simple representative model examples with increasing levels of complexity in the Appendix A. A more thorough discussion can be found in [25,26] (and references therein).

In this paper, we set out to analyze a different type of structural degeneracy/unidentifiability present in the family of $\mathsf{\Lambda }$-$\mathsf{\Omega }$ models [67,68,69,70,87]. We found that $\mathsf{\Lambda }$-$\mathsf{\Omega }$ models have infinitely many parameter sets that produce identical stable oscillations, except possible for a phase shift (reflecting the initial phase), but this parameters are not identifiable combinations of unidentifiable parameters as is the case in structural degeneracy. The number of model parameters in the $\mathsf{\Lambda }$-$\mathsf{\Omega }$ models is minimal in the sense that each one controls a different aspect of the model dynamics and the dynamic complexity of the system would be reduced by reducing the number of parameters. The level sets and the phase-plane diagrams can be computed analytically as well as the degenerate solutions. This adds clarity to the problem, leaving out other possible causes for unidentifiability (e.g., numerical, algorithmic), and facilitates the analytical investigation of the underlying compensation mechanisms.

Mathematically, degeneracy can be described in terms of level sets in parameter space for a given attribute (e.g., frequency, amplitude, duty cycle). An attribute level set is the set of all combinations of parameter values that produce solutions with a constant attribute [32]. In general, for a given model, level sets for different attributes may or may not coincide or be constrained by one another. For the $\mathsf{\Lambda }$-$\mathsf{\Omega }$ models, the level sets for all possible attributes coincide.

The family of $\mathsf{\Lambda }$-$\mathsf{\Omega }$ models produces a rich dynamic behavior, including monostable oscillations, bistability between oscillations and equilibria, Hopf bifurcations, and saddle-node bifurcations of limit cycles. The classical $\mathsf{\Lambda }$-$\mathsf{\Omega }$ models show sinusoidal oscillations, while the extended $\mathsf{\Lambda }$-$\mathsf{\Omega }$ models exhibit more complex waveforms, including square-wave (relaxation-type) and spike-like oscillations. Therefore, the $\mathsf{\Lambda }$-$\mathsf{\Omega }$ models serve as canonical models for degeneracy/unidentifiability for a large variety of realistic oscillatory processes. In addition, these canonical models can be used to investigate the consequences of degeneracy in single cells for their response to external perturbations [87] and the dynamics of networks in which these cells are embedded.

Degeneracy of oscillatory patterns has been the focus of various studies in the context of neuronal systems both experimentally and theoretically [32,33,64,88,89,90,91,92,93,94]. Degeneracy is an inherent property of dynamical systems and is believed to arise from compensation mechanisms that generate balances among the participating processes (e.g., neuronal ionic currents) that control the dynamics [32,33]. Because the dynamic mechanisms of generation of biological oscillations involve the interplay of positive and negative feedback effects operating at different time scales, unidentifiability is likely to be present in the majority, if not all biological oscillatory systems. Therefore, the $\mathsf{\Lambda }$-$\mathsf{\Omega }$ models can serve as a first, guiding step to understand the mechanisms leading to degeneracy in biological oscillators and biological oscillatory networks [67,68].

The full degeneracy present in the $\mathsf{\Lambda }$-$\mathsf{\Omega }$ models results from the symmetries present in the model, which are captured by the phase-plane diagrams. We found that a systematic break of these symmetries for the $\lambda$-$\omega$ models leading to the phase-plane diagrams characteristic of more realistic models (e.g., relaxation oscillators) preserve the degeneracies of the attributes, although the level sets for different attributes do not necessarily coincide. Therefore, the $\lambda$-$\omega$ models serve as reference models to investigate the mechanisms of generation of attribute level sets for biological oscillators as perturbations to the corresponding $\lambda$-$\omega$ models that exhibits full degeneracy. We hypothesize that this is the case for the whole family of $\mathsf{\Lambda }$-$\mathsf{\Omega }$ models. More research is needed to test this hypothesis and clarify these issues.

In addition to biological oscillators, our results predict that the phenomenon of degeneracy may arise in other oscillatory systems in physics and chemistry since the $\lambda$-$\omega$ models are real-valued special cases of the complex Ginzburg–Landau equation [71,72,73,74,95,96]. Establishing these ideas requires additional research.

Unidentifiability in parameter estimation is a fundamental problem arising in the field of data science in connection to the models needed to analyze the available data. Degeneracy, the other side of the coin, is also a biological fact reported in biological experiments [64,65,66] (see also [32,33]) and is expected to be pervasive in oscillatory systems. While on one hand degeneracy reflects the lack of information necessary to understand the underlying mechanism in term of the participating process, on the other hand, degeneracy has been proposed to endow systems with functional flexibility [97], thus adding a new dimension to the investigation of degenerate/unidentifiable systems. The family of $\mathsf{\Lambda }$-$\mathsf{\Omega }$ models serves as a framework to investigate these issues and is expected to have implications for the data-driven discovery of non-linear dynamic equations [19,20]. Furthermore, the family of $\mathsf{\Lambda }$-$\mathsf{\Omega }$ models serves as the canonical models to calibrate parameter estimation algorithms. In contrast to the classical unidentifiable models, the family of $\mathsf{\Lambda }$-$\mathsf{\Omega }$ models is potentially identifiable if one has enough information about the transient behavior. However, as we pointed out above, this is information is not available in realistic situations. Our results using noisy ground truth data where transients are activated by noise were not useful to resolve the unidentifiability.

The family of $\mathsf{\Lambda }$-$\mathsf{\Omega }$ models also serves as a framework for the systematic investigation of degeneracy in dynamic models and the interplay between structural unidentifiability and the various forms of unidentifiability raised above resulting for lack of experimental/observational access to the state variables. The development of this framework and the implications for the investigation of degeneracy and identifiability in more realistic models require more research.